      program lc                                                        
c  Main program for computing light and radial velocity curves, and     
c      line profiles                                                    
c                                                                       
c  Version of February 2, 1999                                          
C                                                                       
C     TO PRINT VELOCITIES IN KM/SEC, SET VUNIT=1.                       
C     TO PRINT NORMALIZED VELOCITIES IN SAME COLUMNS, SET VUNIT EQUAL TO
C     DESIRED VELOCITY UNIT IN KM/SEC.                                  
C                                                                       
C     PARAMETER PSHIFT IS DEFINED AS THE PHASE AT WHICH PRIMARY         
C     CONJUNCTION (STAR 1 AWAY FROM OBSERVER) WOULD OCCUR IF THE        
C     ARGUMENT OF PERIASTRON WERE 90 DEGREES. SINCE THE NOMINAL VALUE   
C     OF THIS QUANTITY IS ZERO, PSHIFT MAY BE USED TO INTRODUCE AN      
C     ARBITRARY PHASE SHIFT.                                            
C                                                                       
      implicit real*8(a-h,o-z)                                          
      dimension rad(4),drdo(4),xtha(4),xfia(4),po(2)                    
      dimension rv(3011),grx(3011),gry(3011),grz(3011),rvq(3011),       
     $grxq(3011),gryq(3011),grzq(3011),slump1(3011),slump2(3011),       
     $fr1(3011),fr2(3011),glump1(3011),glump2(3011),xx1(3011),          
     $xx2(3011),yy1(3011),yy2(3011),zz1(3011),zz2(3011),grv1(3011),     
     $grv2(3011),rftemp(3011),rf1(3011),rf2(3011),csbt1(3011),          
     $csbt2(3011),gmag1(3011),gmag2(3011),hld(3200),snfi(6400),         
     $csfi(6400),tld(6400),snth(260),csth(260),theta(520),rho(520),     
     $aa(20),bb(20),mmsave(124)                                         
      dimension fbin1(2000),fbin2(2000),delv1(2000),delv2(2000),        
     $count1(2000),count2(2000),delwl1(2000),delwl2(2000),resf1(2000),  
     $resf2(2000),wl1(2000),wl2(2000),dvks1(100),dvks2(100),wll1(100),  
     $wll2(100),tau1(100),tau2(100),ewid1(100),ewid2(100),depth1(100),  
     $depth2(100),hbarw1(100),hbarw2(100),taug(2000)                    
      DIMENSION XLAT(2,100),xlong(2,100)                                
      dimension xcl(100),ycl(100),zcl(100),rcl(100),op1(100),fcl(100),  
     $dens(100),encl(100),edens(100),xmue(100),yskp(14000),zskp(14000)  
      COMMON /FLVAR/ PSHIFT,DP,DS,EF,EFC,ECOS,perr0,PHPER,PCONJ,        
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc,qfacd       
      COMMON /DPDX/ DPDX1,DPDX2,PHSV,PCSV                               
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC     
      COMMON /KFAC/ KFF1,KFF2,kfo1,kfo2                                 
      COMMON /INVAR/ KH,IPB,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,ifsmv2, 
     $icor1,icor2,ld,ncl,jdphs,ipc                                      
      COMMON /SPOTS/SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG    
     $(2,100),RADSP(2,100),temsp(2,100),xlng(2,100),kks(2,100),         
     $Lspot(2,100)                                                      
      common /cld/ acm,opsf                                             
      common /ardot/ dperdt,hjd,hjd0,perr                               
      common /prof2/ du1,du2,du3,du4,binw1,binw2,sc1,sc2,sl1,sl2,       
     $clight                                                            
      common /inprof/ in1min,in1max,in2min,in2max,mpage,nl1,nl2         
      common /ipro/ nbins,nl,inmax,inmin,nf1,nf2                        
      COMMON /NSPT/ NSP1,NSP2,IFAT1,IFAT2                               
      data xtha(1),xtha(2),xtha(3),xtha(4),xfia(1),xfia(2),xfia(3),     
     $xfia(4)/0.d0,1.570796d0,1.570796d0,1.570796d0,                    
     $0.d0,0.d0,1.5707963d0,3.14159365d0/                               
  205 format('**********************************************************
     $************')                                                    
   79 format(6x,'JD',17x,'Phase     light 1     light 2     (1+2+3)    n
     $orm lite   dist      mag+K')                                      
   96 FORMAT(6x,'JD',13x,'Phase',5x,'r1pol',6x,'r1pt',5x,'r1sid',5x,'r1b
     $ak',5x,'r2pol',5x,'r2pt',6x,'r2sid',5x,'r2bak')                   
  296 format(f14.6,f13.5,8f10.5)                                        
   45 FORMAT(6x,'JD',14x,'Phase     V Rad 1     V Rad 2      del V1     
     $ del V2   V1 km/s      V2 km/s')                                  
   93 format(f14.6,f13.5,4f12.6,2d13.4)                                 
   47 FORMAT(2x,'wv lth',8x,'L1',9x,'L2',7x,'x1',6x,'x2',6x,'y1',6x,    
     $'y2',6x,'el3     opsf      m zero   factor')                      
   48 FORMAT('  ecc',4x,'s-m axis',6x,'F1',9x,'F2',7x,'Vgam',7x,        
     $'Incl',6x,'g1',6x,'g2  Nspot1 Nspot 2')                           
   54 FORMAT(2x,'T1',6x,'T2',5x,'Alb 1  Alb 2',4x,'Pot 1',8x,'Pot 2',   
     $8x,'M2/M1',2x,'x1(bolo) x2(bolo) y1(bolo) y2(bolo)')              
   33 FORMAT(I4,I5,I6,I6,I6,I4,f13.6,d14.5,f9.5,f10.2,d16.4)            
   74 FORMAT(' DIMENSIONLESS RADIAL VELOCITIES CONTAIN FACTOR P/(2PI*A)'
     $)                                                                 
   46 FORMAT('GRID1/4    GRID2/4',2X,'POLAR SBR 1',3X,'POLAR SBR 2'     
     $,3X,'SURF. AREA 1',2X,'SURF. AREA 2',3X,'PERI. PH.',3X,'CONJ. PH. 
     $')                                                                
   50 FORMAT(40H PRIMARY COMPONENT EXCEEDS CRITICAL LOBE)               
   51 FORMAT(42H SECONDARY COMPONENT EXCEEDS CRITICAL LOBE)             
   42 FORMAT(1H )                                                       
   41 FORMAT('star',4X,'r pole',5X,'deriv',5X,'r point',5X,'deriv',     
     $6X,'r side',6X,'deriv',5X,'r back',6X,'deriv')                    
   80 FORMAT(1X,F6.4,1X,F7.5)                                           
    2 FORMAT(F6.5,F10.4,2F10.4,F10.4,f9.3,2f7.3)                        
    5 FORMAT(F6.5,F11.4,2F11.4,F11.4,F10.3,2f8.3,i5,i7)                 
    6 FORMAT(2(F7.4,1X),2f7.3,2d12.5,f10.5,4F7.3)                       
    8 FORMAT(f7.4,f8.4,2F7.3,2d13.5,f10.5,f8.3,f9.3,f9.3,f9.3)          
    3 FORMAT(f15.6,F15.5,4F12.8,F10.5,f10.4)                            
    1 FORMAT(4I2,2I3,f13.6,d12.5,f7.5,F8.2)                             
    4 FORMAT(F9.6,2F10.5,4F7.3,F8.4,d10.4,F8.3,F8.4)                    
   34 FORMAT(F9.6,2F11.5,4f8.3,F9.4,d11.4,F9.3,F9.4)                    
   49 FORMAT(' PROGRAM SHOULD NOT BE USED IN MODE 1 OR 3 WITH NON-ZERO E
     $CCENTRICITY')                                                     
   10 FORMAT('MODE   IPB  IFAT1 IFAT2  N1  N2',4x,'Arg. Per',7x,'dPerdt 
     $',4x,'Th e',4x,'V UNIT(km/s)    V FAC')                           
  148 format('   mpage  nref   mref   ifsmv1   ifsmv2   icor1   icor2   
     $ld')                                                              
  171 format('JDPHS',5x,'J.D. zero',7x,'Period',11x,'dPdt',             
     $6x,'Ph. shift',3x,'fract. sd.',2x,'noise',5x,'seed')              
  244 format('Note: The light curve output contains simulated observa', 
     $'tional scatter, as requested,')                                  
  245 format('with standard deviation',f9.5,' of light at the reference'
     $,' phase.')                                                       
  149 format(i6,2i7,i8,i9,i9,i8,i6)                                     
  170 format(i3,f17.6,d18.10,d14.6,f10.4,d13.4,i5,f13.0)                
   40 FORMAT(I3,8F11.5)                                                 
   94 FORMAT(i6,i11,4F14.6,F12.5,F12.5)                                 
   84 FORMAT(1X,I4,4F12.5)                                              
   85 FORMAT(4f9.5)                                                     
   83 FORMAT(1X,'STAR  CO-LATITUDE  LONGITUDE  SPOT RADIUS  TEMP. FACTOR
     $')                                                                
  150 format(' Star     M/Msun   (Mean Radius)/Rsun   M Bol    Log g (cg
     $s)')                                                              
  250 format(4x,I1,4x,f6.2,11x,f7.2,6x,f6.2,5x,f5.2)                    
  350 format(' Primary star exceeds outer contact surface')             
  351 format(' Secondary star exceeds outer contact surface')           
   22 format(8(i1,1x))                                                  
  649 format(i1,f15.6,d15.10,d13.6,f10.4,d10.4,i2,f11.0)                
   63 format(3f9.4,f7.4,d11.4,f9.4,d11.3,f9.4,f7.3)                     
   64 format(3f10.4,f9.4,d12.4,f10.4,d12.4,f9.4,f9.3,d12.4)             
   69 format('      xcl       ycl       zcl      rcl       op1         f
     $cl        ne       mu e      encl     dens')                      
 2048 format(d11.5,f9.4,f9.2,i3)                                        
 2049 format(i3,d14.5,f18.2,f20.2,i14)                                  
  907 format(6x,'del v',6x,'del wl (mic.)',7x,'wl',9x,'profile',6x,'res 
     $flux')                                                            
  903 format(6f14.7)                                                    
   92 format('Phase =',f14.6)                                           
  142 format('star',4x,'bin width (microns)',3x,'continuum scale',4x,'co
     $ntinuum slope',2x,'nfine')                                        
  167 format(30x,'star',i2)                                             
  138 format(f9.6,d12.5,f10.5,i5)                                       
  152 format(f20.6,d23.5,17x,f13.5,i6)                                  
  157 format('star ',i1,'   line wavelength',4x,'equivalent width (micro
     $ns)',5x,'rect. line depth',2x,'kks')                              
  217 format(f14.6,f15.6,f13.6,4f10.4)                                  
  218 format(f14.6,f16.6,f14.6,4f10.4)                                  
  219 format(5x,'JD start',9x,'JD stop',6x,'JD incr',4x,                
     $'Ph start',2x,'Ph. stop',3x,'Ph incr',3x,'Ph norm')               
      ot=1.d0/3.d0                                                      
      KH=17                                                             
      pi=dacos(-1.d0)                                                   
      clight=2.99791d5                                                  
      en0=6.0254d23                                                     
      nn=100                                                            
      gau=0.d0                                                          
      DO 1000 IT=1,1000                                                 
      read(5,22) mpage,nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld           
      if(mpage.eq.9) stop                                               
      read(5,649) jdphs,hjd0,period,dpdt,pshift,stdev,noise,seed        
      read(5,217) hjdst,hjdsp,hjdin,phstrt,phstop,phin,phn              
      READ(5,1) MODE,IPB,IFAT1,IFAT2,N1,N2,perr0,dperdt,the,VUNIT       
      READ(5,2) E,A,F1,F2,VGA,XINCL,GR1,GR2                             
      read(5,6) tavh,tavc,alb1,alb2,poth,potc,rm,xbol1,xbol2,ybol1,     
     $ybol2                                                             
      READ(5,4)WL,HLUM,CLUM,XH,xc,yh,yc,EL3,opsf,ZERO,FACTOR            
      acm=6.960d10*a                                                    
      nf1=1                                                             
      nf2=1                                                             
      if(mpage.ne.3) goto 897                                           
      colam=clight/wl                                                   
      read(5,2048) binwm1,sc1,sl1,nf1                                   
      binw1=colam*binwm1                                                
      do 86 iln=1,100                                                   
      read(5,138) wll1(iln),ewid1(iln),depth1(iln),kks(1,iln)           
      if(wll1(iln).lt.0.d0) goto 89                                     
      tau1(iln)=-dlog(1.d0-depth1(iln))                                 
      hbarw1(iln)=0.d0                                                  
      if(depth1(iln).ne.0.d0) hbarw1(iln)=.5d0*clight*ewid1(iln)/       
     $(wll1(iln)*dabs(depth1(iln)))                                     
      nl1=iln                                                           
   86 continue                                                          
   89 continue                                                          
      read(5,2048) binwm2,sc2,sl2,nf2                                   
      binw2=colam*binwm2                                                
      do 99 iln=1,100                                                   
      read(5,138) wll2(iln),ewid2(iln),depth2(iln),kks(2,iln)           
      if(wll2(iln).lt.0.d0) goto 91                                     
      tau2(iln)=-dlog(1.d0-depth2(iln))                                 
      hbarw2(iln)=0.d0                                                  
      if(depth2(iln).ne.0.d0) hbarw2(iln)=.5d0*clight*ewid2(iln)/       
     $(wll2(iln)*dabs(depth2(iln)))                                     
      nl2=iln                                                           
   99 continue                                                          
   91 continue                                                          
      do 622 iln=1,nl1                                                  
      flam=(wll1(iln)/wl)**2                                            
  622 dvks1(iln)=clight*(flam-1.d0)/(flam+1.d0)                         
      do 623 iln=1,nl2                                                  
      flam=(wll2(iln)/wl)**2                                            
  623 dvks2(iln)=clight*(flam-1.d0)/(flam+1.d0)                         
  897 continue                                                          
      NSP1=0                                                            
      NSP2=0                                                            
      DO 88 KP=1,2                                                      
      DO 87 I=1,100                                                     
      READ(5,85)XLAT(KP,I),XLONG(KP,I),RADSP(KP,I),TEMSP(KP,I)          
      xlng(kp,i)=xlong(kp,i)                                            
      IF(XLAT(KP,I).GE.200.d0) GOTO 88                                  
      SINLAT(KP,I)=dsin(XLAT(KP,I))                                     
      COSLAT(KP,I)=dcos(XLAT(KP,I))                                     
      SINLNG(KP,I)=dsin(XLONG(KP,I))                                    
      COSLNG(KP,I)=dcos(XLONG(KP,I))                                    
      IF(KP.EQ.1)NSP1=NSP1+1                                            
   87 IF(KP.EQ.2)NSP2=NSP2+1                                            
   88 CONTINUE                                                          
      ncl=0                                                             
      do 62 i=1,100                                                     
      read(5,63) xcl(i),ycl(i),zcl(i),rcl(i),op1(i),fcl(i),edens(i),    
     $xmue(i),encl(i)                                                   
      if(xcl(i).gt.100.d0) goto 66                                      
      ncl=ncl+1                                                         
      dens(i)=edens(i)*xmue(i)/en0                                      
   62 continue                                                          
   66 dint1=pi*(1.d0-xbol1/3.d0)                                        
      dint2=pi*(1.d0-xbol2/3.d0)                                        
      if(ld.eq.2) DINT1=dint1+PI*2.d0*ybol1/9.d0                        
      if(ld.eq.2) DINT2=dint2+PI*2.d0*ybol2/9.d0                        
      if(ld.eq.3) dint1=dint1-.2d0*pi*ybol1                             
      if(ld.eq.3) dint2=dint2-.2d0*pi*ybol2                             
      NSTOT=NSP1+NSP2                                                   
      NP1=N1+1                                                          
      NP2=N1+N2+2                                                       
      IRTE=0                                                            
      IRVOL1=0                                                          
      IRVOL2=0                                                          
      CALL SINCOS(1,N1,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)                   
      CALL SINCOS(2,N2,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)                   
      hjd=hjd0                                                          
      CALL modlog(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD, 
     $rm,poth,potc,gr1,gr2,alb1,alb2,n1,n2,f1,f2,mod,xincl,the,mode,    
     $snth,csth,snfi,csfi,grv1,grv2,xx1,yy1,zz1,xx2,yy2,zz2,glump1,     
     $glump2,csbt1,csbt2,gmag1,gmag2)                                   
      CALL VOLUME(VOL1,RM,POTH,DP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,        
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,  
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,      
     $GMAG1,GMAG2,GR1,1)                                                
      CALL VOLUME(VOL2,RM,POTC,DP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,        
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,  
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,      
     $GMAG1,GMAG2,GR2,1)                                                
      if(e.eq.0.d0) goto 117                                            
      DAP=1.d0+E                                                        
      P1AP=POTH-2.d0*E*RM/(1.d0-E*E)                                    
      VL1=VOL1                                                          
      CALL VOLUME(VL1,RM,P1AP,DAP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,        
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,  
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,      
     $GMAG1,GMAG2,GR1,2)                                                
      DPDX1=(POTH-P1AP)*(1.d0-E*E)*.5d0/E                               
      P2AP=POTC-2.d0*E/(1.d0-E*E)                                       
      VL2=VOL2                                                          
      CALL VOLUME(VL2,RM,P2AP,DAP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,        
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,  
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,      
     $GMAG1,GMAG2,GR2,2)                                                
      DPDX2=(POTC-P2AP)*(1.d0-E*E)*.5d0/E                               
  117 CONTINUE                                                          
      PHSV=POTH                                                         
      PCSV=POTC                                                         
      IF(E.EQ.0.d0) GOTO 61                                             
      IF(MOD.EQ.1) WRITE(6,49)                                          
   61 CONTINUE                                                          
      CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,    
     $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POTH,POTC,N1,N2,F1,F2,D,HLUM        
     $,clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1, 
     $ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phn,rm,xincl,   
     $hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,   
     $zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,csbt2,gmag1,   
     $gmag2,fbin1,fbin2,delv1,delv2,count1,count2,delwl1,delwl2,        
     $resf1,resf2,wl1,wl2,dvks1,dvks2,tau1,tau2,hbarw1,hbarw2,          
     $xcl,ycl,zcl,rcl,op1,fcl,dens,encl,edens,taug,yskp,zskp,mode,1)    
      KH=0                                                              
      if(kfo1.eq.0) goto 380                                            
      write(6,350)                                                      
      goto 381                                                          
  380 IF(KFF1.EQ.1) WRITE(6,50)                                         
  381 if(kfo2.eq.0) goto 382                                            
      write(6,351)                                                      
      goto 383                                                          
  382 IF(KFF2.EQ.1) WRITE(6,51)                                         
  383 IF((KFF1+KFF2+kfo1+kfo2).GT.0) WRITE(6,42)                        
      write(6,148)                                                      
      write(6,149) mpage,nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld         
      write(6,42)                                                       
      write(6,171)                                                      
      write(6,170) jdphs,hjd0,period,dpdt,pshift,stdev,noise,seed       
      write(6,42)                                                       
      write(6,219)                                                      
      write(6,218) hjdst,hjdsp,hjdin,phstrt,phstop,phin,phn             
      write(6,42)                                                       
      WRITE(6,10)                                                       
      WRITE(6,33)MODE,IPB,IFAT1,IFAT2,N1,N2,perr0,dperdt,the,VUNIT,vfac 
      WRITE(6,42)                                                       
      WRITE(6,48)                                                       
      WRITE(6,5)E,A,F1,F2,VGA,XINCL,GR1,GR2,NSP1,NSP2                   
      WRITE(6,42)                                                       
      WRITE(6,54)                                                       
      WRITE(6,8)TAVH,TAVC,ALB1,ALB2,PHSV,PCSV,rm,XBOL1,xbol2,ybol1,     
     $ybol2                                                             
      WRITE(6,42)                                                       
      WRITE(6,47)                                                       
      WRITE(6,34)WL,HLUM,CLUM,XH,XC,yh,yc,el3,opsf,ZERO,FACTOR          
      ns1=1                                                             
      ns2=2                                                             
      if(mpage.ne.3) goto 174                                           
      write(6,42)                                                       
      write(6,142)                                                      
      write(6,2049) ns1,binwm1,sc1,sl1,nf1                              
      write(6,2049) ns2,binwm2,sc2,sl2,nf2                              
      write(6,42)                                                       
      write(6,157) ns1                                                  
      do 155 iln=1,nl1                                                  
  155 write(6,152) wll1(iln),ewid1(iln),depth1(iln),kks(1,iln)          
      write(6,42)                                                       
      write(6,157) ns2                                                  
      do 151 iln=1,nl2                                                  
  151 write(6,152) wll2(iln),ewid2(iln),depth2(iln),kks(2,iln)          
  174 continue                                                          
      write(6,42)                                                       
      WRITE(6,42)                                                       
      IF(NSTOT.GT.0) WRITE(6,83)                                        
      DO 188 KP=1,2                                                     
      IF((NSP1+KP-1).EQ.0) GOTO 188                                     
      IF((NSP2+(KP-2)**2).EQ.0) GOTO 188                                
      NSPOT=NSP1                                                        
      IF(KP.EQ.2) NSPOT=NSP2                                            
      DO 187 I=1,NSPOT                                                  
  187 WRITE(6,84)KP,XLAT(KP,I),XLONG(KP,I),RADSP(KP,I),TEMSP(KP,I)      
  188 WRITE(6,42)                                                       
      if(ncl.eq.0) goto 67                                              
      write(6,69)                                                       
      do 68 i=1,ncl                                                     
   68 write(6,64) xcl(i),ycl(i),zcl(i),rcl(i),op1(i),fcl(i),edens(i),   
     $xmue(i),encl(i),dens(i)                                           
      write(6,42)                                                       
   67 continue                                                          
      write(6,150)                                                      
      rr1=.6203505d0*vol1**ot                                           
      rr2=.6203505d0*vol2**ot                                           
      tav1=10000.d0*tavh                                                
      tav2=10000.d0*tavc                                                
      call mlrg(a,period,rm,rr1,rr2,tav1,tav2,sms1,sms2,sr1,sr2,        
     $bolm1,bolm2,xlg1,xlg2)                                            
      write(6,250) ns1,sms1,sr1,bolm1,xlg1                              
      write(6,250) ns2,sms2,sr2,bolm2,xlg2                              
      write(6,42)                                                       
      WRITE(6,46)                                                       
      WRITE(6,94) MMSAVE(NP1),MMSAVE(NP2),SBRH,SBRC,SM1,SM2,PHPERI,     
     $PCONJ                                                             
      WRITE(6,42)                                                       
      if(stdev.eq.0.d0.or.mpage.ne.1) goto 246                          
      write(6,244)                                                      
      write(6,245) stdev                                                
  246 continue                                                          
      WRITE(6,42)                                                       
      ALL=HOT+COOL+EL3                                                  
      IF(MODE.EQ.-1) ALL=COOL+EL3                                       
      if(mpage.eq.1) write(6,79)                                        
      if(mpage.eq.2) write(6,45)                                        
      if(mpage.eq.4) write(6,96)                                        
      LL1=MMSAVE(N1)+1                                                  
      NPP2=NP2-1                                                        
      LL2=MMSAVE(NPP2)+1                                                
      LLL1=MMSAVE(NP1)                                                  
      LLL2=MMSAVE(NP2)                                                  
      LLLL1=(LL1+LLL1)/2                                                
      LLLL2=(LL2+LLL2)/2                                                
      POTH=PHSV                                                         
      POTC=PCSV                                                         
      PO(1)=POTH                                                        
      PO(2)=POTC                                                        
      IF(E.EQ.0.d0) IRVOL1=1                                            
      IF(E.EQ.0.d0) IRVOL2=1                                            
      IF(E.EQ.0.d0) IRTE=1                                              
      start=hjdst                                                       
      stop=hjdsp                                                        
      step=hjdin                                                        
      if(jdphs.ne.2) goto 887                                           
      start=phstrt                                                      
      stop=phstop                                                       
      step=phin                                                         
  887 continue                                                          
      do 20 phjd=start,stop,step                                        
      hjdi=phjd                                                         
      phasi=phjd                                                        
      call jdph(hjdi,phasi,hjd0,period,dpdt,hjdo,phaso)                 
      hjd=hjdi                                                          
      phas=phasi                                                        
      if(jdphs.ne.1) hjd=hjdo                                           
      if(jdphs.ne.2) phas=phaso                                         
      CALL modlog(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD, 
     $rm,poth,potc,gr1,gr2,alb1,alb2,n1,n2,f1,f2,mod,xincl,the,mode,    
     $snth,csth,snfi,csfi,grv1,grv2,xx1,yy1,zz1,xx2,yy2,zz2,glump1,     
     $glump2,csbt1,csbt2,gmag1,gmag2)                                   
      CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,    
     $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POTH,POTC,N1,N2,F1,F2,D,hlum,       
     $clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1,  
     $ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phas,rm,xincl,  
     $hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,   
     $zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,csbt2,gmag1,   
     $gmag2,fbin1,fbin2,delv1,delv2,count1,count2,delwl1,delwl2,resf1,  
     $resf2,wl1,wl2,dvks1,dvks2,tau1,tau2,hbarw1,hbarw2,                
     $xcl,ycl,zcl,rcl,op1,fcl,dens,encl,edens,taug,yskp,zskp,mode,0)    
  128 format('HJD = ',f14.5,'    Phase = ',f14.5)                       
  131 format(3x,'Y Sky Coordinate',4x,'Z Sky Coordinate')               
  130 format(f16.6,f20.6)                                               
      if(mpage.ne.5) goto 127                                           
      write(6,42)                                                       
      write(6,42)                                                       
      write(6,128) hjd,phas                                             
      write(6,42)                                                       
      write(6,131)                                                      
      do 129 imp=1,ipc                                                  
      write(6,130) yskp(imp),zskp(imp)                                  
  129 continue                                                          
      goto 20                                                           
  127 continue                                                          
      HTT=HOT                                                           
      IF(MODE.EQ.-1) HTT=0.d0                                           
      TOTAL=HTT+COOL+EL3                                                
      TOTALL=TOTAL/ALL                                                  
      TOT=TOTALL*FACTOR                                                 
      if(stdev.le.0.d0) goto 348                                        
      call rangau(seed,nn,stdev,gau)                                    
      ranf=1.d0+gau*dsqrt(totall**noise)                                
      total=total*ranf                                                  
      tot=tot*ranf                                                      
      totall=totall*ranf                                                
  348 continue                                                          
      SMAGG=-1.085736d0*dlog(TOTALL)+ZERO                               
      if(mpage.eq.1) write(6,3) hjd,phas,htt,cool,total,tot,d,smagg     
      if(mpage.eq.2) write(6,93) hjd,phas,vsum1,vsum2,vra1,vra2,vkm1,   
     $vkm2                                                              
      if(mpage.ne.3) goto 81                                            
      write(6,92) phas                                                  
      write(6,42)                                                       
      write(6,167) ns1                                                  
      write(6,907)                                                      
      do 906 i=in1min,in1max                                            
  906 write(6,903) delv1(i),delwl1(i),wl1(i),fbin1(i),resf1(i)          
      write(6,42)                                                       
      write(6,167) ns2                                                  
      write(6,907)                                                      
      do 908 i=in2min,in2max                                            
  908 write(6,903) delv2(i),delwl2(i),wl2(i),fbin2(i),resf2(i)          
      write(6,42)                                                       
      write(6,205)                                                      
      write(6,42)                                                       
      write(6,42)                                                       
   81 continue                                                          
      if(mpage.eq.4) write(6,296) hjd,phas,rv(1),rv(ll1),rv(llll1),     
     $rv(lll1),rvq(1),rvq(ll2),rvq(llll2),rvq(lll2)                     
   20 CONTINUE                                                          
      if(mpage.eq.5) stop                                               
      WRITE(6,42)                                                       
      WRITE(6,41)                                                       
      WRITE(6,42)                                                       
      do 119 ii=1,2                                                     
      gt1=dfloat(2-ii)                                                  
      gt2=dfloat(ii-1)                                                  
      f=f1*gt1+f2*gt2                                                   
      do 118 i=1,4                                                      
      call romq(po(ii),rm,f,dp,e,xtha(i),xfia(i),rad(i),drdo(i),        
     $drdq,dodq,ii,mode)                                                
  118 continue                                                          
      write(6,40) ii,rad(1),drdo(1),rad(2),drdo(2),rad(3),drdo(3),      
     $rad(4),drdo(4)                                                    
  119 continue                                                          
      if(mpage.ne.2) goto 1000                                          
      WRITE(6,42)                                                       
      WRITE(6,74)                                                       
 1000 CONTINUE                                                          
      STOP                                                              
      END                                                               
      Subroutine light(phs,xincl,xh,xc,yh,yc,n1,n2,sumhot,sumkul,rv,grx,
     $gry,grz,rvq,grxq,gryq,grzq,mmsave,theta,rho,aa,bb,slump1,slump2,  
     $somhot,somkul,d,wl,snth,csth,snfi,csfi,tld,gmag1,gmag2,fbin1,     
     $fbin2,delv1,delv2,count1,count2,delwl1,delwl2,resf1,resf2,wl1,    
     $wl2,dvks1,dvks2,tau1,tau2,hbarw1,hbarw2,xcl,ycl,zcl,rcl,op1,fcl,  
     $edens,encl,dens,taug,yskp,zskp,ifphn)                             
c   Version of October 10, 2000                                       
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),THETA(*),RHO(*),AA(*),BB(*)       
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*),tld(*),gmag1(*),gmag2(*)
      dimension xcl(*),ycl(*),zcl(*),rcl(*),op1(*),fcl(*),dens(*),      
     $encl(*),edens(*),yskp(*),zskp(*)                                  
      dimension fbin1(*),fbin2(*),delv1(*),delv2(*),count1(*),count2(*),
     $delwl1(*),delwl2(*),resf1(*),resf2(*),wl1(*),wl2(*),dvks1(*),     
     $dvks2(*),tau1(*),tau2(*),hbarw1(*),hbarw2(*),taug(*)              
      COMMON X1                                                         
      COMMON /KFAC/ KFF1,KFF2,kfo1,kfo2                                 
      COMMON /NSPT/ NSP1,NSP2,IFAT1,IFAT2                               
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm    
     $,ifs1dm,ifs2dm,icr1dm,icr2dm,ld,ncl,jdphs,ipc                     
      common /flvar/ du2,du3,du4,du5,du6,du7,du8,du9,du10,du11,         
     $du12,du13,du14,du15,du16,du17,vunit,vfvu,du20,qfacd               
      common /prof2/ vo1,vo2,ff1,ff2,binw1,binw2,sc1,sc2,sl1,sl2,       
     $clight                                                            
      common /cld/ acm,opsf                                             
      common /inprof/ in1min,in1max,in2min,in2max,mpage,nl1,nl2         
      common /setest/ sefac                                             
      common /flpro/ vksf,binc,binw,difp,deldel,renfsq                  
      common /ipro/ nbins,nl,inmax,inmin,nf1,nf2                        
      COMMON /SPOTS/ SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG   
     $(2,100),RAD(2,100),TEMSP(2,100),xlng(2,100),kks(2,100),           
     $Lspot(2,100)                                                      
      pi=dacos(-1.d0)                                                   
      twopi=pi+pi                                                       
      pih=.5d0*pi                                                       
      dtr=pi/180.d0                                                     
      kstp=4                                                            
      cirf=.002d0                                                       
      if(ifphn.eq.1) goto 16                                            
      if(mpage.ne.3) goto 16                                            
      nbins=2000                                                        
      binc1=.5d0*dfloat(nbins)                                          
      binc2=binc1                                                       
      in1max=0                                                          
      in2max=0                                                          
      in1min=30000                                                      
      in2min=30000                                                      
      marm1=10                                                          
      marp1=10                                                          
      marm2=10                                                          
      marp2=10                                                          
      clfac=clight/vunit                                                
      do 916 i=1,nbins                                                  
      fbin1(i)=0.d0                                                     
      fbin2(i)=0.d0                                                     
      count1(i)=0.d0                                                    
      count2(i)=0.d0                                                    
      delv1(i)=0.d0                                                     
      delv2(i)=0.d0                                                     
  916 continue                                                          
   16 continue                                                          
      PHA=PHS*twopi                                                     
      K=6                                                               
      KK=K+1                                                            
      XLUMP=14384.d0/WL                                                 
      XINC=XINCL*dtr                                                    
      L=1                                                               
      TEST=(PHS-.5d0)**2                                                
      TESTS=(TEST-.071525d0)**2                                         
      SINI=dsin(XINC)                                                   
      COSPH=dcos(PHA)                                                   
      SINPH=dsin(PHA)                                                   
      SINSQ=SINPH**2                                                    
      COSI=dcos(XINC)                                                   
      NP1=N1+1                                                          
      NP2=N1+N2+2                                                       
      LLL1=MMSAVE(NP1)                                                  
      LLL2=MMSAVE(NP2)                                                  
      NPP2=NP2-1                                                        
      LL1=MMSAVE(N1)+1                                                  
      LL2=MMSAVE(NPP2)+1                                                
      LLLL1=(LL1+LLL1)/2                                                
      LLLL2=(LL2+LLL2)/2                                                
      SINSQE=0.d0                                                       
      IF(SINI.GT.0.d0) SINSQE=((1.10d0*(RV(LLL1)+RVQ(LLL2))/D)**2       
     $-cosi**2)/SINI**2                                                 
      CICP=COSI*COSPH                                                   
      CISP=COSI*SINPH                                                   
      XLOS=COSPH*SINI                                                   
      YLOS=-SINPH*SINI                                                  
      ZLOS=COSI                                                         
      SUM=0.d0                                                          
      SOM=0.d0                                                          
      IF(TEST.LE..0625d0) GOTO 18                                       
      COMP=-1.d0                                                        
      CMP=1.d0                                                          
      COMPP=1.d0                                                        
      KOMP=2                                                            
      nl=nl2                                                            
      ffc=ff2                                                           
      voc=vo2*vfvu                                                      
      NSPOT=NSP2                                                        
      IFAT=IFAT2                                                        
      CMPP=0.d0                                                         
      X=XC                                                              
      y=yc                                                              
      EN=N2                                                             
      NPH=N2                                                            
      NP=2*N2                                                           
      nf=nf2                                                            
      GOTO 28                                                           
   18 X=XH                                                              
      y=yh                                                              
      COMP=1.d0                                                         
      KOMP=1                                                            
      nl=nl1                                                            
      ffc=ff1                                                           
      voc=vo1*vfvu                                                      
      NSPOT=NSP1                                                        
      IFAT=IFAT1                                                        
      CMP=0.d0                                                          
      COMPP=-1.d0                                                       
      CMPP=1.d0                                                         
      EN=N1                                                             
      NPH=N1                                                            
      NP=2*N1                                                           
      nf=nf1                                                            
   28 DELTH=pih/EN                                                      
      enf=dfloat(nf)                                                    
      renfsq=1.d0/(enf*enf)                                             
      nfm1=nf-1                                                         
      r2nfdt=0.5d0*delth/enf                                            
      vfvuff=vfvu*ffc                                                   
      AR=CMPP*RV(LLLL1)+CMP*RVQ(LLLL2)                                  
      BR=CMPP*RV(1)+CMP*RVQ(1)                                          
      ASQ=AR*AR                                                         
      BSQ=BR*BR                                                         
      AB=AR*BR                                                          
      absq=ab*ab                                                        
      ASBS=ASQ-BSQ                                                      
      KF=(2-KOMP)*KFF1+(KOMP-1)*KFF2                                    
      CMPPD=CMPP*D                                                      
      CMPD=CMP*D                                                        
      NPP=NP+1                                                          
      TEMF=1.d0                                                         
      ipc=0                                                             
      DO 36 I=1,NP                                                      
      IF(I.GT.NPH)GOTO 54                                               
      UPDOWN=1.d0                                                       
      IK=I                                                              
      GOTO 55                                                           
   54 UPDOWN=-1.d0                                                      
      IK=NPP-I                                                          
   55 CONTINUE                                                          
      IPN1=IK+(KOMP-1)*N1                                               
      SINTH=SNTH(IPN1)                                                  
      COSTH=CSTH(IPN1)*UPDOWN                                           
      tanth=sinth/costh                                                 
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      XM=MM                                                             
      MH=MM                                                             
      MM=2*MM                                                           
      DELFI=pi/XM                                                       
      r2nfdf=.5d0/enf                                                   
      deldel=delth*delfi                                                
      IP=(KOMP-1)*NP1+IK                                                
      IY=MMSAVE(IP)+1                                                   
      IF(TEST.LE..0625d0)GOTO 19                                        
      GX=GRXQ(IY)                                                       
      GY=-GRYQ(IY)                                                      
      GZ=UPDOWN*GRZQ(IY)                                                
      grmag=gmag2(iy)                                                   
      GOTO 29                                                           
   19 GX=GRX(IY)                                                        
      GY=-GRY(IY)                                                       
      GZ=UPDOWN*GRZ(IY)                                                 
      grmag=gmag1(iy)                                                   
   29 COSSAV=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                            
      SUMJ=0.d0                                                         
      SOMJ=0.d0                                                         
      MPP=MM+1                                                          
      IY=IY-1                                                           
      DO 26 J=1,MM                                                      
      IF(J.GT.MH) GOTO 58                                               
      RTLEFT=1.d0                                                       
      JK=J                                                              
      GOTO 59                                                           
   58 RTLEFT=-1.d0                                                      
      JK=MPP-J                                                          
   59 CONTINUE                                                          
      IX=IY+JK                                                          
      IS=IX+(KOMP-1)*LLL1                                               
      SINFI=SNFI(IS)*RTLEFT                                             
      COSFI=CSFI(IS)                                                    
      STSF=SINTH*SINFI                                                  
      STCF=SINTH*COSFI                                                  
      IF(TEST.LE..0625d0)GOTO 39                                        
      IF(RVQ(IX).EQ.-1.d0) GOTO 26                                      
      GX=GRXQ(IX)                                                       
      GY=RTLEFT*GRYQ(IX)                                                
      GZ=UPDOWN*GRZQ(IX)                                                
      R=RVQ(IX)                                                         
      grmag=gmag2(ix)                                                   
      GOTO 49                                                           
   39 IF(RV(IX).EQ.-1.d0) GOTO 26                                       
      GX=GRX(IX)                                                        
      GY=RTLEFT*GRY(IX)                                                 
      GZ=UPDOWN*GRZ(IX)                                                 
      R=RV(IX)                                                          
      grmag=gmag1(ix)                                                   
   49 COSGAM=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                            
      ZZ=R*COSTH                                                        
      YY=R*COMP*STSF                                                    
      XX=CMPD+COMP*STCF*R                                               
      if(mpage.ne.5) goto 174                                           
      if(cosgam.gt.0.d0) goto 174                                       
      ipc=ipc+1                                                         
      yskp(ipc)=(xx-qfacd)*sinph+yy*cosph                               
      zskp(ipc)=(-xx+qfacd)*cicp+yy*cisp+zz*sini                        
      if(nspot.eq.0) goto 174                                           
      call spot(komp,nspot,sinth,costh,sinfi,cosfi,temf)                
      if(temf.eq.1.d0) goto 174                                         
      yskr=yskp(ipc)                                                    
      zskr=zskp(ipc)                                                    
      kstp=4                                                            
      cirf=.002d0                                                       
      stp=twopi/dfloat(kstp)                                            
      do 179 ang=stp,twopi,stp                                          
      ipc=ipc+1                                                         
      yskp(ipc)=yskr+dsin(ang)*cirf                                     
      zskp(ipc)=zskr+dcos(ang)*cirf                                     
  179 continue                                                          
  174 continue                                                          
      if(sinsq.gt.sinsqe) goto 27                                       
      IF(TESTS.LT.2.2562d-3) GOTO 170                                   
      IF((STCF*R).GT.(sefac*(CMP+COMP*X1))) GOTO 129                    
  170 PROD=COSSAV*COSGAM                                                
      IF(PROD.GT.0.d0) GOTO 22                                          
      COSSAV=-COSSAV                                                    
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                 
      ZSKY=-XX*CICP+yy*CISP+ZZ*SINI+CMPD*CICP                           
      RHO(L)=dsqrt(YSKY**2+ZSKY**2)                                     
      THETA(L)=dasin(ZSKY/RHO(L))                                       
      IF(YSKY.LT.0.d0) GOTO 92                                          
      THETA(L)=twopi+THETA(L)                                           
      GOTO 93                                                           
   92 THETA(L)=pi-THETA(L)                                              
   93 IF (THETA(L).GE.twopi) THETA(L)=THETA(L)                          
     $-twopi                                                            
      L=L+1                                                             
      GOTO 27                                                           
   22 COSSAV=COSGAM                                                     
      GOTO 27                                                           
  129 COSSAV=COSGAM                                                     
      IF(KF.LE.0) GOTO 27                                               
      ZZ=R*COSTH                                                        
      YY=R*COMP*STSF                                                    
      XX=CMPD+COMP*STCF*R                                               
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                 
      ZSKY=-XX*CICP+YY*CISP+ZZ*SINI+CMPD*CICP                           
      rptsq=YSKY**2+ZSKY**2                                             
      rtstsq=absq/(BSQ+ASBS*(ZSKY**2/rptsq))                            
      IF(rptsq.LE.rtstsq) GOTO 26                                       
   27 IF(COSGAM.GE.0.d0) GOTO 26                                        
      COSGAM=-COSGAM                                                    
      DARKEN=1.d0-X+X*COSGAM                                            
      if(ld.ne.2) goto 141                                              
      if(cosgam.eq.0.d0) goto 141                                       
      darken=darken-y*cosgam*dlog(cosgam)                               
      goto 147                                                          
  141 continue                                                          
      if(ld.eq.3) darken=darken-y*(1.d0-dsqrt(cosgam))                  
  147 if(darken.lt.0.d0) darken=0.d0                                    
      ATRATN=1.d0                                                       
      ATRATO=1.d0                                                       
      CORFAC=1.d0                                                       
      do 923 jn=1,nl                                                    
      Lspot(komp,jn)=0                                                  
  923 if(kks(komp,jn).eq.0) Lspot(komp,jn)=1                            
      IF(NSPOT.EQ.0) GOTO 640                                           
      CALL SPOT(KOMP,NSPOT,SINTH,COSTH,SINFI,COSFI,TEMF)                
      IF(TEMF.EQ.1.d0) GOTO 640                                         
      TSP=TLD(IS)*TEMF                                                  
      IF(IFAT.EQ.0) GOTO 941                                            
      CALL ATM(TLD(IS),WL,ATRATO)                                       
      CALL ATM(TSP,WL,ATRATN)                                           
  941 CORFAC=ATRATN*(dexp(XLUMP/TLD(IS))-1.d0)/(ATRATO*                 
     $(dexp(xlump/tsp)-1.d0))                                           
  640 CONTINUE                                                          
      rit=1.d0                                                          
      if(ncl.eq.0) goto 818                                             
      do 815 icl=1,ncl                                                  
      opsfcl=opsf*fcl(icl)                                              
      call cloud(xlos,ylos,zlos,xx,yy,zz,xcl(icl),ycl(icl),zcl(icl),    
     $rcl(icl),wl,op1(icl),opsfcl,edens(icl),acm,encl(icl),cmpd,        
     $ri,dx,dens(icl),tau)                                              
      rit=rit*ri                                                        
  815 continue                                                          
  818 continue                                                          
      DIF=rit*COSGAM*DARKEN*CORFAC*(CMP*SLUMP2(IX)+CMPP*SLUMP1(IX))     
      v=-r*(STCF*YLOS-stsf*XLOS)*COMP                                   
      if(ifphn.eq.1) goto 423                                           
      if(mpage.ne.3) goto 423                                           
      vflump=vfvuff*r*comp*costh                                        
      vcks=v*vfvuff                                                     
      vks=vcks+voc                                                      
      vksf=vks                                                          
      dvdr=vcks/r                                                       
      dvdth=vcks/tanth                                                  
      dvdfib=vfvuff*r*comp*(sinfi*ylos+cosfi*xlos)                      
c     dvdfic=dvdfib*sinth                                               
      difp=dif*deldel*renfsq                                            
c  dvdth and dvdfi (below) each need another term involving dr/d(theta) 
c    or dr/d(fi), that I will put in later. There will be a small loss  
c    of accuracy for distorted stars without those terms. See notes.    
      if(komp.eq.2) goto 422                                            
      binc=binc1                                                        
      binw=binw1                                                        
      do 1045 ifn=-nfm1,nfm1,2                                          
      dthf=dfloat(ifn)*r2nfdt                                           
      dvdfi=dvdfib*(sinth+costh*dthf)                                   
      do 1046 jfn=-nfm1,nfm1,2                                          
      if(nf.eq.1) goto 1047                                             
      dfif=dfloat(jfn)*r2nfdf*delfi                                     
      dvdth=-vflump*((cosfi-sinfi*dfif)*ylos-(sinfi+cosfi*dfif)*xlos)   
      dlr=0.d0                                                          
      vksf=vks+dvdr*dlr+dvdth*dthf+dvdfi*dfif                           
 1047 call linpro(komp,dvks1,hbarw1,tau1,count1,taug,fbin1,delv1)       
      if(inmin.lt.in1min) in1min=inmin                                  
      if(inmax.gt.in1max) in1max=inmax                                  
 1046 continue                                                          
 1045 continue                                                          
      goto 423                                                          
  422 continue                                                          
      binc=binc2                                                        
      binw=binw2                                                        
      do 1145 ifn=-nfm1,nfm1,2                                          
      dthf=dfloat(ifn)*r2nfdt                                           
      dvdfi=dvdfib*(sinth+costh*dthf)                                   
      do 1146 jfn=-nfm1,nfm1,2                                          
      if(nf.eq.1) goto 1147                                             
      dfif=dfloat(jfn)*r2nfdf*delfi                                     
      dvdth=-vflump*((cosfi-sinfi*dfif)*ylos-(sinfi+cosfi*dfif)*xlos)   
      dlr=0.d0                                                          
      vksf=vks+dvdr*dlr+dvdth*dthf+dvdfi*dfif                           
      ffi=dacos(cosfi)                                                  
      if(sinfi.lt.0.d0) ffi=twopi-ffi                                   
 1147 call linpro(komp,dvks2,hbarw2,tau2,count2,taug,fbin2,delv2)       
      if(inmin.lt.in2min) in2min=inmin                                  
      if(inmax.gt.in2max) in2max=inmax                                  
 1146 continue                                                          
 1145 continue                                                          
  423 continue                                                          
      DIFF=DIF*V                                                        
      SOMJ=SOMJ+DIFF                                                    
   45 SUMJ=SUMJ+DIF                                                     
   26 CONTINUE                                                          
      SOMJ=SOMJ*DELFI                                                   
      SUMJ=SUMJ*DELFI                                                   
      SOM=SOM+SOMJ                                                      
   36 SUM=SUM+SUMJ                                                      
      IF(SINSQ.GE.SINSQE) GOTO 75                                       
      L=L-1                                                             
      LK=k                                                              
      if(L.lt.14) LK=L/2-1                                              
      CALL fourls(theta,rho,L,LK,aa,bb)                                 
   75 IF(TEST.LE..0625d0) GOTO 118                                      
      SUMKUL=SUM*DELTH                                                  
      SOMKUL=SOM*DELTH                                                  
      X=XH                                                              
      y=yh                                                              
      KOMP=1                                                            
      nl=nl1                                                            
      ffc=ff1                                                           
      voc=vo1*vfvu                                                      
      NSPOT=NSP1                                                        
      IFAT=IFAT1                                                        
      EN=N1                                                             
      SAFTY=2.6d0*RV(LLL1)/EN                                           
      RMAX=RVQ(LLL2)+SAFTY                                              
      RMIN=RVQ(1)-SAFTY                                                 
      NPH=N1                                                            
      NP=2*N1                                                           
      nf=nf1                                                            
      GOTO 128                                                          
  118 X=XC                                                              
      y=yc                                                              
      KOMP=2                                                            
      nl=nl2                                                            
      ffc=ff2                                                           
      voc=vo2*vfvu                                                      
      NSPOT=NSP2                                                        
      IFAT=IFAT2                                                        
      SUMHOT=SUM*DELTH                                                  
      SOMHOT=SOM*DELTH                                                  
      if(inmax.gt.in1max) in1max=inmax                                  
      if(inmin.lt.in1min) in1min=inmin                                  
      EN=N2                                                             
      SAFTY=2.6d0*RVQ(LLL2)/EN                                          
      RMAX=RV(LLL1)+SAFTY                                               
      RMIN=RV(1)-SAFTY                                                  
      NPH=N2                                                            
      NP=2*N2                                                           
      nf=nf2                                                            
  128 DELTH=pih/EN                                                      
      enf=dfloat(nf)                                                    
      nfm1=nf-1                                                         
      renfsq=1.d0/(enf*enf)                                             
      r2nfdt=.5d0*delth/enf                                             
      vfvuff=vfvu*ffc                                                   
      SOM=0.d0                                                          
      SUM=0.d0                                                          
      NPP=NP+1                                                          
      TEMF=1.d0                                                         
      inmin=30000                                                       
      inmax=0                                                           
      DO 136 I=1,NP                                                     
      IF(I.GT.NPH) GOTO 154                                             
      UPDOWN=1.d0                                                       
      IK=I                                                              
      GOTO 155                                                          
  154 UPDOWN=-1.d0                                                      
      IK=NPP-I                                                          
  155 CONTINUE                                                          
      IPN1=IK+(KOMP-1)*N1                                               
      SINTH=SNTH(IPN1)                                                  
      COSTH=CSTH(IPN1)*UPDOWN                                           
      tanth=sinth/costh                                                 
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      XM=MM                                                             
      MH=MM                                                             
      MM=2*MM                                                           
      DELFI=pi/XM                                                       
      deldel=delth*delfi                                                
      SOMJ=0.d0                                                         
      SUMJ=0.d0                                                         
      SIGN=0.d0                                                         
      DRHO=1.d0                                                         
      MPP=MM+1                                                          
      DO 126 J=1,MM                                                     
      IF(J.GT.MH) GOTO 158                                              
      RTLEFT=1.d0                                                       
      JK=J                                                              
      GOTO 159                                                          
  158 RTLEFT=-1.d0                                                      
      JK=MPP-J                                                          
  159 CONTINUE                                                          
      IP=(KOMP-1)*NP1+IK                                                
      IX=MMSAVE(IP)+JK                                                  
      IS=IX+LLL1*(KOMP-1)                                               
      SINFI=SNFI(IS)*RTLEFT                                             
      COSFI=CSFI(IS)                                                    
      STSF=SINTH*SINFI                                                  
      STCF=SINTH*COSFI                                                  
      IF(TEST.LE..0625d0)GOTO 139                                       
      IF(RV(IX).EQ.-1.d0) GOTO 126                                      
      GX=GRX(IX)                                                        
      GY=RTLEFT*GRY(IX)                                                 
      GZ=UPDOWN*GRZ(IX)                                                 
      R=RV(IX)                                                          
      grmag=gmag1(ix)                                                   
      GOTO 149                                                          
  139 IF(RVQ(IX).EQ.-1.d0) GOTO 126                                     
      GX=GRXQ(IX)                                                       
      GY=RTLEFT*GRYQ(IX)                                                
      GZ=UPDOWN*GRZQ(IX)                                                
      R=RVQ(IX)                                                         
      grmag=gmag2(ix)                                                   
  149 COSGAM=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                            
      IF(COSGAM.LT.0.d0) GOTO 104                                       
      SIGN=0.d0                                                         
      OLSIGN=0.d0                                                       
      GOTO 126                                                          
  104 COSGAM=-COSGAM                                                    
      ZZ=R*COSTH                                                        
      YY=R*COMPP*STSF                                                   
      XX=CMPPD+COMPP*STCF*R                                             
      DARKEN=1.d0-X+X*COSGAM                                            
      if(ld.ne.2) goto 142                                              
      if(cosgam.eq.0.d0) goto 142                                       
      darken=darken-y*cosgam*dlog(cosgam)                               
      goto 148                                                          
  142 continue                                                          
      if(ld.eq.3) darken=darken-y*(1.d0-dsqrt(cosgam))                  
  148 if(darken.lt.0.d0) darken=0.d0                                    
      OLDIF=DIF                                                         
      ATRATN=1.d0                                                       
      ATRATO=1.d0                                                       
      CORFAC=1.d0                                                       
      do 823 jn=1,nl                                                    
      Lspot(komp,jn)=0                                                  
  823 if(kks(komp,jn).eq.0) Lspot(komp,jn)=1                            
      IF(NSPOT.EQ.0) GOTO 660                                           
      CALL SPOT(KOMP,NSPOT,SINTH,COSTH,SINFI,COSFI,TEMF)                
      IF(TEMF.EQ.1.d0) GOTO 660                                         
      TSP=TLD(IS)*TEMF                                                  
      IF(IFAT.EQ.0) GOTO 661                                            
      CALL ATM(TLD(IS),WL,ATRATO)                                       
      CALL ATM(TSP,WL,ATRATN)                                           
  661 CORFAC=ATRATN*(dexp(XLUMP/TLD(IS))-1.d0)/(ATRATO*                 
     $(dexp(xlump/tsp)-1.d0))                                           
  660 CONTINUE                                                          
      rit=1.d0                                                          
      if(ncl.eq.0) goto 718                                             
      do 715 icl=1,ncl                                                  
      opsfcl=opsf*fcl(icl)                                              
      call cloud(xlos,ylos,zlos,xx,yy,zz,xcl(icl),ycl(icl),zcl(icl),    
     $rcl(icl),wl,op1(icl),opsfcl,edens(icl),acm,encl(icl),cmppd,       
     $ri,dx,dens(icl),tau)                                              
      rit=rit*ri                                                        
  715 continue                                                          
  718 continue                                                          
      DIF=rit*COSGAM*DARKEN*CORFAC*(CMPP*SLUMP2(IX)+CMP*SLUMP1(IX))     
      v=R*(STCF*YLOS-STSF*XLOS)*COMP                                    
      DIFF=DIF*V                                                        
      IF(SINSQ.GT.SINSQE) GOTO 63                                       
      OLSIGN=SIGN                                                       
      OLDRHO=DRHO                                                       
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                 
      ZSKY=-XX*CICP+yy*CISP+ZZ*SINI+CMPD*CICP                           
      RRHO=dsqrt(ysky*ysky+zsky*zsky)                                   
      IF(RRHO.GT.RMAX)GOTO 63                                           
      IF(RRHO.LT.RMIN)GOTO 126                                          
      THET=dasin(ZSKY/RRHO)                                             
      IF(YSKY.LT.0.d0) GOTO 192                                         
      THET=twopi+THET                                                   
      GOTO 193                                                          
  192 THET=pi-THET                                                      
  193 IF(THET.GE.twopi) THET=THET-twopi                                 
      RHHO=0.d0                                                         
      DO 52 N=1,KK                                                      
      ENNN=N-1                                                          
      ENTHET=ENNN*THET                                                  
   52 RHHO=RHHO+AA(N)*dcos(ENTHET)+BB(N)*dsin(ENTHET)                   
      SIGN=1.d0                                                         
      IF(RRHO.LE.RHHO) sign=-1.d0                                       
      if(mpage.eq.3) goto 861                                           
      DRHO=dabs(RRHO-RHHO)                                              
      IF((SIGN*OLSIGN).GE.0.d0) GOTO 60                                 
      SUMDR=DRHO+OLDRHO                                                 
      FACT=-(.5d0-DRHO/SUMDR)                                           
      IF(FACT.LT.0.d0) GOTO 198                                         
      RDIF=OLDIF                                                        
      GOTO 199                                                          
  198 RDIF=DIF                                                          
  199 CORR=FACT*RDIF*SIGN                                               
      CORRR=CORR*V                                                      
      SUMJ=SUMJ+CORR                                                    
      SOMJ=SOMJ+CORRR                                                   
   60 IF(SIGN.LT.0.d0) GOTO 126                                         
   63 SUMJ=SUMJ+DIF                                                     
      SOMJ=SOMJ+DIFF                                                    
      if(mpage.ne.5) goto 127                                           
      ipc=ipc+1                                                         
      yskp(ipc)=(xx-qfacd)*sinph+yy*cosph                               
      zskp(ipc)=(-xx+qfacd)*cicp+yy*cisp+zz*sini                        
      if(nspot.eq.0) goto 126                                           
      call spot(komp,nspot,sinth,costh,sinfi,cosfi,temf)                
      if(temf.eq.1.d0) goto 126                                         
      yskr=yskp(ipc)                                                    
      zskr=zskp(ipc)                                                    
      stp=twopi/dfloat(kstp)                                            
      do 189 ang=stp,twopi,stp                                          
      ipc=ipc+1                                                         
      yskp(ipc)=yskr+dsin(ang)*cirf                                     
      zskp(ipc)=zskr+dcos(ang)*cirf                                     
  189 continue                                                          
      goto 126                                                          
  127 continue                                                          
      if(mpage.ne.3) goto 126                                           
      if(ifphn.eq.1) goto 126                                           
  861 vflump=vfvuff*r*comp*costh                                        
      vcks=v*vfvuff                                                     
      vks=vcks+voc                                                      
      vksf=vks                                                          
      dvdr=vcks/r                                                       
      dvdth=vcks/tanth                                                  
      dvdfib=vfvuff*r*comp*(sinfi*ylos+cosfi*xlos)                      
      difp=dif*deldel*renfsq                                            
      if(komp.eq.2) goto 452                                            
      binc=binc1                                                        
      binw=binw1                                                        
      do 1245 ifn=-nfm1,nfm1,2                                          
      dthf=dfloat(ifn)*r2nfdt                                           
      snthl=costh*dthf                                                  
      zz=r*(costh-sinth*dthf)                                           
      dvdfi=dvdfib*(sinth+costh*dthf)                                   
      do 1246 jfn=-nfm1,nfm1,2                                          
      if(nf.eq.1) goto 1247                                             
      dfif=dfloat(jfn)*r2nfdf*delfi                                     
      dlr=0.d0                                                          
      xx=cmppd+compp*r*snthl*(cosfi-sinfi*dfif)                         
      yy=r*compp*snthl*(sinfi+cosfi*dfif)                               
      ysky=(xx-cmpd)*sinph+yy*cosph                                     
      zsky=(cmpd-xx)*cicp+yy*cisp+zz*sini                               
      rrho=dsqrt(ysky*ysky+zsky*zsky)                                   
      if(rrho.lt.rhho) goto 1246                                        
      dvdth=-vflump*((cosfi-sinfi*dfif)*ylos-(sinfi+cosfi*dfif)*xlos)   
      vksf=vks+dvdr*dlr+dvdth*dthf+dvdfi*dfif                           
 1247 call linpro(komp,dvks1,hbarw1,tau1,count1,taug,fbin1,delv1)       
      if(inmax.gt.in1max) in1max=inmax                                  
      if(inmin.lt.in1min) in1min=inmin                                  
 1246 continue                                                          
 1245 continue                                                          
      goto 126                                                          
  452 continue                                                          
      binc=binc2                                                        
      binw=binw2                                                        
      do 1445 ifn=-nfm1,nfm1,2                                          
      dthf=dfloat(ifn)*r2nfdt                                           
      snthl=costh*dthf                                                  
      zz=r*(costh-sinth*dthf)                                           
      dvdfi=dvdfib*(sinth+costh*dthf)                                   
      do 1446 jfn=-nfm1,nfm1,2                                          
      if(nf.eq.1) goto 1447                                             
      dfif=dfloat(jfn)*r2nfdf*delfi                                     
      dvdth=-vflump*((cosfi-sinfi*dfif)*ylos-(sinfi+cosfi*dfif)*xlos)   
      dlr=0.d0                                                          
      xx=cmppd+compp*r*snthl*(cosfi-sinfi*dfif)                         
      yy=r*compp*snthl*(sinfi+cosfi*dfif)                               
      ysky=(xx-cmpd)*sinph+yy*cosph                                     
      zsky=(cmpd-xx)*cicp+yy*cisp+zz*sini                               
      rrho=dsqrt(ysky*ysky+zsky*zsky)                                   
      if(rrho.lt.rhho) goto 1446                                        
      vksf=vks+dvdr*dlr+dvdth*dthf+dvdfi*dfif                           
 1447 call linpro(komp,dvks2,hbarw2,tau2,count2,taug,fbin2,delv2)       
      if(inmax.gt.in2max) in2max=inmax                                  
      if(inmin.lt.in2min) in2min=inmin                                  
 1446 continue                                                          
 1445 continue                                                          
  126 CONTINUE                                                          
      SOMJ=SOMJ*DELFI                                                   
      SUMJ=SUMJ*DELFI                                                   
      SOM=SOM+SOMJ                                                      
  136 SUM=SUM+SUMJ                                                      
      if(mpage.eq.5) return                                             
      IF(TEST.LE..0625d0) GOTO 120                                      
      SOMHOT=SOM*DELTH                                                  
      SUMHOT=SUM*DELTH                                                  
      GOTO 121                                                          
  120 SUMKUL=SUM*DELTH                                                  
      SOMKUL=SOM*DELTH                                                  
  121 continue                                                          
      if(ifphn.eq.1) return                                             
      if(mpage.ne.3) return                                             
      in1min=in1min-marm1                                               
      in1max=in1max+marp1                                               
      in2min=in2min-marm2                                               
      in2max=in2max+marp2                                               
      if(nl1.eq.0) goto 3115                                            
      do 2912 i=in1min,in1max                                           
      fbin1(i)=1.d0-fbin1(i)/sumhot                                       
      if(count1(i).eq.0.d0) goto 2918                                   
      delv1(i)=delv1(i)/count1(i)                                       
      goto 2919                                                         
 2918 delv1(i)=binw1*(dfloat(i)-binc1)                                  
 2919 vdc=delv1(i)/clfac                                                
      vfac=dsqrt((1.d0+vdc)/(1.d0-vdc))                                 
      delwl1(i)=wl*(vfac-1.d0)                                          
      wl1(i)=wl*vfac                                                    
      resf1(i)=(sl1*delwl1(i)+sc1)*fbin1(i)                             
 2912 continue                                                          
 3115 if(nl2.eq.0) return                                               
      do 2914 i=in2min,in2max                                           
      fbin2(i)=1.d0-fbin2(i)/sumkul                                       
      if(count2(i).eq.0.d0) goto 2917                                   
      delv2(i)=delv2(i)/count2(i)                                       
      goto 2920                                                         
 2917 delv2(i)=binw2*(dfloat(i)-binc2)                                  
 2920 vdc=delv2(i)/clfac                                                
      vfac=dsqrt((1.d0+vdc)/(1.d0-vdc))                                 
      delwl2(i)=wl*(vfac-1.d0)                                          
      wl2(i)=wl*vfac                                                    
      resf2(i)=(sl2*delwl2(i)+sc2)*fbin2(i)                             
 2914 continue                                                          
      return                                                            
      END                                                               
      SUBROUTINE SINCOS (KOMP,N,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)          
c  Version of November 9, 1995                                          
      implicit real*8 (a-h,o-z)                                         
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*),MMSAVE(*)               
      IP=(KOMP-1)*(N1+1)+1                                              
      IQ=IP-1                                                           
      IS=0                                                              
      IF(KOMP.EQ.2) IS=MMSAVE(IQ)                                       
      MMSAVE(IP)=0                                                      
      EN=N                                                              
      DO 8 I=1,N                                                        
      EYE=I                                                             
      EYE=EYE-.5d0                                                      
      TH=1.570796326794897d0*EYE/EN                                     
      IPN1=I+N1*(KOMP-1)                                                
      SNTH(IPN1)=dsin(TH)                                               
      CSTH(IPN1)=dcos(TH)                                               
      EM=SNTH(IPN1)*EN*1.3d0                                            
      MM=EM+1.d0                                                        
      XM=MM                                                             
      IP=(KOMP-1)*(N1+1)+I+1                                            
      IQ=IP-1                                                           
      MMSAVE(IP)=MMSAVE(IQ)+MM                                          
      DO 8 J=1,MM                                                       
      IS=IS+1                                                           
      XJ=J                                                              
      FI=3.141592653589793d0*(XJ-.5d0)/XM                               
      CSFI(IS)=dcos(FI)                                                 
      SNFI(IS)=dsin(FI)                                                 
    8 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE SURFAS(RMASS,POTENT,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,      
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,FF,D,SNTH,CSTH,SNFI,CSFI,GRV1,  
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,     
     $GMAG2,GREXP)                                                      
c  Version of October 5, 1995                                           
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*
     $),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*)  
     $,GRV1(*),GRV2(*),XX1(*),YY1(*),ZZ1(*),XX2(*),YY2(*),ZZ2(*),GLUMP1 
     $(*),GLUMP2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                 
      common /radi/ R1H,RLH,R1C,RLC                                     
      COMMON X1                                                         
      COMMON /ECCEN/e,adum,perdum,vgadum,sindum,vfdum,vfadum,vgmdum,    
     $v1dum,v2dum,ifcdum                                                
      DSQ=D*D                                                           
      RMAS=RMASS                                                        
      IF(KOMP.EQ.2) RMAS=1.d0/RMASS                                     
      RF=FF**2                                                          
      RTEST=0.d0                                                        
      IP=(KOMP-1)*(N1+1)+1                                              
      IQ=IP-1                                                           
      IS=0                                                              
      ISX=(KOMP-1)*MMSAVE(IQ)                                           
      MMSAVE(IP)=0                                                      
      KFLAG=0                                                           
      CALL ELLONE (FF,D,RMAS,X1,OMEGA,XL2,OM2)                          
      IF(KOMP.EQ.2) OMEGA=RMASS*OMEGA+.5d0*(1.d0-RMASS)                 
      X2=X1                                                             
      IF(KOMP.EQ.2) X1=1.d0-X1                                          
      IF(E.NE.0.d0) GOTO 43                                             
      IF(POTENT.LT.OMEGA) CALL NEKMIN(RMASS,POTENT,X1,ZZ)               
      IF(POTENT.LT.OMEGA) X2=1.d0-X1                                    
   43 COMP=dfloat(3-2*KOMP)                                             
      CMP=dfloat(KOMP-1)                                                
      CMPD=CMP*D                                                        
      TESTER=CMPD+COMP*X1                                               
      RM1=RMASS+1.d0                                                    
      RMS=RMASS                                                         
      RM1S=RM1                                                          
      IF(KOMP.NE.2) GOTO 15                                             
      POT=POTENT/RMASS+.5d0*(RMASS-1.d0)/RMASS                          
      RM=1.d0/RMASS                                                     
      RM1=RM+1.d0                                                       
      GOTO 20                                                           
   15 POT=POTENT                                                        
      RM=RMASS                                                          
   20 EN=N                                                              
      RSAVE=1.d0/POT                                                    
      DO 8 I=1,N                                                        
      IF(I.NE.2) GOTO 82                                                
      IF(KOMP.EQ.1) RTEST=.3d0*RV(1)                                    
      IF(KOMP.EQ.2) RTEST=.3d0*RVQ(1)                                   
   82 CONTINUE                                                          
      IPN1=I+N1*(KOMP-1)                                                
      SINTH=SNTH(IPN1)                                                  
      XNU=CSTH(IPN1)                                                    
      XNUSQ=XNU**2                                                      
      EM=SINTH*EN*1.3d0                                                 
      XLUMP=1.d0-XNUSQ                                                  
      MM=EM+1.d0                                                        
      DO 8 J=1,MM                                                       
      KOUNT=0                                                           
      IS=IS+1                                                           
      ISX=ISX+1                                                         
      DELR=0.d0                                                         
      COSFI=CSFI(ISX)                                                   
      XMU=SNFI(ISX)*SINTH                                               
      XLAM=SINTH*COSFI                                                  
      R=RSAVE                                                           
   14 R=R+DELR                                                          
      KOUNT=KOUNT+1                                                     
      IF(KOUNT.LT.20) GOTO 70                                           
      KFLAG=1                                                           
      R=-1.d0                                                           
      GOTO 86                                                           
   70 RSQ=R*R                                                           
      PAR=DSQ-2.d0*XLAM*R*D+RSQ                                         
      RPAR=dsqrt(PAR)                                                   
      OM=1.d0/R+RM*((1.d0/RPAR)-XLAM*R/DSQ)+RM1*.5d0*RSQ*XLUMP*RF       
      DOMR=1.d0/(RF*RM1*XLUMP*R-1.d0/RSQ-(RM*(R-XLAM*D))/(PAR*RPAR)     
     $-rm*xlam/dsq)                                                     
      DELR=(POT-OM)*DOMR                                                
      ABDELR=dabs(DELR)                                                 
      IF(ABDELR.GT..00001d0) GOTO 14                                    
      ABR=dabs(R)                                                       
      IF(R.GT.RTEST) GOTO 74                                            
      KFLAG=1                                                           
      R=-1.d0                                                           
      IF(KOMP.EQ.2) GOTO 98                                             
      GOTO 97                                                           
   74 IF(ABR.LT.TESTER) RSAVE=R                                         
      Z=R*XNU                                                           
      Y=COMP*R*XMU                                                      
      X2T=ABR*XLAM                                                      
      X=CMPD+COMP*X2T                                                   
      IF(KOMP.EQ.2) GOTO 62                                             
      IF(X.LT.X1) GOTO 65                                               
      KFLAG=1                                                           
      R=-1.d0                                                           
      GOTO 97                                                           
   62 IF(X2T.LT.X2) GOTO 65                                             
      KFLAG=1                                                           
      R=-1.d0                                                           
      GOTO 98                                                           
   65 SUMSQ=Y**2+Z**2                                                   
      PAR1=X**2+SUMSQ                                                   
      RPAR1=dsqrt(PAR1)                                                 
      XNUM1=1.d0/(PAR1*RPAR1)                                           
      XL=D-X                                                            
      PAR2=XL**2+SUMSQ                                                  
      RPAR2=dsqrt(PAR2)                                                 
      XNUM2=1.d0/(PAR2*RPAR2)                                           
      OMZ=-Z*(XNUM1+RMS*XNUM2)                                          
      OMY=Y*(RM1S*RF-XNUM1-RMS*XNUM2)                                   
      OMX=RMS*XL*XNUM2-X*XNUM1+RM1S*X*RF-RMS/DSQ                        
      IF(KOMP.EQ.2) OMX=RMS*XL*XNUM2-X*XNUM1-RM1S*XL*RF+1.d0/DSQ        
      GRMAG=dsqrt(OMX*OMX+OMY*OMY+OMZ*OMZ)                              
      IF(IS.EQ.1) GRPOLE=GRMAG                                          
      GRAV=(GRMAG/GRPOLE)**GREXP                                        
      A=COMP*XLAM*OMX                                                   
      B=COMP*XMU*OMY                                                    
      C=XNU*OMZ                                                         
      COSBET=-(A+B+C)/GRMAG                                             
      IF(COSBET.LT..7d0) COSBET=.7d0                                    
   86 IF(KOMP.EQ.2) GOTO 98                                             
   97 RV(IS)=R                                                          
      GRX(IS)=OMX                                                       
      GRY(IS)=OMY                                                       
      GRZ(IS)=OMZ                                                       
      GMAG1(IS)=dsqrt(OMX*OMX+OMY*OMY+OMZ*OMZ)                          
      FR1(IS)=1.d0                                                      
      GLUMP1(IS)=R*R*SINTH/COSBET                                       
      GRV1(IS)=GRAV                                                     
      XX1(IS)=X                                                         
      YY1(IS)=Y                                                         
      ZZ1(IS)=Z                                                         
      CSBT1(IS)=COSBET                                                  
      GOTO 8                                                            
   98 RVQ(IS)=R                                                         
      GRXQ(IS)=OMX                                                      
      GRYQ(IS)=OMY                                                      
      GRZQ(IS)=OMZ                                                      
      GMAG2(IS)=dsqrt(OMX*OMX+OMY*OMY+OMZ*OMZ)                          
      FR2(IS)=1.d0                                                      
      GLUMP2(IS)=R*R*SINTH/COSBET                                       
      GRV2(IS)=GRAV                                                     
      XX2(IS)=X                                                         
      YY2(IS)=Y                                                         
      ZZ2(IS)=Z                                                         
      CSBT2(IS)=COSBET                                                  
    8 CONTINUE                                                          
      IF(KFLAG.EQ.0) GOTO 53                                            
      ISS=IS-1                                                          
      IF(KOMP.NE.1) GOTO 50                                             
      CALL RING(RMASS,POTENT,1,N,FR1,HLD,R1H,RLH)                       
      DO 55 I=1,ISS                                                     
      IPL=I+1                                                           
      IF(RV(I).GE.0.d0)GOTO 55                                          
      FR1(IPL)=FR1(IPL)+FR1(I)                                          
      FR1(I)=0.d0                                                       
   55 CONTINUE                                                          
   53 IF(KOMP.EQ.2) GOTO 54                                             
      IS=0                                                              
      DO 208 I=1,N                                                      
      IPN1=I+N1*(KOMP-1)                                                
      EM=SNTH(IPN1)*EN*1.3d0                                            
      MM=EM+1.d0                                                        
      DO 208 J=1,MM                                                     
      IS=IS+1                                                           
      GLUMP1(IS)=FR1(IS)*GLUMP1(IS)                                     
  208 CONTINUE                                                          
      RETURN                                                            
   50 CALL RING(RMASS,POTENT,2,N,FR2,HLD,R1C,RLC)                       
      DO 56 I=1,IS                                                      
      IPL=I+1                                                           
      IF(RVQ(I).GE.0.d0) GOTO 56                                        
      FR2(IPL)=FR2(IPL)+FR2(I)                                          
      FR2(I)=0.d0                                                       
   56 CONTINUE                                                          
   54 CONTINUE                                                          
      IS=0                                                              
      DO 108 I=1,N                                                      
      IPN1=I+N1*(KOMP-1)                                                
      EM=SNTH(IPN1)*EN*1.3d0                                            
      MM=EM+1.d0                                                        
      DO 108 J=1,MM                                                     
      IS=IS+1                                                           
      GLUMP2(IS)=FR2(IS)*GLUMP2(IS)                                     
  108 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,  
     $HLD,SLUMP1,SLUMP2,THETA,RHO,AA,BB,PHSV,PCSV,N1,N2,F1,F2,d,hlum,   
     $clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1   
     $,ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phas,rm,       
     $xincl,hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,     
     $yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,       
     $csbt2,gmag1,gmag2,fbin1,fbin2,delv1,delv2,count1,count2,          
     $delwl1,delwl2,resf1,resf2,wl1,wl2,dvks1,dvks2,tau1,tau2,hbarw1,   
     $hbarw2,xcl,ycl,zcl,rcl,op1,fcl,dens,encl,edens,taug,yskp,         
     $zskp,mode,ifphn)                                                  
c  Version of September 11, 1998                                        
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),      
     $GRZQ(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),SLUMP1(*),SLUMP2(*),       
     $THETA(*),RHO(*),AA(*),BB(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*),TLD(*)
     $,GLUMP1(*),GLUMP2(*),XX1(*),XX2(*),YY1(*),YY2(*),ZZ1(*),ZZ2(*),   
     $GRV1(*),GRV2(*),RFTEMP(*),RF1(*),RF2(*),CSBT1(*),CSBT2(*)         
     $,GMAG1(*),GMAG2(*)                                                
      dimension fbin1(*),fbin2(*),delv1(*),delv2(*),count1(*),count2(*),
     $delwl1(*),delwl2(*),resf1(*),resf2(*),wl1(*),wl2(*),dvks1(*),     
     $dvks2(*),tau1(*),tau2(*),hbarw1(*),hbarw2(*),taug(*)              
      dimension xcl(*),ycl(*),zcl(*),rcl(*),op1(*),fcl(*),dens(*),      
     $edens(*),encl(*),yskp(*),zskp(*)                                  
      COMMON /INVAR/ KH,IPBDUM,IRTE,NREF,IRVOL1,irvol2,mref,ifsmv1,     
     $ifsmv2,icor1,icor2,ld,ncl,jdphs,ipc                               
      COMMON /FLVAR/ PSHIFT,DP,DS,EF,EFC,ECOS,perr0,PHPER,PCONJ,        
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc,qfacd       
      common /nspt/ nsp1,nsp2,ifdum1,ifdum2                             
      common /ardot/ dperdt,hjd,hjd0,perr                               
      common /spots/ snlat(2,100),cslat(2,100),snlng(2,100),            
     $cslng(2,100),rdsp(2,100),tmsp(2,100),xlng(2,100),kks(2,100),      
     $Lspot(2,100)                                                      
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC     
      common /prof2/ vo1,vo2,ff1,ff2,du1,du2,du3,du4,du5,du6,du7        
      pi=3.141592653589793d0                                            
      twopi=pi+pi                                                       
      ff1=f1                                                            
      ff2=f2                                                            
      qfac1=1.d0/(1.d0+rm)                                              
      qfac=rm*qfac1                                                     
      MOD=(MODE-2)**2                                                   
      IF(MOD.EQ.1) XC=XH                                                
      if(mod.eq.1) yc=yh                                                
      PSFT=PHAS-PHPERI                                                  
   29 if(PSFT.GT.1.d0) PSFT=PSFT-1.d0                                   
      if(psft.gt.1.d0) goto 29                                          
   30 if(PSFT.LT.0.d0) PSFT=PSFT+1.d0                                   
      if(psft.lt.0.d0) goto 30                                          
      XMEAN=PSFT*twopi                                                  
      tr=xmean                                                          
      do 60 kp=1,2                                                      
      nsp=nsp1*(2-kp)+nsp2*(kp-1)                                       
      ff=f1*dfloat(2-kp)+f2*dfloat(kp-1)                                
      ifsmv=ifsmv1*(2-kp)+ifsmv2*(kp-1)                                 
      if(ifsmv.eq.0) goto 60                                            
      do 61 i=1,nsp                                                     
      xlg=xlng(kp,i)+twopi*ff*(phas-pconj)-(tr-trc)                     
      snlng(kp,i)=dsin(xlg)                                             
      cslng(kp,i)=dcos(xlg)                                             
   61 continue                                                          
   60 continue                                                          
      if(e.ne.0.d0) call KEPLER(XMEAN,E,DUM,TR)                         
      U=TR+PERR                                                         
      COSU=dcos(U)                                                      
      GPHA=U*.1591549d0-.25d0                                           
   40 if(GPHA.lt.0.d0) GPHA=GPHA+1.d0                                   
      if(gpha.lt.0.d0) goto 40                                          
   50 if(GPHA.GE.1.d0) GPHA=GPHA-1.d0                                   
      if(gpha.ge.1.d0) goto 50                                          
      D=EF/(1.d0+E*dcos(TR))                                            
      qfacd=qfac*d                                                      
      IF(IRTE.EQ.1) GOTO 19                                             
      CALL LCR(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SLUM
     $P1,SLUMP2,RM,PHSV,PCSV,N1,N2,F1,F2,D,HLUM,CLUM,xh,xc,yh,yc,gr1,gr2
     $,wl,SM1,SM2,TPOLH,TPOLC,SBRH,SBRC,IFAT1,IFAT2,TAVH,TAVC,alb1,alb2,
     $xbol1,xbol2,ybol1,ybol2,vol1,vol2,snth,csth,snfi,csfi,tld,glump1, 
     $glump2,xx1,xx2,yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,csbt1,csbt2, 
     $rftemp,rf1,rf2,gmag1,gmag2,mode,ifc)                              
   19 CONTINUE                                                          
      VO1=qfac*SINI*(ECOS+COSU)/EFC+VGAM                                
      VO2=-qfac1*SINI*(ECOS+COSU)/EFC+VGAM                              
      call light(gpha,xincl,xh,xc,yh,yc,n1,n2,hot,cool,rv,grx,gry,grz,  
     $rvq,grxq,gryq,grzq,mmsave,theta,rho,aa,bb,slump1,slump2,somhot,   
     $somkul,d,wl,snth,csth,snfi,csfi,tld,gmag1,gmag2,fbin1,fbin2,      
     $delv1,delv2,count1,count2,delwl1,delwl2,resf1,resf2,wl1,wl2,dvks1,
     $dvks2,tau1,tau2,hbarw1,hbarw2,xcl,ycl,zcl,rcl,op1,fcl,edens,      
     $encl,dens,taug,yskp,zskp,ifphn)                                   
      VRA1=0.d0                                                         
      VRA2=0.d0                                                         
      IF(HOT.GT.0.d0) VRA1=F1*SOMHOT/HOT                                
      IF(COOL.GT.0.d0) VRA2=F2*SOMKUL/COOL                              
      vsum1=vo1                                                         
      vsum2=vo2                                                         
      if(icor1.eq.1) vsum1=vo1+vra1                                     
      if(icor2.eq.1) vsum2=vo2+vra2                                     
      VKM1=VSUM1*VFVU                                                   
      VKM2=VSUM2*VFVU                                                   
      RETURN                                                            
      END                                                               
      SUBROUTINE DURA(F,XINCL,RM,D,THE,OMEG,R)                          
c  Version of May 19, 1996                                              
C                                                                       
C     PARAMETER 'THE' IS THE SEMI-DURATION OF X-RAY ECLIPSE, AND SHOULD 
C     BE IN CIRCULAR MEASURE.                                           
      IMPLICIT REAL*8(A-H,O-Z)                                          
      DELX=0.D0                                                         
      FSQ=F*F                                                           
      RMD=1.d0/RM                                                       
      RMD1=RMD+1.D0                                                     
      XINC=.017453293d0*XINCL                                           
      TH=6.2831853071795865d0*THE                                       
      CI=DCOS(XINC)                                                     
      SI=DSIN(XINC)                                                     
      DSQ=D*D                                                           
      ST=DSIN(TH)                                                       
      CT=DCOS(TH)                                                       
      COTI=CI/SI                                                        
      TT=ST/CT                                                          
      C1=CT*SI                                                          
      C2=TT*ST*SI                                                       
      C3=C1+C2                                                          
      C4=COTI*CI/CT                                                     
      C5=C3+C4                                                          
      C6=C2+C4                                                          
      C7=(ST*ST+COTI*COTI)/CT**2                                        
      X=D*(SI*SI*ST*ST+CI*CI)+.00001D0                                  
   15 X=X+DELX                                                          
      PAR=X*X+C7*(D-X)**2                                               
      RPAR=DSQRT(PAR)                                                   
      PAR32=PAR*RPAR                                                    
      PAR52=PAR*PAR32                                                   
      FC=(C6*D-C5*X)/PAR32+C1**3*C5*RMD/(D-X)**2+C3*FSQ*RMD1*X-C2*FSQ*D*
     $RMD1-C1*RMD/DSQ                                                   
      DFCDX=(-C5*PAR-3.D0*(C6*D-C5*X)*((1.D0+C7)*X-C7*D))/PAR52+2.D0*C1 
     $**3*C5*RMD/(D-X)**3+C3*FSQ*RMD1                                   
      DELX=-FC/DFCDX                                                    
      ABDELX=DABS(DELX)                                                 
      IF(ABDELX.GT.1.d-7) GOTO 15                                       
      Y=-(D-X)*TT                                                       
      Z=-(D-X)*COTI/CT                                                  
      YZ2=Y*Y+Z*Z                                                       
      OMEG=1.D0/DSQRT(X*X+YZ2)+RMD/DSQRT((D-X)**2+YZ2)+.5D0*RMD1*FSQ*   
     $(X*X+Y*Y)-RMD*X/DSQ                                               
      OMEG=RM*OMEG+.5d0*(1.d0-RM)                                       
      R=DSQRT(X*X+YZ2)                                                  
      RETURN                                                            
      END                                                               
      SUBROUTINE LCR(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HL
     $D,SLUMP1,SLUMP2,RM,POTH,POTC,N1,N2,F1,F2,D,HLUM,CLUM,xh,xc,yh,yc, 
     $GR1,GR2,WL,SM1,SM2,TPOLH,TPOLC,SBRH,SBRC,IFAT1,IFAT2,TAVH,TAVC,   
     $alb1,alb2,xbol1,xbol2,ybol1,ybol2,vol1,vol2,snth,csth,snfi,csfi,  
     $tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,  
     $csbt1,csbt2,rftemp,rf1,rf2,gmag1,gmag2,mode,ifc)                  
c  Version of September 14, 1998                                        
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),     
     $CSTH(*),SNFI(*),CSFI(*),TLD(*),GLUMP1(*),GLUMP2(*),XX1(*),XX2(*)  
     $,YY1(*),YY2(*),ZZ1(*),ZZ2(*),GRV1(*),GRV2(*),RFTEMP(*),RF1(*),    
     $RF2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                        
      COMMON /DPDX/ DPDX1,DPDX2,PHSV,PCSV                               
      COMMON /ECCEN/ E,dum1,dum2,dum3,dum4,dum5,dum6,dum7,dum8,dum9,idum
      COMMON /SUMM/ SUMM1,SUMM2                                         
      COMMON /INVAR/ KHDUM,IPB,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,     
     $ifsmv2,icor1,icor2,ld,ncl,jdphs,ipc                               
      XLUMP=1.4384d0/WL                                                 
      VL1=VOL1                                                          
      VL2=VOL2                                                          
      DP=1.d0-E                                                         
      IF(IRVOL1.EQ.1) GOTO 88                                           
      CALL VOLUME(VL1,RM,POTH,DP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,GRXQ,    
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM1,SM1,GRV1,  
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2
     $,GR1,1)                                                           
      IF(E.EQ.0.d0) GOTO 88                                             
      POTHD=PHSV                                                        
      IF(IFC.EQ.2) POTHD=PHSV+DPDX1*(1.d0/D-1.d0/(1.d0-E))              
      CALL VOLUME(VL1,RM,POTHD,D,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,GRXQ,    
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM1,SM1,GRV1,  
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2
     $,GR1,IFC)                                                         
   88 CONTINUE                                                          
      IF(IRVOL2.EQ.1) GOTO 86                                           
      CALL VOLUME(VL2,RM,POTC,DP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ,    
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM2,SM2,GRV1,  
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2
     $,GR2,1)                                                           
      IF(E.EQ.0.d0) GOTO 86                                             
      POTCD=PCSV                                                        
      IF(IFC.EQ.2) POTCD=PCSV+DPDX2*(1.d0/D-1.d0/(1.d0-E))              
      CALL VOLUME(VL2,RM,POTCD,D,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ,    
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM2,SM2,GRV1,  
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2
     $,GR2,IFC)                                                         
   86 CONTINUE                                                          
      TPOLH=TAVH*dsqrt(dsqrt(SM1/SUMM1))                                
      TPOLC=TAVC*dsqrt(dsqrt(SM2/SUMM2))                                
      g1=gmag1(1)                                                       
      g2=gmag2(1)                                                       
      IF(MODE.EQ.1)TPOLC=TPOLH*dsqrt(dsqrt((G2/G1)**GR1))               
      IF(MODE.EQ.1)TAVC=TPOLC/dsqrt(dsqrt(SM2/SUMM2))                   
      RATH=1.d0                                                         
      RATC=1.d0                                                         
      tph=10000.d0*tpolh                                                
      tpc=10000.d0*tpolc                                                
      IF(IFAT1.NE.0) CALL ATM(tph,WL,RATH)                              
      IF(IFAT2.NE.0) CALL ATM(tpc,WL,RATC)                              
      RATCH=RATC/RATH                                                   
      call lum(hlum,xh,yh,wl,tpolh,n1,n1,1,sbrh,rv,rvq,glump1,          
     $glump2,grv1,grv2,mmsave,summ1d,fr1,sm1d,ifat1,vold,rm,poth,       
     $f1,d,snth)                                                        
      dhfac=9.d0-3.d0*xh                                                
      dcfac=9.d0-3.d0*xc                                                
      if(ld.eq.2) dhfac=dhfac+2.d0*yh                                   
      if(ld.eq.2) dcfac=dcfac+2.d0*yc                                   
      if(ld.eq.3) dhfac=dhfac-1.8d0*yh                                  
      if(ld.eq.3) dcfac=dcfac-1.8d0*yc                                  
      sbrc=ratch*sbrh*dhfac*(dexp(xlump/tpolh)-1.d0)/(dcfac*            
     $(dexp(xlump/tpolc)-1.d0))                                         
      call lum(clum,xc,yc,wl,tpolc,n2,n1,2,sbrt,rv,rvq,glump1,          
     $glump2,grv1,grv2,mmsave,summ2d,fr2,sm2d,ifat2,vold,rm,potc,       
     $f2,d,snth)                                                        
      IF(IPB.EQ.1) SBRC=SBRT                                            
      IF(MODE.GT.0)CLUM=CLUM*SBRC/SBRT                                  
      IF(MODE.LE.0)SBRC=SBRT                                            
      if(mref.eq.2) goto 30                                             
      ratlum=hlum/clum                                                  
      call bolo(ratlum,tavh,tavc,wl,ifat1,ifat2,ratbol)                 
      rb=1.d0/ratbol                                                    
      call olump(rv,grx,gry,grz,rvq,grxq,gryq,grzq,slump1,slump2,mmsave 
     $,gr1,alb1,rb,wl,tpolh,sbrh,summ1,n1,n2,1,ifat1,xc,yc,d,snth       
     $,csth,snfi,csfi,tld,glump1,glump2)                                
      rb=ratbol                                                         
      call olump(rv,grx,gry,grz,rvq,grxq,gryq,grzq,slump1,slump2,mmsave 
     $,gr2,alb2,rb,wl,tpolc,sbrc,summ2,n1,n2,2,ifat2,xh,yh,d,snth       
     $,csth,snfi,csfi,tld,glump1,glump2)                                
      return                                                            
   30 continue                                                          
      sbr1b=tpolh**4/dint1                                              
      sbr2b=tpolc**4/dint2                                              
      LT=N1+1                                                           
      IMAX1=MMSAVE(LT)                                                  
      DO 80 I=1,IMAX1                                                   
      RFTEMP(I)=1.d0                                                    
   80 RF1(I)=1.d0                                                       
      LT=N1+N2+2                                                        
      IMAX2=MMSAVE(LT)                                                  
      DO 81 I=1,IMAX2                                                   
   81 RF2(I)=1.d0                                                       
      DO 93 NR=1,NREF                                                   
      CALL LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,MMSAVE,        
     $alb1,wl,tpolh,sbrh,n1,n2,1,ifat1,fr1,snth,                        
     $tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,xbol2,ybol2,grv1,       
     $grv2,sbr1b,sbr2b,rftemp,rf2,gmag1,gmag2,dint1)                    
      CALL LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,MMSAVE,        
     $ALB2,WL,TPOLC,SBRC,N1,N2,2,IFAT2,fr2,snth,                        
     $tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,xbol1,ybol1,            
     $grv1,grv2,sbr1b,sbr2b,rf2,rf1,gmag1,gmag2,dint2)                  
      DO 70 I=1,IMAX1                                                   
   70 RF1(I)=RFTEMP(I)                                                  
   93 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE VOLUME(V,Q,P,D,FF,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,        
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM,SM,    
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1 
     $,GMAG2,GREXP,IFC)                                                 
