      SUBROUTINE SINCOS (KOMP,N,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)             19410
C     VERSION OF March 10, 1992. DESTROY EARLIER VERSIONS.                 19420
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*),MMSAVE(*)                  19430
      IP=(KOMP-1)*(N1+1)+1                                                 19440
      IQ=IP-1                                                              19450
      IS=0                                                                 19460
      IF(KOMP.EQ.2) IS=MMSAVE(IQ)                                          19470
      MMSAVE(IP)=0                                                         19480
      EN=N                                                                 19490
      DO 8 I=1,N                                                           19500
      EYE=I                                                                19510
      EYE=EYE-.5                                                           19520
      TH=1.570796*EYE/EN                                                   19530
      IPN1=I+N1*(KOMP-1)                                                   19540
      SNTH(IPN1)=SIN(TH)                                                   19550
      CSTH(IPN1)=COS(TH)                                                   19560
      EM=SNTH(IPN1)*EN*1.3                                                 19570
      MM=EM+1.                                                             19580
      XM=MM                                                                19590
      IP=(KOMP-1)*(N1+1)+I+1                                               19600
      IQ=IP-1                                                              19610
      MMSAVE(IP)=MMSAVE(IQ)+MM                                             19620
      DO 8 J=1,MM                                                          19630
      IS=IS+1                                                              19640
      XJ=J                                                                 19650
      FI=3.141593*(XJ-.5)/XM                                               19660
      CSFI(IS)=COS(FI)                                                     19670
      SNFI(IS)=SIN(FI)                                                     19680
    8 CONTINUE                                                             19690
      RETURN                                                               19700
      END                                                                  19710
      SUBROUTINE SURFAS(RMASS,POTENT,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,         19720
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,FF,D,SNTH,CSTH,SNFI,CSFI,GRV1,     19730
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,        19740
     $GMAG2,GREXP)                                                         19750
c  Experimental version of April 9, 1992                                   19760
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*   19770
     $),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*)     19780
     $,GRV1(*),GRV2(*),XX1(*),YY1(*),ZZ1(*),XX2(*),YY2(*),ZZ2(*),GLUMP1    19790
     $(*),GLUMP2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                    19800
      common /radi/ R1H,RLH,R1C,RLC                                        19810
      COMMON X1                                                            19820
      COMMON /ECCEN/E                                                      21820
      DSQ=D*D                                                              19840
      RMAS=RMASS                                                           19850
      IF(KOMP.EQ.2) RMAS=1./RMASS                                          19860
      RF=FF**2                                                             19870
      RTEST=0.                                                             19880
      IP=(KOMP-1)*(N1+1)+1                                                 19890
      IQ=IP-1                                                              19900
      IS=0                                                                 19910
      ISX=(KOMP-1)*MMSAVE(IQ)                                              19920
      MMSAVE(IP)=0                                                         19930
      KFLAG=0                                                              19940
      CALL ELLONE (FF,D,RMAS,X1,OMEGA,XL2,OM2)                             19950
      IF(KOMP.EQ.2) OMEGA=RMASS*OMEGA+.5*(1.-RMASS)                        19960
      X2=X1                                                                19970
      IF(KOMP.EQ.2) X1=1.-X1                                               19980
      IF(E.NE.0.) GOTO 43                                                  19990
      IF(POTENT.LT.OMEGA) CALL NEKMIN(RMASS,POTENT,X1,ZZ)                  20000
      IF(POTENT.LT.OMEGA) X2=1.-X1                                         20010
   43 COMP=3-2*KOMP                                                        20020
      CMP=KOMP-1                                                           20030
      CMPD=CMP*D                                                           20040
      TESTER=CMPD+COMP*X1                                                  20050
      RM1=RMASS+1.                                                         20060
      RMS=RMASS                                                            20070
      RM1S=RM1                                                             20080
      IF(KOMP.NE.2) GOTO 15                                                20090
      POT=POTENT/RMASS+.5*(RMASS-1.)/RMASS                                 20100
      RM=1./RMASS                                                          20110
      RM1=RM+1.                                                            20120
      GOTO 20                                                              20130
   15 POT=POTENT                                                           20140
      RM=RMASS                                                             20150
   20 EN=N                                                                 20160
      RSAVE=1./POT                                                         20170
      DO 8 I=1,N                                                           20180
      IF(I.NE.2) GOTO 82                                                   20190
      IF(KOMP.EQ.1) RTEST=.3*RV(1)                                         20200
      IF(KOMP.EQ.2) RTEST=.3*RVQ(1)                                        20210
   82 CONTINUE                                                             20220
      IPN1=I+N1*(KOMP-1)                                                   20230
      SINTH=SNTH(IPN1)                                                     20240
      XNU=CSTH(IPN1)                                                       20250
      XNUSQ=XNU**2                                                         20260
      EM=SINTH*EN*1.3                                                      20270
      XLUMP=1.-XNUSQ                                                       20280
      MM=EM+1.                                                             20290
      DO 8 J=1,MM                                                          20300
      KOUNT=0                                                              20310
      IS=IS+1                                                              20320
      ISX=ISX+1                                                            20330
      DELR=0.                                                              20340
      COSFI=CSFI(ISX)                                                      20350
      XMU=SNFI(ISX)*SINTH                                                  20360
      XLAM=SINTH*COSFI                                                     20370
      R=RSAVE                                                              20380
   14 R=R+DELR                                                             20390
      KOUNT=KOUNT+1                                                        20400
      IF(KOUNT.LT.20) GOTO 70                                              20410
      KFLAG=1                                                              20420
      R=-1.                                                                20430
      GOTO 86                                                              20440
   70 RSQ=R*R                                                              20450
      PAR=DSQ-2.*XLAM*R*D+RSQ                                              20460
      RPAR=SQRT(PAR)                                                       20470
      OM=1./R+RM*((1./RPAR)-XLAM*R/DSQ)+RM1*.5*RSQ*XLUMP*RF                20480
      DOMR=1./(RF*RM1*XLUMP*R-1./RSQ-(RM*(R-XLAM*D))/(PAR*RPAR)-RM*XLAM/   20490
     $DSQ)                                                                 20500
      DELR=(POT-OM)*DOMR                                                   20510
      ABDELR=ABS(DELR)                                                     20520
      IF(ABDELR.GT..00001) GOTO 14                                         20530
      ABR=ABS(R)                                                           20540
      IF(R.GT.RTEST) GOTO 74                                               20550
      KFLAG=1                                                              20560
      R=-1.                                                                20570
      IF(KOMP.EQ.2) GOTO 98                                                20580
      GOTO 97                                                              20590
   74 IF(ABR.LT.TESTER) RSAVE=R                                            20600
      Z=R*XNU                                                              20610
      Y=COMP*R*XMU                                                         20620
      X2T=ABR*XLAM                                                         20630
      X=CMPD+COMP*X2T                                                      20640
      IF(KOMP.EQ.2) GOTO 62                                                20650
      IF(X.LT.X1) GOTO 65                                                  20660
      KFLAG=1                                                              20670
      R=-1.                                                                20680
      GOTO 97                                                              20690
   62 IF(X2T.LT.X2) GOTO 65                                                20700
      KFLAG=1                                                              20710
      R=-1.                                                                20720
      GOTO 98                                                              20730
   65 SUMSQ=Y**2+Z**2                                                      20740
      PAR1=X**2+SUMSQ                                                      20750
      RPAR1=SQRT(PAR1)                                                     20760
      XNUM1=1./(PAR1*RPAR1)                                                20770
      XL=D-X                                                               20780
      PAR2=XL**2+SUMSQ                                                     20790
      RPAR2=SQRT(PAR2)                                                     20800
      XNUM2=1./(PAR2*RPAR2)                                                20810
      OMZ=-Z*(XNUM1+RMS*XNUM2)                                             20820
      OMY=Y*(RM1S*RF-XNUM1-RMS*XNUM2)                                      20830
      OMX=RMS*XL*XNUM2-X*XNUM1+RM1S*X*RF-RMS/DSQ                           20840
      IF(KOMP.EQ.2)OMX=RMS*XL*XNUM2-X*XNUM1-RM1S*XL*RF+1./DSQ              20850
      GRMAG=SQRT(OMX*OMX+OMY*OMY+OMZ*OMZ)                                  20860
      IF(IS.EQ.1) GRPOLE=GRMAG                                             20870
      GRAV=(GRMAG/GRPOLE)**GREXP                                           20880
      A=COMP*XLAM*OMX                                                      20890
      B=COMP*XMU*OMY                                                       20900
      C=XNU*OMZ                                                            20910
      COSBET=-(A+B+C)/GRMAG                                                20920
      IF(COSBET.LT..7) COSBET=.7                                           20930
   86 IF(KOMP.EQ.2) GOTO 98                                                20940
   97 RV(IS)=R                                                             20950
      GRX(IS)=OMX                                                          20960
      GRY(IS)=OMY                                                          20970
      GRZ(IS)=OMZ                                                          20980
      GMAG1(IS)=SQRT(OMX*OMX+OMY*OMY+OMZ*OMZ)                              20990
      FR1(IS)=1.                                                           21000
      GLUMP1(IS)=R*R*SINTH/COSBET                                          21010
      GRV1(IS)=GRAV                                                        21020
      XX1(IS)=X                                                            21030
      YY1(IS)=Y                                                            21040
      ZZ1(IS)=Z                                                            21050
      CSBT1(IS)=COSBET                                                     21060
      GOTO 8                                                               21070
   98 RVQ(IS)=R                                                            21080
      GRXQ(IS)=OMX                                                         21090
      GRYQ(IS)=OMY                                                         21100
      GRZQ(IS)=OMZ                                                         21110
      GMAG2(IS)=SQRT(OMX*OMX+OMY*OMY+OMZ*OMZ)                              21120
      FR2(IS)=1.                                                           21130
      GLUMP2(IS)=R*R*SINTH/COSBET                                          21140
      GRV2(IS)=GRAV                                                        21150
      XX2(IS)=X                                                            21160
      YY2(IS)=Y                                                            21170
      ZZ2(IS)=Z                                                            21180
      CSBT2(IS)=COSBET                                                     21190
    8 CONTINUE                                                             21200
      IF(KFLAG.EQ.0) GOTO 53                                               21210
      ISS=IS-1                                                             21220
      IF(KOMP.NE.1) GOTO 50                                                21230
      CALL RING(RMASS,POTENT,1,N,FR1,HLD,R1H,RLH)                          21240
      DO 55 I=1,ISS                                                        21250
      IPL=I+1                                                              21260
      IF(RV(I).GE.0.)GOTO 55                                               21270
      FR1(IPL)=FR1(IPL)+FR1(I)                                             21280
      FR1(I)=0.                                                            21290
   55 CONTINUE                                                             21300
   53 IF(KOMP.EQ.2) GOTO 54                                                21310
      IS=0                                                                 21320
      DO 208 I=1,N                                                         21330
      IPN1=I+N1*(KOMP-1)                                                   21340
      EM=SNTH(IPN1)*EN*1.3                                                 21350
      MM=EM+1.                                                             21360
      DO 208 J=1,MM                                                        21370
      IS=IS+1                                                              21380
      GLUMP1(IS)=FR1(IS)*GLUMP1(IS)                                        21390
  208 CONTINUE                                                             21400
      RETURN                                                               21410
   50 CALL RING(RMASS,POTENT,2,N,FR2,HLD,R1C,RLC)                          21420
      DO 56 I=1,IS                                                         21430
      IPL=I+1                                                              21440
      IF(RVQ(I).GE.0.) GOTO 56                                             21450
      FR2(IPL)=FR2(IPL)+FR2(I)                                             21460
      FR2(I)=0.                                                            21470
   56 CONTINUE                                                             21480
   54 CONTINUE                                                             21490
      IS=0                                                                 21500
      DO 108 I=1,N                                                         21510
      IPN1=I+N1*(KOMP-1)                                                   21520
      EM=SNTH(IPN1)*EN*1.3                                                 21530
      MM=EM+1.                                                             21540
      DO 108 J=1,MM                                                        21550
      IS=IS+1                                                              21560
      GLUMP2(IS)=FR2(IS)*GLUMP2(IS)                                        21570
  108 CONTINUE                                                             21580
      RETURN                                                               21590
      END                                                                  21600
      SUBROUTINE BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,     21610
     $HLD,SLUMP1,SLUMP2,THETA,RHO,AA,BB,PHSV,PCSV,N1,N2,F1,F2,d,hlum,      21620
     $clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1      21630
     $,ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phas,rm,          21640
     $xincl,hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,        21650
     $yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,          21660
     $csbt2,gmag1,gmag2,mode)                                              21670
c  Experimental version of March 31, 1992                                  21680
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),         21690
     $GRZQ(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),SLUMP1(*),SLUMP2(*),          21700
     $THETA(*),RHO(*),AA(*),BB(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*),TLD(*)   21710
     $,GLUMP1(*),GLUMP2(*),XX1(*),XX2(*),YY1(*),YY2(*),ZZ1(*),ZZ2(*),      21720
     $GRV1(*),GRV2(*),RFTEMP(*),RF1(*),RF2(*),CSBT1(*),CSBT2(*)            21730
     $,GMAG1(*),GMAG2(*)                                                   21740
      COMMON /INVAR/ KH,IPBDUM,IRTE,NREF,IRVOL1,irvol2,mref,ifsmv1,        21750
     $ifsmv2,icor1,icor2,ld                                                21760
      COMMON /FLVAR/ PER,PSHIFT,DP,DS,EF,EFC,ECOS,PERR,PHPER,PCONJ,        21770
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc                21780
      common /nspt/ nsp1,nsp2,ifdum1,ifdum2                                21790
      common /spots/ snlat(2,100),cslat(2,100),snlng(2,100),               21800
     $cslng(2,100),rdsp(2,100),tmsp(2,100),xlng(2,100)                     21810
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC        21820
      pi=3.141593                                                          21830
      twopi=pi+pi                                                          21840
      MOD=(MODE-2)**2                                                      21850
      IF(MOD.EQ.1) XC=XH                                                   21860
      if(mod.eq.1) yc=yh                                                   21870
      PSFT=PHAS-PHPERI                                                     21880
   29 if(PSFT.GT.1.) PSFT=PSFT-1.                                          21890
      if(psft.gt.1.) goto 29                                               21900
   30 if(PSFT.LT.0.) PSFT=PSFT+1.                                          21910
      if(psft.lt.0.) goto 30                                               21920
      XMEAN=PSFT*6.283185                                                  21930
      tr=xmean                                                             21940
      do 60 kp=1,2                                                         21950
      nsp=nsp1*(2-kp)+nsp2*(kp-1)                                          21960
      ff=f1*(2-kp)+f2*(kp-1)                                               21970
      ifsmv=ifsmv1*(2-kp)+ifsmv2*(kp-1)                                    21980
      if(ifsmv.eq.0) goto 60                                               21990
      do 61 i=1,nsp                                                        22000
      xlg=xlng(kp,i)+twopi*ff*(phas-pconj)-(tr-trc)                        22010
      snlng(kp,i)=sin(xlg)                                                 22020
      cslng(kp,i)=cos(xlg)                                                 22030
   61 continue                                                             22040
   60 continue                                                             22050
      if(e.ne.0.) call KEPLER(XMEAN,E,DUM,TR)                              22060
      U=TR+PERR                                                            22070
      COSU=COS(U)                                                          22080
      GPHA=U*.1591549-.25                                                  22090
   40 if(GPHA.lt.0.) GPHA=GPHA+1.                                          22100
      if(gpha.lt.0.) goto 40                                               22110
   50 if(GPHA.GE.1.) GPHA=GPHA-1.                                          22120
      if(gpha.ge.1.) goto 50                                               22130
      D=EF/(1.+E*COS(TR))                                                  22140
      IF(IRTE.EQ.1) GOTO 19                                                22150
      CALL LCR(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SLUM   22160
     $P1,SLUMP2,RM,PHSV,PCSV,N1,N2,F1,F2,D,HLUM,CLUM,xh,xc,yh,yc,gr1,gr2   22170
     $,wl,SM1,SM2,TPOLH,TPOLC,SBRH,SBRC,IFAT1,IFAT2,TAVH,TAVC,alb1,alb2,   22180
     $xbol1,xbol2,ybol1,ybol2,vol1,vol2,snth,csth,snfi,csfi,tld,glump1,    22190
     $glump2,xx1,xx2,yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,csbt1,csbt2,    22200
     $rftemp,rf1,rf2,gmag1,gmag2,mode,ifc)                                 22210
   19 CONTINUE                                                             22220
      VO1=RM*SINI*(ECOS+COSU)/((1.+RM)*EFC)+VGAM                           22230
      VO2=-SINI*(ECOS+COSU)/((1.+RM)*EFC)+VGAM                             22240
      call light(gpha,xincl,xh,xc,yh,yc,n1,n2,hot,cool,rv,grx,gry,grz,     22250
     $rvq,grxq,gryq,grzq,mmsave,theta,rho,aa,bb,slump1,slump2,somhot,      22260
     $somkul,d,wl,snth,csth,snfi,csfi,tld,gmag1,gmag2)                     22270
      VRA1=0.                                                              22280
      VRA2=0.                                                              22290
      IF(HOT.GT.0.) VRA1=F1*SOMHOT/HOT                                     22300
      IF(COOL.GT.0.) VRA2=F2*SOMKUL/COOL                                   22310
      vsum1=vo1                                                            22320
      vsum2=vo2                                                            22330
      if(icor1.eq.1) vsum1=vo1+vra1                                        22340
      if(icor2.eq.1) vsum2=vo2+vra2                                        22350
      VKM1=VSUM1*VFVU                                                      22360
      VKM2=VSUM2*VFVU                                                      22370
      RETURN                                                               22380
      END                                                                  22390
      SUBROUTINE DURA(F,XINCL,RM,DD,THE,OMEG,R)                            22400
