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