c  Version of September 14, 1998                                        
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*
     $),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*)  
     $,GRV1(*),GRV2(*),GLUMP1(*),GLUMP2(*),XX1(*),YY1(*),ZZ1(*),XX2(*), 
     $YY2(*),ZZ2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                 
      DP=1.d-3*P                                                        
      IF (IFC.EQ.1) DP=0.d0                                             
      TOL=1.d-6*P**2                                                    
      DELP=0.d0                                                         
      KNTR=0                                                            
   16 P=P+DELP                                                          
      KNTR=KNTR+1                                                       
      IF(KNTR.GE.20) TOL=TOL+TOL                                        
      PS=P                                                              
      DO 17 I=1,IFC                                                     
      P=PS                                                              
      IF(I.EQ.1) P=P+DP                                                 
      CALL SURFAS(Q,P,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,      
     $MMSAVE,FR1,FR2,HLD,FF,D,SNTH,CSTH,SNFI,CSFI,GRV1,GRV2,XX1,YY1,ZZ1,
     $XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2,GREXP)          
      IF(KOMP.EQ.2) GOTO 14                                             
      call lum(1.d0,1.d0,0.d0,1.d0,1.d0,n,n1,1,sbrd,rv,rvq,glump1,      
     $glump2,grv1,grv2,mmsave,summ,fr1,sm,0,vol,q,p,ff,d,snth)          
      GOTO 15                                                           
   14 call lum(1.d0,1.d0,0.d0,1.d0,1.d0,n,n1,2,sbrd,rv,rvq,glump1,      
     $glump2,grv1,grv2,mmsave,summ,fr2,sm,0,vol,q,p,ff,d,snth)          
   15 CONTINUE                                                          
      IF(I.EQ.1) VOLS=VOL                                               
      VOL2=VOLS                                                         
   17 VOL1=VOL                                                          
      IF(IFC.EQ.1) V=VOL                                                
      IF(IFC.EQ.1) RETURN                                               
      DPDV=DP/(VOL2-VOL1)                                               
      DELP=(V-VOL1)*DPDV                                                
      ABDELP=dabs(DELP)                                                 
      IF(ABDELP.GT.TOL) GOTO 16                                         
      P=PS                                                              
      RETURN                                                            
      END                                                               
      SUBROUTINE ATM(TEMP,WAVE,RATIO)                                   