C     VERSION OF 4/10/78. DESTROY EARLIER VERSIONS.                        22410
C                                                                          22420
C     PARAMETER 'THE' IS THE SEMI-DURATION OF X-RAY ECLIPSE, AND SHOULD    22430
C     BE IN CIRCULAR MEASURE.                                              22440
      IMPLICIT REAL*8(A-H,O-Z)                                             22450
      REAL F,XINCL,RM,THE,OMEG,R,DD                                        22460
      DELX=0.D0                                                            22470
      FSQ=F*F                                                              22480
      RMD=1./RM                                                            22490
      RMD1=RMD+1.D0                                                        22500
      XINC=.0174533*XINCL                                                  22510
      TH=6.283185*THE                                                      22520
      CI=DCOS(XINC)                                                        22530
      SI=DSIN(XINC)                                                        22540
      D=DD                                                                 22550
      DSQ=D*D                                                              22560
      ST=DSIN(TH)                                                          22570
      CT=DCOS(TH)                                                          22580
      COTI=CI/SI                                                           22590
      TT=ST/CT                                                             22600
      C1=CT*SI                                                             22610
      C2=TT*ST*SI                                                          22620
      C3=C1+C2                                                             22630
      C4=COTI*CI/CT                                                        22640
      C5=C3+C4                                                             22650
      C6=C2+C4                                                             22660
      C7=(ST*ST+COTI*COTI)/CT**2                                           22670
      X=D*(SI*SI*ST*ST+CI*CI)+.00001D0                                     22680
   15 X=X+DELX                                                             22690
      PAR=X*X+C7*(D-X)**2                                                  22700
      RPAR=DSQRT(PAR)                                                      22710
      PAR32=PAR*RPAR                                                       22720
      PAR52=PAR*PAR32                                                      22730
      FC=(C6*D-C5*X)/PAR32+C1**3*C5*RMD/(D-X)**2+C3*FSQ*RMD1*X-C2*FSQ*D*   22740
     $RMD1-C1*RMD/DSQ                                                      22750
      DFCDX=(-C5*PAR-3.D0*(C6*D-C5*X)*((1.D0+C7)*X-C7*D))/PAR52+2.D0*C1    22760
     $**3*C5*RMD/(D-X)**3+C3*FSQ*RMD1                                      22770
      DELX=-FC/DFCDX                                                       22780
      ABDELX=DABS(DELX)                                                    22790
      IF(ABDELX.GT..00001D0) GOTO 15                                       22800
      Y=-(D-X)*TT                                                          22810
      Z=-(D-X)*COTI/CT                                                     22820
      YZ2=Y*Y+Z*Z                                                          22830
      OMEG=1.D0/DSQRT(X*X+YZ2)+RMD/DSQRT((D-X)**2+YZ2)+.5D0*RMD1*FSQ*      22840
     $(X*X+Y*Y)-RMD*X/DSQ                                                  22850
      OMEG=RM*OMEG+.5*(1.-RM)                                              22860
      R=DSQRT(X*X+YZ2)                                                     22870
      RETURN                                                               22880
      END                                                                  22890
      SUBROUTINE LCR(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HL   22900
     $D,SLUMP1,SLUMP2,RM,POTH,POTC,N1,N2,F1,F2,D,HLUM,CLUM,xh,xc,yh,yc,    22910
     $GR1,GR2,WL,SM1,SM2,TPOLH,TPOLC,SBRH,SBRC,IFAT1,IFAT2,TAVH,TAVC,      22920
     $alb1,alb2,xbol1,xbol2,ybol1,ybol2,vol1,vol2,snth,csth,snfi,csfi,     22930
     $tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,dint1,dint2,grv1,grv2,     22940
     $csbt1,csbt2,rftemp,rf1,rf2,gmag1,gmag2,mode,ifc)                     22950
c  Experimental version of March 31, 1992                                  22960
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*   22970
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),        22980
     $CSTH(*),SNFI(*),CSFI(*),TLD(*),GLUMP1(*),GLUMP2(*),XX1(*),XX2(*)     22990
     $,YY1(*),YY2(*),ZZ1(*),ZZ2(*),GRV1(*),GRV2(*),RFTEMP(*),RF1(*),       23000
     $RF2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                           23010
      COMMON /DPDX/ DPDX1,DPDX2,PHSV,PCSV                                  23020
      COMMON /ECCEN/ E                                                     21820
      COMMON /SUMM/ SUMM1,SUMM2                                            23040
      COMMON /INVAR/ KHDUM,IPB,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,        23050
     $ifsmv2,icor1,icor2,ld                                                23060
      XLUMP=1.4384/WL                                                      23070
      VL1=VOL1                                                             23080
      VL2=VOL2                                                             23090
      DP=1.-E                                                              23100
      IF(IRVOL1.EQ.1) GOTO 88                                              23110
      CALL VOLUME(VL1,RM,POTH,DP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,GRXQ,       23120
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM1,SM1,GRV1,     23130
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2   23140
     $,GR1,1)                                                              23150
      IF(E.EQ.0.) GOTO 88                                                  23160
      POTHD=PHSV                                                           23170
      IF(IFC.EQ.2) POTHD=PHSV+DPDX1*(1./D-1./(1.-E))                       23180
      CALL VOLUME(VL1,RM,POTHD,D,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,GRXQ,       23190
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM1,SM1,GRV1,     23200
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2   23210
     $,GR1,IFC)                                                            23220
   88 CONTINUE                                                             23230
      IF(IRVOL2.EQ.1) GOTO 86                                              23240
      CALL VOLUME(VL2,RM,POTC,DP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ,       23250
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM2,SM2,GRV1,     23260
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2   23270
     $,GR2,1)                                                              23280
      IF(E.EQ.0.) GOTO 86                                                  23290
      POTCD=PCSV                                                           23300
      IF(IFC.EQ.2) POTCD=PCSV+DPDX2*(1./D-1./(1.-E))                       23310
      CALL VOLUME(VL2,RM,POTCD,D,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ,       23320
     $GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM2,SM2,GRV1,     23330
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2   23340
     $,GR2,IFC)                                                            23350
   86 CONTINUE                                                             23360
      TPOLH=TAVH*SQRT(SQRT(SM1/SUMM1))                                     23370
      TPOLC=TAVC*SQRT(SQRT(SM2/SUMM2))                                     23380
      g1=gmag1(1)                                                          23390
      g2=gmag2(1)                                                          23400
      IF(MODE.EQ.1)TPOLC=TPOLH*SQRT(SQRT((G2/G1)**GR1))                    23410
      IF(MODE.EQ.1)TAVC=TPOLC/SQRT(SQRT(SM2/SUMM2))                        23420
      RATH=1.                                                              23430
      RATC=1.                                                              23440
      tph=10000.*tpolh                                                     23450
      tpc=10000.*tpolc                                                     23460
      IF(IFAT1.NE.0) CALL ATM(tph,WL,RATH)                                 23470
      IF(IFAT2.NE.0) CALL ATM(tpc,WL,RATC)                                 23480
      RATCH=RATC/RATH                                                      23490
      call lum(hlum,xh,yh,gr1,wl,tpolh,n1,n1,1,sbrh,rv,rvq,glump1,         23500
     $glump2,grv1,grv2,mmsave,summ1d,fr1,sm1d,ifat1,vold,rm,poth,          23510
     $f1,d,hld,snth)                                                       23520
      sbrc=ratch*sbrh*(9.-3.*xh+2.*yh)*(exp(xlump/tpolh)-1.)/              23530
     $((9.-3.*xc+2.*yc)*(exp(xlump/tpolc)-1.))                             23540
      call lum(clum,xc,yc,gr2,wl,tpolc,n2,n1,2,sbrt,rv,rvq,glump1,         23550
     $glump2,grv1,grv2,mmsave,summ2d,fr2,sm2d,ifat2,vold,rm,potc,          23560
     $f2,d,hld,snth)                                                       23570
      IF(IPB.EQ.1) SBRC=SBRT                                               23580
      IF(MODE.GT.0)CLUM=CLUM*SBRC/SBRT                                     23590
      IF(MODE.LE.0)SBRC=SBRT                                               23600
      if(mref.eq.2) goto 30                                                23610
      ratlum=hlum/clum                                                     23620
      call bolo(ratlum,tavh,tavc,wl,ifat1,ifat2,ratbol)                    23630
      rb=1./ratbol                                                         23640
      call olump(rv,grx,gry,grz,rvq,grxq,gryq,grzq,slump1,slump2,mmsave    23650
     $,gr1,alb1,rb,wl,tpolh,sbrh,summ1,n1,n2,1,ifat1,fr1,xc,yc,d,snth      23660
     $,csth,snfi,csfi,tld,glump1,glump2)                                   23670
      rb=ratbol                                                            23680
      call olump(rv,grx,gry,grz,rvq,grxq,gryq,grzq,slump1,slump2,mmsave    23690
     $,gr2,alb2,rb,wl,tpolc,sbrc,summ2,n1,n2,2,ifat2,fr2,xh,yh,d,snth      23700
     $,csth,snfi,csfi,tld,glump1,glump2)                                   23710
      return                                                               23720
   30 continue                                                             23730
      sbr1b=tpolh**4/dint1                                                 23740
      sbr2b=tpolc**4/dint2                                                 23750
      LT=N1+1                                                              23760
      IMAX1=MMSAVE(LT)                                                     23770
      DO 80 I=1,IMAX1                                                      23780
      RFTEMP(I)=1.                                                         23790
   80 RF1(I)=1.                                                            23800
      LT=N1+N2+2                                                           23810
      IMAX2=MMSAVE(LT)                                                     23820
      DO 81 I=1,IMAX2                                                      23830
   81 RF2(I)=1.                                                            23840
      DO 93 NR=1,NREF                                                      23850
      CALL LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,MMSAVE,           23860
     $gr1,alb1,wl,tpolh,sbrh,summ1,n1,n2,1,ifat1,fr1,snth,csth,snfi,       23870
     $csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,xbol2,ybol2,grv1,     23880
     $grv2,sbr1b,sbr2b,rftemp,rf2,gmag1,gmag2,dint1)                       23890
      CALL LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,MMSAVE,           23900
     $GR2,ALB2,WL,TPOLC,SBRC,SUMM2,N1,N2,2,IFAT2,fr2,snth,csth,snfi,       23910
     $csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,zz1,zz2,xbol1,ybol1,          23920
     $grv1,grv2,sbr1b,sbr2b,rf2,rf1,gmag1,gmag2,dint2)                     23930
      DO 70 I=1,IMAX1                                                      23940
   70 RF1(I)=RFTEMP(I)                                                     23950
   93 CONTINUE                                                             23960
      RETURN                                                               23970
      END                                                                  23980
      SUBROUTINE VOLUME(V,Q,P,D,FF,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,           23990
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMM,SM,       24000
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1    24010
     $,GMAG2,GREXP,IFC)                                                    24020
c  Experimental version of April 9, 1992                                   24030
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*   24040
     $),MMSAVE(*),FR1(*),FR2(*),HLD(*),SNTH(*),CSTH(*),SNFI(*),CSFI(*)     24050
     $,GRV1(*),GRV2(*),GLUMP1(*),GLUMP2(*),XX1(*),YY1(*),ZZ1(*),XX2(*),    24060
     $YY2(*),ZZ2(*),CSBT1(*),CSBT2(*),GMAG1(*),GMAG2(*)                    24070
      DP=1.E-3*P                                                           24080
      IF (IFC.EQ.1) DP=0.                                                  24090
      TOL=1.E-6*P**2                                                       24100
      DELP=0.                                                              24110
      KNTR=0                                                               24120
   16 P=P+DELP                                                             24130
      KNTR=KNTR+1                                                          24140
      IF(KNTR.GE.20) TOL=TOL+TOL                                           24150
      PS=P                                                                 24160
      DO 17 I=1,IFC                                                        24170
      P=PS                                                                 24180
      IF(I.EQ.1) P=P+DP                                                    24190
      CALL SURFAS(Q,P,N,N1,KOMP,RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,         24200
     $MMSAVE,FR1,FR2,HLD,FF,D,SNTH,CSTH,SNFI,CSFI,GRV1,GRV2,XX1,YY1,ZZ1,   24210
     $XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,GMAG2,GREXP)             24220
      IF(KOMP.EQ.2) GOTO 14                                                24230
      call lum(1.,1.,0.,grexp,1.,1.,n,n1,1,sbrd,rv,rvq,glump1,glump2,      24240
     $grv1,grv2,mmsave,summ,fr1,sm,0,vol,q,p,ff,d,hld,snth)                24250
      GOTO 15                                                              24260
   14 call lum(1.,1.,0.,grexp,1.,1.,n,n1,2,sbrd,rv,rvq,glump1,glump2,      24270
     $grv1,grv2,mmsave,summ,fr2,sm,0,vol,q,p,ff,d,hld,snth)                24280
   15 CONTINUE                                                             24290
      IF(I.EQ.1) VOLS=VOL                                                  24300
      VOL2=VOLS                                                            24310
   17 VOL1=VOL                                                             24320
      IF(IFC.EQ.1) V=VOL                                                   24330
      IF(IFC.EQ.1) RETURN                                                  24340
      DPDV=DP/(VOL2-VOL1)                                                  24350
      DELP=(V-VOL1)*DPDV                                                   24360
      ABDELP=ABS(DELP)                                                     24370
      IF(ABDELP.GT.TOL) GOTO 16                                            24380
      P=PS                                                                 24390
      RETURN                                                               24400
      END                                                                  24410
      SUBROUTINE ATM(TEMP,WAVE,RATIO)                                      24420
C     VERSION OF OCTOBER 30, 1987. DESTROY EARLIER VERSIONS.               24430
C     THIS SUBROUTINE COMPUTES THE RATIO OF STELLAR ATMOSPHERE TO          24440
C     PLANCKIAN FLUX FOR MAIN SEQUENCE STARS.                              24450
C     WARNING: LIMITS OF APPLICABILITY OF THIS SUBROUTINE AS FOLLOWS.      24460
C     IN TEMP RANGE 4,000 TO 25,000 K, WAVE LENGTH LIMITS .0912 TO         24470
C     6.5647 MICRONS. IN TEMP RANGE 25,000 TO 50,000 K, WAVE LENGTH        24480
C     LIMITS .0912 TO 2.2794 MICRONS. SUBROUTINE SHOULD NOT BE USED        24490
C     OUTSIDE THESE RANGES UNDER ANY CIRCUMSTANCES. WARNING MESSAGES       24500
C     ARE PRINTED IF RANGES EXCEEDED.                                      24510
      DIMENSION A(500),AA(196),AB(182),AC(105)                             24520
      EQUIVALENCE (A(1),AA(1)),(A(197),AB(1)),(A(379),AC(1))               24530
      DATA AA/.907883,-1.66066,.691436,4*0.0,.280961,-.447328,.130532,     24540
     $4*0.0,.142757,-.248947,.0620680,4*0.0,1.20223,-2.33426,1.07589,      24550
     $4*0.0,.403049,-.586521,.171767,4*0.0,.161172,-.260420,.0640562,      24560
     $4*0.0,1.78103,-3.64231,1.70609,4*0.0,.667290,-1.06532,.341809,       24570
     $4*0.0,.181770,-.345798,.0830478,4*0.0,-.329302E+2,.598843E+3,        24580
     $-.418910E+4,.146768E+5,-.274273E+5,.26093E+5,-.992741E+4,            24590
     $.442164,-1.17612,.558862,4*0.0,-.0671895,-.114919,.0244139,4*0.0,    24600
     $-.207863E+3,.243857E+4,-.115507E+5,.283797E+5,-.38233E+5,            24610
     $.268213E+5,-.76661E+4,.44286,-.984141,.391429,4*0.0,-.0128738,       24620
     $-.165689,.0419392,4*0.0,-.245193E+3,.287031E+4,-.135822E+5,.333606   24630
     $E+5,-.449465E+5,.315396E+5,-.901805E+4,.542751,-1.13157,.450211,     24640
     $4*0.0,.0371818,-.212137,.0551231,4*0.0,-.276230E+3,.322251E+4,       24650
     $-.1522E+5,.373523E+5,-.503142E+5,.353117E+5,-.101001E+5,             24660
     $.685486,-1.3322,.527152,4*0.0,.099154,-.266633,.0700431,4*0.0,       24670
     $-.280198E+3,.324071E+4,-.152149E+5,.371972E+5,-.499884E+5,           24680
     $.350354E+5,-.100136E+5,.912472,-1.64954,.647666,4*0.0,.192641,       24690
     $-.345623,.0898134,4*0.0,-.649967,.332897E+2,-.395009E+3,             24700
     $.229755E+4,-.701922E+4,.105395E+5,-.609668E+4,.66853,-1.55813,       24710
     $.773398,4*0.0,.202931,-.391874,.117945,4*0.0,-1.25062,.246184E+2,    24720
     $-.934153E+2,-.358292E+2,.745529E+3,-.133152E+4,.738503E+3/           24730
      DATA AB/.779381,-1.72595,.842554,4*0.0,.269815,-.455823,.135618,     24740
     $4*0.0,-1.86875,-.233737E+1,.435333E+3,-.330973E+4,.100351E+5,        24750
     $-.137955E+5,.714497E+4,.234313,-.196784,-.142740,4*0.0,.341826,      24760
     $-.513333,.149367,4*0.0,                                              24770
     $        -1.97002,-.28814E+2,.778254E+3,-.498095E+4,.139559E+5,       24780
     $-.182625E+5,.913243E+4,.607074,-.994799,.290486,4*0.0,.360036,       24790
     $-.520072,.147628,4*0.0,-1.99412,-.300984E+2,.735863E+3,-.448795E+4   24800
     $,.121397E+5,-.154572E+5,.756405E+4,.6119,-.976964,.278779,4*0.0,     24810
     $.346578,-.487175,.136367,4*0.0,-.227562E+1,-.734396E+1,              24820
     $.250026E+3,-.110376E+4,.179596E+4,-.947930E+3,-.903513E+2,           24830
     $.650197,-1.01937,.297332,4*0.0,.295735,-.382176,.09522,4*0.0,        24840
     $-.242829E+1,.26082E+2,-495.932,.400884E+4,-.135251E+5,.201766E+5,    24850
     $-.110779E+5,.515008,-.840408,.268146,4*0.0,.0567424,-.0326554,       24860
     $-.0248079,4*0.0,-1.98073,5.89063,.923565E+2,-.180606E+2,-.324872E+   24870
     $4,.106725E+5,-.981050E+4,.696006,-.158377E+1,.840335,4*0.0,          24880
     $-.0372941,.113422,-.0827454,4*0.0,-2.20129,2.57431,.954215E+2,       24890
     $.140538E+3,-.358564E+4,.10078E+5,-.843473E+4,.877022,-2.11887,       24900
     $.123868E+1,4*0.0,-.157214,.322907,-.162776,4*0.0,-2.28223,           24910
     $1.88918,.906366E+2,.182422E+3,-.344699E+4,.915905E+4,-.736667E+4,    24920
     $.422139,-.673966,.13635,4*0.0,-.269755,.515384,-.234871,4*0.0/       24930
      DATA AC/-1.72294,.111136E+2,.546943E+2,-.26662E+3,-.827417E+3,       24940
     $.463518E+4,-.501896E+4,.750032,-1.87636,1.12717,4*0.0,-.327529,      24950
     $.664143,-.307411,4*0.0,-1.82965,7.19243,.333359E+2,-.189727E+2,      24960
     $-.567324E+3,.71019E+3,.553124E+3,.968157,-.249486E+1,1.57956,4*0.0   24970
     $,-.413511,.836881,-.379774,4*0.0,-2.12746,5.72613,.206067E+2,        24980
     $.103493E+3,-.379424E+3,-.115088E+4,.285070E+4,1.30615,-3.49355,      24990
     $2.31489,4*0.0,-.551384,.109752E+1,-.485262,4*0.0,-2.44968,           25000
     $5.38777,.918011E+1,.156549E+3,-.162539E+3,-.191797E+4,.316283E+4,    25010
     $.184856E+1,-.512211E+1,.348495E+1,4*0.0,-.860859,.163724E+1,         25020
     $-.699413,4*0.0,-2.36238,6.98545,.340574E+2,.139450E+3,               25030
     $-.680221E+3,-.155574E+4,.462879E+4,.18498E+1,-.579597E+1,            25040
     $.43628E+1,4*0.0,-.111754E+1,.214057E+1,-.924189,4*0.0/               25050
   51 FORMAT(' TEMPERATURE RANGE EXCEEDED IN SUBROUTINE ATM')              25060
   52 FORMAT(' WAVE LENGTH RANGE EXCEEDED IN SUBROUTINE ATM')              25070
      IF (TEMP.GT.50000.) WRITE(6,51)                                      25080
      IF (TEMP.LT.4000.) WRITE (6,51)                                      25090
      IF (WAVE.LT..0912) WRITE(6,52)                                       25100
      IF(TEMP.LE.25000.) GOTO 53                                           25110
      IF(WAVE.GT.2.2794) WRITE(6,52)                                       25120
      GOTO 54                                                              25130
   53 IF(WAVE.GT.6.5647) WRITE(6,52)                                       25140
   54 CONTINUE                                                             25150
      IF (WAVE.GT..3647) GOTO 15                                           25160
      NWAVE=0                                                              25170
      GOTO 30                                                              25180
   15 IF (WAVE.GT..8206) GOTO 25                                           25190
      NWAVE=1                                                              25200
      GOTO 30                                                              25210
   25 NWAVE=2                                                              25220
   30 IF (TEMP.LT.40000.) GOTO 300                                         25230
      FRACT=(TEMP-40000.)/10000.                                           25240
      DIV1=.0203                                                           25250
      DIV2=.0280                                                           25260
      NUMBER=0                                                             25270
      GOTO 500                                                             25280
  300 IF (TEMP.LT.30000.) GOTO 305                                         25290
      FRACT=(TEMP-30000.)/10000.                                           25300
      DIV1=.0280                                                           25310
      DIV2=.0302                                                           25320
      NUMBER=1                                                             25330
      GOTO 500                                                             25340
  305 IF (TEMP.LT.25000.) GOTO 310                                         25350
      FRACT=(TEMP-25000.)/5000.                                            25360
      DIV1=.0302                                                           25370
      DIV2=.0701                                                           25380
      NUMBER=2                                                             25390
      GOTO 500                                                             25400
  310 IF (TEMP.LT.20000.) GOTO 315                                         25410
      FRACT=(TEMP-20000.)/5000.                                            25420
      DIV1=.0701                                                           25430
      DIV2=.0504                                                           25440
      NUMBER=3                                                             25450
      GOTO 500                                                             25460
  315 IF (TEMP.LT.18000.) GOTO 320                                         25470
      FRACT=(TEMP-18000.)/2000.                                            25480
      DIV1=.0504                                                           25490
      DIV2=.0504                                                           25500
      NUMBER=4                                                             25510
      GOTO 500                                                             25520
  320 IF (TEMP.LT.16000.) GOTO 325                                         25530
      FRACT=(TEMP-16000.)/2000.                                            25540
      DIV1=.0504                                                           25550
      DIV2=.0504                                                           25560
      NUMBER=5                                                             25570
      GOTO 500                                                             25580
  325 IF (TEMP.LT.14000.) GOTO 330                                         25590
      FRACT=(TEMP-14000.)/2000.                                            25600
      DIV1=.0504                                                           25610
      DIV2=.0504                                                           25620
      NUMBER=6                                                             25630
      GOTO 500                                                             25640
  330 IF (TEMP.LT.12000.) GOTO 335                                         25650
      FRACT=(TEMP-12000.)/2000.                                            25660
      DIV1=.0504                                                           25670
      DIV2=.1000                                                           25680
      NUMBER=7                                                             25690
      GOTO 500                                                             25700
  335 IF (TEMP.LT.11000.) GOTO 340                                         25710
      FRACT=(TEMP-11000.)/1000.                                            25720
      DIV1=.1000                                                           25730
      DIV2=.1000                                                           25740
      NUMBER=8                                                             25750
      GOTO 500                                                             25760
  340 IF (TEMP.LT.10000.) GOTO 345                                         25770
      FRACT=(TEMP-10000.)/1000.                                            25780
      DIV1=.1000                                                           25790
      DIV2=.0912                                                           25800
      NUMBER=9                                                             25810
      GOTO 500                                                             25820
  345 IF (TEMP.LT.9500.) GOTO 350                                          25830
      FRACT=(TEMP-9500.)/500.                                              25840
      DIV1=.0912                                                           25850
      DIV2=.0912                                                           25860
      NUMBER=10                                                            25870
      GOTO 500                                                             25880
  350 IF (TEMP.LT.9000.) GOTO 355                                          25890
      FRACT=(TEMP-9000.)/500.                                              25900
      DIV1=.0912                                                           25910
      DIV2=.0912                                                           25920
      NUMBER=11                                                            25930
      GOTO 500                                                             25940
  355 IF (TEMP.LT.8500.) GOTO 360                                          25950
      FRACT=(TEMP-8500.)/500.                                              25960
      DIV1=.0912                                                           25970
      DIV2=.0912                                                           25980
      NUMBER=12                                                            25990
      GOTO 500                                                             26000
  360 IF (TEMP.LT.8000.) GOTO 365                                          26010
      FRACT=(TEMP-8000.)/500.                                              26020
      DIV1=.0912                                                           26030
      DIV2=.0912                                                           26040
      NUMBER=13                                                            26050
      GOTO 500                                                             26060
  365 IF (TEMP.LT.7500.) GOTO 370                                          26070
      FRACT=(TEMP-7500.)/500.                                              26080
      DIV1=.0912                                                           26090
      DIV2=.1239                                                           26100
      NUMBER=14                                                            26110
      GOTO 500                                                             26120
  370 IF (TEMP.LT.7000.) GOTO 375                                          26130
      FRACT=(TEMP-7000.)/500.                                              26140
      DIV1=.1239                                                           26150
      DIV2=.1239                                                           26160
      NUMBER=15                                                            26170
      GOTO 500                                                             26180
  375 IF (TEMP.LT.6500.) GOTO 380                                          26190
      FRACT=(TEMP-6500.)/500.                                              26200
      DIV1=.1239                                                           26210
      DIV2=.1239                                                           26220
      NUMBER=16                                                            26230
      GOTO 500                                                             26240
  380 IF (TEMP.LT.6000.) GOTO 385                                          26250
      FRACT=(TEMP-6000.)/500.                                              26260
      DIV1=.1239                                                           26270
      DIV2=.1444                                                           26280
      NUMBER=17                                                            26290
      GOTO 500                                                             26300
  385 IF (TEMP.LT.5500.) GOTO 390                                          26310
      FRACT=(TEMP-5500.)/500.                                              26320
      DIV1=.1444                                                           26330
      DIV2=.1444                                                           26340
      NUMBER=18                                                            26350
      GOTO 500                                                             26360
  390 IF (TEMP.LT.5000.) GOTO 395                                          26370
      FRACT=(TEMP-5000.)/500.                                              26380
      DIV1=.1444                                                           26390
      DIV2=.1444                                                           26400
      NUMBER=19                                                            26410
      GOTO 500                                                             26420
  395 IF (TEMP.LT.4500.) GOTO 400                                          26430
      FRACT=(TEMP-4500.)/500.                                              26440
      DIV1=.1444                                                           26450
      DIV2=.1444                                                           26460
      NUMBER=20                                                            26470
      GOTO 500                                                             26480
  400 FRACT=(TEMP-4000.)/500.                                              26490
      DIV1=.1444                                                           26500
      DIV2=.1623                                                           26510
      NUMBER=21                                                            26520
  500 M=NUMBER*21+NWAVE*7+1                                                26530
      N=(NUMBER+1)*21+NWAVE*7+1                                            26540
      SUM1=A(M)                                                            26550
      SUM2=A(N)                                                            26560
      ALWAV1=ALOG10(WAVE/DIV1)                                             26570
      ALWAV2=ALOG10(WAVE/DIV2)                                             26580
      DO 20 I=1,6                                                          26590
      M1=M+I                                                               26600
      SUM1=SUM1+A(M1)*ALWAV1**I                                            26610
      N1=N+I                                                               26620
      SUM2=SUM2+A(N1)*ALWAV2**I                                            26630
   20 CONTINUE                                                             26640
      ALRATI=SUM2*(1.-FRACT)+SUM1*FRACT                                    26650
      RATIO=10.**ALRATI                                                    26660
      RETURN                                                               26670
      END                                                                  26680
      SUBROUTINE FOUR (L,K,THETA,RHO,AA,BB)                                26690