c  Version of October 6, 1995                                           
C     THIS SUBROUTINE COMPUTES THE RATIO OF STELLAR ATMOSPHERE TO       
C     PLANCKIAN FLUX FOR MAIN SEQUENCE STARS.                           
C     WARNING: LIMITS OF APPLICABILITY OF THIS SUBROUTINE AS FOLLOWS.   
C     IN TEMP RANGE 4,000 TO 25,000 K, WAVE LENGTH LIMITS .0912 TO      
C     6.5647 MICRONS. IN TEMP RANGE 25,000 TO 50,000 K, WAVE LENGTH     
C     LIMITS .0912 TO 2.2794 MICRONS. SUBROUTINE SHOULD NOT BE USED     
C     OUTSIDE THESE RANGES UNDER ANY CIRCUMSTANCES. WARNING MESSAGES    
C     ARE PRINTED IF RANGES EXCEEDED.                                   
      implicit real*8 (a-h,o-z)                                         
      DIMENSION A(500),AA(196),AB(182),AC(105)                          
      EQUIVALENCE (A(1),AA(1)),(A(197),AB(1)),(A(379),AC(1))            
      DATA AA/.907883,-1.66066,.691436,4*0.0,.280961,-.447328,.130532,  
     $4*0.0,.142757,-.248947,.0620680,4*0.0,1.20223,-2.33426,1.07589,   
     $4*0.0,.403049,-.586521,.171767,4*0.0,.161172,-.260420,.0640562,   
     $4*0.0,1.78103,-3.64231,1.70609,4*0.0,.667290,-1.06532,.341809,    
     $4*0.0,.181770,-.345798,.0830478,4*0.0,-.329302E+2,.598843E+3,     
     $-.418910E+4,.146768E+5,-.274273E+5,.26093E+5,-.992741E+4,         
     $.442164,-1.17612,.558862,4*0.0,-.0671895,-.114919,.0244139,4*0.0, 
     $-.207863E+3,.243857E+4,-.115507E+5,.283797E+5,-.38233E+5,         
     $.268213E+5,-.76661E+4,.44286,-.984141,.391429,4*0.0,-.0128738,    
     $-.165689,.0419392,4*0.0,-.245193E+3,.287031E+4,-.135822E+5,.333606
     $E+5,-.449465E+5,.315396E+5,-.901805E+4,.542751,-1.13157,.450211,  
     $4*0.0,.0371818,-.212137,.0551231,4*0.0,-.276230E+3,.322251E+4,    
     $-.1522E+5,.373523E+5,-.503142E+5,.353117E+5,-.101001E+5,          
     $.685486,-1.3322,.527152,4*0.0,.099154,-.266633,.0700431,4*0.0,    
     $-.280198E+3,.324071E+4,-.152149E+5,.371972E+5,-.499884E+5,        
     $.350354E+5,-.100136E+5,.912472,-1.64954,.647666,4*0.0,.192641,    
     $-.345623,.0898134,4*0.0,-.649967,.332897E+2,-.395009E+3,          
     $.229755E+4,-.701922E+4,.105395E+5,-.609668E+4,.66853,-1.55813,    
     $.773398,4*0.0,.202931,-.391874,.117945,4*0.0,-1.25062,.246184E+2, 
     $-.934153E+2,-.358292E+2,.745529E+3,-.133152E+4,.738503E+3/        
      DATA AB/.779381,-1.72595,.842554,4*0.0,.269815,-.455823,.135618,  
     $4*0.0,-1.86875,-.233737E+1,.435333E+3,-.330973E+4,.100351E+5,     
     $-.137955E+5,.714497E+4,.234313,-.196784,-.142740,4*0.0,.341826,   
     $-.513333,.149367,4*0.0,                                           
     $        -1.97002,-.28814E+2,.778254E+3,-.498095E+4,.139559E+5,    
     $-.182625E+5,.913243E+4,.607074,-.994799,.290486,4*0.0,.360036,    
     $-.520072,.147628,4*0.0,-1.99412,-.300984E+2,.735863E+3,-.448795E+4
     $,.121397E+5,-.154572E+5,.756405E+4,.6119,-.976964,.278779,4*0.0,  
     $.346578,-.487175,.136367,4*0.0,-.227562E+1,-.734396E+1,           
     $.250026E+3,-.110376E+4,.179596E+4,-.947930E+3,-.903513E+2,        
     $.650197,-1.01937,.297332,4*0.0,.295735,-.382176,.09522,4*0.0,     
     $-.242829E+1,.26082E+2,-495.932,.400884E+4,-.135251E+5,.201766E+5, 
     $-.110779E+5,.515008,-.840408,.268146,4*0.0,.0567424,-.0326554,    
     $-.0248079,4*0.0,-1.98073,5.89063,.923565E+2,-.180606E+2,-.324872E+
     $4,.106725E+5,-.981050E+4,.696006,-.158377E+1,.840335,4*0.0,       
     $-.0372941,.113422,-.0827454,4*0.0,-2.20129,2.57431,.954215E+2,    
     $.140538E+3,-.358564E+4,.10078E+5,-.843473E+4,.877022,-2.11887,    
     $.123868E+1,4*0.0,-.157214,.322907,-.162776,4*0.0,-2.28223,        
     $1.88918,.906366E+2,.182422E+3,-.344699E+4,.915905E+4,-.736667E+4, 
     $.422139,-.673966,.13635,4*0.0,-.269755,.515384,-.234871,4*0.0/    
      DATA AC/-1.72294,.111136E+2,.546943E+2,-.26662E+3,-.827417E+3,    
     $.463518E+4,-.501896E+4,.750032,-1.87636,1.12717,4*0.0,-.327529,   
     $.664143,-.307411,4*0.0,-1.82965,7.19243,.333359E+2,-.189727E+2,   
     $-.567324E+3,.71019E+3,.553124E+3,.968157,-.249486E+1,1.57956,4*0.0
     $,-.413511,.836881,-.379774,4*0.0,-2.12746,5.72613,.206067E+2,     
     $.103493E+3,-.379424E+3,-.115088E+4,.285070E+4,1.30615,-3.49355,   
     $2.31489,4*0.0,-.551384,.109752E+1,-.485262,4*0.0,-2.44968,        
     $5.38777,.918011E+1,.156549E+3,-.162539E+3,-.191797E+4,.316283E+4, 
     $.184856E+1,-.512211E+1,.348495E+1,4*0.0,-.860859,.163724E+1,      
     $-.699413,4*0.0,-2.36238,6.98545,.340574E+2,.139450E+3,            
     $-.680221E+3,-.155574E+4,.462879E+4,.18498E+1,-.579597E+1,         
     $.43628E+1,4*0.0,-.111754E+1,.214057E+1,-.924189,4*0.0/            
   51 FORMAT(' TEMPERATURE RANGE EXCEEDED IN SUBROUTINE ATM')           
   52 FORMAT(' WAVE LENGTH RANGE EXCEEDED IN SUBROUTINE ATM')           
      IF (TEMP.GT.50000.d0) WRITE(6,51)                                 
      IF (TEMP.LT.4000.d0) WRITE (6,51)                                 
      IF (WAVE.LT..0912d0) WRITE(6,52)                                  
      IF(TEMP.LE.25000.d0) GOTO 53                                      
      IF(WAVE.GT.2.2794d0) WRITE(6,52)                                  
      GOTO 54                                                           
   53 IF(WAVE.GT.6.5647d0) WRITE(6,52)                                  
   54 CONTINUE                                                          
      IF (WAVE.GT..3647d0) GOTO 15                                      
      NWAVE=0                                                           
      GOTO 30                                                           
   15 IF (WAVE.GT..8206d0) GOTO 25                                      
      NWAVE=1                                                           
      GOTO 30                                                           
   25 NWAVE=2                                                           
   30 IF (TEMP.LT.40000.d0) GOTO 300                                    
      FRACT=(TEMP-40000.d0)/10000.d0                                    
      DIV1=.0203d0                                                      
      DIV2=.0280d0                                                      
      NUMBER=0                                                          
      GOTO 500                                                          
  300 IF (TEMP.LT.30000.d0) GOTO 305                                    
      FRACT=(TEMP-30000.d0)/10000.d0                                    
      DIV1=.0280d0                                                      
      DIV2=.0302d0                                                      
      NUMBER=1                                                          
      GOTO 500                                                          
  305 IF (TEMP.LT.25000.d0) GOTO 310                                    
      FRACT=(TEMP-25000.d0)/5000.d0                                     
      DIV1=.0302d0                                                      
      DIV2=.0701d0                                                      
      NUMBER=2                                                          
      GOTO 500                                                          
  310 IF (TEMP.LT.20000.d0) GOTO 315                                    
      FRACT=(TEMP-20000.d0)/5000.d0                                     
      DIV1=.0701d0                                                      
      DIV2=.0504d0                                                      
      NUMBER=3                                                          
      GOTO 500                                                          
  315 IF (TEMP.LT.18000.d0) GOTO 320                                    
      FRACT=(TEMP-18000.d0)/2000.d0                                     
      DIV1=.0504d0                                                      
      DIV2=.0504d0                                                      
      NUMBER=4                                                          
      GOTO 500                                                          
  320 IF (TEMP.LT.16000.d0) GOTO 325                                    
      FRACT=(TEMP-16000.d0)/2000.d0                                     
      DIV1=.0504d0                                                      
      DIV2=.0504d0                                                      
      NUMBER=5                                                          
      GOTO 500                                                          
  325 IF (TEMP.LT.14000.d0) GOTO 330                                    
      FRACT=(TEMP-14000.d0)/2000.d0                                     
      DIV1=.0504d0                                                      
      DIV2=.0504d0                                                      
      NUMBER=6                                                          
      GOTO 500                                                          
  330 IF (TEMP.LT.12000.d0) GOTO 335                                    
      FRACT=(TEMP-12000.d0)/2000.d0                                     
      DIV1=.0504d0                                                      
      DIV2=.1000d0                                                      
      NUMBER=7                                                          
      GOTO 500                                                          
  335 IF (TEMP.LT.11000.d0) GOTO 340                                    
      FRACT=(TEMP-11000.d0)/1000.d0                                     
      DIV1=.1000d0                                                      
      DIV2=.1000d0                                                      
      NUMBER=8                                                          
      GOTO 500                                                          
  340 IF (TEMP.LT.10000.d0) GOTO 345                                    
      FRACT=(TEMP-10000.d0)/1000.d0                                     
      DIV1=.1000d0                                                      
      DIV2=.0912d0                                                      
      NUMBER=9                                                          
      GOTO 500                                                          
  345 IF (TEMP.LT.9500.d0) GOTO 350                                     
      FRACT=(TEMP-9500.d0)/500.d0                                       
      DIV1=.0912d0                                                      
      DIV2=.0912d0                                                      
      NUMBER=10                                                         
      GOTO 500                                                          
  350 IF (TEMP.LT.9000.d0) GOTO 355                                     
      FRACT=(TEMP-9000.d0)/500.d0                                       
      DIV1=.0912d0                                                      
      DIV2=.0912d0                                                      
      NUMBER=11                                                         
      GOTO 500                                                          
  355 IF (TEMP.LT.8500.d0) GOTO 360                                     
      FRACT=(TEMP-8500.d0)/500.d0                                       
      DIV1=.0912d0                                                      
      DIV2=.0912d0                                                      
      NUMBER=12                                                         
      GOTO 500                                                          
  360 IF (TEMP.LT.8000.d0) GOTO 365                                     
      FRACT=(TEMP-8000.d0)/500.d0                                       
      DIV1=.0912d0                                                      
      DIV2=.0912d0                                                      
      NUMBER=13                                                         
      GOTO 500                                                          
  365 IF (TEMP.LT.7500.d0) GOTO 370                                     
      FRACT=(TEMP-7500.d0)/500.d0                                       
      DIV1=.0912d0                                                      
      DIV2=.1239d0                                                      
      NUMBER=14                                                         
      GOTO 500                                                          
  370 IF (TEMP.LT.7000.d0) GOTO 375                                     
      FRACT=(TEMP-7000.d0)/500.d0                                       
      DIV1=.1239d0                                                      
      DIV2=.1239d0                                                      
      NUMBER=15                                                         
      GOTO 500                                                          
  375 IF (TEMP.LT.6500.d0) GOTO 380                                     
      FRACT=(TEMP-6500.d0)/500.d0                                       
      DIV1=.1239d0                                                      
      DIV2=.1239d0                                                      
      NUMBER=16                                                         
      GOTO 500                                                          
  380 IF (TEMP.LT.6000.d0) GOTO 385                                     
      FRACT=(TEMP-6000.d0)/500.d0                                       
      DIV1=.1239d0                                                      
      DIV2=.1444d0                                                      
      NUMBER=17                                                         
      GOTO 500                                                          
  385 IF (TEMP.LT.5500.d0) GOTO 390                                     
      FRACT=(TEMP-5500.d0)/500.d0                                       
      DIV1=.1444d0                                                      
      DIV2=.1444d0                                                      
      NUMBER=18                                                         
      GOTO 500                                                          
  390 IF (TEMP.LT.5000.d0) GOTO 395                                     
      FRACT=(TEMP-5000.d0)/500.d0                                       
      DIV1=.1444d0                                                      
      DIV2=.1444d0                                                      
      NUMBER=19                                                         
      GOTO 500                                                          
  395 IF (TEMP.LT.4500.d0) GOTO 400                                     
      FRACT=(TEMP-4500.d0)/500.d0                                       
      DIV1=.1444d0                                                      
      DIV2=.1444d0                                                      
      NUMBER=20                                                         
      GOTO 500                                                          
  400 FRACT=(TEMP-4000.d0)/500.d0                                       
      DIV1=.1444d0                                                      
      DIV2=.1623d0                                                      
      NUMBER=21                                                         
  500 M=NUMBER*21+NWAVE*7+1                                             
      N=(NUMBER+1)*21+NWAVE*7+1                                         
      SUM1=A(M)                                                         
      SUM2=A(N)                                                         
      ALWAV1=dlog10(WAVE/DIV1)                                          
      ALWAV2=dlog10(WAVE/DIV2)                                          
      DO 20 I=1,6                                                       
      M1=M+I                                                            
      SUM1=SUM1+A(M1)*ALWAV1**I                                         
      N1=N+I                                                            
      SUM2=SUM2+A(N1)*ALWAV2**I                                         
   20 CONTINUE                                                          
      ALRATI=SUM2*(1.d0-FRACT)+SUM1*FRACT                               
      RATIO=10.d0**ALRATI                                               
      RETURN                                                            
      END                                                               
      subroutine fourls(th,ro,nobs,nth,aa,bb)                           
      implicit real*8(a-h,o-z)                                          
c   version of September 14, 1998                                       
c                                                                       
c    Input integer nth is the largest Fourier term fitted (e.g.         
c       for nth=6, terms up to sine & cosine of 6 theta are             
c       evaluated).                                                     
c    This subroutine can handle nth only up to 6. Additional            
c      programming is needed for larger values.                         
c                                                                       
      dimension aa(*),bb(*),th(*),ro(*),obs(5000),ll(14),mm(14),        
     $cn(196),cl(14),out(14)                                            
      mpl=nth+1                                                         
      ml=mpl+nth                                                        
      jjmax=ml*ml                                                       
      nobsml=nobs*ml                                                    
      nobmpl=nobs*mpl                                                   
      do 90 i=1,nobs                                                    
      obs(i)=1.d0                                                       
      iz=nobsml+i                                                       
      obs(iz)=ro(i)                                                     
      if(nth.eq.0) goto 90                                              
      ic=i+nobs                                                         
      is=i+nobmpl                                                       
      sint=dsin(th(i))                                                  
      cost=dcos(th(i))                                                  
      obs(ic)=cost                                                      
      obs(is)=sint                                                      
      if(nth.eq.1) goto 90                                              
      ic=ic+nobs                                                        
      is=is+nobs                                                        
      sncs=sint*cost                                                    
      cs2=cost*cost                                                     
      obs(ic)=cs2+cs2-1.d0                                              
      obs(is)=sncs+sncs                                                 
      if(nth.eq.2) goto 90                                              
      ic=ic+nobs                                                        
      is=is+nobs                                                        
      sn3=sint*sint*sint                                                
      cs3=cs2*cost                                                      
      obs(ic)=4.d0*cs3-3.d0*cost                                        
      obs(is)=3.d0*sint-4.d0*sn3                                        
      if(nth.eq.3) goto 90                                              
      ic=ic+nobs                                                        
      is=is+nobs                                                        
      cs4=cs2*cs2                                                       
      obs(ic)=8.d0*(cs4-cs2)+1.d0                                       
      obs(is)=4.d0*(2.d0*cs3*sint-sncs)                                 
      if(nth.eq.4) goto 90                                              
      ic=ic+nobs                                                        
      is=is+nobs                                                        
      cs5=cs3*cs2                                                       
      obs(ic)=16.d0*cs5-20.d0*cs3+5.d0*cost                             
      obs(is)=16.d0*sn3*sint*sint-20.d0*sn3+5.d0*sint                   
      if(nth.eq.5) goto 90                                              
      ic=ic+nobs                                                        
      is=is+nobs                                                        
      obs(ic)=32.d0*cs3*cs3-48.d0*cs4+18.d0*cs2-1.d0                    
      obs(is)=32.d0*sint*(cs5-cs3)+6.d0*sncs                            
   90 continue                                                          
      do 20 jj=1,jjmax                                                  
   20 cn(jj)=0.d0                                                       
      do 21 j=1,ml                                                      
   21 cl(j)=0.d0                                                        
      do 24 nob=1,nobs                                                  
      iii=nob+nobsml                                                    
      do 23 k=1,ml                                                      
      do 23 i=1,ml                                                      
      ii=nob+nobs*(i-1)                                                 
      kk=nob+nobs*(k-1)                                                 
      j=i+(k-1)*ml                                                      
   23 cn(j)=cn(j)+obs(ii)*obs(kk)                                       
      do 24 i=1,ml                                                      
      ii=nob+nobs*(i-1)                                                 
   24 cl(i)=cl(i)+obs(iii)*obs(ii)                                      
      call dminv(cn,ml,d,ll,mm)                                         
      call dgmprd(cn,cl,out,ml,ml,1)                                    
      do 51 i=2,mpl                                                     
      aa(i)=out(i)                                                      
      ipl=i+nth                                                         
   51 bb(i)=out(ipl)                                                    
      aa(1)=out(1)                                                      
      bb(1)=0.d0                                                        
      return                                                            
      end                                                               
      SUBROUTINE ELLONE(FF,dd,rm,xl1,OM1,XL2,OM2)                       