c   Version of April 9, 1992                                               26700
      DIMENSION THETA(*),RHO(*),AA(*),BB(*)                                26710
      KK=K+1                                                               26720
      EL=L                                                                 26730
      DEL=2./EL                                                            26740
      KST=2                                                                26750
      FAC=6.2831852/EL                                                     26760
      DO 65 N=1,KK                                                         26770
      AA(N)=0.                                                             26780
   65 BB(N)=0.                                                             26790
      THETA(L)=3.1415927*(2.-1./EL)                                        26800
      DO 29 J=1,L                                                          26810
      XJ=J-1                                                               26820
      TH=FAC*XJ                                                            26830
      KT=KST                                                               26840
      DO 50 M=KT,L                                                         26850
      IF(TH.LT.THETA(M)) GOTO 22                                           26860
      GOTO 50                                                              26870
   22 MP=M-1                                                               26880
      P=(TH-THETA(MP))/(THETA(M)-THETA(MP))                                26890
      RRO=(1.-P)*RHO(MP)+P*RHO(M)                                          26900
      KST=M                                                                26910
      GOTO 25                                                              26920
   50 CONTINUE                                                             26930
   25 DO 29 N=1,KK                                                         26940
      EN=N-1                                                               26950
      ENTHET=EN*TH                                                         26960
      AA(N)=AA(N)+RRO*COS(ENTHET)*DEL                                      26970
   29 BB(N)=BB(N)+RRO*SIN(ENTHET)*DEL                                      26980
      AA(1)=AA(1)*.5                                                       26990
      RETURN                                                               27000
      END                                                                  27010
      SUBROUTINE SORTER (L,THETA,RHO)                                      27020
c  Version of April 9, 1992                                                27030
      DIMENSION THETA(*),RHO(*)                                            27040
      M=1                                                                  27050
      DO 1 I=1,L                                                           27060
      M=M+M                                                                27070
      IF((M-L).GT.0) GOTO 2                                                27080
    1 CONTINUE                                                             27090
    2 M=(M-1)/2                                                            27100
    3 IF(M.EQ.0) RETURN                                                    27110
      NN=L-M                                                               27120
      DO 5 J=1,NN                                                          27130
      DO 5 I=1,J,M                                                         27140
      IL=J+1-I                                                             27150
      IML=IL+M                                                             27160
      IF((THETA(IML)-THETA(IL)).GE.0.) GOTO 5                              27170
      X=THETA(IML)                                                         27180
      Y=RHO(IML)                                                           27190
      THETA(IML)=THETA(IL)                                                 27200
      RHO(IML)=RHO(IL)                                                     27210
      THETA(IL)=X                                                          27220
      RHO(IL)=Y                                                            27230
    5 CONTINUE                                                             27240
      M=M/2                                                                27250
      GOTO 3                                                               27260
      END                                                                  27270
      SUBROUTINE ELLONE(FFF,DD,RMSP,XL1,OM1,XL2,OM2)                       27280
C     VERSION OF 4/8/82. DESTROY EARLIER VERSIONS.                         27290
C     XL2 AND OM2 VALUES ASSUME SYNCHRONOUS ROTATION AND CIRCULAR ORBIT.   27300
C     THEY ARE NOT NEEDED FOR NON-SYNCHRONOUS OR NON-CIRCULAR CASES.       27310
      IMPLICIT REAL*8(A-H,O-Z)                                             27320
      REAL RMSP,XL1,OM1,XL2,OM2,FFF,DD                                     27330
      RMASS=RMSP                                                           27340
      FF=FFF                                                               27350
      D=DD                                                                 27360
      XL=.5D0*D                                                            27370
      DO 5 I=1,2                                                           27380
      RFAC=FF**2                                                           27390
      IF(I.EQ.2) RFAC=1.D0                                                 27400
      IF(I.EQ.2) D=1.D0                                                    27410
      DSQ=D*D                                                              27420
      DELXL=0.D0                                                           27430
      RM1=RMASS+1.D0                                                       27440
   88 XL=XL+DELXL                                                          27450
      XSQ=XL*XL                                                            27460
      P=(D-XL)**2                                                          27470
      RP=DABS(D-XL)                                                        27480
      PRP=P*RP                                                             27490
      F=RFAC*RM1*XL-1.D0/XSQ-RMASS*(XL-D)/PRP-RMASS/DSQ                    27500
      DXLDF=1.D0/(RFAC*RM1+2.D0/(XSQ*XL)+2.D0*RMASS/PRP)                   27510
      DELXL=-F*DXLDF                                                       27520
      ABDEL=DABS(DELXL)                                                    27530
      IF(ABDEL.GT..000001D0) GOTO 88                                       27540
      IF(I.EQ.2) GOTO 8                                                    27550
      XL1=XL                                                               27560
      OM1=1.D0/XL+RMASS*((1.D0/RP)-XL/DSQ)+RM1*.5D0*XSQ*RFAC               27570
      IF(RMSP.GT.1.)RMASS=1.D0/RMASS                                       27580
      XMU3=RMASS/(3.D0*(RMASS+1.D0))                                       27590
      XMU3CR=XMU3**.33333333333D0                                          27600
    5 XL=1.D0+XMU3CR+XMU3CR*XMU3CR/3.D0+XMU3/9.D0                          27610
    8 IF(RMSP.GT.1.)XL=D-XL                                                27620
      RMASS=RMSP                                                           27630
      XL2=XL                                                               27640
      OM2=1.D0/DABS(XL)+RMASS*((1.D0/DSQRT(1.D0-XL-XL+XL*XL))-XL)+RM1*     27650
     $.5D0*XL*XL                                                           27660
      RETURN                                                               27670
      END                                                                  27680
      SUBROUTINE RING(QSP,OMSP,KOMP,L,FR,HLD,R1SP,RLSP)                    27690
c  Version of April 9, 1992                                                27700
      IMPLICIT REAL*8(A-H,O-Z)                                             27710
      REAL QSP,OMSP,FR,HLD,Z,XLSP,OM1SP,XL2SP,OM2SP,EMSP,EYSP,THESP,ELSP   27720
     $,R1SP,RLSP                                                           27730
      DIMENSION RAD(100),THET(100),AA(3),BB(3),FI(150),THA(150),FR(*),     27740
     $HLD(*)                                                               27750
      IA=1                                                                 27760
      IX=0                                                                 27770
      LR=L+1                                                               27780
      DO 92 I=1,LR                                                         27790
      THA(I)=0.D0                                                          27800
   92 FI(I)=-.1D0                                                          27810
      Q=QSP                                                                27820
      OMEGA=OMSP                                                           27830
      K=3                                                                  27840
      EL=L                                                                 27850
      ELSP=L                                                               27860
      DEL=2.D0/EL                                                          27870
      CALL ELLONE(1.,1.,QSP,XLSP,OM1SP,XL2SP,OM2SP)                        27880
      CALL NEKMIN(QSP,OMSP,XLSP,Z)                                         27890
      XL=XLSP                                                              27900
      OM1=OM1SP                                                            27910
      XL2=XL2SP                                                            27920
      OM2=OM2SP                                                            27930
      QQ=Q                                                                 27940
      XLSQ=XL*XL                                                           27950
      IF(Q.GT.1.D0) QQ=1.D0/Q                                              27960
      RMAX=DEXP(.345D0*DLOG(QQ)-1.125D0)                                   27970
      R=RMAX*(OM1-OMEGA)/(OM1-OM2)                                         27980
      DO 22 IT=1,L                                                         27990
      EYT=IT                                                               28000
      TH=EYT*1.5707963D0/EL                                                28010
      COSQ=DCOS(TH)**2                                                     28020
      DELR=0.D0                                                            28030
   14 R=DABS(R+DELR)                                                       28040
      RSQ=R*R                                                              28050
      X2R2=XLSQ+RSQ                                                        28060
      RX2R2=DSQRT(X2R2)                                                    28070
      XM2R2=(XL-1.D0)**2+RSQ                                               28080
      RXM2R2=DSQRT(XM2R2)                                                  28090
      OM=1.D0/RX2R2+Q*(1.D0/RXM2R2-XL)+.5D0*(Q+1.D0)*(XLSQ+RSQ*COSQ)       28100
      DOMDR=-R/(X2R2*RX2R2)-Q*R/(XM2R2*RXM2R2)+(Q+1.D0)*COSQ*R             28110
      DELR=(OMEGA-OM)/DOMDR                                                28120
      ABDELR=DABS(DELR)                                                    28130
      IF(ABDELR.GT..00001D0) GOTO 14                                       28140
      RAD(IT)=R                                                            28150
   22 THET(IT)=TH*4.D0                                                     28160
      R1=RAD(1)                                                            28170
      RL=RAD(L)                                                            28180
      R1SP=R1                                                              28190
      RLSP=RL                                                              28200
      R90SQ=RL*RL                                                          28210
      DO 18 IJ=1,L                                                         28220
      EYJ=IJ                                                               28230
      RAD(IJ)=RAD(IJ)-(RL-R1)*(EYJ-1.D0)/EL                                28240
   18 CONTINUE                                                             28250
      DO 65 N=1,K                                                          28260
      AA(N)=0.D0                                                           28270
   65 BB(N)=0.D0                                                           28280
      DO 29 J=1,L                                                          28290
      DO 29 N=1,K                                                          28300
      EN=N-1                                                               28310
      ENTHET=EN*THET(J)                                                    28320
      AA(N)=AA(N)+RAD(J)*DCOS(ENTHET)*DEL                                  28330
   29 BB(N)=BB(N)+RAD(J)*DSIN(ENTHET)*DEL                                  28340
      AA(1)=AA(1)*.5D0                                                     28350
      IF(KOMP.EQ.2) XL=1.-XLSP                                             28360
      XLSQ=XL*XL                                                           28370
      DIS=RL/XL-.0005D0                                                    28380
      DO 42 IR=1,L                                                         28390
      LL=IR-1                                                              28400
      EY=L+1-IR                                                            28410
      THA(IR)=1.5707963D0*EY/EL                                            28420
      IF(THA(IR).LT.1.57079D0) GOTO 82                                     28430
      COT=0.D0                                                             28440
      GOTO 83                                                              28450
   82 COT=1.D0/DTAN(THA(IR))                                               28460
   83 IF(COT.GE.DIS) GOTO 50                                               28470
      COSSQ=DCOS(THA(IR))**2                                               28480
      A0=AA(1)                                                             28490
      A1=AA(2)                                                             28500
      A2=AA(3)                                                             28510
      B1=BB(2)                                                             28520
      B2=BB(3)                                                             28530
      DELSIN=0.D0                                                          28540
      KNTR=0                                                               28550
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                                28560
   88 SINTH=SINTH+DELSIN                                                   28570
      KNTR=KNTR+1                                                          28580
      IF(SINTH.GT.1.D0)SINTH=1.D0/SINTH                                    28590
      CSQ=1.D0-SINTH*SINTH                                                 28600
      COSTH=DSQRT(CSQ)                                                     28610
      SINSQ=SINTH*SINTH                                                    28620
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                            28630
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                        28640
      C4SQ=COS4*COS4                                                       28650
      SINCOS=SIN4*COS4                                                     28660
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS            28670
      ARC=DARSIN(SINTH)                                                    28680
      RR=RRR+(RL-R1)*(2.D0*ARC/3.1415926536D0-1.D0/EL)                     28690
      IF(KNTR.GT.30) GOTO 42                                               28700
      P=RR*SINTH                                                           28710
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2   28720
     $)*COSTH+(RL+RL-R1-R1)/(3.1415926536D0*COSTH)                         28730
      DPDSIN=RR+SINTH*DRDSIN                                               28740
      F=P*P/COSSQ-RR*RR-XLSQ                                               28750
      DFDSIN=(P+P)*DPDSIN/COSSQ-(RR+RR)*DRDSIN                             28760
      DELSIN=-F/DFDSIN                                                     28770
      ABDEL=DABS(DELSIN)                                                   28780
      IF(ABDEL.GT..00001D0)  GOTO 88                                       28790
   42 FI(IR)=DATAN(RR*COSTH/XL)                                            28800
   50 LL1=LL+1                                                             28810
      DELTH=1.5707963D0/EL                                                 28820
      DO 75 I=1,L                                                          28830
      EY=L+1-I                                                             28840
      EY=EY-.5D0                                                           28850
      THE=1.5707963D0*EY/EL                                                28860
      SNTH=DSIN(THE)                                                       28870
      EYSP=L+1-I                                                           28880
      EYSP=EYSP-.5                                                         28890
      THESP=THE                                                            28900
      CSTH=DCOS(THE)                                                       28910
      EMSP=SIN(THESP)*ELSP*1.3                                             28920
      MM=EMSP+1.                                                           28930
      XM=MM                                                                28940
      DELFI=3.1415926D0/XM                                                 28950
      HDELFI=1.5707963D0/XM                                                28960
      DO 75 J=1,MM                                                         28970
      IX=IX+1                                                              28980
      IF(I.LE.LL1) GOTO 43                                                 28990
      HLD(IX)=1.                                                           29000
      GOTO 75                                                              29010
   43 XJ=MM+1-J                                                            29020
      FE=3.1415926D0*(XJ-.5D0)/XM                                          29030
      PH2=FE+HDELFI                                                        29040
      PHB=PH2                                                              29050
      IF(FI(I).GT.(FE-HDELFI)) GOTO 51                                     29060
      HLD(IX)=1.                                                           29070
      GOTO 75                                                              29080
   51 IPL=I+1                                                              29090
      IF(FI(IPL).GT.0.D0) GOTO 66                                          29100
      RR=A0+A1-A2+(RL-R1)*(1.D0-1.D0/EL)                                   29110
      PH1=DELFI*(XJ-1.D0)                                                  29120
      TH1=DATAN(XL/RR)                                                     29130
      GOTO 56                                                              29140
   66 IF(FI(IPL).LT.(FE+HDELFI)) GOTO 52                                   29150
      HLD(IX)=0.                                                           29160
      GOTO 75                                                              29170
   52 IF(FI(IPL).LT.(FE-HDELFI)) GOTO 53                                   29180
      PH1=FI(IPL)                                                          29190
      TH1=THA(IPL)                                                         29200
      GOTO 56                                                              29210
   53 DELSIN=0.D0                                                          29220
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                                29230
      TANFE=DTAN(FE-HDELFI)                                                29240
   77 SINTH=SINTH+DELSIN                                                   29250
      IF(SINTH.GT.1.D0) SINTH=1.D0/SINTH                                   29260
      SINSQ=SINTH*SINTH                                                    29270
      CSQ=1.D0-SINSQ                                                       29280
      COSTH=DSQRT(CSQ)                                                     29290
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                            29300
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                        29310
      C4SQ=COS4*COS4                                                       29320
      SINCOS=SIN4*COS4                                                     29330
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS            29340
      ARC=DARSIN(SINTH)                                                    29350
      RR=RRR+(RL-R1)*(2.D0*ARC/3.1415926536D0-1.D0/EL)                     29360
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2   29370
     $)*COSTH+(RL+RL-R1-R1)/(3.1415926536D0*COSTH)                         29380
      F=RR*COSTH-XL*TANFE                                                  29390
      DFDSIN=COSTH*DRDSIN-RR*SINTH/COSTH                                   29400
      DELSIN=-F/DFDSIN                                                     29410
      ABDEL=DABS(DELSIN)                                                   29420
      IF(ABDEL.GT..00001D0)  GOTO 77                                       29430
      PH1=FE-HDELFI                                                        29440
      TH1=DATAN(XL/(RR*SINTH*DCOS(PH1)))                                   29450
   56 IF(FI(I).GT.(FE+HDELFI)) GOTO 57                                     29460
      PHB=FI(I)                                                            29470
      TH2=THA(I)                                                           29480
      GOTO 60                                                              29490
   57 DELSIN=0.D0                                                          29500
      SINTH=DSQRT(COSSQ*(XLSQ+R90SQ)/R90SQ)                                29510
      TANFE=DTAN(FE+HDELFI)                                                29520
   78 SINTH=SINTH+DELSIN                                                   29530
      IF(SINTH.GT.1.D0) SINTH=1.D0/SINTH                                   29540
      SINSQ=SINTH*SINTH                                                    29550
      CSQ=1.D0-SINSQ                                                       29560
      COSTH=DSQRT(CSQ)                                                     29570
      SIN4=8.D0*COSTH**3*SINTH-4.D0*COSTH*SINTH                            29580
      COS4=8.D0*CSQ*(CSQ-1.D0)+1.D0                                        29590
      C4SQ=COS4*COS4                                                       29600
      SINCOS=SIN4*COS4                                                     29610
      RRR=A0+A1*COS4+A2*(C4SQ+C4SQ-1.D0)+B1*SIN4+(B2+B2)*SINCOS            29620
      ARC=DARSIN(SINTH)                                                    29630
      RR=RRR+(RL-R1)*(2.D0*ARC/3.1415926536D0-1.D0/EL)                     29640
      DRDSIN=-A1*SINTH/COSTH-4.D0*A2*SINTH+B1-(B2+B2)*SINSQ/COSTH+(B2+B2   29650
     $)*COSTH+(RL+RL-R1-R1)/(3.1415926536D0*COSTH)                         29660
      F=RR*COSTH-XL*TANFE                                                  29670
      DFDSIN=COSTH*DRDSIN-RR*SINTH/COSTH                                   29680
      DELSIN=-F/DFDSIN                                                     29690
      ABDEL=DABS(DELSIN)                                                   29700
      IF(ABDEL.GT..00001D0)  GOTO 78                                       29710
      TH2=DATAN(XL/(RR*SINTH*DCOS(PH2)))                                   29720
   60 CTHT=DCOS(THA(IPL))                                                  29730
      CTH1=DCOS(TH1)                                                       29740
      CTH2=DCOS(TH2)                                                       29750
      STH1=DSIN(TH1)                                                       29760
      STH2=DSIN(TH2)                                                       29770
      DTH=TH2-TH1                                                          29780
      DCTH=CTH1-CTH2                                                       29790
      OMDP=PH2*DCTH-.5D0*(PH1*STH1+PHB*STH2)*DTH                           29800
      OMP=DELFI*(CTHT-CTH1)                                                29810
      OMN=OMP+OMDP                                                         29820
      HLD(IX)=OMN/(DELTH*DELFI*SNTH)                                       29830
   75 CONTINUE                                                             29840
      DO 94 JB=1,IX                                                        29850
      JA=IX+1-JB                                                           29860
   94 FR(JB)=HLD(JA)                                                       29870
      RETURN                                                               29880
      END                                                                  29890
      SUBROUTINE MODLOG(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2   29900
     $,HLD,RM,POTH,POTC,GR1,GR2,ALB1,ALB2,N1,N2,F1,F2,MOD,XINCL,THE,       29910
     $MODE,SNTH,CSTH,SNFI,CSFI,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,GLUMP1    29920
     $,GLUMP2,CSBT1,CSBT2,GMAG1,GMAG2)                                     29930