c  Version of October 6, 1995                                           
C     XL2 AND OM2 VALUES ASSUME SYNCHRONOUS ROTATION AND CIRCULAR ORBIT.
C     THEY ARE NOT NEEDED FOR NON-SYNCHRONOUS OR NON-CIRCULAR CASES.    
      IMPLICIT REAL*8(A-H,O-Z)                                          
      rmass=rm                                                          
      d=dd                                                              
      XL=.5D0*D                                                         
      DO 5 I=1,2                                                        
      RFAC=ff*ff                                                        
      IF(I.EQ.2) RFAC=1.D0                                              
      IF(I.EQ.2) D=1.D0                                                 
      DSQ=D*D                                                           
      DELXL=0.D0                                                        
      RM1=RMASS+1.D0                                                    
   88 XL=XL+DELXL                                                       
      XSQ=XL*XL                                                         
      P=(D-XL)**2                                                       
      RP=DABS(D-XL)                                                     
      PRP=P*RP                                                          
      F=RFAC*RM1*XL-1.D0/XSQ-RMASS*(XL-D)/PRP-RMASS/DSQ                 
      DXLDF=1.D0/(RFAC*RM1+2.D0/(XSQ*XL)+2.D0*RMASS/PRP)                
      DELXL=-F*DXLDF                                                    
      ABDEL=DABS(DELXL)                                                 
      IF(ABDEL.GT..000001D0) GOTO 88                                    
      IF(I.EQ.2) GOTO 8                                                 
      xl1=xl                                                            
      OM1=1.D0/XL+RMASS*((1.D0/RP)-XL/DSQ)+RM1*.5D0*XSQ*RFAC            
      IF(rm.GT.1.d0) RMASS=1.D0/RMASS                                   
      XMU3=RMASS/(3.D0*(RMASS+1.D0))                                    
      XMU3CR=XMU3**.33333333333D0                                       
    5 XL=1.D0+XMU3CR+XMU3CR*XMU3CR/3.D0+XMU3/9.D0                       
    8 IF(rm.GT.1.d0) XL=D-XL                                            
      RMASS=rm                                                          
      XL2=XL                                                            
      OM2=1.D0/DABS(XL)+RMASS*((1.D0/DSQRT(1.D0-XL-XL+XL*XL))-XL)+RM1*  
     $.5D0*XL*XL                                                        
      RETURN                                                            
      END                                                               
      SUBROUTINE RING(Q,OM,KOMP,L,FR,HLD,R1,RL)                         
c   Version of September 14, 1998                                       
      IMPLICIT REAL*8(A-H,O-Z)                                          
      DIMENSION RAD(100),THET(100),AA(3),BB(3),FI(150),THA(150),FR(*),  
     $HLD(*)                                                            
      IX=0                                                              
      LR=L+1                                                            
      DO 92 I=1,LR                                                      
      THA(I)=0.D0                                                       
   92 FI(I)=-.1D0                                                       
      OMEGA=OM                                                          
      K=3                                                               
      EL=dfloat(L)                                                      
      DEL=2.D0/EL                                                       
      CALL ELLONE(1.d0,1.d0,Q,xlsv,OM1,XL2,OM2)                         
      CALL NEKMIN(Q,OM,xlsv,Z)                                          
      XL=xlsv                                                           
      QQ=Q                                                              
      XLSQ=XL*XL                                                        
      IF(Q.GT.1.D0) QQ=1.D0/Q                                           
      RMAX=DEXP(.345D0*DLOG(QQ)-1.125D0)                                
      R=RMAX*(OM1-OMEGA)/(OM1-OM2)                                      
      DO 22 IT=1,L                                                      
      EYT=dfloat(IT)                                                    
      TH=EYT*1.570796326794897d0/EL                                     
      COSQ=DCOS(TH)**2                                                  
      DELR=0.D0                                                         
   14 R=DABS(R+DELR)                                                    
      RSQ=R*R                                                           
      X2R2=XLSQ+RSQ                                                     
      RX2R2=DSQRT(X2R2)                                                 
      XM2R2=(XL-1.D0)**2+RSQ                                            
      RXM2R2=DSQRT(XM2R2)                                               
      OM=1.D0/RX2R2+Q*(1.D0/RXM2R2-XL)+.5D0*(Q+1.D0)*(XLSQ+RSQ*COSQ)    
      DOMDR=-R/(X2R2*RX2R2)-Q*R/(XM2R2*RXM2R2)+(Q+1.D0)*COSQ*R          
      DELR=(OMEGA-OM)/DOMDR                                             
      ABDELR=DABS(DELR)                                                 
      IF(ABDELR.GT..00001D0) GOTO 14                                    
      RAD(IT)=R                                                         
   22 THET(IT)=TH*4.D0                                                  
      R1=RAD(1)                                                         
      RL=RAD(L)                                                         
      R90SQ=RL*RL                                                       
      DO 18 IJ=1,L                                                      
      EYJ=IJ                                                            
      RAD(IJ)=RAD(IJ)-(RL-R1)*(EYJ-1.D0)/EL                             
   18 CONTINUE                                                          
      DO 65 N=1,K                                                       
      AA(N)=0.D0                                                        
   65 BB(N)=0.D0                                                        
      DO 29 J=1,L                                                       
      DO 29 N=1,K                                                       
      EN=N-1                                                            
      ENTHET=EN*THET(J)                                                 
      AA(N)=AA(N)+RAD(J)*DCOS(ENTHET)*DEL                               
   29 BB(N)=BB(N)+RAD(J)*DSIN(ENTHET)*DEL                               
      AA(1)=AA(1)*.5D0                                                  
      IF(KOMP.EQ.2) XL=1.d0-xlsv                                        
      XLSQ=XL*XL                                                        
      DIS=RL/XL-.0005D0                                                 
      DO 42 IR=1,L                                                      
      LL=IR-1                                                           
      EY=dfloat(L+1-IR)                                                 
      THA(IR)=1.570796326794897D0*EY/EL                                 
      IF(THA(IR).LT.1.570796326794897D0) GOTO 82                        
      COT=0.D0                                                          
      GOTO 83                                                           
   82 COT=1.D0/DTAN(THA(IR))                                            
   83 IF(COT.GE.DIS) GOTO 50                                            
      COSSQ=DCOS(THA(IR))**2                                            
      A0=AA(1)                                                          
      A1=AA(2)                                                          
      A2=AA(3)                                                          
      B1=BB(2)                                                          
      B2=BB(3)                                                          
      DELSIN=0.D0                                                       
      KNTR=0                                                            
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                             
   88 SINTH=SINTH+DELSIN                                                
      KNTR=KNTR+1                                                       
      IF(SINTH.GT.1.D0) SINTH=1.D0/SINTH                                
      CSQ=1.D0-SINTH*SINTH                                              
      COSTH=DSQRT(CSQ)                                                  
      SINSQ=SINTH*SINTH                                                 
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                         
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                     
      C4SQ=COS4*COS4                                                    
      SINCOS=SIN4*COS4                                                  
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS         
      ARC=dasin(SINTH)                                                  
      RR=RRR+(RL-R1)*(2.D0*ARC/3.141592653589793D0-1.D0/EL)             
      IF(KNTR.GT.30) GOTO 42                                            
      P=RR*SINTH                                                        
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2
     $)*COSTH+(RL+RL-R1-R1)/(3.141592653589793D0*COSTH)                 
      DPDSIN=RR+SINTH*DRDSIN                                            
      F=P*P/COSSQ-RR*RR-XLSQ                                            
      DFDSIN=(P+P)*DPDSIN/COSSQ-(RR+RR)*DRDSIN                          
      DELSIN=-F/DFDSIN                                                  
      ABDEL=DABS(DELSIN)                                                
      IF(ABDEL.GT..00001D0)  GOTO 88                                    
   42 FI(IR)=DATAN(RR*COSTH/XL)                                         
   50 LL1=LL+1                                                          
      DELTH=1.570796326794897D0/EL                                      
      DO 75 I=1,L                                                       
      EY=dfloat(L+1-I)-.5d0                                             
      THE=1.570796326794897D0*EY/EL                                     
      SNTH=DSIN(THE)                                                    
      EM=dsin(THE)*EL*1.3d0                                             
      MM=EM+1.d0                                                        
      XM=dfloat(MM)                                                     
      DELFI=3.141592653589793D0/XM                                      
      HDELFI=1.570796326794897D0/XM                                     
      DO 75 J=1,MM                                                      
      IX=IX+1                                                           
      IF(I.LE.LL1) GOTO 43                                              
      HLD(IX)=1.d0                                                      
      GOTO 75                                                           
   43 XJ=MM+1-J                                                         
      FE=3.141592653589793D0*(XJ-.5D0)/XM                               
      PH2=FE+HDELFI                                                     
      PHB=PH2                                                           
      IF(FI(I).GT.(FE-HDELFI)) GOTO 51                                  
      HLD(IX)=1.d0                                                      
      GOTO 75                                                           
   51 IPL=I+1                                                           
      IF(FI(IPL).GT.0.D0) GOTO 66                                       
      RR=A0+A1-A2+(RL-R1)*(1.D0-1.D0/EL)                                
      PH1=DELFI*(XJ-1.D0)                                               
      TH1=DATAN(XL/RR)                                                  
      GOTO 56                                                           
   66 IF(FI(IPL).LT.(FE+HDELFI)) GOTO 52                                
      HLD(IX)=0.d0                                                      
      GOTO 75                                                           
   52 IF(FI(IPL).LT.(FE-HDELFI)) GOTO 53                                
      PH1=FI(IPL)                                                       
      TH1=THA(IPL)                                                      
      GOTO 56                                                           
   53 DELSIN=0.D0                                                       
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                             
      TANFE=DTAN(FE-HDELFI)                                             
   77 SINTH=SINTH+DELSIN                                                
      IF(SINTH.GT.1.D0) SINTH=1.D0/SINTH                                
      SINSQ=SINTH*SINTH                                                 
      CSQ=1.D0-SINSQ                                                    
      COSTH=DSQRT(CSQ)                                                  
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                         
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                     
      C4SQ=COS4*COS4                                                    
      SINCOS=SIN4*COS4                                                  
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS         
      ARC=dasin(SINTH)                                                  
      RR=RRR+(RL-R1)*(2.D0*ARC/3.141592653589793D0-1.D0/EL)             
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2
     $)*COSTH+(RL+RL-R1-R1)/(3.141592653589793D0*COSTH)                 
      F=RR*COSTH-XL*TANFE                                               
      DFDSIN=COSTH*DRDSIN-RR*SINTH/COSTH                                
      DELSIN=-F/DFDSIN                                                  
      ABDEL=DABS(DELSIN)                                                
      IF(ABDEL.GT..00001D0)  GOTO 77                                    
      PH1=FE-HDELFI                                                     
      TH1=DATAN(XL/(RR*SINTH*DCOS(PH1)))                                
   56 IF(FI(I).GT.(FE+HDELFI)) GOTO 57                                  
      PHB=FI(I)                                                         
      TH2=THA(I)                                                        
      GOTO 60                                                           
   57 DELSIN=0.D0                                                       
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                             
      TANFE=DTAN(FE+HDELFI)                                             
   78 SINTH=SINTH+DELSIN                                                
      IF(SINTH.GT.1.D0) SINTH=1.D0/SINTH                                
      SINSQ=SINTH*SINTH                                                 
      CSQ=1.D0-SINSQ                                                    
      COSTH=DSQRT(CSQ)                                                  
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                         
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                     
      C4SQ=COS4*COS4                                                    
      SINCOS=SIN4*COS4                                                  
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS         
      ARC=dasin(SINTH)                                                  
      RR=RRR+(RL-R1)*(2.D0*ARC/3.141592653589793D0-1.D0/EL)             
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2
     $)*COSTH+(RL+RL-R1-R1)/(3.141592653589793D0*COSTH)                 
      F=RR*COSTH-XL*TANFE                                               
      DFDSIN=COSTH*DRDSIN-RR*SINTH/COSTH                                
      DELSIN=-F/DFDSIN                                                  
      ABDEL=DABS(DELSIN)                                                
      IF(ABDEL.GT..00001D0)  GOTO 78                                    
      TH2=DATAN(XL/(RR*SINTH*DCOS(PH2)))                                
   60 CTHT=DCOS(THA(IPL))                                               
      CTH1=DCOS(TH1)                                                    
      CTH2=DCOS(TH2)                                                    
      STH1=DSIN(TH1)                                                    
      STH2=DSIN(TH2)                                                    
      DTH=TH2-TH1                                                       
      DCTH=CTH1-CTH2                                                    
      OMDP=PH2*DCTH-.5D0*(PH1*STH1+PHB*STH2)*DTH                        
      OMP=DELFI*(CTHT-CTH1)                                             
      OMN=OMP+OMDP                                                      
      HLD(IX)=OMN/(DELTH*DELFI*SNTH)                                    
   75 CONTINUE                                                          
      DO 94 JB=1,IX                                                     
      JA=IX+1-JB                                                        
   94 FR(JB)=HLD(JA)                                                    
      RETURN                                                            
      END                                                               
      SUBROUTINE MODLOG(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2
     $,HLD,RM,POTH,POTC,GR1,GR2,ALB1,ALB2,N1,N2,F1,F2,MOD,XINCL,THE,    
     $MODE,SNTH,CSTH,SNFI,CSFI,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,GLUMP1 
     $,GLUMP2,CSBT1,CSBT2,GMAG1,GMAG2)                                  