c    Experimental version of March 31, 1992.                               29940
      DOUBLE PRECISION RES,DRR,ANS,DUMM                                    29950
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ     29960
     $(*),MMSAVE(*),FR1(*),FR2(*),HLD(*),GRV1(*),GRV2(*),XX1(*),YY1(*),    29970
     $ZZ1(*),XX2(*),YY2(*),ZZ2(*),GLUMP1(*),GLUMP2(*),CSBT1(*),CSBT2(*)    29980
     $,GMAG1(*),GMAG2(*)                                                   29990
      DIMENSION DRR(4),RES(2),ANS(2),LX(2),MX(2)                           30000
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*)                            30010
      common /kfac/ kff1,kff2,kfo1,kfo2                                    30020
      COMMON /FLVAR/ PER,PSHIFT,DP,DS,EF,EFC,ECOS,PERR,PHPER,PCONJ,        30030
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc                30040
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC        30050
      COMMON /INVAR/ KH,IPBDUM,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,        30060
     $ifsmv2,icor1,icor2,ld                                                30070
   95 FORMAT(' WARNING: ALTHOUGH COMPONENT 2 DOES NOT EXCEED ITS LIMITIN   30080
     $G LOBE AT THE END OF ECLIPSE, IT DOES EXCEED THE LOBE AT PERIASTRO   30090
     $N')                                                                  30100
   99 FORMAT(' SPECIFIED ECLIPSE DURATION INCONSISTENT WITH OTHER PARAME   30110
     $TERS')                                                               30120
      DP=1.-E                                                              30130
      DS=DP*DP                                                             30140
      MOD=(MODE-2)**2                                                      30150
      IF(MOD.EQ.1) GR2=GR1                                                 30160
      IF(MOD.EQ.1) ALB2=ALB1                                               30170
      IF(MOD.EQ.1) POTC=POTH                                               30180
      MD4=(MODE-5)**2                                                      30190
      MD5=(2*MODE-11)**2                                                   30200
      call ellone(f1,dp,rm,xl1,po1cr,xl2,omo1)                             30210
      RMR=1./RM                                                            30220
      CALL ELLONE(F2,DP,RMR,XL1,po2c,XL2,omo2)                             30230
      po2cr=rm*po2c+(1.-rm)*.5                                             30240
      if(md4.eq.1) poth=po1cr                                              30250
      if(md5.eq.1) potc=po2cr                                              30260
      kff1=0                                                               30270
      kff2=0                                                               30280
      if(poth.lt.po1cr) kff1=1                                             30290
      if(potc.lt.po2cr) kff2=1                                             30300
      kfo1=0                                                               30310
      kfo2=0                                                               30320
      if(e.ne.0.) goto 100                                                 30330
      if(f1.ne.1.) goto 105                                                30340
      if(poth.lt.omo1) kfo1=1                                              30350
  105 if(f2.ne.1.) goto 100                                                30360
      if(potc.lt.omo1) kfo2=1                                              30370
  100 continue                                                             30380
      SINI=SIN(.0174533*XINCL)                                             30390
      VF=50.61455/PERIOD                                                   30400
      VFAC=VF*A                                                            30410
      VGAM=VGA*VUNIT/VFAC                                                  30420
      VFVU=VFAC/VUNIT                                                      30430
      IFC=2                                                                30440
      IF(E.NE.0.) GOTO 60                                                  30450
      PERR=1.570796                                                        30460
      IFC=1                                                                30470
   60 CONTINUE                                                             30480
      TRC=1.570796-PERR                                                    30490
   39 if(TRC.LT.0.) TRC=TRC+6.283185                                       30500
      if(trc.lt.0.) goto 39                                                30510
   40 if(trc.ge.6.283185) trc=trc-6.283185                                 30520
      if(trc.ge.6.283185) goto 40                                          30530
      HTRC=.5*TRC                                                          30540
      IF(ABS(1.570796-HTRC).LT.7.E-6) GOTO 101                             30550
      IF(ABS(4.712389-HTRC).LT.7.E-6) GOTO 101                             30560
      ECAN=2.*ATAN(SQRT((1.-E)/(1.+E))*TAN(HTRC))                          30570
      GOTO 103                                                             30580
  101 ECAN=3.141593                                                        30590
  103 XMC=ECAN-E*SIN(ECAN)                                                 30600
      IF(XMC.LT.0.) XMC=XMC+6.283185                                       30610
      PHPER=1.-XMC/6.283185                                                30620
      PCONJ=(XMC+PERR)/6.283185-.25+PSHIFT                                 30630
   38 if(pconj.ge.1.) pconj=pconj-1.                                       30640
      if(pconj.ge.1.) goto 38                                              30650
   41 if(pconj.lt.0.) pconj=pconj+1.                                       30660
      if(pconj.lt.0.) goto 41                                              30670
      PHPERI=PHPER+PCONJ                                                   30680
      EF=1.-E*E                                                            30690
      EFC=SQRT(EF)                                                         30700
      ECOS=E*COS(PERR)                                                     30710
      IF(MODE.NE.-1) RETURN                                                30720
      if(kh.eq.17) goto 241                                                30730
      if((kh-12)**2.eq.1) goto 241                                         30740
      if((kh-12)**2.eq.4) goto 241                                         30750
      IF((KH-11)**2.LE.1) GOTO 241                                         30760
      IF((2*KH-41)**2.EQ.81) GOTO 241                                      30770
      RETURN                                                               30780
  241 CONTINUE                                                             30790
      EFCC=SQRT((1.-E)/(1.+E))                                             30800
      THER=THE*6.283185                                                    30810
      DELTR=.001                                                           30820
      DTR1=0.                                                              30830
      DTR2=0.                                                              30840
      VOLTOL=5.E-6                                                         30850
      DXMTOL=5.E-6                                                         30860
      TR0=1.570796-PERR                                                    30870
      HTR0=.5*TR0                                                          30880
      IF((1.570796-ABS(HTR0)).LT.7.E-6) GOTO 201                           30890
      IF((4.712389-ABS(HTR0)).LT.7.E-6) GOTO 201                           30900
      ECAN0=2.*ATAN(SQRT((1.-E)/(1.+E))*TAN(HTR0))                         30910
      GOTO 203                                                             30920
  201 ECAN0=3.141593                                                       30930
  203 XM0=ECAN0-E*SIN(ECAN0)                                               30940
      XM1=XM0-THER*(1.-.2*E)                                               30950
      XM2=XM0+THER*(1.-.2*E)                                               30960
      CALL KEPLER(XM1,E,DUM,TRR1)                                          30970
      CALL KEPLER(XM2,E,DUM,TRR2)                                          30980
  160 TRR1=TRR1+DTR1                                                       30990
      TRR2=TRR2+DTR2                                                       31000
      DO 161 IB=1,3                                                        31010
      TR1=TRR1                                                             31020
      TR2=TRR2                                                             31030
      IF(IB.EQ.2) TR1=TRR1+DELTR                                           31040
      IF(IB.EQ.3) TR2=TRR2+DELTR                                           31050
      IF(TR1.GT.TR0) TR0=TR0+6.283185                                      31060
      IF(TR0.GT.TR2) TR2=TR2+6.283185                                      31070
      DS1=EF/(1.+E*COS(TR1))                                               31080
      DS2=EF/(1.+E*COS(TR2))                                               31090
      RD1=DS1*SIN(ABS(TR1-TR0))                                            31100
      RD2=DS2*SIN(ABS(TR2-TR0))                                            31110
      TRE1=(TR0-TR1)/6.283185                                              31120
      TRE2=(TR2-TR0)/6.283185                                              31130
      CALL DURA(F2,XINCL,RM,DS1,TRE1,POTR,RA)                              31140
      CALL VOLUME(VS1,RM,POTR,DS1,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ       31150
     $,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,GRV1,    31160
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,        31170
     $GMAG2,GR1,1)                                                         31180
      CALL DURA(F2,XINCL,RM,DS2,TRE2,POTR,RA)                              31190
      CALL VOLUME(VS2,RM,POTR,DS2,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,GRXQ       31200
     $,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,GRV1,    31210
     $GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,GMAG1,        31220
     $GMAG2,GR2,1)                                                         31230
      IF(IB.NE.1) GOTO 185                                                 31240
      ECAN1=2.*ATAN(SQRT((1.-E)/(1.+E))*TAN(.5*TR1))                       31250
      ECAN2=2.*ATAN(SQRT((1.-E)/(1.+E))*TAN(.5*TR2))                       31260
      POTC=POTR                                                            31270
      DTHE=DS2                                                             31280
      DVOL=VS2-VS1                                                         31290
      XM1=ECAN1-E*SIN(ECAN1)                                               31300
      XM2=ECAN2-E*SIN(ECAN2)                                               31310
      IF(XM1.LT.0.) XM1=XM1+6.283185                                       31320
      IF(XM2.LT.0.) XM2=XM2+6.283185                                       31330
      DXM=XM2-XM1-2.*THER                                                  31340
      DDMDN1=-EFCC*(1.-E*COS(ECAN1))*COS(.5*ECAN1)**2/COS(.5*TR1)**2       31350
      DDMDN2=EFCC*(1.-E*COS(ECAN2))*COS(.5*ECAN2)**2/COS(.5*TR2)**2        31360
  185 CONTINUE                                                             31370
      IF(IB.NE.2) GOTO 162                                                 31380
      DRR(1)=(VS2-VS1-DVOL)/DELTR                                          31390
      DRR(2)=DDMDN1                                                        31400
  162 CONTINUE                                                             31410
      IF(IB.NE.3) GOTO 161                                                 31420
      DRR(3)=(VS2-VS1-DVOL)/DELTR                                          31430
      DRR(4)=DDMDN2                                                        31440
  161 CONTINUE                                                             31450
      RES(1)=-DVOL                                                         31460
      RES(2)=-DXM                                                          31470
      CALL DMINV(DRR,2,DUMM,LX,MX)                                         31480
      CALL DGMPRD(DRR,RES,ANS,2,2,1)                                       31490
      DTR1=ANS(1)                                                          31500
      DTR2=ANS(2)                                                          31510
      IF(ABS(DTR1).GT.VOLTOL) GOTO 160                                     31520
      IF(ABS(DTR2).GT.DXMTOL) GOTO 160                                     31530
      POTH=9999.99                                                         31540
      RMR=1./RM                                                            31550
      CALL ELLONE(F2,DTHE,RMR,XLA,OM1,XL2,OM2)                             31560
      OM1=RM*OM1+(1.-RM)*.5                                                31570
      IF(POTC.LT.OM1) GOTO 22                                              31580
      IF(RA.LE.XLA) GOTO 28                                                31590
   22 WRITE(6,99)                                                          31600
      RETURN                                                               31610
   28 CONTINUE                                                             31620
      IF(E.NE.0.) CALL VOLUME(VTHE,RM,POTC,DTHE,F2,N2,N1,2,RV,GRX,         31630
     $GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,   31640
     $SUMMD,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,      31650
     $GLUMP2,GMAG1,GMAG2,GR2,1)                                            31660
      IF(E.NE.0.) CALL VOLUME(VTHE,RM,POTC,DP,F2,N2,N1,2,RV,GRX,           31670
     $GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,   31680
     $SUMMD,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,      31690
     $GLUMP2,GMAG1,GMAG2,GR2,2)                                            31700
      CALL ELLONE(F2,DP,RMR,XLD,OMP,XL2,OM2)                               31710
      OMP=RM*OMP+(1.-RM)*.5                                                31720
      IF(POTC.LT.OMP) WRITE(6,95)                                          31730
      RETURN                                                               31740
      END                                                                  31750
      SUBROUTINE KEPLER(XMSP,ECSP,ECANSP,TRSP)                             31760
C     ORIGINAL VERSION OF 3/23/78.                                         31770
      IMPLICIT REAL*8(A-H,O-Z)                                             31780
      REAL XMSP,ECSP,ECANSP,TRSP                                           31790
      TOL=.00000001D0                                                      31800
      XM=XMSP                                                              31810
      EC=ECSP                                                              31820
      DLECAN=0.D0                                                          31830
      ECAN=XM                                                              31840
   18 ECAN=ECAN+DLECAN                                                     31850
      XMC=ECAN-EC*DSIN(ECAN)                                               31860
      DEDM=1.D0/(1.D0-EC*DCOS(ECAN))                                       31870
      DLECAN=(XM-XMC)*DEDM                                                 31880
      ABDLEC=DABS(DLECAN)                                                  31890
      IF(ABDLEC.GT.TOL) GOTO 18                                            31900
      TRSP=2.D0*DATAN(DSQRT((1.D0+EC)/(1.D0-EC))*DTAN(.5D0*ECAN))          31910
      IF(TRSP.LT.0.) TRSP=TRSP+6.283185                                    31920
      ECANSP=ECAN                                                          31930
      RETURN                                                               31940
      END                                                                  31950
      SUBROUTINE DGMPRD(A,B,R,N,M,L)                                       31960
c  Version of April 9, 1992                                                31970
      DIMENSION A(*),B(*),R(*)                                             31980
      DOUBLE PRECISION A,B,R                                               31990
      IR=0                                                                 32000
      IK=-M                                                                32010
      DO 10 K=1,L                                                          32020
      IK=IK+M                                                              32030
      DO 10 J=1,N                                                          32040
      IR=IR+1                                                              32050
      JI=J-N                                                               32060
      IB=IK                                                                32070
      R(IR)=0.D0                                                           32080
      DO 10 I=1,M                                                          32090
      JI=JI+N                                                              32100
      IB=IB+1                                                              32110
   10 R(IR)=R(IR)+A(JI)*B(IB)                                              32120
      RETURN                                                               32130
      END                                                                  32140
      SUBROUTINE DMINV(A,N,D,L,M)                                          32150
c  Version of April 9, 1992                                                32160
      DIMENSION A(*),L(*),M(*)                                             32170
      DOUBLE PRECISION A,D,BIGA,HOLD                                       32180
      D=1.D0                                                               32190
      NK=-N                                                                32200
      DO 80 K=1,N                                                          32210
      NK=NK+N                                                              32220
      L(K)=K                                                               32230
      M(K)=K                                                               32240
      KK=NK+K                                                              32250
      BIGA=A(KK)                                                           32260
      DO 20 J=K,N                                                          32270
      IZ=N*(J-1)                                                           32280
      DO 20 I=K,N                                                          32290
      IJ=IZ+I                                                              32300
   10 IF(DABS(BIGA).GE.DABS(A(IJ))) GOTO 20                                32310
      BIGA=A(IJ)                                                           32320
      L(K)=I                                                               32330
      M(K)=J                                                               32340
   20 CONTINUE                                                             32350
      J=L(K)                                                               32360
      IF(J.LE.K) GOTO 35                                                   32370
      KI=K-N                                                               32380
      DO 30 I=1,N                                                          32390
      KI=KI+N                                                              32400
      HOLD=-A(KI)                                                          32410
      JI=KI-K+J                                                            32420
      A(KI)=A(JI)                                                          32430
   30 A(JI) =HOLD                                                          32440
   35 I=M(K)                                                               32450
      IF(I.LE.K) GOTO 45                                                   32460
      JP=N*(I-1)                                                           32470
      DO 40 J=1,N                                                          32480
      JK=NK+J                                                              32490
      JI=JP+J                                                              32500
      HOLD=-A(JK)                                                          32510
      A(JK)=A(JI)                                                          32520
   40 A(JI) =HOLD                                                          32530
   45 IF(BIGA.NE.0.D0) GOTO 48                                             32540
      D=0.D0                                                               32550
      RETURN                                                               32560
   48 DO 55 I=1,N                                                          32570
      IF(I.EQ.K) GOTO 55                                                   32580
      IK=NK+I                                                              32590
      A(IK)=A(IK)/(-BIGA)                                                  32600
   55 CONTINUE                                                             32610
      DO 65 I=1,N                                                          32620
      IK=NK+I                                                              32630
      HOLD=A(IK)                                                           32640
      IJ=I-N                                                               32650
      DO 65 J=1,N                                                          32660
      IJ=IJ+N                                                              32670
      IF(I.EQ.K) GOTO 65                                                   32680
      IF(J.EQ.K) GOTO 65                                                   32690
      KJ=IJ-I+K                                                            32700
      A(IJ)=HOLD*A(KJ)+A(IJ)                                               32710
   65 CONTINUE                                                             32720
      KJ=K-N                                                               32730
      DO 75 J=1,N                                                          32740
      KJ=KJ+N                                                              32750
      IF(J.EQ.K) GOTO 75                                                   32760
      A(KJ)=A(KJ)/BIGA                                                     32770
   75 CONTINUE                                                             32780
      D=D*BIGA                                                             32790
      A(KK)=1.D0/BIGA                                                      32800
   80 CONTINUE                                                             32810
      K=N                                                                  32820
  100 K=(K-1)                                                              32830
      IF(K.LE.0) RETURN                                                    32840
      I=L(K)                                                               32850
      IF(I.LE.K) GOTO 120                                                  32860
      JQ=N*(K-1)                                                           32870
      JR=N*(I-1)                                                           32880
      DO 110 J=1,N                                                         32890
      JK=JQ+J                                                              32900
      HOLD=A(JK)                                                           32910
      JI=JR+J                                                              32920
      A(JK)=-A(JI)                                                         32930
  110 A(JI) =HOLD                                                          32940
  120 J=M(K)                                                               32950
      IF(J.LE.K) GOTO 100                                                  32960
  125 KI=K-N                                                               32970
      DO 130 I=1,N                                                         32980
      KI=KI+N                                                              32990
      HOLD=A(KI)                                                           33000
      JI=KI-K+J                                                            33010
      A(KI)=-A(JI)                                                         33020
  130 A(JI) =HOLD                                                          33030
      GO TO 100                                                            33040
      END                                                                  33050
      SUBROUTINE NEKMIN(RMSP,OMEGSP,XSP,ZSP)                               33060
C     VERSION OF 2/19/75. DESTROY EARLIER VERSIONS.                        33070
      IMPLICIT REAL*8(A-H,O-Z)                                             33080
      REAL RMSP,OMEGSP,XSP,ZSP                                             33090
      DIMENSION DN(4),EN(2),OUT(2),LL(2),MM(2)                             33100
      RM=RMSP                                                              33110
      OMEG=OMEGSP                                                          33120
      X=XSP                                                                33130
      Z=ZSP                                                                33140
      Z=.05                                                                33150
   15 P1=X*X+Z*Z                                                           33160
      RP1=DSQRT(P1)                                                        33170
      P115=P1*RP1                                                          33180
      P2=(1.-X)**2+Z*Z                                                     33190
      RP2=DSQRT(P2)                                                        33200
      P215=P2*RP2                                                          33210
      DODZ=-Z/P115-RM*Z/P215                                               33220
      OM=1./RP1+RM/RP2+(1.+RM)*.5*X*X-RM*X                                 33230
      DELOM=OMEG-OM                                                        33240
      DELZ=DELOM/DODZ                                                      33250
      Z=DABS(Z+DELZ)                                                       33260
      ABDELZ=DABS(DELZ)                                                    33270
      IF(ABDELZ.GT..00001) GOTO 15                                         33280
   16 P1=X*X+Z*Z                                                           33290
      RP1=DSQRT(P1)                                                        33300
      P115=P1*RP1                                                          33310
      P125=P1*P115                                                         33320
      P2=(1.-X)**2+Z*Z                                                     33330
      RP2=DSQRT(P2)                                                        33340
      P215=P2*RP2                                                          33350
      P225=P2*P215                                                         33360
      DN(1)=-X/P115+RM*(1.-X)/P215+(1.+RM)*X-RM                            33370
      DN(2)=(3.*X*X-P1)/P125+(3.*RM*(1.-X)**2-RM*((1.-X)**2+Z*Z))/P225     33380
     $+(RM+1.)                                                             33390
      DN(3)=-Z/P115-RM*Z/P215                                              33400
      DN(4)=3.*X*Z/P125-3.*RM*Z*(1.-X)/P225                                33410
      OME=1./RP1+RM/RP2+(1.+RM)*.5*X*X-RM*X                                33420
      EN(1)=OMEG-OME                                                       33430
      EN(2)=-DN(1)                                                         33440
      CALL DMINV(DN,2,D,LL,MM)                                             33450
      CALL DGMPRD(DN,EN,OUT,2,2,1)                                         33460
      DT1=OUT(1)                                                           33470
      DT2=OUT(2)                                                           33480
      ABDX=DABS(DT1)                                                       33490
      X=X+DT1                                                              33500
      ABDZ=DABS(DT2)                                                       33510
      Z=Z+DT2                                                              33520
      IF(ABDX.GT..000001) GOTO 16                                          33530
      IF(ABDZ.GT..000001) GOTO 16                                          33540
      RMSP=RM                                                              33550
      OMEGSP=OMEG                                                          33560
      XSP=X                                                                33570
      ZSP=Z                                                                33580
      RETURN                                                               33590
      END                                                                  33600
      SUBROUTINE SQUARE(OBS,NOBS,ML,OUT,PE,D,CN,CNN,SS,CL,LL,MM)           33610
c  Version of April 9, 1992                                                33620
      DIMENSION OBS(*),OUT(*),PE(*),CN(*),CL(*),LL(*),MM(*),CNN(*)         33630
      DOUBLE PRECISION CN,D,CL,OUT,OBSII,OBSKK,CAY,S,CLL,OBSQQ,CNN         33640
      S=0.D0                                                               33650
      CLL=0.D0                                                             33660
      CAY=NOBS-ML                                                          33670
      JMAX=ML*ML                                                           33680
      DO 20 J=1,JMAX                                                       33690
   20 CN(J)=0.D0                                                           33700
      DO 21 J=1,ML                                                         33710
   21 CL(J)=0.D0                                                           33720
      DO 25 NOB=1,NOBS                                                     33730
      III=NOB+NOBS*ML                                                      33740
      OBSQQ=OBS(III)                                                       33750
      DO 23 K=1,ML                                                         33760
      DO 23 I=1,ML                                                         33770
      II=NOB+NOBS*(I-1)                                                    33780
      KK=NOB+NOBS*(K-1)                                                    33790
      J=I+(K-1)*ML                                                         33800
      OBSII=OBS(II)                                                        33810
      OBSKK=OBS(KK)                                                        33820
   23 CN(J)=CN(J)+OBSII*OBSKK                                              33830
      DO 24 I=1,ML                                                         33840
      II=NOB+NOBS*(I-1)                                                    33850
      OBSII=OBS(II)                                                        33860
   24 CL(I)=CL(I)+OBSQQ*OBSII                                              33870
   25 CLL=CLL+OBSQQ*OBSQQ                                                  33880
      DO 50 J=1,JMAX                                                       33890
   50 CNN(J)=CN(J)                                                         33900
      CALL DMINV(CNN,ML,D,LL,MM)                                           33910
      CALL DGMPRD(CNN,CL,OUT,ML,ML,1)                                      33920
      DO 26 I=1,ML                                                         33930
   26 S=S+CL(I)*OUT(I)                                                     33940
      S=CLL-S                                                              33950
      SS=S                                                                 33960
      SIGSQ=S/CAY                                                          33970
      DO 27 J=1,ML                                                         33980
      JJ=J*ML+J-ML                                                         33990
      CNJJ=CNN(JJ)                                                         34000
      ARG=SIGSQ*CNJJ                                                       34010
   27 PE(J)=.6745*SQRT(ARG)                                                34020
      RETURN                                                               34030
      END                                                                  34040
      Subroutine lum (xlum,x,y,grexp,wl,tpoll,n,n1,nstar,sbr,rv,rvq,       34050
     $glump1,glump2,grv1,grv2,mmsave,summ,fr,sm,ifat,vol,rm,om,f,d,        34060
     $hld,snth)                                                            34070
c   Experimental version of March 31, 1992                                 34080
      dimension rv(*),rvq(*),mmsave(*),fr(*),hld(*),snth(*),glump1(*),     34090
     $glump2(*),grv1(*),grv2(*)                                            34100
      common /radi/ R1H,RLH,R1C,RLC                                        34110
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm,      34120
     $is1dm,is2dm,ic1dm,ic2dm,ld                                           34130
      TPOLE=10000.*TPOLL                                                   34140
      KR=0                                                                 34150
      RAT=1.                                                               34160
      RATPOL=1.                                                            34170
      IF(IFAT.NE.0) CALL ATM(TPOLE,WL,RATPOL)                              34180
      XLUMP=14384./WL                                                      34190
      COMP=3-2*NSTAR                                                       34200
      POLFAC=EXP(XLUMP/TPOLE)-1.                                           34210
      EN=N                                                                 34220
      DELTH=1.570796/EN                                                    34230
      IQ=(NSTAR-1)*(N1+1)                                                  34240
      IS=(NSTAR-1)*MMSAVE(IQ)                                              34250
      SUM=0.                                                               34260
      SUMM=0.                                                              34270
      SM=0.                                                                34280
      VOL=0.                                                               34290
      DO 36 I=1,N                                                          34300
      IPN1=I+N1*(NSTAR-1)                                                  34310
      SINTH=SNTH(IPN1)                                                     34320
      EM=SINTH*EN*1.3                                                      34330
      MM=EM+1.                                                             34340
      XM=MM                                                                34350
      DELFI=3.141593/XM                                                    34360
      DFST=DELFI*SINTH                                                     34370
      SUMJ=0.                                                              34380
      SUMMJ=0.                                                             34390
      SMJ=0.                                                               34400
      VOLJ=0.                                                              34410
      DO 26 J=1,MM                                                         34420
      IS=IS+1                                                              34430
      IP=(NSTAR-1)*(N1+1)+I                                                34440
      IX=MMSAVE(IP)+J                                                      34450
      IF(NSTAR.EQ.1) GOTO 39                                               34460
      IF(RVQ(IX).EQ.-1.) GOTO 25                                           34470
      R=RVQ(IX)                                                            34480
      GOTO 49                                                              34490
   39 IF(RV(IX).EQ.-1.) GOTO 25                                            34500
      R=RV(IX)                                                             34510
   49 grav=(2-nstar)*grv1(ix)+(nstar-1)*grv2(ix)                           34520
      TLOCAL=TPOLE*SQRT(SQRT(GRAV))                                        34530
      IF(IFAT.NE.0) CALL ATM(TLOCAL,WL,RAT)                                34540
      GRAVM=RAT*POLFAC/(EXP(XLUMP/TLOCAL)-1.)                              34550
      di=(2-nstar)*glump1(ix)+(nstar-1)*glump2(ix)                         34560
      DIF=DI*GRAVM                                                         34570
      DIFF=DI*GRAV                                                         34580
      SMJ=SMJ+DI                                                           34590
      SUMJ=SUMJ+DIF                                                        34600
      SUMMJ=SUMMJ+DIFF                                                     34610
      VOLJ=VOLJ+R*R*R*FR(IX)                                               34620
      GOTO 26                                                              34630
   25 KR=1                                                                 34640
   26 CONTINUE                                                             34650
      SMJ=SMJ*DELFI                                                        34660
      SUMJ=SUMJ*DELFI                                                      34670
      SUMMJ=SUMMJ*DELFI                                                    34680
      SM=SM+SMJ                                                            34690
      SUMM=SUMM+SUMMJ                                                      34700
      VOL=VOL+VOLJ*DFST                                                    34710
   36 SUM=SUM+SUMJ                                                         34720
      if(ld.eq.1) darkin=3.141593*(1.-x/3.)                                34730
      if(ld.eq.2) darkin=3.141593*(1.-x/3.+2.*y/9.)                        34740
      SBR=.25*RATPOL*XLUM/(SUM*DELTH*DARKIN)                               34750
      SM=SM*DELTH*4.                                                       34760
      SUMM=SUMM*DELTH*4.                                                   34770
      VOL=VOL*1.333333*DELTH                                               34780
      IF(KR.EQ.0) RETURN                                                   34790
      CALL ELLONE(F,D,RM,XL1,OMD,XLD,OMD)                                  34800
      CALL NEKMIN(RM,OM,XL1,ZD)                                            34810
      IF(NSTAR.EQ.2) XL1=D-XL1                                             34820
      R1=(2-nstar)*R1H+(nstar-1)*R1C                                       34830
      RL=(2-nstar)*RLH+(nstar-1)*RLC                                       34840
      VOL=VOL+1.047198*XL1*R1*RL                                           34850
      RETURN                                                               34860
      END                                                                  34870
      SUBROUTINE LUMP(GRX,GRY,GRZ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,            34880
     $MMSAVE,GREXP,ALB,WL,TPOLL,SBR,SUMM,N1,N2,KOMP,IFAT,fr,snth,          34890
     $CSTH,SNFI,CSFI,TLD,GLUMP1,GLUMP2,XX1,XX2,YY1,YY2,ZZ1,ZZ2,xbol,ybol   34900
     $,GRV1,GRV2,SBR1B,SBR2B,RF,RFO,GMAG1,GMAG2,DINT)                      34910
c   Experimental version of March 31, 1992                                 34920
      DIMENSION GRX(*),GRY(*),GRZ(*),GRXQ(*),GRYQ(*),grzq(*),              34930
     $SLUMP1(*),SLUMP2(*),MMSAVE(*),FR(*),SNTH(*),CSTH(*),                 34940
     $SNFI(*),CSFI(*),TLD(*),GLUMP1(*),GLUMP2(*),XX1(*),XX2(*),YY1(*)      34950
     $,YY2(*),ZZ1(*),ZZ2(*),GRV1(*),GRV2(*),RF(*),RFO(*),                  34960
     $GMAG1(*),GMAG2(*)                                                    34970
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm       34980
     $,ifs1dm,ifs2dm,icr1dm,icr2dm,ld                                      34990
      RTOL=1.                                                              35000
      IQ=(KOMP-1)*(N1+1)                                                   35010
      IS=(KOMP-1)*MMSAVE(IQ)                                               35020
      ATRAT=1.                                                             35030
      RATPOL=1.                                                            35040
      PI=3.141593                                                          35050
      PIH=.5D0*PI                                                          35060
      TPOLE=10000.*TPOLL                                                   35070
      IF(IFAT.NE.0) CALL ATM(tpole,WL,RATPOL)                              35080
      CMPP=2-KOMP                                                          35090
      COMPP=2*KOMP-3                                                       35100
      COMP=-COMPP                                                          35110
      CMP=KOMP-1                                                           35120
      N=(2-KOMP)*N1+(KOMP-1)*N2                                            35130
      NO=(2-KOMP)*N2+(KOMP-1)*N1                                           35140
      NOD=2*NO                                                             35150
      EN=N                                                                 35160
      ENO=NO                                                               35170
      DELTHO=PIH/ENO                                                       35180
      XLUMP=14384./WL                                                      35190
      POLFAC=EXP(XLUMP/TPOLE)-1.                                           35200
      CNST=ALB*DELTHO*SBR2B/(DINT*SBR1B)                                   35210
      IF(KOMP.EQ.2)CNST=ALB*DELTHO*SBR1B/(DINT*SBR2B)                      35220
      DO 191 I=1,N                                                         35230
      IPN1=I+N1*(KOMP-1)                                                   35240
      SINTH=SNTH(IPN1)                                                     35250
      EM=SINTH*EN*1.3                                                      35260
      MM=EM+1.                                                             35270
      IP=(KOMP-1)*(N1+1)+I                                                 35280
      IY=MMSAVE(IP)                                                        35290
      DO 193 J=1,MM                                                        35300
      IX=IY+J                                                              35310
      SUM=0.                                                               35320
      IF(FR(IX).EQ.0.) GOTO 194                                            35330
      DO 190 IOTH=1,NOD                                                    35340
      IOTHS=IOTH                                                           35350
      IF(IOTH.GT.NO) IOTHS=NOD-IOTH+1                                      35360
      IPNO=IOTHS+N1*(2-KOMP)                                               35370
      SINTHO=SNTH(IPNO)                                                    35380
      EMO=SINTHO*ENO*1.3                                                   35390
      MMO=EMO+1.                                                           35400
      MMOD=2*MMO                                                           35410
      IPO=(2-KOMP)*(N1+1)+IOTHS                                            35420
      IYO=MMSAVE(IPO)                                                      35430
      XMO=MMO                                                              35440
      DELFIO=PI/XMO                                                        35450
      DO 190 JOFI=1,MMOD                                                   35460
      JOFU=JOFI                                                            35470
      IF(JOFI.GT.MMO) JOFU=MMOD-JOFI+1                                     35480
      IXO=IYO+JOFU                                                         35490
      IX1=IX                                                               35500
      IX2=IXO                                                              35510
      IF(KOMP.EQ.1) GOTO 200                                               35520
      IF(GLUMP1(IXO).EQ.0.) GOTO 181                                       35530
      IX1=IXO                                                              35540
      IX2=IX                                                               35550
      GOTO 201                                                             35560
  200 CONTINUE                                                             35570
      IF(GLUMP2(IXO).EQ.0.) GOTO 179                                       35580
  201 RTL1=1.                                                              35590
      RTL2=1.                                                              35600
      UPD1=1.                                                              35610
      UPD2=1.                                                              35620
      IF(KOMP.EQ.2) GOTO 22                                                35630
      IF(JOFI.GT.MMO) RTL2=-1.                                             35640
      IF(IOTH.GT.NO) UPD2=-1.                                              35650
      GOTO 23                                                              35660
   22 IF(JOFI.GT.MMO) RTL1=-1.                                             35670
      IF(IOTH.GT.NO) UPD1=-1.                                              35680
   23 CONTINUE                                                             35690
      GX2=GRXQ(IX2)                                                        35700
      GY2=GRYQ(IX2)*RTL2                                                   35710
      GZ2=GRZQ(IX2)*UPD2                                                   35720
      X1C=XX1(IX1)                                                         35730
      X2C=XX2(IX2)                                                         35740
      Y1C=YY1(IX1)*RTL1                                                    35750
      Y2C=YY2(IX2)*RTL2                                                    35760
      Z1C=ZZ1(IX1)*UPD1                                                    35770
      Z2C=ZZ2(IX2)*UPD2                                                    35780
      DX=(X2C-X1C)*COMP                                                    35790
      DY=(Y2C-Y1C)*COMP                                                    35800
      DZ=(Z2C-Z1C)*COMP                                                    35810
      DLRSQ=DX*DX+DY*DY+DZ*DZ                                              35820
      CSNUM2=(DX*GX2+DY*GY2+DZ*GZ2)*COMPP                                  35830
      IF(CSNUM2.GE.0.) GOTO 190                                            35840
      GX1=GRX(IX1)                                                         35850
      GY1=GRY(IX1)*RTL1                                                    35860
      GZ1=GRZ(IX1)*UPD1                                                    35870
      CSNUM1=(DX*GX1+DY*GY1+DZ*GZ1)*COMP                                   35880
      IF(CSNUM1.GE.0.) GOTO 190                                            35890
      DMAG=SQRT(DLRSQ)                                                     35900
      CSGM1=-CSNUM1/(DMAG*GMAG1(IX1))                                      35910
      CSGM2=-CSNUM2/(DMAG*GMAG2(IX2))                                      35920
      IF(KOMP.EQ.2) GOTO 181                                               35930
      DGAM2=1.-XBOL+XBOL*CSGM2                                             35940
      if(ld.eq.1) goto 179                                                 35950
      if(csgm2.eq.0.) goto 179                                             35960
      dgam2=dgam2-ybol*csgm2*alog(csgm2)                                   35970
      if(dgam2.lt.0.) dgam2=0.                                             35980
  179 DSUM=GRV2(IXO)*GLUMP2(IXO)*RFO(IXO)*CSGM1*CSGM2*DGAM2/DLRSQ          35990
      GOTO 182                                                             36000
  181 DGAM1=1.-XBOL+XBOL*CSGM1                                             36010
      if(ld.eq.1) goto 184                                                 36020
      if(csgm1.eq.0.) goto 184                                             36030
      dgam1=dgam1-ybol*csgm1*alog(csgm1)                                   36040
      if(dgam1.lt.0.) dgam1=0.                                             36050
  184 DSUM=GRV1(IXO)*GLUMP1(IXO)*RFO(IXO)*CSGM2*CSGM1*DGAM1/DLRSQ          36060
  182 CONTINUE                                                             36070
      SUM=SUM+DSUM*DELFIO                                                  36080
  190 CONTINUE                                                             36090
  194 RF(IX)=(CNST*SUM/(CMPP*GRV1(IX)+CMP*GRV2(IX)))+1.                    36100
      IF(RF(IX).LE.RTOL) GOTO 191                                          36110
  193 CONTINUE                                                             36120
  191 CONTINUE                                                             36130
      DO 8 I=1,N                                                           36140
      IPN1=I+N1*(KOMP-1)                                                   36150
      SINTH=SNTH(IPN1)                                                     36160
      EM=SINTH*EN*1.3                                                      36170
      MM=EM+1.                                                             36180
      XM=MM                                                                36190
      DELFI=PI/XM                                                          36200
      IP=(KOMP-1)*(N1+1)+I                                                 36210
      IY=MMSAVE(IP)                                                        36220
      DO 8 J=1,MM                                                          36230
      IS=IS+1                                                              36240
      IX=IY+J                                                              36250
      IF(FR(IX).EQ.0.) GOTO 8                                              36260
      IF(KOMP.EQ.1) GRV=GRV1(IX)                                           36270
      IF(KOMP.EQ.2) GRV=GRV2(IX)                                           36280
      TNEW=TPOLE*SQRT(SQRT(GRV*RF(IX)))                                    36290
      TLD(IS)=TNEW                                                         36300
      IF(IFAT.EQ.0) GOTO 18                                                36310
      CALL ATM(TNEW,WL,RATNEW)                                             36320
      ATRAT=RATNEW/RATPOL                                                  36330
   18 GRREFL=POLFAC*ATRAT/(EXP(XLUMP/TNEW)-1.)                             36340
   37 IF(KOMP.EQ.1) GOTO 77                                                36350
      slump2(ix)=glump2(ix)*grrefl*sbr                                     36360
      GOTO 8                                                               36370
   77 slump1(ix)=glump1(ix)*grrefl*sbr                                     36380
    8 CONTINUE                                                             36390
      RETURN                                                               36400
      END                                                                  36410
      FUNCTION ARSIN(X)                                                    36420
      ARSIN=ASIN(X)                                                        36430
      RETURN                                                               36440
      END                                                                  36450
      FUNCTION ARCOS(X)                                                    36460
      ARCOS=ACOS(X)                                                        36470
      RETURN                                                               36480
      END                                                                  36490
      FUNCTION DARSIN(X)                                                   36500
      IMPLICIT REAL*8 (A-H,O-Z)                                            36510
      DARSIN=DASIN(X)                                                      36520
      RETURN                                                               36530
      END                                                                  36540
      SUBROUTINE ROMQSP(OME,Q,F,D,EC,TH,FI,R,DRDO,DRDQ,DODQ,KOMP,MODE)     36550
c  Version of May 27, 1991. Destroy earlier versions.                      36560
      PI=3.141593                                                          36570
      theq=1.570796                                                        36580
      MOD46=(MODE-5)**2                                                    36590
      MOD56=(2*MODE-11)**2                                                 36600
      DQ=1.E-4*Q                                                           36610
      QP=Q+DQ                                                              36620
      TOL=1.E-6                                                            36630
C     TH, FI SHOULD BE IN RADIANS.                                         36640
      XLAM=SIN(TH)*COS(FI)                                                 36650
      XNU=COS(TH)                                                          36660
      XNUSQ=1.-XNU*XNU                                                     36670
      RMA=Q                                                                36680
      QF=1.                                                                36690
      DP=1.-EC                                                             36700
      QFM=1.                                                               36710
      IF(KOMP.NE.2) GOTO 23                                                36720
      RMA=1./Q                                                             36730
      QF=1./Q                                                              36740
      QFM=-1./Q**2                                                         36750
   23 CONTINUE                                                             36760
      CALL ELLONE(F,DP,RMA,X,OMEG,XLD,OMD)                                 36770
      OM2SAV=OMEG                                                          36780
      RMAP=QP                                                              36790
      IF(KOMP.NE.2) GOTO 92                                                36800
      OMEG=OMEG*Q+(1.-Q)*.5                                                36810
      IF(MOD56.EQ.1) OME=OMEG                                              36820
      RMAP=1./QP                                                           36830
      GOTO 93                                                              36840
   92 CONTINUE                                                             36850
      IF(MOD46.EQ.1) OME=OMEG                                              36860
   93 CONTINUE                                                             36870
      POT=OME                                                              36880
      IF(KOMP.EQ.2) POT=OME/Q+.5*(Q-1.)/Q                                  36890
      CALL ELLONE(F,DP,RMAP,XP,OMP,XLD,OMD)                                36900
      DODQ=(OMP-OM2SAV)/DQ                                                 36910
      RM1=RMA+1.                                                           36920
      DS=D*D                                                               36930
      RF=F*F                                                               36940
      R=1./POT                                                             36950
      KOUNT=0                                                              36960
      DELR=0.                                                              36970
      IF(FI.NE.0.) GOTO 85                                                 36980
      IF(TH.NE.THEQ) GOTO 85                                               36990
      IF(MODE.EQ.6) GOTO 115                                               37000
      IF(MODE.NE.4) GOTO 80                                                37010
      IF(KOMP.EQ.1) GOTO 115                                               37020
      GOTO 85                                                              37030
   80 IF(MODE.NE.5) GOTO 85                                                37040
      IF(KOMP.EQ.2) GOTO 115                                               37050
   85 CONTINUE                                                             37060
   14 R=R+DELR                                                             37070
      KOUNT=KOUNT+1                                                        37080
      IF(KOUNT.LT.20) GOTO 70                                              37090
  217 DOMR=-1.E25                                                          37100
      R=-1.                                                                37110
      GOTO 116                                                             37120
   70 RSQ=R*R                                                              37130
      PAR=DS-2.*XLAM*R*D+RSQ                                               37140
      RPAR=SQRT(PAR)                                                       37150
      OM=1./R+RMA*(1./RPAR-XLAM*R/DS)+RM1*.5*RSQ*XNUSQ*RF                  37160
      DOMR=1./(RF*RM1*XNUSQ*R-1./RSQ-(RMA*(R-XLAM*D))/(PAR*RPAR)-          37170
     $RMA*XLAM/DS)                                                         37180
      DELR=(POT-OM)*DOMR                                                   37190
      ABDELR=ABS(DELR)                                                     37200
      IF(ABDELR.GT.TOL) GOTO 14                                            37210
      DOMRSV=DOMR                                                          37220
      IF(R.GE.1.) GOTO 217                                                 37230
      IF(FI.NE.0.) GO TO 116                                               37240
      IF(TH.NE.THEQ)GO TO 116                                              37250
      IF(OME-OMEG.LT.0) GOTO 217                                           37260
      IF(OME-OMEG.EQ.0) GOTO 114                                           37260
      IF(OME-OMEG.GT.0) GOTO 116                                           37260
  114 DOMR=1.E15                                                           37270
  115 R=X                                                                  37280
  116 DRDO=DOMR*QF                                                         37290
  118 CONTINUE                                                             37300
      DRDQ=(1./RPAR-R*XLAM/DS+.5*RF*RSQ*XNUSQ)/(1./RSQ+RMA*((1.            37310
     $/(PAR*RPAR))*(R-XLAM*D)+XLAM/DS)-RF*XNUSQ*RM1*R)                     37320
      DRDQ=DRDQ*QFM                                                        37330
      IF(MODE.EQ.6) GOTO 215                                               37340
      IF(MODE.NE.4) GOTO 180                                               37350
      IF(KOMP.EQ.1) GOTO 215                                               37360
      RETURN                                                               37370
  180 IF(MODE.NE.5) RETURN                                                 37380
      IF(KOMP.EQ.2) GOTO 215                                               37390
      RETURN                                                               37400
  215 IF(FI.NE.0.) GOTO 230                                                37410
      IF(TH.NE.THEQ) GOTO 230                                              37420
      DRDQ=(XP-X)/DQ                                                       37430
      RETURN                                                               37440
  230 DRDQ=DRDQ+DOMRSV*DODQ                                                37450
      RETURN                                                               37460
      END                                                                  37470
      SUBROUTINE SPOT(KOMP,N,SINTH,COSTH,SINFI,COSFI,TEMF)                 37480
C                                                                          37490
c   If a surface point is in more than one spot, this subroutine           37500
c      adopts the product of the spot temperature factors.                 37510
C                                                                          37520
c   "Latitudes" here are actually co-latitudes (running from 0 at one      37530
c      pole to 180 deg. at the other.                                      37540
C                                                                          37550
c   Version of May 20, 1991                                                37560
C                                                                          37570
      COMMON /SPOTS/ SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG      37580
     $(2,100),RAD(2,100),TEMSP(2,100),DUMMY(2,100)                         37590
      TEMF=1.                                                              37600
      DO 15 I=1,N                                                          37610
      COSDFI=COSFI*COSLNG(KOMP,I)+SINFI*SINLNG(KOMP,I)                     37620
      S=ARCOS(COSTH*COSLAT(KOMP,I)+SINTH*SINLAT(KOMP,I)*COSDFI)*57.29578   37630
      IF(S.GT.RAD(KOMP,I)) GOTO 15                                         37640
      TEMF=TEMF*TEMSP(KOMP,I)                                              37650
   15 CONTINUE                                                             37660
      RETURN                                                               37670
      END                                                                  37680
      SUBROUTINE OLUMP(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,SLUMP1,SLUMP2,    37690
     $MMSAVE,GREXP,ALB,RB,WL,TPOLL,SBR,SUMM,N1,N2,KOMP,IFAT,FR,x,y,D,      37700
     $SNTH,CSTH,SNFI,CSFI,tld,glump1,glump2)                               37710
c   Experimental version of March 31, 1992.                                37720
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*   37730
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),FR(*),F(3),W(3),SNTH(*),CSTH(*),     37740
     $SNFI(*),CSFI(*),tld(*),glump1(*),glump2(*)                           37750
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm,      37760
     $ifs1dm,ifs2dm,icr1dm,icr2dm,ld                                       37770
      IQ=(KOMP-1)*(N1+1)                                                   37780
      IS=(KOMP-1)*MMSAVE(IQ)                                               37790
      FP=7.957747E-2                                                       37800
      ATRAT=1.                                                             37810
      RATPOL=1.                                                            37820
      RAT=1.                                                               37830
      PI=3.141593                                                          37840
      PIH=1.570796                                                         37850
      PI32=4.712389                                                        37860
      F(1)=.1127017                                                        37870
      F(2)=.5                                                              37880
      F(3)=.8872983                                                        37890
      W(1)=.2777778                                                        37900
      W(2)=.4444444                                                        37910
      W(3)=.2777778                                                        37920
      TPOLE=10000.*TPOLL                                                   37930
      IF(IFAT.NE.0) CALL ATM(tpole,WL,RATPOL)                              37940
      CMPP=2-KOMP                                                          37950
      COMPP=2*KOMP-3                                                       37960
      COMP=-COMPP                                                          37970
      CMP=KOMP-1                                                           37980
      CMPD=CMP*D                                                           37990
      CMPPD=CMPP*D                                                         38000
      N=(2-KOMP)*N1+(KOMP-1)*N2                                            38010
      ENN=(15.+X)*(1.+GREXP)/(15.-5.*X)                                    38020
      NP=N1+1+(2-KOMP)*(N2+1)                                              38030
      NPP=N1*(KOMP-1)+(NP-1)*(2-KOMP)                                      38040
      LL=MMSAVE(NPP)+1                                                     38050
      LLL=MMSAVE(NP)                                                       38060
      LLLL=(LL+LLL)/2                                                      38070
      AR=RV(LLL)*CMP+RVQ(LLL)*CMPP                                         38080
      BR=RV(LLLL)*CMP+RVQ(LLLL)*CMPP                                       38090
      CR=RV(1)*CMP+RVQ(1)*CMPP                                             38100
      BOA=BR/AR                                                            38110
      BOAL=1.-BOA*BOA                                                      38120
      BOC2=(BR/CR)**2                                                      38130
      CC=1./(1.-.25*ENN*(1.-BOA**2)*(.9675-.3008*BOA))                     38140
      HCN=.5*CC*ENN                                                        38150
      DF=1.-X/3.                                                           38160
      if(ld.eq.2) df=1.-x/3.+2.*y/9.                                       38170
      XLUMP=14384./WL                                                      38180
      POLFAC=EXP(XLUMP/TPOLE)-1.                                           38190
      IF(KOMP.EQ.1) GOTO 82                                                38200
      GRPOLE=SQRT(GRXQ(1)**2+GRYQ(1)**2+GRZQ(1)**2)                        38210
      GOTO 83                                                              38220
   82 GRPOLE=SQRT(GRX(1)**2+GRY(1)**2+GRZ(1)**2)                           38230
   83 EN=N                                                                 38240
      DO 8 I=1,N                                                           38250
      IPN1=I+N1*(KOMP-1)                                                   38260
      SINTH=SNTH(IPN1)                                                     38270
      COSTH=CSTH(IPN1)                                                     38280
      EM=SINTH*EN*1.3                                                      38290
      MM=EM+1.                                                             38300
      XM=MM                                                                38310
      DELFI=3.141593/XM                                                    38320
      IP=(KOMP-1)*(N1+1)+I                                                 38330
      IY=MMSAVE(IP)                                                        38340
      DO 8 J=1,MM                                                          38350
      IS=IS+1                                                              38360
      STCF=SINTH*CSFI(IS)                                                  38370
      STSF=SINTH*SNFI(IS)                                                  38380
      IX=IY+J                                                              38390
      IF(KOMP.EQ.1) GOTO 39                                                38400
      IF(RVQ(IX).EQ.-1.) GOTO 8                                            38410
      GX=GRXQ(IX)                                                          38420
      GY=GRYQ(IX)                                                          38430
      GZ=GRZQ(IX)                                                          38440
      R=RVQ(IX)                                                            38450
      GOTO 49                                                              38460
   39 IF(RV(IX).EQ.-1.)GOTO 8                                              38470
      GX=GRX(IX)                                                           38480
      GY=GRY(IX)                                                           38490
      GZ=GRZ(IX)                                                           38500
      R=RV(IX)                                                             38510
   49 GRMAG=SQRT(GX*GX+GY*GY+GZ*GZ)                                        38520
      ZZ=R*COSTH                                                           38530
      YY=R*COMP*STSF                                                       38540
      XX=CMPD+COMP*STCF*R                                                  38550
      XXREF=(CMPPD+COMPP*XX)*COMPP                                         38560
      GRAV=(GRMAG/GRPOLE)**GREXP                                           38570
      TLOCAL=TPOLE*SQRT(SQRT(GRAV))                                        38580
      DIST=SQRT(XXREF*XXREF+YY*YY+ZZ*ZZ)                                   38590
      RMX=ARSIN(.5*(BR+CR)/DIST)                                           38600
      XCOS=XXREF/DIST                                                      38610
      YCOS=YY/DIST                                                         38620
      ZCOS=ZZ/DIST                                                         38630
      COSINE=(XCOS*GX+YCOS*GY+ZCOS*GZ)/GRMAG                               38640
      RC=PIH-ARCOS(COSINE)                                                 38650
      AH=RC/RMX                                                            38660
      RP=ABS(AH)                                                           38670
      IF(AH.LE..99999) GOTO 22                                             38680
      P=1.                                                                 38690
      GOTO 16                                                              38700
   22 IF(AH.GE.-.99999) GOTO 24                                            38710
      ALBEP=0.                                                             38720
      GOTO 19                                                              38730
   24 SUM=0.                                                               38740
      FIST=ARSIN(RP)                                                       38750
      FII=PIH-FIST                                                         38760
      DO 15 IT=1,3                                                         38770
      FE=FII*F(IT)+FIST                                                    38780
      PAR=1.-(RP/SIN(FE))**2                                               38790
      RPAR=SQRT(PAR)                                                       38800
      SUM=PAR*RPAR*W(IT)+SUM                                               38810
   15 CONTINUE                                                             38820
      FTRI=(1.-X)*RP*SQRT(1.-RP**2)+.6666667*X*FII-.6666667*X*SUM*FII      38830
      FSEC=(PIH+FIST)*DF                                                   38840
      P=(FTRI+FSEC)/(PI*DF)                                                38850
      IF(COSINE.LT.0.) P=1.-P                                              38860
      RTF=SQRT(1.-AH**2)                                                   38870
      DENO=PI32-3.*(AH*RTF+ARSIN(AH))                                      38880
      IF(DENO.NE.0.) GOTO 117                                              38890
      ABAR=1.                                                              38900
      GOTO 116                                                             38910
  117 ABAR=2.*RTF**3/DENO                                                  38920
  116 COSINE=COS(PIH-RMX*ABAR)                                             38930
   16 COSQ=1./(1.+(YY/XXREF)**2)                                           38940
      COT2=(ZZ/XXREF)**2                                                   38950
      Z=BOAL/(1.+BOC2*COT2)                                                38960
      E=CC-HCN*COSQ*Z                                                      38970
      ALBEP=ALB*E*P                                                        38980
   19 IF(COSINE.LE.0.) ALBEP=0.                                            38990
      TNEW=TLOCAL*SQRT(SQRT(1.+(FP*SUMM/(DIST*DIST*GRAV))*COSINE*RB*       39000
     $ALBEP))                                                              39010
      TLD(IS)=TNEW                                                         39020
      IF(IFAT.EQ.0) GOTO 18                                                39030
      CALL ATM(TNEW,WL,RATNEW)                                             39040
      ATRAT=RATNEW/RATPOL                                                  39050
   18 GRREFL=POLFAC*ATRAT/(EXP(XLUMP/TNEW)-1.)                             39060
   37 IF(KOMP.EQ.1) GOTO 77                                                39070
      slump2(ix)=glump2(ix)*grrefl*sbr                                     39080
      GOTO 8                                                               39090
   77 slump1(ix)=glump1(ix)*grrefl*sbr                                     39100
    8 CONTINUE                                                             39110
      RETURN                                                               39120
      END                                                                  39130
      SUBROUTINE BOLO (RATLUM,T1,T2,WL,IFAT1,IFAT2,RATBOL)                 39140
C     VERSION OF 3/29/77. DESTROY EARLIER VERSIONS.                        39150
C     BOLO USES HARRIS CALIBRATION FROM T=3970 TO 5800, MORTON-ADAMS       39160
C     FROM T=5800 TO 37500, AND BLACK BODY LAWS BELOW 3970 AND ABOVE       39170
C     37500.                                                               39180
      VT1=1.                                                               39190
      VT2=1.                                                               39200
      WLT1=1.                                                              39210
      WLT2=1.                                                              39220
      TP1=10000.*T1                                                        39230
      TP2=10000.*T2                                                        39240
      IF(IFAT1.EQ.0) GOTO 95                                               39250
      CALL ATM(TP1,.5500,VT1)                                              39260
      CALL ATM(TP1,WL,WLT1)                                                39270
   95 IF(IFAT2.EQ.0) GOTO 90                                               39280
      CALL ATM(TP2,.5500,VT2)                                              39290
      CALL ATM(TP2,WL,WLT2)                                                39300
   90 RAT12=VT1*WLT2/(VT2*WLT1)                                            39310
      RATV=RAT12*RATLUM*(EXP(1.4384/(WL*T1))-1.)*(EXP(1.4384/(.5500*T2))   39320
     $-1.)/((EXP(1.4384/(.5500*T1))-1.)*(EXP(1.4384/(WL*T2))-1.))          39330
      DELMV=1.0857362*ALOG(RATV)                                           39340
      X1=ALOG10(T1)+.14                                                    39350
      X2=ALOG10(T2)+.14                                                    39360
      IF(T1.GT.3.75) GOTO 22                                               39370
      IF(T1.LT..397) GOTO 22                                               39380
      IF(T1.GT.1.07) GOTO 18                                               39390
      BC1=.0004964+X1*(.184549+X1*(-7.25068+X1*(-26.4889+X1*(-198.069+X1   39400
     $*(-73.801+652.95*X1)))))                                             39410
      GOTO 19                                                              39420
   18 BC1=.868864+X1*(-14.8735+X1*(70.1318+X1*(-227.706+X1*(386.7+X1*(     39430
     $-339.447+121.8*X1)))))                                               39440
      GOTO 19                                                              39450
   22 BC1=-2.5*ALOG10((T1/.6700)**4*(EXP(1.4384/(.5450*T1))-1.)/(EXP(      39460
     $1.4384/.36515)-1.))                                                  39470
   19 IF(T2.GT.3.75) GOTO 32                                               39480
      IF(T2.LT..397) GOTO 32                                               39490
      IF(T2.GT.1.07) GOTO 28                                               39500
      BC2=.0004964+X2*(.184549+X2*(-7.25068+X2*(-26.4889+X2*(-198.069+X2   39510
     $*(-73.801+652.95*X2)))))                                             39520
      GOTO 29                                                              39530
   28 BC2=.868864+X2*(-14.8735+X2*(70.1318+X2*(-227.706+X2*(386.7+X2*(     39540
     $-339.447+121.8*X2)))))                                               39550
      GOTO 29                                                              39560
   32 BC2=-2.5*ALOG10((T2/.6700)**4*(EXP(1.4384/(.5450*T2))-1.)/(EXP(      39570
     $1.4384/.36515)-1.))                                                  39580
   29 DELMB=DELMV+BC2-BC1                                                  39590
      RATBOL=10.**(.4*DELMB)                                               39600
      RETURN                                                               39610
      END                                                                  39620
      subroutine mlrg(a,p,q,r1,r2,t1,t2,sm1,sm2,sr1,sr2,bolm1,             39630
     $bolm2,xlg1,xlg2)                                                     39640
c  This subroutine computes absolute dimensions and other quantities       39650
c  for the stars of a binary star system.                                  39660
c  a = orbital semi-major axis, the sum of the two a's for the two         39670
c  stars. The unit is a solar radius.                                      39680
c  r1,r2 = relative mean (equivalent sphere) radii for stars 1 and 2. Th   39690
c  unit is the orbital semimajor axis.                                     39700
c  p = orbit period in days.                                               39710
c  q = mass ratio, m2/m1.                                                  39720
c  t1,t2= flux-weighted mean surface temperatures for stars 1 and 2,in K   39730
c  sm1,sm2= masses of stars 1 and 2 in solar units.                        39740
c  sr1,sr2= mean radii of stars 1 and 2 in solar units.                    39750
c  bolm1, bolm2= absolute bolometric magnitudes of stars 1, 2.             39760
c  xlg1, xlg2= log (base 10) of mean surface acceleration (effective gra   39770
c  for stars 1 and 2.                                                      39780
c                                                                          39790
      G=6.668e-8                                                           39800
      tsun=5800.                                                           39810
      rsunau=214.8                                                         39820
      sunmas=1.991e33                                                      39830
      sunrad=6.960e10                                                      39840
      gmr=G*sunmas/sunrad**2                                               39850
      sunmb=4.77                                                           39860
      sr1=r1*a                                                             39870
      sr2=r2*a                                                             39880
      yrsid=365.2564                                                       39890
      tmass=(a/rsunau)**3/(p/yrsid)**2                                     39900
      sm1=tmass/(1.+q)                                                     39910
      sm2=tmass*q/(1.+q)                                                   39920
      bol1=(t1/tsun)**4*sr1**2                                             39930
      bol2=(t2/tsun)**4*sr2**2                                             39940
      bolm1=sunmb-2.5*alog10(bol1)                                         39950
      bolm2=sunmb-2.5*alog10(bol2)                                         39960
      xlg1=alog10(gmr*sm1/sr1**2)                                          39970
      xlg2=alog10(gmr*sm2/sr2**2)                                          39980
      return                                                               39990
      end                                                                  40000