c    Version of July 2, 1999                                       
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ  
     $(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),GRV1(*),GRV2(*),XX1(*),YY1(*), 
     $ZZ1(*),XX2(*),YY2(*),ZZ2(*),GLUMP1(*),GLUMP2(*),CSBT1(*),CSBT2(*) 
     $,GMAG1(*),GMAG2(*)                                                
      DIMENSION DRR(4),RES(2),ANS(2),LX(2),MX(2)                        
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*)                         
      common /kfac/ kff1,kff2,kfo1,kfo2                                 
      common /setest/ sefac                                             
      common /ardot/ dperdt,hjd,hjd0,perr                               
      COMMON /FLVAR/ PSHIFT,DP,DS,EF,EFC,ECOS,perr0,PHPER,PCONJ,        
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc,qfacd       
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC     
      COMMON /INVAR/ KH,IPBDUM,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,     
     $ifsmv2,icor1,icor2,ld,ncl,jdphs,ipc                               
   95 FORMAT(' WARNING: ALTHOUGH COMPONENT 2 DOES NOT EXCEED ITS LIMITIN
     $G LOBE AT THE END OF ECLIPSE, IT DOES EXCEED THE LOBE AT PERIASTRO
     $N')                                                               
   99 FORMAT(' SPECIFIED ECLIPSE DURATION INCONSISTENT WITH OTHER PARAME
     $TERS')                                                            
      perr=perr0+dperdt*(hjd-hjd0)                       
      DP=1.d0-E                                                         
      DS=DP*DP                                                          
      MOD=(MODE-2)**2                                                   
      IF(MOD.EQ.1) GR2=GR1                                              
      IF(MOD.EQ.1) ALB2=ALB1                                            
      IF(MOD.EQ.1) POTC=POTH                                            
      MD4=(MODE-5)**2                                                   
      MD5=(2*MODE-11)**2                                                
      call ellone(f1,dp,rm,xl1,po1cr,xl2,omo1)                          
      sefac=.8712d0                                                     
      doc=(po1cr-poth)/(po1cr-omo1)                                     
      if(doc.gt.0.d0) sefac=.201d0*doc*doc-.386d0*doc+.8712d0           
      RMR=1.d0/RM                                                       
      CALL ELLONE(F2,DP,RMR,XL1,po2c,XL2,omo2)                          
      po2cr=rm*po2c+(1.d0-rm)*.5d0                                      
      if(md4.eq.1) poth=po1cr                                           
      if(md5.eq.1) potc=po2cr                                           
      kff1=0                                                            
      kff2=0                                                            
      if(poth.lt.po1cr) kff1=1                                          
      if(potc.lt.po2cr) kff2=1                                          
      kfo1=0                                                            
      kfo2=0                                                            
      if(e.ne.0.d0) goto 100                                            
      if(f1.ne.1.d0) goto 105                                           
      if(poth.lt.omo1) kfo1=1                                           
  105 if(f2.ne.1.d0) goto 100                                           
      if(potc.lt.omo1) kfo2=1                                           
  100 continue                                                          
      SINI=dsin(.017453292519943d0*XINCL)                               
      VF=50.61455d0/PERIOD                                              
      VFAC=VF*A                                                         
      VGAM=VGA*VUNIT/VFAC                                               
      VFVU=VFAC/VUNIT                                                   
      IFC=2                                                             
      IF(E.NE.0.d0) GOTO 60                                             
      perr=1.570796326794897d0                                          
      IFC=1                                                             
   60 CONTINUE                                                          
      TRC=1.570796326794897d0-perr                                      
   39 if(TRC.LT.0.d0) TRC=TRC+6.283185307179586d0                       
      if(trc.lt.0.d0) goto 39                                           
   40 if(trc.ge.6.283185307179586d0) trc=trc-6.283185307179586d0        
      if(trc.ge.6.283185307179586d0) goto 40                            
      HTRC=.5d0*TRC                                                     
      IF(dabs(1.570796326794897d0-HTRC).LT.7.d-6) GOTO 101              
      IF(dabs(4.712388980384690d0-HTRC).LT.7.d-6) GOTO 101              
      ECAN=2.d0*datan(dsqrt((1.d0-E)/(1.d0+E))*dtan(HTRC))              
      GOTO 103                                                          
  101 ECAN=3.141592653589793d0                                          
  103 XMC=ECAN-E*dsin(ECAN)                                             
      IF(XMC.LT.0.d0) XMC=XMC+6.283185307179586d0                       
      PHPER=1.d0-XMC/6.283185307179586d0                                
      PCONJ=(XMC+perr)/6.283185307179586d0-.25d0+PSHIFT                 
   38 if(pconj.ge.1.d0) pconj=pconj-1.d0                                
      if(pconj.ge.1.d0) goto 38                                         
   41 if(pconj.lt.0.d0) pconj=pconj+1.d0                                
      if(pconj.lt.0.d0) goto 41                                         
      PHPERI=PHPER+PCONJ                                                
      EF=1.d0-E*E                                                       
      EFC=dsqrt(EF)                                                     
      ECOS=E*dcos(perr)                                                 
      IF(MODE.NE.-1) RETURN                                             
      if(kh.eq.17) goto 241                                             
      if((kh-12)**2.eq.1) goto 241                                      
      if((kh-12)**2.eq.4) goto 241                                      
      IF((KH-11)**2.LE.1) GOTO 241                                      
      IF((2*KH-41)**2.EQ.81) GOTO 241                                   
      RETURN                                                            
  241 CONTINUE                                                          
      EFCC=dsqrt((1.d0-E)/(1.d0+E))                                     
      THER=THE*6.283185307179586d0                                      
      DELTR=.001d0                                                      
      DTR1=0.d0                                                         
      DTR2=0.d0                                                         
      VOLTOL=5.d-6                                                      
      DXMTOL=5.d-6                                                      
      TR0=1.570796326794897d0-perr                                      
      HTR0=.5d0*TR0                                                     
      IF((1.570796326794897d0-dabs(HTR0)).LT.7.d-6) GOTO 201            
      IF((4.712388980384690d0-dabs(HTR0)).LT.7.d-6) GOTO 201            
      ECAN0=2.d0*datan(dsqrt((1.d0-E)/(1.d0+E))*dtan(HTR0))             
      GOTO 203                                                          
  201 ECAN0=3.141592653589793d0                                         
  203 XM0=ECAN0-E*dsin(ECAN0)                                           
      XM1=XM0-THER*(1.d0-.2d0*E)                                        
      XM2=XM0+THER*(1.d0-.2d0*E)                                        
      CALL KEPLER(XM1,E,DUM,TRR1)                                       
      CALL KEPLER(XM2,E,DUM,TRR2)                                       
  160 TRR1=TRR1+DTR1                                                    
      TRR2=TRR2+DTR2                                                    
      DO 161 IB=1,3                                                     
      TR1=TRR1                                                          
      TR2=TRR2                                                          
      IF(IB.EQ.2) TR1=TRR1+DELTR                                        
      IF(IB.EQ.3) TR2=TRR2+DELTR                                        
      IF(TR1.GT.TR0) TR0=TR0+6.283185307179586d0                        
      IF(TR0.GT.TR2) TR2=TR2+6.283185307179586d0                        
      DS1=EF/(1.d0+E*dcos(TR1))                                         
      DS2=EF/(1.d0+E*dcos(TR2))                                         
      TRE1=(TR0-TR1)/6.283185307179586d0                                
      TRE2=(TR2-TR0)/6.283185307179586d0                                
      CALL DURA(F2,XINCL,RM,DS1,TRE1,POTR,RA)                           
      CALL VOLUME(VS1,RM,POTR,DS1,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ    
     $,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,GRV1, 
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,     
     $GMAG2,GR1,1)                                                      
      CALL DURA(F2,XINCL,RM,DS2,TRE2,POTR,RA)                           
      CALL VOLUME(VS2,RM,POTR,DS2,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ    
     $,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,GRV1, 
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,     
     $GMAG2,GR2,1)                                                      
      IF(IB.NE.1) GOTO 185                                              
      ECAN1=2.d0*datan(dsqrt((1.d0-E)/(1.d0+E))*dtan(.5d0*TR1))         
      ECAN2=2.d0*datan(dsqrt((1.d0-E)/(1.d0+E))*dtan(.5d0*TR2))         
      POTC=POTR                                                         
      DTHE=DS2                                                          
      DVOL=VS2-VS1                                                      
      XM1=ECAN1-E*dsin(ECAN1)                                           
      XM2=ECAN2-E*dsin(ECAN2)                                           
      IF(XM1.LT.0.d0) XM1=XM1+6.283185307179586d0                       
      IF(XM2.LT.0.d0) XM2=XM2+6.283185307179586d0                       
      DXM=XM2-XM1-2.d0*THER                                             
      DDMDN1=-EFCC*(1.d0-E*dcos(ECAN1))*dcos(.5d0*ECAN1)**2/            
     $dcos(.5d0*tr1)**2                                                 
      DDMDN2=EFCC*(1.d0-E*dcos(ECAN2))*dcos(.5d0*ECAN2)**2/             
     $dcos(.5d0*tr2)**2                                                 
  185 CONTINUE                                                          
      IF(IB.NE.2) GOTO 162                                              
      DRR(1)=(VS2-VS1-DVOL)/DELTR                                       
      DRR(2)=DDMDN1                                                     
  162 CONTINUE                                                          
      IF(IB.NE.3) GOTO 161                                              
      DRR(3)=(VS2-VS1-DVOL)/DELTR                                       
      DRR(4)=DDMDN2                                                     
  161 CONTINUE                                                          
      RES(1)=-DVOL                                                      
      RES(2)=-DXM                                                       
      CALL DMINV(DRR,2,DUMM,LX,MX)                                      
      CALL DGMPRD(DRR,RES,ANS,2,2,1)                                    
      DTR1=ANS(1)                                                       
      DTR2=ANS(2)                                                       
      IF(dabs(DTR1).GT.VOLTOL) GOTO 160                                 
      IF(dabs(DTR2).GT.DXMTOL) GOTO 160                                 
      POTH=9999.99d0                                                    
      RMR=1.d0/RM                                                       
      CALL ELLONE(F2,DTHE,RMR,XLA,OM1,XL2,OM2)                          
      OM1=RM*OM1+(1.d0-RM)*.5d0                                         
      IF(POTC.LT.OM1) GOTO 22                                           
      IF(RA.LE.XLA) GOTO 28                                             
   22 WRITE(6,99)                                                       
      RETURN                                                            
   28 CONTINUE                                                          
      IF(E.NE.0.d0) CALL VOLUME(VTHE,RM,POTC,DTHE,F2,N2,N1,2,RV,GRX,    
     $GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,
     $SUMMD,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,   
     $GLUMP2,GMAG1,GMAG2,GR2,1)                                         
      IF(E.NE.0.d0) CALL VOLUME(VTHE,RM,POTC,DP,F2,N2,N1,2,RV,GRX,      
     $GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,
     $SUMMD,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,   
     $GLUMP2,GMAG1,GMAG2,GR2,2)                                         
      CALL ELLONE(F2,DP,RMR,XLD,OMP,XL2,OM2)                            
      OMP=RM*OMP+(1.d0-RM)*.5d0                                         
      IF(POTC.LT.OMP) WRITE(6,95)                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE KEPLER(XM,EC,ECAN,TR)                                  
c  Version of October 6, 1995                                           
      IMPLICIT REAL*8(A-H,O-Z)                                          
      TOL=1.d-8                                                         
      DLECAN=0.D0                                                       
      ECAN=XM                                                           
   18 ECAN=ECAN+DLECAN                                                  
      XMC=ECAN-EC*DSIN(ECAN)                                            
      DEDM=1.D0/(1.D0-EC*DCOS(ECAN))                                    
      DLECAN=(XM-XMC)*DEDM                                              
      ABDLEC=DABS(DLECAN)                                               
      IF(ABDLEC.GT.TOL) GOTO 18                                         
      TR=2.D0*DATAN(DSQRT((1.D0+EC)/(1.D0-EC))*DTAN(.5D0*ECAN))         
      IF(TR.LT.0.) TR=TR+6.2831853071795865d0                           
      RETURN                                                            
      END                                                               
      SUBROUTINE DGMPRD(A,B,R,N,M,L)                                    
c  Version of April 9, 1992                                             
      DIMENSION A(*),B(*),R(*)                                          
      DOUBLE PRECISION A,B,R                                            
      IR=0                                                              
      IK=-M                                                             
      DO 10 K=1,L                                                       
      IK=IK+M                                                           
      DO 10 J=1,N                                                       
      IR=IR+1                                                           
      JI=J-N                                                            
      IB=IK                                                             
      R(IR)=0.D0                                                        
      DO 10 I=1,M                                                       
      JI=JI+N                                                           
      IB=IB+1                                                           
   10 R(IR)=R(IR)+A(JI)*B(IB)                                           
      RETURN                                                            
      END                                                               
      SUBROUTINE DMINV(A,N,D,L,M)                                       
c  Version of April 9, 1992                                             
      DIMENSION A(*),L(*),M(*)                                          
      DOUBLE PRECISION A,D,BIGA,HOLD                                    
      D=1.D0                                                            
      NK=-N                                                             
      DO 80 K=1,N                                                       
      NK=NK+N                                                           
      L(K)=K                                                            
      M(K)=K                                                            
      KK=NK+K                                                           
      BIGA=A(KK)                                                        
      DO 20 J=K,N                                                       
      IZ=N*(J-1)                                                        
      DO 20 I=K,N                                                       
      IJ=IZ+I                                                           
   10 IF(DABS(BIGA).GE.DABS(A(IJ))) GOTO 20                             
      BIGA=A(IJ)                                                        
      L(K)=I                                                            
      M(K)=J                                                            
   20 CONTINUE                                                          
      J=L(K)                                                            
      IF(J.LE.K) GOTO 35                                                
      KI=K-N                                                            
      DO 30 I=1,N                                                       
      KI=KI+N                                                           
      HOLD=-A(KI)                                                       
      JI=KI-K+J                                                         
      A(KI)=A(JI)                                                       
   30 A(JI) =HOLD                                                       
   35 I=M(K)                                                            
      IF(I.LE.K) GOTO 45                                                
      JP=N*(I-1)                                                        
      DO 40 J=1,N                                                       
      JK=NK+J                                                           
      JI=JP+J                                                           
      HOLD=-A(JK)                                                       
      A(JK)=A(JI)                                                       
   40 A(JI) =HOLD                                                       
   45 IF(BIGA.NE.0.D0) GOTO 48                                          
      D=0.D0                                                            
      RETURN                                                            
   48 DO 55 I=1,N                                                       
      IF(I.EQ.K) GOTO 55                                                
      IK=NK+I                                                           
      A(IK)=A(IK)/(-BIGA)                                               
   55 CONTINUE                                                          
      DO 65 I=1,N                                                       
      IK=NK+I                                                           
      HOLD=A(IK)                                                        
      IJ=I-N                                                            
      DO 65 J=1,N                                                       
      IJ=IJ+N                                                           
      IF(I.EQ.K) GOTO 65                                                
      IF(J.EQ.K) GOTO 65                                                
      KJ=IJ-I+K                                                         
      A(IJ)=HOLD*A(KJ)+A(IJ)                                            
   65 CONTINUE                                                          
      KJ=K-N                                                            
      DO 75 J=1,N                                                       
      KJ=KJ+N                                                           
      IF(J.EQ.K) GOTO 75                                                
      A(KJ)=A(KJ)/BIGA                                                  
   75 CONTINUE                                                          
      D=D*BIGA                                                          
      A(KK)=1.D0/BIGA                                                   
   80 CONTINUE                                                          
      K=N                                                               
  100 K=(K-1)                                                           
      IF(K.LE.0) RETURN                                                 
      I=L(K)                                                            
      IF(I.LE.K) GOTO 120                                               
      JQ=N*(K-1)                                                        
      JR=N*(I-1)                                                        
      DO 110 J=1,N                                                      
      JK=JQ+J                                                           
      HOLD=A(JK)                                                        
      JI=JR+J                                                           
      A(JK)=-A(JI)                                                      
  110 A(JI) =HOLD                                                       
  120 J=M(K)                                                            
      IF(J.LE.K) GOTO 100                                               
  125 KI=K-N                                                            
      DO 130 I=1,N                                                      
      KI=KI+N                                                           
      HOLD=A(KI)                                                        
      JI=KI-K+J                                                         
      A(KI)=-A(JI)                                                      
  130 A(JI) =HOLD                                                       
      GO TO 100                                                         
      END                                                               
      SUBROUTINE NEKMIN(RM,OMEG,X,Z)                                    
c  Version of October 9, 1995                                           
      IMPLICIT REAL*8(A-H,O-Z)                                          
      DIMENSION DN(4),EN(2),OUT(2),LL(2),MM(2)                          
      Z=.05d0                                                           
   15 P1=X*X+Z*Z                                                        
      RP1=DSQRT(P1)                                                     
      P115=P1*RP1                                                       
      P2=(1.d0-X)**2+Z*Z                                                
      RP2=DSQRT(P2)                                                     
      P215=P2*RP2                                                       
      DODZ=-Z/P115-RM*Z/P215                                            
      OM=1.d0/RP1+RM/RP2+(1.d0+RM)*.5d0*X*X-RM*X                        
      DELOM=OMEG-OM                                                     
      DELZ=DELOM/DODZ                                                   
      Z=DABS(Z+DELZ)                                                    
      ABDELZ=DABS(DELZ)                                                 
      IF(ABDELZ.GT..00001d0) GOTO 15                                    
   16 P1=X*X+Z*Z                                                        
      RP1=DSQRT(P1)                                                     
      P115=P1*RP1                                                       
      P125=P1*P115                                                      
      P2=(1.d0-X)**2+Z*Z                                                
      RP2=DSQRT(P2)                                                     
      P215=P2*RP2                                                       
      P225=P2*P215                                                      
      DN(1)=-X/P115+RM*(1.d0-X)/P215+(1.d0+RM)*X-RM                     
      DN(2)=(3.d0*X*X-P1)/P125+(3.d0*RM*(1.d0-X)**2-RM*((1.d0-X)**2     
     $+z*z))/p225+(RM+1.d0)                                             
      DN(3)=-Z/P115-RM*Z/P215                                           
      DN(4)=3.d0*X*Z/P125-3.d0*RM*Z*(1.d0-X)/P225                       
      OME=1.d0/RP1+RM/RP2+(1.d0+RM)*.5d0*X*X-RM*X                       
      EN(1)=OMEG-OME                                                    
      EN(2)=-DN(1)                                                      
      CALL DMINV(DN,2,D,LL,MM)                                          
      CALL DGMPRD(DN,EN,OUT,2,2,1)                                      
      DT1=OUT(1)                                                        
      DT2=OUT(2)                                                        
      ABDX=DABS(DT1)                                                    
      X=X+DT1                                                           
      ABDZ=DABS(DT2)                                                    
      Z=Z+DT2                                                           
      IF(ABDX.GT.1.d-8) GOTO 16                                         
      IF(ABDZ.GT.1.d-8) GOTO 16                                         
      RETURN                                                            
      END                                                               
      Subroutine lum (xlum,x,y,wl,tpoll,n,n1,nstar,sbr,rv,rvq,glump1,   
     $glump2,grv1,grv2,mmsave,summ,fr,sm,ifat,vol,rm,om,f,d,snth)       
c   Version of September 14, 1998                                       
      implicit real*8 (a-h,o-z)                                         
      dimension rv(*),rvq(*),mmsave(*),fr(*),snth(*),glump1(*),         
     $glump2(*),grv1(*),grv2(*)                                         
      common /radi/ R1H,RLH,R1C,RLC                                     
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm,   
     $is1dm,is2dm,ic1dm,ic2dm,ld,ncl,jdphs,ipc                          
      TPOLE=10000.d0*TPOLL                                              
      KR=0                                                              
      RAT=1.d0                                                          
      RATPOL=1.d0                                                       
      IF(IFAT.NE.0) CALL ATM(TPOLE,WL,RATPOL)                           
      XLUMP=14384.d0/WL                                                 
      POLFAC=dexp(XLUMP/TPOLE)-1.d0                                     
      EN=dfloat(N)                                                      
      DELTH=1.570796326794897d0/EN                                      
      SUM=0.d0                                                          
      SUMM=0.d0                                                         
      SM=0.d0                                                           
      VOL=0.d0                                                          
      DO 36 I=1,N                                                       
      IPN1=I+N1*(NSTAR-1)                                               
      SINTH=SNTH(IPN1)                                                  
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      XM=dfloat(MM)                                                     
      DELFI=3.141592653589793d0/XM                                      
      DFST=DELFI*SINTH                                                  
      SUMJ=0.d0                                                         
      SUMMJ=0.d0                                                        
      SMJ=0.d0                                                          
      VOLJ=0.d0                                                         
      DO 26 J=1,MM                                                      
      IP=(NSTAR-1)*(N1+1)+I                                             
      IX=MMSAVE(IP)+J                                                   
      IF(NSTAR.EQ.1) GOTO 39                                            
      IF(RVQ(IX).EQ.-1.d0) GOTO 25                                      
      R=RVQ(IX)                                                         
      GOTO 49                                                           
   39 IF(RV(IX).EQ.-1.d0) GOTO 25                                       
      R=RV(IX)                                                          
   49 grav=dfloat(2-nstar)*grv1(ix)+dfloat(nstar-1)*grv2(ix)            
      TLOCAL=TPOLE*dsqrt(dsqrt(GRAV))                                   
      IF(IFAT.NE.0) CALL ATM(TLOCAL,WL,RAT)                             
      GRAVM=RAT*POLFAC/(dexp(XLUMP/TLOCAL)-1.d0)                        
      di=dfloat(2-nstar)*glump1(ix)+dfloat(nstar-1)*glump2(ix)          
      DIF=DI*GRAVM                                                      
      DIFF=DI*GRAV                                                      
      SMJ=SMJ+DI                                                        
      SUMJ=SUMJ+DIF                                                     
      SUMMJ=SUMMJ+DIFF                                                  
      VOLJ=VOLJ+R*R*R*FR(IX)                                            
      GOTO 26                                                           
   25 KR=1                                                              
   26 CONTINUE                                                          
      SMJ=SMJ*DELFI                                                     
      SUMJ=SUMJ*DELFI                                                   
      SUMMJ=SUMMJ*DELFI                                                 
      SM=SM+SMJ                                                         
      SUMM=SUMM+SUMMJ                                                   
      VOL=VOL+VOLJ*DFST                                                 
   36 SUM=SUM+SUMJ                                                      
      darkin=3.141592653589793d0*(1.d0-x/3.d0)                          
      if(ld.eq.2) darkin=darkin+.6981317d0*y                            
      if(ld.eq.3) darkin=darkin-.6283185d0*y                            
      SBR=.25d0*RATPOL*XLUM/(SUM*DELTH*DARKIN)                          
      SM=SM*DELTH*4.d0                                                  
      SUMM=SUMM*DELTH*4.d0                                              
      VOL=VOL*1.3333333333333d0*DELTH                                   
      IF(KR.EQ.0) RETURN                                                
      CALL ELLONE(F,D,RM,XL1,OMD,XLD,OMD)                               
      CALL NEKMIN(RM,OM,XL1,ZD)                                         
      IF(NSTAR.EQ.2) XL1=D-XL1                                          
      R1=dfloat(2-nstar)*R1H+dfloat(nstar-1)*R1C                        
      RL=dfloat(2-nstar)*RLH+dfloat(nstar-1)*RLC                        
      VOL=VOL+1.047198d0*XL1*R1*RL                                      
      RETURN                                                            
      END                                                               
      SUBROUTINE LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,         
     $MMSAVE,ALB,WL,TPOLL,SBR,N1,N2,KOMP,IFAT,fr,snth,                  
     $TLD,GLUMP1,GLUMP2,XX1,XX2,YY1,YY2,ZZ1,ZZ2,xbol,ybol               
     $,GRV1,GRV2,SBR1B,SBR2B,RF,RFO,GMAG1,GMAG2,DINT)                   
c   Version of September 14, 1998                                       
      implicit real*8 (a-h,o-z)                                         
      DIMENSION GRX(*),GRY(*),GRZ(*),GRXQ(*),GRYQ(*),grzq(*),           
     $SLUMP1(*),SLUMP2(*),MMSAVE(*),FR(*),SNTH(*),                      
     $TLD(*),GLUMP1(*),GLUMP2(*),XX1(*),XX2(*),YY1(*)                   
     $,YY2(*),ZZ1(*),ZZ2(*),GRV1(*),GRV2(*),RF(*),RFO(*),               
     $GMAG1(*),GMAG2(*)                                                 
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm    
     $,ifs1dm,ifs2dm,icr1dm,icr2dm,ld,ncl,jdphs,ipc                     
      IQ=(KOMP-1)*(N1+1)                                                
      IS=(KOMP-1)*MMSAVE(IQ)                                            
      ATRAT=1.d0                                                        
      RATPOL=1.d0                                                       
      PI=3.141592653589793d0                                            
      PIH=.5d0*PI                                                       
      TPOLE=10000.d0*TPOLL                                              
      IF(IFAT.NE.0) CALL ATM(tpole,WL,RATPOL)                           
      CMPP=dfloat(2-KOMP)                                               
      COMPP=dfloat(2*KOMP-3)                                            
      COMP=-COMPP                                                       
      CMP=dfloat(KOMP-1)                                                
      N=(2-KOMP)*N1+(KOMP-1)*N2                                         
      NO=(2-KOMP)*N2+(KOMP-1)*N1                                        
      NOD=2*NO                                                          
      EN=dfloat(N)                                                      
      ENO=dfloat(NO)                                                    
      DELTHO=PIH/ENO                                                    
      XLUMP=14384.d0/WL                                                 
      POLFAC=dexp(XLUMP/TPOLE)-1.d0                                     
      CNST=ALB*DELTHO*SBR2B/(DINT*SBR1B)                                
      IF(KOMP.EQ.2) CNST=ALB*DELTHO*SBR1B/(DINT*SBR2B)                  
      DO 191 I=1,N                                                      
      IPN1=I+N1*(KOMP-1)                                                
      SINTH=SNTH(IPN1)                                                  
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      IP=(KOMP-1)*(N1+1)+I                                              
      IY=MMSAVE(IP)                                                     
      DO 193 J=1,MM                                                     
      IX=IY+J                                                           
      SUM=0.d0                                                          
      IF(FR(IX).EQ.0.d0) GOTO 193                                       
      DO 190 IOTH=1,NOD                                                 
      IOTHS=IOTH                                                        
      IF(IOTH.GT.NO) IOTHS=NOD-IOTH+1                                   
      IPNO=IOTHS+N1*(2-KOMP)                                            
      SINTHO=SNTH(IPNO)                                                 
      EMO=SINTHO*ENO*1.3d0                                              
      MMO=EMO+1.d0                                                      
      MMOD=2*MMO                                                        
      IPO=(2-KOMP)*(N1+1)+IOTHS                                         
      IYO=MMSAVE(IPO)                                                   
      XMO=MMO                                                           
      DELFIO=PI/XMO                                                     
      DO 190 JOFI=1,MMOD                                                
      JOFU=JOFI                                                         
      IF(JOFI.GT.MMO) JOFU=MMOD-JOFI+1                                  
      IXO=IYO+JOFU                                                      
      IX1=IX                                                            
      IX2=IXO                                                           
      IF(KOMP.EQ.1) GOTO 200                                            
      IF(GLUMP1(IXO).EQ.0.d0) GOTO 184                                  
      IX1=IXO                                                           
      IX2=IX                                                            
      GOTO 201                                                          
  200 CONTINUE                                                          
      IF(GLUMP2(IXO).EQ.0.d0) GOTO 179                                  
  201 RTL1=1.d0                                                         
      RTL2=1.d0                                                         
      UPD1=1.d0                                                         
      UPD2=1.d0                                                         
      IF(KOMP.EQ.2) GOTO 22                                             
      IF(JOFI.GT.MMO) RTL2=-1.d0                                        
      IF(IOTH.GT.NO) UPD2=-1.d0                                         
      GOTO 23                                                           
   22 IF(JOFI.GT.MMO) RTL1=-1.d0                                        
      IF(IOTH.GT.NO) UPD1=-1.d0                                         
   23 CONTINUE                                                          
      GX2=GRXQ(IX2)                                                     
      GY2=GRYQ(IX2)*RTL2                                                
      GZ2=GRZQ(IX2)*UPD2                                                
      X1C=XX1(IX1)                                                      
      X2C=XX2(IX2)                                                      
      Y1C=YY1(IX1)*RTL1                                                 
      Y2C=YY2(IX2)*RTL2                                                 
      Z1C=ZZ1(IX1)*UPD1                                                 
      Z2C=ZZ2(IX2)*UPD2                                                 
      DX=(X2C-X1C)*COMP                                                 
      DY=(Y2C-Y1C)*COMP                                                 
      DZ=(Z2C-Z1C)*COMP                                                 
      DLRSQ=DX*DX+DY*DY+DZ*DZ                                           
      CSNUM2=(DX*GX2+DY*GY2+DZ*GZ2)*COMPP                               
      IF(CSNUM2.GE.0.d0) GOTO 190                                       
      GX1=GRX(IX1)                                                      
      GY1=GRY(IX1)*RTL1                                                 
      GZ1=GRZ(IX1)*UPD1                                                 
      CSNUM1=(DX*GX1+DY*GY1+DZ*GZ1)*COMP                                
      IF(CSNUM1.GE.0.d0) GOTO 190                                       
      DMAG=dsqrt(DLRSQ)                                                 
      CSGM1=-CSNUM1/(DMAG*GMAG1(IX1))                                   
      CSGM2=-CSNUM2/(DMAG*GMAG2(IX2))                                   
      IF(KOMP.EQ.2) GOTO 181                                            
      DGAM2=1.d0-XBOL+XBOL*CSGM2                                        
      if(ld.ne.2) goto 179                                              
      if(csgm2.eq.0.d0) goto 179                                        
      dgam2=dgam2-ybol*csgm2*dlog(csgm2)                                
      goto 147                                                          
  179 continue                                                          
      if(ld.eq.3) dgam2=dgam2-ybol*(1.d0-dsqrt(csgm2))                  
  147 if(dgam2.lt.0.d0) dgam2=0.d0                                      
      DSUM=GRV2(IXO)*GLUMP2(IXO)*RFO(IXO)*CSGM1*CSGM2*DGAM2/DLRSQ       
      GOTO 182                                                          
  181 DGAM1=1.d0-XBOL+XBOL*CSGM1                                        
      if(ld.ne.2) goto 184                                              
      if(csgm1.eq.0.d0) goto 184                                        
      dgam1=dgam1-ybol*csgm1*dlog(csgm1)                                
      goto 148                                                          
  184 continue                                                          
      if(ld.eq.3) dgam1=dgam1-ybol*(1.d0-dsqrt(csgm1))                  
  148 if(dgam1.lt.0.d0) dgam1=0.d0                                      
      DSUM=GRV1(IXO)*GLUMP1(IXO)*RFO(IXO)*CSGM2*CSGM1*DGAM1/DLRSQ       
  182 CONTINUE                                                          
      SUM=SUM+DSUM*DELFIO                                               
  190 CONTINUE                                                          
      RF(IX)=(CNST*SUM/(CMPP*GRV1(IX)+CMP*GRV2(IX)))+1.d0               
  193 CONTINUE                                                          
  191 CONTINUE                                                          
      DO 8 I=1,N                                                        
      IPN1=I+N1*(KOMP-1)                                                
      SINTH=SNTH(IPN1)                                                  
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      IP=(KOMP-1)*(N1+1)+I                                              
      IY=MMSAVE(IP)                                                     
      DO 8 J=1,MM                                                       
      IS=IS+1                                                           
      IX=IY+J                                                           
      IF(FR(IX).EQ.0.d0) GOTO 8                                         
      IF(KOMP.EQ.1) GRV=GRV1(IX)                                        
      IF(KOMP.EQ.2) GRV=GRV2(IX)                                        
      TNEW=TPOLE*dsqrt(dsqrt(GRV*RF(IX)))                               
      TLD(IS)=TNEW                                                      
      IF(IFAT.EQ.0) GOTO 18                                             
      CALL ATM(TNEW,WL,RATNEW)                                          
      ATRAT=RATNEW/RATPOL                                               
   18 GRREFL=POLFAC*ATRAT/(dexp(XLUMP/TNEW)-1.d0)                       
   37 IF(KOMP.EQ.1) GOTO 77                                             
      slump2(ix)=glump2(ix)*grrefl*sbr                                  
      GOTO 8                                                            
   77 slump1(ix)=glump1(ix)*grrefl*sbr                                  
    8 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE ROMQ(OME,Q,F,D,EC,TH,FI,R,DRDO,DRDQ,DODQ,KOMP,MODE)    
c  Version of September 14, 1998                                        
      implicit real*8 (a-h,o-z)                                         
      theq=1.570796326794897d0                                          
      MOD46=(MODE-5)**2                                                 
      MOD56=(2*MODE-11)**2                                              
      modkom=mode*(komp+komp-3)                                         
      DQ=1.d-4*Q                                                        
      QP=Q+DQ                                                           
      TOL=5.d-6                                                         
C     TH, FI SHOULD BE IN RADIANS.                                      
      sinth=dsin(th)                                                    
      XNUSQ=sinth*sinth                                                 
      XLAM=sinth*dcos(FI)                                               
      RMA=Q                                                             
      QF=1.d0                                                           
      DP=1.d0-EC                                                        
      QFM=1.d0                                                          
      IF(KOMP.NE.2) GOTO 23                                             
      RMA=1.d0/Q                                                        
      QF=1.d0/Q                                                         
      QFM=-1.d0/Q**2                                                    
   23 CONTINUE                                                          
      CALL ELLONE(F,DP,RMA,X,OMEG,XLD,OMD)                              
      OM2SAV=OMEG                                                       
      RMAP=QP                                                           
      IF(KOMP.NE.2) GOTO 92                                             
      OMEG=OMEG*Q+(1.d0-Q)*.5d0                                         
      IF(MOD56.EQ.1) OME=OMEG                                           
      RMAP=1.d0/QP                                                      
      GOTO 93                                                           
   92 CONTINUE                                                          
      IF(MOD46.EQ.1) OME=OMEG                                           
   93 CONTINUE                                                          
      POT=OME                                                           
      IF(KOMP.EQ.2) POT=OME/Q+.5d0*(Q-1.d0)/Q                           
      CALL ELLONE(F,DP,RMAP,XP,OMP,XLD,OMD)                             
      DODQ=(OMP-OM2SAV)/DQ                                              
      RM1=RMA+1.d0                                                      
      DS=D*D                                                            
      RF=F*F                                                            
      R=1.d0/POT                                                        
      KOUNT=0                                                           
      DELR=0.d0                                                         
      IF(FI.NE.0.d0) GOTO 85                                            
      IF(TH.NE.THEQ) GOTO 85                                            
      IF(MODE.EQ.6) GOTO 114                                            
      IF(MODE.NE.4) GOTO 80                                             
      IF(KOMP.EQ.1) GOTO 114                                            
      GOTO 85                                                           
   80 IF(MODE.NE.5) GOTO 85                                             
      IF(KOMP.EQ.2) GOTO 114                                            
   85 CONTINUE                                                          
   14 R=R+DELR                                                          
      KOUNT=KOUNT+1                                                     
      IF(KOUNT.LT.20) GOTO 70                                           
  217 if(mode.eq.6) goto 114                                            
      if(modkom.eq.-4) goto 114                                         
      if(modkom.eq.5) goto 114                                          
      DOMR=-1.d15                                                       
      R=-1.d0                                                           
      GOTO 116                                                          
   70 RSQ=R*R                                                           
      PAR=DS-2.d0*XLAM*R*D+RSQ                                          
      RPAR=dsqrt(PAR)                                                   
      OM=1.d0/R+RMA*(1.d0/RPAR-XLAM*R/DS)+RM1*.5d0*RSQ*XNUSQ*RF         
      DOMR=1.d0/(RF*RM1*XNUSQ*R-1.d0/RSQ-(RMA*(R-XLAM*D))/(PAR*RPAR)-   
     $RMA*XLAM/DS)                                                      
      DELR=(POT-OM)*DOMR                                                
      ABDELR=dabs(DELR)                                                 
      IF(ABDELR.GT.TOL) GOTO 14                                         
      DOMRSV=DOMR                                                       
      IF(R.GE.1.d0) GOTO 217                                            
      IF(FI.NE.0.d0) GO TO 116                                          
      IF(TH.NE.THEQ)GO TO 116                                           
      IF(OME-OMEG) 217,114,116                                          
  114 DOMR=1.d15                                                        
      R=X                                                               
      goto 118                                                          
  116 DRDQ=(1.d0/RPAR-R*XLAM/DS+.5d0*RF*RSQ*XNUSQ)/(1.d0/RSQ+RMA*       
     $((1.d0/(PAR*RPAR))*(R-XLAM*D)+XLAM/DS)-RF*XNUSQ*RM1*R)            
      DRDQ=DRDQ*QFM                                                     
  118 drdo=domr*qf                                                      
      IF(MODE.EQ.6) GOTO 215                                            
      IF(MODE.NE.4) GOTO 180                                            
      IF(KOMP.EQ.1) GOTO 215                                            
      RETURN                                                            
  180 IF(MODE.NE.5) RETURN                                              
      IF(KOMP.EQ.2) GOTO 215                                            
      RETURN                                                            
  215 IF(FI.NE.0.d0) GOTO 230                                           
      IF(TH.NE.THEQ) GOTO 230                                           
      DRDQ=(XP-X)/DQ                                                    
      RETURN                                                            
  230 DRDQ=DRDQ+DOMRSV*DODQ                                             
      RETURN                                                            
      END                                                               
      SUBROUTINE SPOT(KOMP,N,SINTH,COSTH,SINFI,COSFI,TEMF)              
C                                                                       
c   If a surface point is in more than one spot, this subroutine        
c      adopts the product of the spot temperature factors.              
C                                                                       
c   "Latitudes" here actually run from 0 at one pole to 180 deg.        
c      at the other.                                                    
C                                                                       
c   Version of February 11, 1998                                        
C                                                                       
      implicit real*8 (a-h,o-z)                                         
      common /inprof/ in1min,in1max,in2min,in2max,mpage,nl1,nl2         
      COMMON /SPOTS/ SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG   
     $(2,100),RAD(2,100),TEMSP(2,100),xlng(2,100),kks(2,100),           
     $Lspot(2,100)                                                      
      TEMF=1.d0                                                         
      nl=(2-komp)*nl1+(komp-1)*nl2                                      
      DO 15 I=1,N                                                       
      do 42 j=1,nl                                                      
   42 if(kks(komp,j).eq.-i) Lspot(komp,j)=Lspot(komp,j)+1               
      COSDFI=COSFI*COSLNG(KOMP,I)+SINFI*SINLNG(KOMP,I)                  
      S=dacos(COSTH*COSLAT(KOMP,I)+SINTH*SINLAT(KOMP,I)*COSDFI)         
      IF(S.GT.RAD(KOMP,I)) GOTO 15                                      
      TEMF=TEMF*TEMSP(KOMP,I)                                           
      if(mpage.ne.3) goto 15                                            
      do 24 j=1,nl                                                      
      kk=kks(komp,j)                                                    
      if(kk.eq.-i) Lspot(komp,j)=0                                      
      if(kk.eq.i) Lspot(komp,j)=Lspot(komp,j)+1                         
   24 continue                                                          
   15 continue                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE OLUMP(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2, 
     $MMSAVE,GREXP,ALB,RB,WL,TPOLL,SBR,SUMM,N1,N2,KOMP,IFAT,x,y,D,      
     $SNTH,CSTH,SNFI,CSFI,tld,glump1,glump2)                            
c   Version of September 14, 1998                                       
      implicit real*8 (a-h,o-z)                                         
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),F(3),W(3),SNTH(*),CSTH(*),        
     $SNFI(*),CSFI(*),tld(*),glump1(*),glump2(*)                        
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm,   
     $ifs1dm,ifs2dm,icr1dm,icr2dm,ld,ncl,jdphs,ipc                      
      IQ=(KOMP-1)*(N1+1)                                                
      IS=(KOMP-1)*MMSAVE(IQ)                                            
      FP=7.957747d-2                                                    
      ATRAT=1.d0                                                        
      RATPOL=1.d0                                                       
      PI=3.141592653589793d0                                            
      PIH=1.570796326794897d0                                           
      PI32=4.712388980384690d0                                          
      F(1)=.1127017d0                                                   
      F(2)=.5d0                                                         
      F(3)=.8872983d0                                                   
      W(1)=.277777777777777d0                                           
      W(2)=.444444444444444d0                                           
      W(3)=.277777777777777d0                                           
      TPOLE=10000.d0*TPOLL                                              
      IF(IFAT.NE.0) CALL ATM(tpole,WL,RATPOL)                           
      CMPP=dfloat(2-KOMP)                                               
      COMPP=dfloat(2*KOMP-3)                                            
      COMP=-COMPP                                                       
      CMP=dfloat(KOMP-1)                                                
      CMPD=CMP*D                                                        
      CMPPD=CMPP*D                                                      
      N=(2-KOMP)*N1+(KOMP-1)*N2                                         
      ENN=(15.d0+X)*(1.d0+GREXP)/(15.d0-5.d0*X)                         
      NP=N1+1+(2-KOMP)*(N2+1)                                           
      NPP=N1*(KOMP-1)+(NP-1)*(2-KOMP)                                   
      LL=MMSAVE(NPP)+1                                                  
      LLL=MMSAVE(NP)                                                    
      LLLL=(LL+LLL)/2                                                   
      AR=RV(LLL)*CMP+RVQ(LLL)*CMPP                                      
      BR=RV(LLLL)*CMP+RVQ(LLLL)*CMPP                                    
      CR=RV(1)*CMP+RVQ(1)*CMPP                                          
      BOA=BR/AR                                                         
      BOAL=1.d0-BOA*BOA                                                 
      BOC2=(BR/CR)**2                                                   
      CC=1.d0/(1.d0-.25d0*ENN*(1.d0-BOA**2)*(.9675d0-.3008d0*BOA))      
      HCN=.5d0*CC*ENN                                                   
      DF=1.d0-X/3.d0                                                    
      if(ld.eq.2) df=df+2.d0*y/9.d0                                     
      if(ld.eq.3) df=df-.2d0*y                                          
      XLUMP=14384.d0/WL                                                 
      POLFAC=dexp(XLUMP/TPOLE)-1.d0                                     
      IF(KOMP.EQ.1) GOTO 82                                             
      GRPOLE=dsqrt(GRXQ(1)**2+GRYQ(1)**2+GRZQ(1)**2)                    
      GOTO 83                                                           
   82 GRPOLE=dsqrt(GRX(1)**2+GRY(1)**2+GRZ(1)**2)                       
   83 EN=dfloat(N)                                                      
      DO 8 I=1,N                                                        
      IPN1=I+N1*(KOMP-1)                                                
      SINTH=SNTH(IPN1)                                                  
      COSTH=CSTH(IPN1)                                                  
      EM=SINTH*EN*1.3d0                                                 
      MM=EM+1.d0                                                        
      IP=(KOMP-1)*(N1+1)+I                                              
      IY=MMSAVE(IP)                                                     
      DO 8 J=1,MM                                                       
      IS=IS+1                                                           
      STCF=SINTH*CSFI(IS)                                               
      STSF=SINTH*SNFI(IS)                                               
      IX=IY+J                                                           
      IF(KOMP.EQ.1) GOTO 39                                             
      IF(RVQ(IX).EQ.-1.d0) GOTO 8                                       
      GX=GRXQ(IX)                                                       
      GY=GRYQ(IX)                                                       
      GZ=GRZQ(IX)                                                       
      R=RVQ(IX)                                                         
      GOTO 49                                                           
   39 IF(RV(IX).EQ.-1.d0)GOTO 8                                         
      GX=GRX(IX)                                                        
      GY=GRY(IX)                                                        
      GZ=GRZ(IX)                                                        
      R=RV(IX)                                                          
   49 GRMAG=dsqrt(GX*GX+GY*GY+GZ*GZ)                                    
      ZZ=R*COSTH                                                        
      YY=R*COMP*STSF                                                    
      XX=CMPD+COMP*STCF*R                                               
      XXREF=(CMPPD+COMPP*XX)*COMPP                                      
      GRAV=(GRMAG/GRPOLE)**GREXP                                        
      TLOCAL=TPOLE*dsqrt(dsqrt(GRAV))                                   
      DIST=dsqrt(XXREF*XXREF+YY*YY+ZZ*ZZ)                               
      RMX=dasin(.5d0*(BR+CR)/DIST)                                      
      XCOS=XXREF/DIST                                                   
      YCOS=YY/DIST                                                      
      ZCOS=ZZ/DIST                                                      
      COSINE=(XCOS*GX+YCOS*GY+ZCOS*GZ)/GRMAG                            
      RC=PIH-dacos(COSINE)                                              
      AH=RC/RMX                                                         
      RP=dabs(AH)                                                       
      IF(AH.LE..99999d0) GOTO 22                                        
      P=1.d0                                                            
      GOTO 16                                                           
   22 IF(AH.GE.-.99999d0) GOTO 24                                       
      ALBEP=0.d0                                                        
      GOTO 19                                                           
   24 SUM=0.d0                                                          
      FIST=dasin(RP)                                                    
      FII=PIH-FIST                                                      
      DO 15 IT=1,3                                                      
      FE=FII*F(IT)+FIST                                                 
      PAR=1.d0-(RP/dsin(FE))**2                                         
      RPAR=dsqrt(PAR)                                                   
      SUM=PAR*RPAR*W(IT)+SUM                                            
   15 CONTINUE                                                          
      FTRI=(1.d0-X)*RP*dsqrt(1.d0-RP**2)+.666666666666666d0*X*FII       
     $-.666666666666667d0*x*sum*fii                                     
      FSEC=(PIH+FIST)*DF                                                
      P=(FTRI+FSEC)/(PI*DF)                                             
      IF(COSINE.LT.0.d0) P=1.d0-P                                       
      RTF=dsqrt(1.d0-AH**2)                                             
      DENO=PI32-3.d0*(AH*RTF+dasin(AH))                                 
      IF(DENO.NE.0.d0) GOTO 117                                         
      ABAR=1.d0                                                         
      GOTO 116                                                          
  117 ABAR=2.d0*RTF**3/DENO                                             
  116 COSINE=dcos(PIH-RMX*ABAR)                                         
   16 COSQ=1.d0/(1.d0+(YY/XXREF)**2)                                    
      COT2=(ZZ/XXREF)**2                                                
      Z=BOAL/(1.d0+BOC2*COT2)                                           
      E=CC-HCN*COSQ*Z                                                   
      ALBEP=ALB*E*P                                                     
   19 IF(COSINE.LE.0.d0) ALBEP=0.d0                                     
      TNEW=TLOCAL*dsqrt(dsqrt(1.d0+(FP*SUMM/(DIST*DIST*GRAV))*          
     $cosine*rb*ALBEP))                                                 
      TLD(IS)=TNEW                                                      
      IF(IFAT.EQ.0) GOTO 18                                             
      CALL ATM(TNEW,WL,RATNEW)                                          
      ATRAT=RATNEW/RATPOL                                               
   18 GRREFL=POLFAC*ATRAT/(dexp(XLUMP/TNEW)-1.d0)                       
   37 IF(KOMP.EQ.1) GOTO 77                                             
      slump2(ix)=glump2(ix)*grrefl*sbr                                  
      GOTO 8                                                            
   77 slump1(ix)=glump1(ix)*grrefl*sbr                                  
    8 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE BOLO (RATLUM,T1,T2,WL,IFAT1,IFAT2,RATBOL)              
c  Version of October 11, 1995                                          
C     BOLO USES HARRIS CALIBRATION FROM T=3970 TO 5800, MORTON-ADAMS    
C     FROM T=5800 TO 37500, AND BLACK BODY LAWS BELOW 3970 AND ABOVE    
C     37500.                                                            
      implicit real*8 (a-h,o-z)                                         
      VT1=1.d0                                                          
      VT2=1.d0                                                          
      WLT1=1.d0                                                         
      WLT2=1.d0                                                         
      TP1=10000.d0*T1                                                   
      TP2=10000.d0*T2                                                   
      IF(IFAT1.EQ.0) GOTO 95                                            
      CALL ATM(TP1,.5500d0,VT1)                                         
      CALL ATM(TP1,WL,WLT1)                                             
   95 IF(IFAT2.EQ.0) GOTO 90                                            
      CALL ATM(TP2,.5500d0,VT2)                                         
      CALL ATM(TP2,WL,WLT2)                                             
   90 RAT12=VT1*WLT2/(VT2*WLT1)                                         
      RATV=RAT12*RATLUM*(dexp(1.4384d0/(WL*T1))-1.d0)*(dexp(1.4384d0/   
     $(.5500d0*t2))-1.d0)/((dexp(1.4384d0/(.5500d0*T1))-1.d0)*          
     $(dexp(1.4384d0/(wl*t2))-1.d0))                                    
      DELMV=1.085736204758130d0*dlog(RATV)                              
      X1=dlog10(T1)+.14d0                                               
      X2=dlog10(T2)+.14d0                                               
      IF(T1.GT.3.75d0) GOTO 22                                          
      IF(T1.LT..397d0) GOTO 22                                          
      IF(T1.GT.1.07d0) GOTO 18                                          
      BC1=.0004964d0+X1*(.184549d0+X1*(-7.25068d0+X1*(-26.4889d0        
     $+x1*(-198.069d0+x1*(-73.801d0+652.95d0*X1)))))                    
      GOTO 19                                                           
   18 BC1=.868864d0+X1*(-14.8735d0+X1*(70.1318d0+X1*(-227.706d0         
     $+x1*(386.7d0+x1*(-339.447d0+121.8d0*X1)))))                       
      GOTO 19                                                           
   22 BC1=-2.5d0*dlog10((T1/.6700d0)**4*(dexp(1.4384d0/(.5450d0*T1))    
     $-1.d0)/(dexp(1.4384d0/.36515d0)-1.d0))                            
   19 IF(T2.GT.3.75d0) GOTO 32                                          
      IF(T2.LT..397d0) GOTO 32                                          
      IF(T2.GT.1.07d0) GOTO 28                                          
      BC2=.0004964d0+X2*(.184549d0+X2*(-7.25068d0+X2*(-26.4889d0        
     $+x2*(-198.069d0+x2*(-73.801d0+652.95d0*X2)))))                    
      GOTO 29                                                           
   28 BC2=.868864d0+X2*(-14.8735d0+X2*(70.1318d0+X2*(-227.706d0         
     $+x2*(386.7d0+x2*(-339.447d0+121.8d0*X2)))))                       
      GOTO 29                                                           
   32 BC2=-2.5d0*dlog10((T2/.6700d0)**4*(dexp(1.4384d0/(.5450d0*T2))    
     $-1.d0)/(dexp(1.4384d0/.36515d0)-1.d0))                            
   29 DELMB=DELMV+BC2-BC1                                               
      RATBOL=10.d0**(.4d0*DELMB)                                        
      RETURN                                                            
      END                                                               
      subroutine mlrg(a,p,q,r1,r2,t1,t2,sm1,sm2,sr1,sr2,bolm1,          
     $bolm2,xlg1,xlg2)                                                  
c  This subroutine computes absolute dimensions and other quantities    
c  for the stars of a binary star system.                               
c  a = orbital semi-major axis, the sum of the two a's for the two      
c  stars. The unit is a solar radius.                                   
c  r1,r2 = relative mean (equivalent sphere) radii for stars 1 and 2. Th
c  unit is the orbital semimajor axis.                                  
c  p = orbit period in days.                                            
c  q = mass ratio, m2/m1.                                               
c  t1,t2= flux-weighted mean surface temperatures for stars 1 and 2,in K
c  sm1,sm2= masses of stars 1 and 2 in solar units.                     
c  sr1,sr2= mean radii of stars 1 and 2 in solar units.                 
c  bolm1, bolm2= absolute bolometric magnitudes of stars 1, 2.          
c  xlg1, xlg2= log (base 10) of mean surface acceleration (effective gra
c  for stars 1 and 2.                                                   
c                                                                       
      implicit real*8 (a-h,o-z)                                         
      G=6.668d-8                                                        
      tsun=5800.d0                                                      
      rsunau=214.8d0                                                    
      sunmas=1.991d33                                                   
      sunrad=6.960d10                                                   
      gmr=G*sunmas/sunrad**2                                            
      sunmb=4.77d0                                                      
      sr1=r1*a                                                          
      sr2=r2*a                                                          
      yrsid=365.2564d0                                                  
      tmass=(a/rsunau)**3/(p/yrsid)**2                                  
      sm1=tmass/(1.d0+q)                                                
      sm2=tmass*q/(1.d0+q)                                              
      bol1=(t1/tsun)**4*sr1**2                                          
      bol2=(t2/tsun)**4*sr2**2                                          
      bolm1=sunmb-2.5d0*dlog10(bol1)                                    
      bolm2=sunmb-2.5d0*dlog10(bol2)                                    
      xlg1=dlog10(gmr*sm1/sr1**2)                                       
      xlg2=dlog10(gmr*sm2/sr2**2)                                       
      return                                                            
      end                                                               
      subroutine cloud(cosa,cosb,cosg,x1,y1,z1,xc,yc,zc,rr,wl,op1,      
     $opsf,edens,acm,en,cmpd,ri,dx,dens,tau)                            
c  Version of September 11, 1998                                        
      implicit real*8 (a-h,o-z)                                         
      dx=0.d0                                                           
      tau=0.d0                                                          
      ri=1.d0                                                           
      sige=.6653d-24                                                    
      dtdxes=sige*edens                                                 
c  cosa can be zero, so an alternate path to the solution is needed     
      dabcoa=dabs(cosa)                                                 
      dabcob=dabs(cosb)                                                 
      if(dabcoa.lt.dabcob) goto 32                                      
      w=cosb/cosa                                                       
      v=cosg/cosa                                                       
      u=y1-yc-w*x1                                                      
      t=z1-zc-v*x1                                                      
      aa=1.d0+w*w+v*v                                                   
      bb=2.d0*(w*u+v*t-xc)                                              
      cc=xc*xc+u*u+t*t-rr*rr                                            
      dubaa=aa+aa                                                       
      dis=bb*bb-4.d0*aa*cc                                              
      if(dis.le.0.d0) return                                            
      sqd=dsqrt(dis)                                                    
      xx1=(-bb+sqd)/dubaa                                               
      xx2=(-bb-sqd)/dubaa                                               
      yy1=w*(xx1-x1)+y1                                                 
      yy2=w*(xx2-x1)+y1                                                 
      zz1=v*(xx1-x1)+z1                                                 
      zz2=v*(xx2-x1)+z1                                                 
      goto 39                                                           
   32 w=cosa/cosb                                                       
      v=cosg/cosb                                                       
      u=x1-xc-w*y1                                                      
      t=z1-zc-v*y1                                                      
      aa=1.d0+w*w+v*v                                                   
      bb=2.d0*(w*u+v*t-yc)                                              
      cc=yc*yc+u*u+t*t-rr*rr                                            
      dubaa=aa+aa                                                       
      dis=bb*bb-4.d0*aa*cc                                              
      if(dis.le.0.d0) return                                            
      sqd=dsqrt(dis)                                                    
      yy1=(-bb+sqd)/dubaa                                               
      yy2=(-bb-sqd)/dubaa                                               
      xx1=w*(yy1-y1)+x1                                                 
      xx2=w*(yy2-y1)+x1                                                 
      zz1=v*(yy1-y1)+z1                                                 
      zz2=v*(yy2-y1)+z1                                                 
   39 dis=bb*bb-4.d0*aa*cc                                              
      if(dis.le.0.d0) return                                            
      sqd=dsqrt(dis)                                                    
c     goto 42                                                           
   42 xs1=(xx1-cmpd)*cosa+yy1*cosb+zz1*cosg                             
      xs2=(xx2-cmpd)*cosa+yy2*cosb+zz2*cosg                             
      xxnear=xx1                                                        
      yynear=yy1                                                        
      zznear=zz1                                                        
      xxfar=xx2                                                         
      yyfar=yy2                                                         
      zzfar=zz2                                                         
      xsnear=xs1                                                        
      xsfar=xs2                                                         
      if(xs1.gt.xs2) goto 38                                            
      xxnear=xx2                                                        
      yynear=yy2                                                        
      zznear=zz2                                                        
      xxfar=xx1                                                         
      yyfar=yy1                                                         
      zzfar=zz1                                                         
      xsnear=xs2                                                        
      xsfar=xs1                                                         
   38 continue                                                          
      xss=(x1-cmpd)*cosa+y1*cosb+z1*cosg                                
      if(xss.ge.xsnear) return                                          
      if(xss.le.xsfar) goto 20                                          
      xxfar=x1                                                          
      yyfar=y1                                                          
      zzfar=z1                                                          
   20 continue                                                          
      dtaudx=dtdxes+(op1*wl**en+opsf)*dens                              
      dx=dsqrt((xxnear-xxfar)**2+(yynear-yyfar)**2+(zznear-zzfar)**2)   
      tau=dx*dtaudx*acm                                                 
      ri=dexp(-tau)                                                     
      return                                                            
      end                                                               
      SUBROUTINE xxsqar(OBS,NOBS,ML,OUT,sd,xlamda,D,CN,CNN,cnc,         
     $clc,ss,cl,ll,mm)                                                  
c  Version of July 9, 1997                                              
      implicit real*8 (a-h,o-z)                                         
      dimension obs(*),out(*),sd(*),cn(*),cnn(*),cnc(*),clc(*),         
     $cl(*),ll(*),mm(*)                                                 
c                                                                       
c  cnc ("cn copy") is the original normal equation matrix               
c  cn is the re-scaled version of cnc                                   
c  cnn comes in as the original n.e. matrix, then is copied             
c    from cn to become the re-scaled n.e.'s, and finally is             
c    inverted by DMINV to become the inverse of the re-scaled           
c    n.e. matrix.                                                       
c                                                                       
      S=0.D0                                                            
      CLL=0.D0                                                          
      CAY=NOBS-ML                                                       
      JMAX=ML*ML                                                        
      DO 20 J=1,JMAX                                                    
   20 CN(J)=0.D0                                                        
      DO 21 J=1,ML                                                      
   21 CL(J)=0.D0                                                        
      DO 25 NOB=1,NOBS                                                  
      III=NOB+NOBS*ML                                                   
      OBSQQ=OBS(III)                                                    
      DO 23 K=1,ML                                                      
      DO 23 I=1,ML                                                      
      II=NOB+NOBS*(I-1)                                                 
      KK=NOB+NOBS*(K-1)                                                 
      J=I+(K-1)*ML                                                      
      OBSII=OBS(II)                                                     
      OBSKK=OBS(KK)                                                     
      CN(J)=CN(J)+OBSII*OBSKK                                           
   23 cnc(j)=cn(j)                                                      
      DO 24 I=1,ML                                                      
      II=NOB+NOBS*(I-1)                                                 
      OBSII=OBS(II)                                                     
   24 CL(I)=CL(I)+OBSQQ*OBSII                                           
   25 CLL=CLL+OBSQQ*OBSQQ                                               
      do 123 k=1,ml                                                     
      do 123 i=1,ml                                                     
      xlf=0.d0                                                          
      if(i.eq.k) xlf=xlamda                                             
      j=i+(k-1)*ml                                                      
      ji=i+(i-1)*ml                                                     
      ki=k+(k-1)*ml                                                     
  123 cn(j)=cn(j)/dsqrt(cnc(ji)*cnc(ki))+xlf                            
      do 124 i=1,ml                                                     
      ji=i+(i-1)*ml                                                     
      clc(i)=cl(i)                                                      
  124 cl(i)=cl(i)/dsqrt(cnc(ji))                                        
      DO 50 J=1,JMAX                                                    
   50 CNN(J)=CN(J)                                                      
      CALL DMINV(CNN,ML,D,LL,MM)                                        
      CALL DGMPRD(CNN,CL,OUT,ML,ML,1)                                   
      do 125 i=1,ml                                                     
      ji=i+(i-1)*ml                                                     
  125 out(i)=out(i)/dsqrt(cnc(ji))                                      
      DO 26 I=1,ML                                                      
   26 S=S+clc(I)*OUT(I)                                                 
      S=CLL-S                                                           
      SS=S                                                              
      SIGSQ=S/CAY                                                       
      DO 27 J=1,ML                                                      
      JJ=J*ML+J-ML                                                      
      CNJJ=CNN(JJ)                                                      
      ARG=SIGSQ*CNJJ                                                    
   27 sd(J)=dsqrt(arg/cnc(jj))                                          
      RETURN                                                            
      END                                                               
      subroutine linpro(komp,dvks,hbarw,tau,count,taug,fbin,delv)       
c  Version of July 14, 1997                                             
      implicit real*8(a-h,o-z)                                          
      dimension dvks(*),hbarw(*),tau(*),count(*),fbin(*),delv(*),taug(*)
      common /flpro/ vks,binc,binw,difp,dum1,dum2                       
      common /ipro/ nbins,nl,inmax,inmin,idum1,idum2                    
      COMMON /SPOTS/ SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG   
     $(2,100),RAD(2,100),TEMSP(2,100),xlng(2,100),kks(2,100),           
     $Lspot(2,100)                                                      
      inmin=30000                                                       
      inmax=0                                                           
      do 83 iln=1,nl                                                    
      if(Lspot(komp,iln).eq.0) goto 83                                  
      vksg=vks+dvks(iln)                                                
      vksp=vksg+hbarw(iln)                                              
      vksm=vksg-hbarw(iln)                                              
      indp=vksp/binw+binc                                               
      indm=vksm/binw+binc                                               
      if(indm.lt.inmin) inmin=indm                                      
      if(indp.gt.inmax) inmax=indp                                      
   83 continue                                                          
      do 82 i=inmin,inmax                                               
   82 taug(i)=0.d0                                                      
      do 84 iln=1,nl                                                    
      if(Lspot(komp,iln).eq.0) goto 84                                  
      vksg=vks+dvks(iln)                                                
      vksp=vksg+hbarw(iln)                                              
      vksm=vksg-hbarw(iln)                                              
      indp=vksp/binw+binc                                               
      indm=vksm/binw+binc                                               
      vks1=(dfloat(indm+1)-binc)*binw                                   
      vks2=(dfloat(indp)-binc)*binw                                     
      fr1=(vks1-vksm)/binw                                              
      fr2=(vksp-vks2)/binw                                              
      taug(indm)=taug(indm)+fr1*tau(iln)                                
      delv(indm)=delv(indm)+fr1*vksm                                    
      count(indm)=count(indm)+fr1                                       
      taug(indp)=taug(indp)+fr2*tau(iln)                                
      delv(indp)=delv(indp)+fr2*vksp                                    
      count(indp)=count(indp)+fr2                                       
      if(indp.ne.indm) goto 28                                          
      taug(indp)=taug(indp)-tau(iln)                                    
      delv(indp)=delv(indp)-.5d0*(vksm+vksp)                            
      count(indp)=count(indp)-1.d0                                      
   28 continue                                                          
      ind=indm                                                          
      idmax=indp-indm-1                                                 
      if(idmax.le.0) goto 84                                            
      do 26 id=1,idmax                                                  
      ind=ind+1                                                         
      vksb=(dfloat(ind)-binc)*binw                                      
      taug(ind)=taug(ind)+tau(iln)                                      
      delv(ind)=delv(ind)+vksb                                          
      count(ind)=count(ind)+1.d0                                        
   26 continue                                                          
   84 continue                                                          
      do 85 i=inmin,inmax                                               
   85 fbin(i)=fbin(i)+(1.d0-dexp(-taug(i)))*difp                        
      return                                                            
      end                                                               
      subroutine jdph(xjdin,phin,t0,p0,dpdt,xjdout,phout)               
c  Version of February 2, 1999                                          
c                                                                       
c  Subroutine jdph computes a phase (phout) based on an input           
c   JD (xjdin), reference epoch (t0), period (p0), and dP/dt (dpdt).    
c   It also computes a JD (xjdout) from an input phase (phin) and the   
c   same ephemeris. So jdph can be used either to get phase from        
c   JD or JD from phase.                                                
c                                                                       
      implicit real*8(a-h,o-z)                                          
      tol=1.d-6                                                         
      abdpdt=dabs(dpdt)                                                 
      deltop=(xjdin-t0)/p0                                              
      fcsq=0.d0                                                         
      fccb=0.d0                                                         
      fc4th=0.d0                                                        
      fc5th=0.d0                                                        
      fc=deltop*dpdt                                                    
      if(dabs(fc).lt.1.d-18) goto 25                                    
      fcsq=fc*fc                                                        
      if(dabs(fcsq).lt.1.d-24) goto 25                                  
      fccb=fc*fcsq                                                      
      if(dabs(fccb).lt.1.d-27) goto 25                                  
      fc4th=fc*fccb                                                     
      if(dabs(fc4th).lt.1.d-28) goto 25                                 
      fc5th=fc*fc4th                                                    
   25 phout=deltop*(1.d0-.5d0*fc+fcsq/3.d0-.25d0*fccb+.2d0*fc4th-       
     $fc5th/6.d0)                                                       
      pddph=dpdt*phin                                                   
      xjdout=p0*phin*(1.d0+.5d0*pddph+pddph**2/6.d0+pddph**3/24.d0      
     $+pddph**4/120.d0+pddph**5/720.d0)+t0                              
      if(abdpdt.lt.tol) return                                          
      phout=dlog(1.d0+deltop*dpdt)/dpdt                                 
      xjdout=(dexp(dpdt*phin)-1.d0)*p0/dpdt+t0                          
      return                                                            
      end                                                               
      subroutine ranuni(sn,smod,sm1p1)                                  
c  Version of February 6, 1997                                          
c                                                                       
c   On each call, subroutine ranuni generates a pseudo-random number,   
c     sm1p1, distributed with uniform probability over the range        
c     -1. to +1.                                                        
c                                                                       
c   The input number sn, from which both output numbers are generated,  
c     should be larger than the modulus 100000008. and smaller          
c     than twice the modulus. The returned number smod will be in       
c     that range and can be used as the input sn on the next call       
c                                                                       
      implicit real*8(a-h,o-z)                                          
      st=23.d0                                                          
      xmod=1.00000001d8                                                 
      smod=st*sn                                                        
      goto 2                                                            
    1 smod=smod-xmod                                                    
    2 if(smod.gt.xmod) goto 1                                           
      sm1p1=(2.d0*smod/xmod-1.d0)                                       
      return                                                            
      end                                                               
      subroutine rangau(smod,nn,sd,gau)                                 
      implicit real*8(a-h,o-z)                                          
c  Version of February 6, 1997                                          
      ffac=0.961d0                                                      
      sfac=ffac*3.d0*sd/(dsqrt(3.d0*dfloat(nn)))                        
      g1=0.d0                                                           
      do 22 i=1,nn                                                      
      sn=smod                                                           
      call ranuni(sn,smod,sm1p1)                                        
      g1=g1+sm1p1                                                       
   22 continue                                                          
      gau=sfac*g1                                                       
      return                                                            
      end                                                               
