program dc c This is the Differential Corrections Main Program. c c Version of October 12, 1998 C C PARAMETER NUMBER 9 IS A, THE RELATIVE ORBITAL SEMI-MAJOR AXIS, IF C SIMULTANEOUS LIGHT AND VELOCITY SOLUTIONS ARE BEING DONE. HOWEVER, C IF ONLY VELOCITY CURVES ARE BEING SOLVED, PARAMETER 9 WILL C EFFECTIVELY BE A*SIN(I), PROVIDED AN INCLINATION OF 90 DEG. IS C ENTERED. CONCEIVABLY,IN SOME RARE SITUATIONS IT MAY BE POSSIBLE TO C FIND A AND I SEPARATELY, USING ONLY VELOCITY CURVES. THIS COULD BE C THE CASE IF THE VELOCITY PROXIMITY EFFECTS ARE IMPORTANT. C C ALTHOUGH THE ARGUMENT OF PERIASTRON IS ENTERED AND PRINTED IN C DEGREES, THE COMPUTED CORRECTION IS IN RADIANS. C C OTHER PROGRAM UNITS: ORBITAL S-M AXIS IN SOLAR RADII (6.96d5 KM), C PERIOD IN DAYS, PHASE IN 2 PI RADIANS, SYSTEMIC VELOCITY AND C THIRD LIGHT IN SAME UNITS AS VELOCITY AND LIGHT OBSERVATIONS, C INCLINATION IN DEGREES, TEMPERATURES IN 10000K., SPOT LATITUDES C IN DEGREES (0=NORTH POLE, 180 DEG.=SOUTH POLE), SPOT LONGITUDES C IN DEGREES (0=LINE OF CENTERS MERIDIAN, INCREASING COUNTER- C CLOCKWISE AS SEEN FROM NORTH POLE TO 360 DEGREES), C SPOT ANGULAR RADII IN DEGREES. SPOT TEMPERATURE FACTOR IS C DIMENSIONLESS. C implicit real*8 (a-h,o-z) dimension phas(3000),flux(3000),wt(3000),br(3000),bl(3000), $obs(60000),hold(60000),cn(2500),cnn(2500),out(50),sd(50),ccl(50), $ll(50),mm(50),hla(25),cla(25),x1a(25),x2a(25),y1a(25),y2a(25), $el3a(25),wla(25),noise(25),sigma(25),knobs(27),mmsavh(124), $mmsavl(124),theta(260),rho(260),aa(20),bb(20),hld(800),v(60), $cnout(2025),snthh(65),csthh(65),snfih(1600),csfih(1600), $snthl(65),csthl(65),snfil(1600),csfil(1600),tldh(1600),tldl(1600) $,opsfa(25),phjd(3000),dfdph(3000),dfdap(3000),para(500) dimension rv(0762),grx(0762),gry(0762),grz(0762),rvq(0762), $grxq(0762),gryq(0762),grzq(0762),slump1(0762),slump2(0762), $srv(0762),sgrx(0762),sgry(0762),sgrz(0762),srvq(0762),sgrxq(0762) $,sgryq(0762),sgrzq(0762),srvl(0762),sgrxl(0762),sgryl(0762), $sgrzl(0762),srvql(0762),sgrxql(0762),sgryql(0762),sgrzql(0762), $slmp1(0762),slmp2(0762),slmp1l(0762),slmp2l(0762),fr1(0762), $fr2(0762),glump1(0762),glump2(0762),grv1(0762),grv2(0762), $xx1(0762),xx2(0762),yy1(0762),yy2(0762),zz1(0762),zz2(0762), $gmag1(0762),gmag2(0762),csbt1(0762),csbt2(0762),rf1(0762), $rf2(0762),rftemp(0762),sxx1(0762),sxx2(0762),syy1(0762), $syy2(0762),szz1(0762),szz2(0762),sgmg1(0762),sgmg2(0762), $sgrv1(0762),sgrv2(0762),sglm1(0762),sglm2(0762),scsb1(0762), $scsb2(0762),srf1(0762),srf2(0762),sglm1l(0762),sglm2l(0762), $sgrv1l(0762),sgrv2l(0762),sxx1l(0762),sxx2l(0762),syy1l(0762), $syy2l(0762),szz1l(0762),szz2l(0762),sgmg1l(0762),sgmg2l(0762), $scsb1l(0762),scsb2l(0762),srf1l(0762),srf2l(0762),yskp(1), $zskp(1) dimension erv(0762),egrx(0762),egry(0762),egrz(0762), $elmp1(0762),eglm1(0762),egrv1(0762),exx1(0762),eyy1(0762), $ezz1(0762),egmg1(0762),ecsb1(0762),erf1(0762),ervq(0762), $egrxq(0762),egryq(0762),egrzq(0762),elmp2(0762),eglm2(0762), $egrv2(0762),exx2(0762),eyy2(0762),ezz2(0762),egmg2(0762), $ecsb2(0762),erf2(0762), $ervl(0762),egrxl(0762),egryl(0762),egrzl(0762),elmp1l(0762), $eglm1l(0762),egrv1l(0762),exx1l(0762),eyy1l(0762),ezz1l(0762), $egmg1l(0762),ecsb1l(0762),erf1l(0762),ervql(0762),egrxql(0762), $egryql(0762),egrzql(0762),elmp2l(0762),eglm2l(0762), $egrv2l(0762),exx2l(0762),eyy2l(0762),ezz2l(0762),egmg2l(0762), $ecsb2l(0762),erf2l(0762),sfr1(0762),sfr1l(0762),efr1(0762), $efr1l(0762),sfr2(0762),sfr2l(0762),efr2(0762),efr2l(0762), $stldh(1600),stldl(1600),etldh(1600),etldl(1600) dimension cnc(2500),clc(50) dimension del(35),keep(36),kep(35),nshift(36),low(35),dlif(35), $xtha(4),xfia(4),arad(4),xlat(2,100),xlong(2,100),radsp(2,100) $,temsp(2,100),po(2),omcr(2) dimension xcl(100),ycl(100),zcl(100),rcl(100),op1(100),fcl(100), $dens(100),encl(100),edens(100),xmue(100) c c The following dimensioned variables are not used by DC. They are c dimensioned only for compatibility with usage of subroutine c LIGHT by program LC. c dimension fbin1(1),fbin2(1),delv1(1),delv2(1),count1(1), $count2(1),delwl1(1),delwl2(1),resf1(1),resf2(1),wl1(1),wl2(1), $dvks1(1),dvks2(1),tau1(1),tau2(1),taug(1),hbarw1(1),hbarw2(1) COMMON /ECCEN/ EC,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC COMMON /DPDX/ DPDX1,DPDX2,POT1,POT2 COMMON /SUMM/ SUMM1,SUMM2 common /ardot/ dperdt,hjd,hjd0,perdum COMMON /FLVAR/ PSH,DP,DSDUM,EF,EFC,ECOS,pert,PHPER,PCONJ, $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc,qfacd COMMON /INVAR/ KH,IPB,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,ifsmv2, $icor1,icor2,ld,ncl,jdphs,ipc 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) COMMON /NSPT/ NSP1,NSP2,IFAT1,IFAT2 common /cld/ acm,opsf common /ipro/ nbins,nl,inmax,inmin,nf1,nf2 DATA ARAD(1),ARAD(2),ARAD(3),ARAD(4)/4HPOLE,5HPOINT,4HSIDE,4HBACK/ 15 FORMAT(1X,16(F11.5)) 16 FORMAT(1X,18(F7.4)) 67 FORMAT(20A4) 99 FORMAT(' SPECIFIED ECLIPSE DURATION INCONSISTENT WITH OTHER PARAME $TERS') 17 FORMAT(1X,22(F6.3)) 19 FORMAT(1X,26(F5.2)) 21 FORMAT(f27.16,f36.16,d22.6) 55 FORMAT(10(3X,d8.1)) 56 FORMAT(10(1X,d7.1)) 20 format(1x,2(4i1,1x),7i1,1x,4(5i1,1x),i1,1x,i1,1x,i1,d10.3) 102 FORMAT('1') 101 FORMAT(' ') 1 FORMAT(I3,I6,I6,I7,I6,I4,I4,I4,f15.6,d13.5,f10.5,f16.3,f14.4) 701 FORMAT(4I2,4I3,f13.6,d12.5,F7.5,F9.3) 2 FORMAT(5(F14.5,F8.4,F5.1)) 85 FORMAT(f10.6,2F10.5,4(1X,F5.3),f8.4,d10.3,i6,d14.5) 18 format(f10.6,2f10.5,4f6.3,f8.4,d10.3,i2,d12.5) 218 FORMAT(F9.6,2F10.5,4F6.3,d10.3,d12.5) 37 FORMAT(1X,11F12.7) 137 FORMAT(1X,F11.7) 715 format(22x,'Input-Output in F Format') 716 format(22x,'Input-Output in D Format') 43 format('No.',2x,'Curve',4x,'Input Param.',8x,'Correction',5x, $'Output Param.',4x,'Standard Deviation') 615 format(i2,i7,4f18.10) 616 format(i2,i7,4d18.10) 138 FORMAT(1X,'SUM OF ABSOLUTE VALUES OF CHECKS IS',1X,D12.6) 181 FORMAT(7X,'NORMAL EQUATIONS') 183 FORMAT (7X,'CORRELATION COEFFICIENTS') 184 FORMAT(7X,'NORMAL EQUATIONS TIMES INVERSE') 185 FORMAT(1X,'CHECK OF COMPUTED DEPENDENT VARIABLES FROM NORMAL EQUAT $IONS') 82 FORMAT(7X,'UNWEIGHTED OBSERVATIONAL EQUATIONS') 83 FORMAT(7X,'WEIGHTED OBSERVATIONAL EQUATIONS') 9 FORMAT(33X,'OBSERVATIONS') 955 FORMAT(3(9x,'phase V rad wt')) 10 FORMAT(3(9x,'phase light wt')) 755 FORMAT(3(9x,'JD V rad wt')) 756 FORMAT(3(9x,'JD light wt')) 149 FORMAT(' PROGRAM SHOULD NOT BE USED IN MODE 1 OR 3 WITH NON-ZERO E $CCENTRICITY') 40 FORMAT(' Sum(W*Res**2) for input values Sum(W*Res**2) pred $icted determinant') 11 format(3x,'Wave L',6x,'L1',8x,'L2 x1 x2 y1 y2 3rd $lt',5x,'opsf NOISE Sigma') 111 FORMAT(3x,'Wave L',6x,'L1',8x,'L2 x1 x2 y1 y2 ops $f',7x,'Sigma') 12 FORMAT('MODE IPB IFAT1 IFAT2 N1 N2 N1L N2L',7x,'Arg Per',5x $,'dperdt',7x,'TH e',8x,'V unit(km/s) V FAC') 206 FORMAT(F6.5,F11.4,2F11.4,F11.4,f10.3,2f8.3,i6,i9) 205 FORMAT(' ecc',3x,'S-M axis',6x,'F1',9x,'F2',8x,'V Gam', $7x,'INCL G1 G2',2x,'Nspot 1',2x,'Nspot 2') 402 FORMAT(' DEL EC DEL PER DEL F1 DEL F2 DEL PHS $ DEL INCL DEL G1 DEL G2 DEL T1 DEL T2') 403 FORMAT(' DEL ALB1 DEL ALB2 DEL POT1 DEL POT2 DEL Q $ DEL L1 DEL L2 DEL X1 DEL X2') 406 FORMAT(' ADJUSTMENT CONTROL INTEGERS; 1 SUPPRESSES ADJUSTMENT, 0 A $LLOWS ADJUSTMENT.') 702 FORMAT(F6.5,F10.4,2F10.4,F10.4,f9.3,2f7.3) 706 FORMAT(F7.4,f8.4,2F7.3,2d12.5,F10.5,4f7.3) 408 FORMAT(2F8.4,2F9.3,2d15.5,F12.5,4f9.3) 705 FORMAT(I1,1X,I1,1X,5I2) 54 FORMAT(' T1 T2',5x,'Alb 1',4x,'Alb 2',9x,'Pot 1',10x, $'Pot 2',7x,'M2/M1 x1(bolo) x2(bolo) y1(bolo) y2(bolo)') 707 FORMAT(' IFVC1 IFVC2 NLC KO KDISK ISYM nppl') 917 format('nref',3x,'mref',3x,'ifsmv1',3x,'ifsmv2',3x,'icor1',3x, $'icor2',3x,'ld') 708 FORMAT(8(4X,I3)) 912 format(i3,i7,i8,i9,i8,i8,i7) 650 FORMAT(20X,'RADII AND RELATED QUANTITIES (FROM INPUT)') 651 FORMAT(5X,'DOM1/DQ',5X,'DOM2/DQ',5X,'OM1-Q CORR.',5X,'OM2-Q CORR.' $,5X,'OM1 S.D.',4X,'OM2 S.D.',4X,'Q S.D.') 652 FORMAT(1X,3F12.6,4X,F12.6,4X,3F12.6) 653 FORMAT(' COMPONENT',11X,'R',9X,'DR/DOM',8X,'DR/DQ',11X,'S.D.') 654 FORMAT(2X,I1,1X,A6,4F14.6) 684 FORMAT(i2,4F13.5) 985 FORMAT(4f9.5) 983 FORMAT(1X,'STAR CO-LATITUDE LONGITUDE SPOT RADIUS TEMP.FACTOR $') 399 FORMAT(' DEL LAT DEL LONG DEL RAD DEL TEMPF DEL LAT $ del LONG del RAD del TEMPF') 60 FORMAT(4I3) 61 FORMAT(1X,4I6) 66 FORMAT(' STAR SPOT STAR SPOT') 166 FORMAT(' SPOTS TO BE ADJUSTED') 440 FORMAT(' AS1=FIRST ADJUSTED SPOT') 441 FORMAT(' AS2=SECOND ADJUSTED SPOT') 405 FORMAT(' ORDER OF PARAMETERS IS AS FOLLOWS:') 1440 FORMAT(' (1) - AS1 LATITUDE') 1441 FORMAT(' (2) - AS1 LONGITUDE') 1442 FORMAT(' (3) - AS1 ANGULAR RADIUS') 1443 FORMAT(' (4) - AS1 TEMPERATURE FACTOR') 1444 FORMAT(' (5) - AS2 LATITUDE') 1445 FORMAT(' (6) - AS2 LONGITUDE') 1446 FORMAT(' (7) - AS2 ANGULAR RADIUS') 1447 FORMAT(' (8) - AS2 TEMPERATURE FACTOR') 1448 FORMAT(' (9) - A=ORBITAL SEMI-MAJOR AXIS') 1449 FORMAT(' (10) - E=ORBITAL ECCENTRICITY') 1450 FORMAT(' (11) - PERR0=ARGUMENT of PERIASTRON at time HJD0') 1451 FORMAT(' (12) - F1=STAR 1 ROTATION PARAMETER') 1452 FORMAT(' (13) - F2=STAR 2 ROTATION PARAMETER') 1453 FORMAT(' (14) - PHASE SHIFT= PHASE OF PRIMARY CONJUNCTION') 1454 FORMAT(' (15) - VGAM=SYSTEMIC RADIAL VELOCITY') 1455 FORMAT(' (16) - INCL=ORBITAL INCLINATION') 1456 FORMAT(' (17) - g1=STAR 1 GRAVITY DARKENING EXPONENT') 1457 FORMAT(' (18) - g2=STAR 2 GRAVITY DARKENING EXPONENT') 1458 FORMAT(' (19) - T1=STAR 1 AVERAGE SURFACE TEMPERATURE') 1459 FORMAT(' (20) - T2=STAR 2 AVERAGE SURFACE TEMPERATURE') 1460 FORMAT(' (21) - ALB1=STAR 1 BOLOMETRIC ALBEDO') 1461 FORMAT(' (22) - ALB2=STAR 2 BOLOMETRIC ALBEDO') 1462 FORMAT(' (23) - POT1=STAR 1 SURFACE POTENTIAL') 1463 FORMAT(' (24) - POT2=STAR 2 SURFACE POTENTIAL') 1464 FORMAT(' (25) - Q=MASS RATIO (STAR 2/STAR 1)') 1470 FORMAT(' (26) - HJD0= Hel. JD reference time') 1471 FORMAT(' (27) - PERIOD= orbital period') 1472 FORMAT(' (28) - DPDT= time derivative of orbital period') 1473 FORMAT(' (29) - DPERDT= time derivative of argument of periastron' $) 1474 FORMAT(' (30) - unused channel reserved for future expansion') 1465 FORMAT(' (31) - L1=STAR 1 RELATIVE MONOCHROMATIC LUMINOSITY') 1466 FORMAT(' (32) - L2=STAR 2 RELATIVE MONOCHROMATIC LUMINOSITY') 1467 FORMAT(' (33) - X1=STAR 1 LIMB DARKENING COEFFICIENT') 1468 FORMAT(' (34) - X2=STAR 2 LIMB DARKENING COEFFICIENT') 1469 FORMAT(' (35) - el3=third light') 119 format(1x,i6,i13,f18.8) 159 format(' Sums of squares of residuals for separate curves, includi $ng only individual weights') 169 format(' Curve No. of obs. Sum of squares') 1063 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') 170 format(i3,f17.6,d18.10,d14.6,f10.4) 649 format(i1,f15.6,d15.10,d13.6,f10.4) 171 format('JDPHS',5x,'J.D. zero',7x,'Period',11x,'dPdt', $6x,'Ph. shift') 911 format(7(i1,1x)) 840 format('Do not try to adjust the ephemeris or any time derivative $parameters','when JDPHS = 2') 839 format('Ordinarily one should not try to adjust both PSHIFT and', $' HJD0. They are perfectly correlated if the period is constant', $' and extremely highly correlated if the period is not constant') tt=2.d0/3.d0 nbins=1 nl=1 inmax=1 inmin=1 nf1=1 nf2=1 toldis=1.d-5 pi=dacos(-1.d0) en0=6.0254d23 XTHA(1)=0.d0 XTHA(2)=.5d0*PI XTHA(3)=.5d0*PI XTHA(4)=.5d0*PI XFIA(1)=0.d0 XFIA(2)=0.d0 XFIA(3)=.5d0*PI XFIA(4)=PI IBEF=0 KH=25 NS=1 NI=0 NY=0 KNOBS(1)=0 WRITE(6,102) WRITE(6,405) WRITE(6,101) WRITE(6,440) WRITE(6,441) WRITE(6,101) WRITE(6,1440) WRITE(6,1441) WRITE(6,1442) WRITE(6,1443) WRITE(6,1444) WRITE(6,1445) WRITE(6,1446) WRITE(6,1447) WRITE(6,1448) WRITE(6,1449) WRITE(6,1450) WRITE(6,1451) WRITE(6,1452) WRITE(6,1453) WRITE(6,1454) WRITE(6,1455) WRITE(6,1456) WRITE(6,1457) WRITE(6,1458) WRITE(6,1459) WRITE(6,1460) WRITE(6,1461) WRITE(6,1462) WRITE(6,1463) WRITE(6,1464) WRITE(6,1470) WRITE(6,1471) WRITE(6,1472) WRITE(6,1473) WRITE(6,1474) WRITE(6,1465) WRITE(6,1466) WRITE(6,1467) WRITE(6,1468) WRITE(6,1469) READ(5,56)(DEL(I),I=1,8) READ(5,56)(DEL(I),I=10,14),(DEL(I),I=16,20) READ(5,56)(DEL(I),I=21,25),(del(i),i=31,34) READ(5,20)(KEP(I),I=1,35),IFDER,IFM,IFR,xlamda READ(5,60) KSPA,NSPA,KSPB,NSPB READ(5,705)IFVC1,IFVC2,NLC,KO,KDISK,ISYM,nppl read(5,911)nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld read(5,649) jdphs,hjd0,period,dpdt,pshift if(jdphs.eq.2.and.kep(26).eq.0) write(6,840) if(jdphs.eq.2.and.kep(27).eq.0) write(6,840) if(jdphs.eq.2.and.kep(28).eq.0) write(6,840) if(kep(14).eq.0.and.kep(26).eq.0) write(6,839) READ(5,701)MODE,IPB,IFAT1,IFAT2,N1,N2,N1L,N2L,perr0,dperdt,THE, $vunit READ(5,702)E,A,F1,F2,VGA,XINCL,GR1,GR2 READ(5,706) TAVH,TAVC,ALB1,ALB2,PHSV,PCSV,RM,xbol1,xbol2,ybol1, $ybol2 acm=6.960d10*a CALL SINCOS(1,N1,N1,SNTHH,CSTHH,SNFIH,CSFIH,MMSAVH) CALL SINCOS(2,N2,N1,SNTHH,CSTHH,SNFIH,CSFIH,MMSAVH) CALL SINCOS(1,N1L,N1L,SNTHL,CSTHL,SNFIL,CSFIL,MMSAVL) CALL SINCOS(2,N2L,N1L,SNTHL,CSTHL,SNFIL,CSFIL,MMSAVL) dint1=pi*(1.d0-xbol1/3.d0) if(ld.eq.2) dint1=dint1+pi*ybol1*2.d0/9.d0 if(ld.eq.3) dint1=dint1-pi*ybol1*.2d0 dint2=pi*(1.d0-xbol2/3.d0) if(ld.eq.2) dint2=dint2+pi*ybol2*2.d0/9.d0 if(ld.eq.3) dint2=dint2-pi*ybol2*.2d0 IS=ISYM+1 KEEP(36)=0 MM1=N1+1 MM2=N1+N2+2 MM3=N1L+1 MM4=N1L+N2L+2 M1H=MMSAVH(MM1) M2H=MMSAVH(MM2) M1L=MMSAVL(MM3) M2L=MMSAVL(MM4) MTLH=M1H+M2H MTLL=M1L+M2L NVC=IFVC1+IFVC2 NLVC=NLC+NVC NVCP=NVC+1 IF(NVC.NE.0) GOTO 288 KEP(9)=1 KEP(15)=1 288 CONTINUE DO 84 I=1,35 KEEP(I)=KEP(I) 84 LOW(I)=1 LOW(1)=0 LOW(2)=0 LOW(3)=0 LOW(5)=0 LOW(6)=0 LOW(7)=0 LOW(10)=0 LOW(11)=0 LOW(12)=0 LOW(13)=0 LOW(14)=0 LOW(16)=0 LOW(23)=0 LOW(24)=0 LOW(25)=0 ifap=1-keep(29) ifphi=1-keep(26)*keep(27)*keep(28) KOSQ=(KO-2)*(KO-2) IF(NVC.EQ.0) GOTO 195 DO 90 I=1,NVC 90 READ(5,218) WLA(I),HLA(I),CLA(I),X1A(I),X2A(I),y1a(i),y2a(i), $opsfa(i),sigma(i) 195 CONTINUE IF(NLVC.EQ.NVC) GOTO 194 DO 190 I=NVCP,NLVC 190 read(5,18)wla(i),hla(i),cla(i),x1a(i),x2a(i),y1a(i),y2a(i), $el3a(i),opsfa(i),noise(i),sigma(i) 194 CONTINUE NSP1=0 NSP2=0 DO 988 KP=1,2 DO 987 I=1,100 READ(5,985)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 988 SNLAT(KP,I)=dsin(XLAT(KP,I)) CSLAT(KP,I)=dcos(XLAT(KP,I)) SNLNG(KP,I)=dsin(XLONG(KP,I)) CSLNG(KP,I)=dcos(XLONG(KP,I)) RDSP(KP,I)=RADSP(KP,I) TMSP(KP,I)=TEMSP(KP,I) IF(KP.EQ.1) NSP1=NSP1+1 987 IF(KP.EQ.2) NSP2=NSP2+1 988 CONTINUE NSTOT=NSP1+NSP2 ncl=0 do 1062 i=1,100 read(5,1063) 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 1066 ncl=ncl+1 dens(i)=edens(i)*xmue(i)/en0 1062 continue 1066 continue para(1)=xlat(kspa,nspa) para(2)=xlong(kspa,nspa) para(3)=radsp(kspa,nspa) para(4)=temsp(kspa,nspa) para(5)=xlat(kspb,nspb) para(6)=xlong(kspb,nspb) para(7)=radsp(kspb,nspb) para(8)=temsp(kspb,nspb) para(9)=a para(10)=e para(11)=perr0 para(12)=f1 para(13)=f2 para(14)=pshift para(15)=vga para(16)=xincl para(17)=gr1 para(18)=gr2 para(19)=tavh para(20)=tavc para(21)=alb1 para(22)=alb2 para(23)=phsv para(24)=pcsv para(25)=rm para(26)=hjd0 para(27)=period para(28)=dpdt para(29)=dperdt para(30)=0.d0 ib=nvc do 191 irx=31,30+nlc ib=ib+1 191 para(irx)=hla(ib) ib=nvc do 186 irx=31+nlc,30+2*nlc ib=ib+1 186 para(irx)=cla(ib) ib=nvc do 187 irx=31+2*nlc,30+3*nlc ib=ib+1 187 para(irx)=x1a(ib) ib=nvc do 188 irx=31+3*nlc,30+4*nlc ib=ib+1 188 para(irx)=x2a(ib) ib=nvc do 189 irx=31+4*nlc,30+5*nlc ib=ib+1 189 para(irx)=el3a(ib) PERT=perr0 EC=E hjd=hjd0 PSH=PSHIFT IRTE=0 IRVOL1=0 IRVOL2=0 CALL MODLOG(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD, $RM,PHSV,PCSV,GR1,GR2,ALB1,ALB2,N1,N2,F1,F2,MOD,XINCL,THE,MODE, $SNTHH,CSTHH,SNFIH,CSFIH,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,GLUMP1, $GLUMP2,CSBT1,CSBT2,GMAG1,GMAG2) call ellone(f1,dp,rm,xld,omcr(1),xldd,omd) rmr=1.d0/rm call ellone(f2,dp,rmr,xld,omcr(2),xldd,omd) omcr(2)=rm*omcr(2)+(1.d0-rm)*.5d0 po(1)=phsv po(2)=pcsv IF(E.EQ.0.d0) GOTO 134 CALL VOLUME(VOL1,RM,PHSV,DP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ, $GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD,SNTHH,CSTHH,SNFIH,CSFIH,SUMMD $,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2 $,GMAG1,GMAG2,GR1,1) CALL VOLUME(VOL2,RM,PCSV,DP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ, $GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD,SNTHH,CSTHH,SNFIH,CSFIH,SUMMD $,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2 $,GMAG1,GMAG2,GR2,1) DAP=1.d0+E P1AP=PHSV-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,MMSAVH,FR1,FR2,HLD,SNTHH,CSTHH,SNFIH,CSFIH,SUMMD $,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2 $,GMAG1,GMAG2,GR1,2) DPDX1=(PHSV-P1AP)*(1.d0-E*E)*.5d0/E P2AP=PCSV-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,MMSAVH,FR1,FR2,HLD,SNTHH,CSTHH,SNFIH,CSFIH,SUMMD $,SMD,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2 $,GMAG1,GMAG2,GR2,2) DPDX2=(PCSV-P2AP)*(1.d0-E*E)*.5d0/E 134 CONTINUE PHP=PHPER POTH=PHSV POTC=PCSV POT1=PHSV POT2=PCSV DO 24 I=1,NLVC opsf=opsfa(i) 24 CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD, $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POTH,POTC,N1,N2,F1,F2,d,hla(i), $cla(i),x1a(i),x2a(i),y1a(i),y2a(i),gr1,gr2,wla(i),sm1,sm2,tpolh, $tpolc,sbrh,sbrc,ifat1,ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2, $ybol1,ybol2,php,rm,xincl,hot,cool,snthh,csthh,snfih,csfih,tldh, $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) DEL(9)=0.d0 DEL(15)=0.d0 del(26)=0.d0 del(27)=0.d0 del(28)=0.d0 del(29)=0.d0 DEL(35)=0.d0 WRITE(6,101) WRITE(6,399) WRITE(6,55)(DEL(I),I=1,8) WRITE(6,101) WRITE(6,402) WRITE(6,55)(DEL(I),I=10,14),(DEL(I),I=16,20) WRITE(6,101) WRITE(6,403) WRITE(6,55)(DEL(I),I=21,25),(del(i),i=31,34) WRITE(6,101) WRITE(6,406) WRITE(6,20)(KEP(I),I=1,35),IFDER,IFM,IFR,xlamda WRITE(6,101) WRITE(6,166) WRITE(6,66) WRITE(6,61) KSPA,NSPA,KSPB,NSPB WRITE(6,101) WRITE(6,707) WRITE(6,708)IFVC1,IFVC2,NLC,KO,KDISK,ISYM,nppl WRITE(6,101) write(6,917) write(6,912)nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld WRITE(6,101) write(6,171) write(6,170) jdphs,hjd0,period,dpdt,pshift WRITE(6,101) WRITE(6,12) WRITE(6,1)MODE,IPB,IFAT1,IFAT2,N1,N2,N1L,N2L,perr0,dperdt,THE, $vunit,vfac WRITE(6,101) WRITE(6,205) WRITE(6,206) E,A,F1,F2,VGA,XINCL,GR1,GR2,nsp1,nsp2 WRITE(6,101) WRITE(6,54) WRITE(6,408) TAVH,TAVC,ALB1,ALB2,PHSV,PCSV,RM,xbol1,xbol2,ybol1, $ybol2 IF(NVC.EQ.0) GOTO 196 WRITE(6,101) WRITE(6,111) DO 91 I=1,NVC 91 WRITE(6,218)WLA(I),HLA(I),CLA(I),X1A(I),X2A(I),y1a(i),y2a(i), $opsfa(i),sigma(i) 196 CONTINUE IF(NLVC.EQ.NVC) GOTO 197 WRITE(6,101) WRITE(6,11) DO 92 I=NVCP,NLVC 92 write(6,85)wla(i),hla(i),cla(i),x1a(i),x2a(i),y1a(i),y2a(i), $el3a(i),opsfa(i),noise(i),sigma(i) 197 CONTINUE WRITE(6,101) IF(NSTOT.GT.0) WRITE(6,983) DO 688 KP=1,2 IF((NSP1+KP-1).EQ.0) GOTO 688 IF((NSP2+(KP-2)**2).EQ.0) GOTO 688 NSPOT=NSP1 IF(KP.EQ.2) NSPOT=NSP2 DO 687 I=1,NSPOT 687 WRITE(6,684)KP,XLAT(KP,I),XLONG(KP,I),RADSP(KP,I),TEMSP(KP,I) 688 WRITE(6,101) if(ncl.eq.0) goto 1067 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,101) 1067 continue WRITE(6,101) WRITE(6,9) DO 75 LCV=1,NLVC WRITE(6,101) IF(LCV.LE.NVC.and.jdphs.eq.2) WRITE(6,955) IF(LCV.GT.NVC.and.jdphs.eq.2) WRITE(6,10) IF(LCV.LE.NVC.and.jdphs.eq.1) WRITE(6,755) IF(LCV.GT.NVC.and.jdphs.eq.1) WRITE(6,756) DO 74 I=NS,7000 ifirst=nppl*(i-1)+NY+1 last=ifirst+nppl-1 READ(5,2) (phjd(in),flux(in),wt(in),in=ifirst,last) WRITE(6,2) (phjd(in),flux(in),wt(in),in=ifirst,last) IF(phjd(ifirst).gt.-10000.d0) GOTO 74 NI=-(phjd(ifirst)+10000.d0) NY=NY+NI NOBS=nppl*(I-NS-1)+NI GOTO 150 74 CONTINUE 150 NS=I-1 LC1=LCV+1 75 KNOBS(LC1)=NOBS+KNOBS(LCV) do 275 ijp=1,knobs(lc1) phas(ijp)=phjd(ijp) if(jdphs.eq.1) call jdph(phjd(ijp),0.d0,hjd0,period,dpdt,xjddum, $phas(ijp)) 275 continue matrix=31 do 427 ima=1,30 427 matrix=matrix-keep(ima) matrix=matrix+nlc*(5-keep(31)-keep(32)-keep(33)-keep(34)-keep(35)) MAT=MATRIX-1 EM=dfloat(MATRIX-15) KTR=.24d0*EM+2.2d0 IF(EM.LE.1.5d0) KTR=1 IF(EM.GT.12.d0) KTR=5 NCOEFF=MATRIX*KNOBS(LC1) NMAT=MAT*KNOBS(LC1) NCOF=NCOEFF DO 63 J=1,NCOF 63 HOLD(J)=0.d0 IF(KOSQ.EQ.1) GOTO 71 DO 416 IOBS=1,NCOEFF 416 OBS(IOBS)=0.d0 KSMAX=37 KSSMAX=37 IF(E.EQ.0.d0) KSMAX=1 IF(E.NE.0.d0) KSSMAX=1 DO 419 KSS=1,KSSMAX DO 417 IB=1,NLVC VC1=0.d0 VC2=0.d0 ELIT=0.d0 IF(IB.GT.NVC) ELIT=1.d0 IF(IB.EQ.IFVC1) VC1=1.d0 IF(IB.EQ.(IFVC2*(1+IFVC1))) VC2=1.d0 IST=KNOBS(IB)+1 IB1=IB+1 ISP=KNOBS(IB1) DO 418 IX=IST,ISP DO 420 KS=1,KSMAX hjd=phjd(ix) dtime=hjd-hjd0 IRTE=0 IRVOL1=0 IRVOL2=0 IF(E.NE.0.d0) GOTO 297 IF(IX.NE.IST) IRTE=1 IF(IX.EQ.1) GOTO 297 IRVOL1=1 IRVOL2=1 297 CONTINUE KSR=KS IF(E.EQ.0.d0) KSR=KSS if(e.ne.0.d0) goto 1110 IF(ISYM.NE.1) GOTO 1110 IF(KSR.EQ.1) GOTO 420 1110 CONTINUE KH=KSR-2 IF(KH.GT.0) GOTO 740 KH=1 GOTO 941 740 CONTINUE if(ifap*kh.eq.11) goto 842 if(ifphi*kh.eq.14) goto 842 IF(KEEP(KH).EQ.1) GOTO 420 842 IF(E.EQ.0.d0) GOTO 889 IF(KSR.LE.2) GOTO 889 IF(KH.LE.9) IRTE=1 IF(KH.LE.9) IRVOL1=1 IF(KH.LE.9) IRVOL2=1 IF(KH.EQ.12) IRVOL2=1 IF(KH.EQ.13) IRVOL1=1 IF(KH.EQ.15) IRTE=1 IF(KH.EQ.15) IRVOL1=1 IF(KH.EQ.15) IRVOL2=1 IF(KH.EQ.16) IRTE=1 IF(KH.EQ.16) IRVOL1=1 IF(KH.EQ.16) IRVOL2=1 IF(KH.EQ.17) IRVOL2=1 IF(KH.EQ.18) IRVOL1=1 IF(KH.EQ.19) IRVOL1=1 IF(KH.EQ.19) IRVOL2=1 IF(KH.EQ.20) IRVOL1=1 IF(KH.EQ.20) IRVOL2=1 IF(KH.EQ.21) IRVOL1=1 IF(KH.EQ.21) IRVOL2=1 IF(KH.EQ.22) IRVOL1=1 IF(KH.EQ.22) IRVOL2=1 IF(KH.EQ.23) IRVOL2=1 IF(KH.EQ.24) IRVOL1=1 IF(KH.GE.31) IRVOL1=1 IF(KH.GE.31) IRVOL2=1 889 CONTINUE LCF=0 IF(KH.GT.30) LCF=IB-NVC KPCT1=0 KPCT2=0 KSP=KH IF(KH.GT.30) KSP=30 IF(KH.LT.2) GOTO 808 DO 804 ICT=1,KSP 804 KPCT1=KPCT1+1-KEEP(ICT) GOTO 809 808 KPCT1=1 809 CONTINUE IF(KH.LT.31) GOTO 806 DO 805 ICT=31,KH 805 KPCT2=KPCT2+1-KEEP(ICT) GOTO 807 806 KPCT2=1 807 CONTINUE II=(KPCT1+NLC*(KPCT2-1)+LCF-1)*KNOBS(LC1)+IX IF(KH.EQ.9) GOTO 300 IF(KH.EQ.15) GOTO 308 if(kh.eq.26) goto 844 if(kh.eq.27) goto 845 if(kh.eq.28) goto 846 if(kh.eq.29) goto 847 IF(KH.EQ.35) GOTO 301 IF(KH.NE.31) GOTO 941 IF(MODE.LE.0) GOTO 941 IF(IPB.EQ.1) GOTO 941 IF(IB.GT.NVC) OBS(II)=(BR(IX)-EL3A(IB))/HLA(IB) GOTO 420 941 CONTINUE DL=DEL(KH) IF(ISYM.EQ.1) DL=.5d0*DEL(KH) SIGN=1.d0 ISS=1 DO 421 IH=1,35 421 DLIF(IH)=0.d0 IF(KSR.LE.2) GOTO 777 ISS=IS DLIF(KH)=1.d0 777 CONTINUE DO 319 IL=1,ISS IF(E.NE.0.d0) GOTO 4011 IF(ISYM.EQ.1) GOTO 4012 GOTO 940 4011 IF(KSR.LE.2) GOTO 940 goto 4014 4012 IF(IL.EQ.2) GOTO 4016 4014 IF(LOW(KH).EQ.1) GOTO 314 VOL1=SVOL1 VOL2=SVOL2 SUMM1=SSUM1 SUMM2=SSUM2 SM1=SSM1 SM2=SSM2 DO 851 IRE=1,MTLH 851 TLDH(IRE)=STLDH(IRE) DO 508 IRE=1,M1H RV(IRE)=SRV(IRE) GRX(IRE)=SGRX(IRE) GRY(IRE)=SGRY(IRE) GRZ(IRE)=SGRZ(IRE) GLUMP1(IRE)=SGLM1(IRE) GRV1(IRE)=SGRV1(IRE) XX1(IRE)=SXX1(IRE) YY1(IRE)=SYY1(IRE) ZZ1(IRE)=SZZ1(IRE) GMAG1(IRE)=SGMG1(IRE) CSBT1(IRE)=SCSB1(IRE) RF1(IRE)=SRF1(IRE) FR1(IRE)=SFR1(IRE) 508 SLUMP1(IRE)=SLMP1(IRE) DO 309 IRE=1,M2H RVQ(IRE)=SRVQ(IRE) GRXQ(IRE)=SGRXQ(IRE) GRYQ(IRE)=SGRYQ(IRE) GRZQ(IRE)=SGRZQ(IRE) GLUMP2(IRE)=SGLM2(IRE) GRV2(IRE)=SGRV2(IRE) XX2(IRE)=SXX2(IRE) YY2(IRE)=SYY2(IRE) ZZ2(IRE)=SZZ2(IRE) GMAG2(IRE)=SGMG2(IRE) CSBT2(IRE)=SCSB2(IRE) RF2(IRE)=SRF2(IRE) FR2(IRE)=SFR2(IRE) 309 SLUMP2(IRE)=SLMP2(IRE) GOTO 940 4016 IF(LOW(KH).EQ.1) GOTO 4018 VOL1=EVOL1 VOL2=EVOL2 SUMM1=ESUM1 SUMM2=ESUM2 SM1=ESM1 SM2=ESM2 DO 852 IRE=1,MTLH 852 TLDH(IRE)=ETLDH(IRE) DO 1508 IRE=1,M1H RV(IRE)=ERV(IRE) GRX(IRE)=EGRX(IRE) GRY(IRE)=EGRY(IRE) GRZ(IRE)=EGRZ(IRE) GLUMP1(IRE)=EGLM1(IRE) GRV1(IRE)=EGRV1(IRE) XX1(IRE)=EXX1(IRE) YY1(IRE)=EYY1(IRE) ZZ1(IRE)=EZZ1(IRE) GMAG1(IRE)=EGMG1(IRE) CSBT1(IRE)=ECSB1(IRE) RF1(IRE)=ERF1(IRE) FR1(IRE)=EFR1(IRE) 1508 SLUMP1(IRE)=ELMP1(IRE) DO 1309 IRE=1,M2H RVQ(IRE)=ERVQ(IRE) GRXQ(IRE)=EGRXQ(IRE) GRYQ(IRE)=EGRYQ(IRE) GRZQ(IRE)=EGRZQ(IRE) GLUMP2(IRE)=EGLM2(IRE) GRV2(IRE)=EGRV2(IRE) XX2(IRE)=EXX2(IRE) YY2(IRE)=EYY2(IRE) ZZ2(IRE)=EZZ2(IRE) GMAG2(IRE)=EGMG2(IRE) CSBT2(IRE)=ECSB2(IRE) RF2(IRE)=ERF2(IRE) FR2(IRE)=EFR2(IRE) 1309 SLUMP2(IRE)=ELMP2(IRE) GOTO 940 4018 CONTINUE VOL1=EVOL1L VOL2=EVOL2L SUMM1=ESUM1L SUMM2=ESUM2L SM1=ESM1L SM2=ESM2L DO 853 IRE=1,MTLL 853 TLDL(IRE)=ETLDL(IRE) DO 310 IRE=1,M1L RV(IRE)=ERVL(IRE) GRX(IRE)=EGRXL(IRE) GRY(IRE)=EGRYL(IRE) GRZ(IRE)=EGRZL(IRE) GLUMP1(IRE)=EGLM1L(IRE) GRV1(IRE)=EGRV1L(IRE) XX1(IRE)=EXX1L(IRE) YY1(IRE)=EYY1L(IRE) ZZ1(IRE)=EZZ1L(IRE) GMAG1(IRE)=EGMG1L(IRE) CSBT1(IRE)=ECSB1L(IRE) RF1(IRE)=ERF1L(IRE) FR1(IRE)=EFR1L(IRE) 310 SLUMP1(IRE)=ELMP1L(IRE) DO 311 IRE=1,M2L RVQ(IRE)=ERVQL(IRE) GRXQ(IRE)=EGRXQL(IRE) GRYQ(IRE)=EGRYQL(IRE) GRZQ(IRE)=EGRZQL(IRE) GLUMP2(IRE)=EGLM2L(IRE) GRV2(IRE)=EGRV2L(IRE) XX2(IRE)=EXX2L(IRE) YY2(IRE)=EYY2L(IRE) ZZ2(IRE)=EZZ2L(IRE) GMAG2(IRE)=EGMG2L(IRE) CSBT2(IRE)=ECSB2L(IRE) RF2(IRE)=ERF2L(IRE) FR2(IRE)=EFR2L(IRE) 311 SLUMP2(IRE)=ELMP2L(IRE) GOTO 940 314 CONTINUE VOL1=SVOL1L VOL2=SVOL2L SUMM1=SSUM1L SUMM2=SSUM2L SM1=SSM1L SM2=SSM2L DO 854 IRE=1,MTLL 854 TLDL(IRE)=STLDL(IRE) DO 1310 IRE=1,M1L RV(IRE)=SRVL(IRE) GRX(IRE)=SGRXL(IRE) GRY(IRE)=SGRYL(IRE) GRZ(IRE)=SGRZL(IRE) GLUMP1(IRE)=SGLM1L(IRE) GRV1(IRE)=SGRV1L(IRE) XX1(IRE)=SXX1L(IRE) YY1(IRE)=SYY1L(IRE) ZZ1(IRE)=SZZ1L(IRE) GMAG1(IRE)=SGMG1L(IRE) CSBT1(IRE)=SCSB1L(IRE) RF1(IRE)=SRF1L(IRE) FR1(IRE)=SFR1L(IRE) 1310 SLUMP1(IRE)=SLMP1L(IRE) DO 1311 IRE=1,M2L RVQ(IRE)=SRVQL(IRE) GRXQ(IRE)=SGRXQL(IRE) GRYQ(IRE)=SGRYQL(IRE) GRZQ(IRE)=SGRZQL(IRE) GLUMP2(IRE)=SGLM2L(IRE) GRV2(IRE)=SGRV2L(IRE) XX2(IRE)=SXX2L(IRE) YY2(IRE)=SYY2L(IRE) ZZ2(IRE)=SZZ2L(IRE) GMAG2(IRE)=SGMG2L(IRE) CSBT2(IRE)=SCSB2L(IRE) RF2(IRE)=SRF2L(IRE) FR2(IRE)=SFR2L(IRE) 1311 SLUMP2(IRE)=SLMP2L(IRE) 940 CONTINUE DELS=DL*SIGN SIGN=-1.d0 IF(NSPA.EQ.0) GOTO 470 xlt=xlat(kspa,nspa)+dels*dlif(1) xlng(kspa,nspa)=xlong(kspa,nspa)+dels*dlif(2) snlat(kspa,nspa)=dsin(xlt) cslat(kspa,nspa)=dcos(xlt) snlng(kspa,nspa)=dsin(xlng(kspa,nspa)) cslng(kspa,nspa)=dcos(xlng(kspa,nspa)) RDSP(KSPA,NSPA)=RADSP(KSPA,NSPA)+DELS*DLIF(3) TMSP(KSPA,NSPA)=TEMSP(KSPA,NSPA)+DELS*DLIF(4) 470 CONTINUE IF(NSPB.EQ.0) GOTO 471 xlt=xlat(kspb,nspb)+dels*dlif(5) xlng(kspb,nspb)=xlong(kspb,nspb)+dels*dlif(6) snlat(kspb,nspb)=dsin(xlt) cslat(kspb,nspb)=dcos(xlt) snlng(kspb,nspb)=dsin(xlng(kspb,nspb)) cslng(kspb,nspb)=dcos(xlng(kspb,nspb)) RDSP(KSPB,NSPB)=RADSP(KSPB,NSPB)+DELS*DLIF(7) TMSP(KSPB,NSPB)=TEMSP(KSPB,NSPB)+DELS*DLIF(8) 471 CONTINUE EC=E+DELS*DLIF(10) PERT=perr0+DELS*DLIF(11) FF1=F1+DELS*DLIF(12) FF2=F2+DELS*DLIF(13) PSH=PSHIFT+DELS*DLIF(14) XINC=XINCL+DELS*DLIF(16) G1=GR1+DELS*DLIF(17) G2=GR2+DELS*DLIF(18) T1=TAVH+DELS*DLIF(19) T2=TAVC+DELS*DLIF(20) A1=ALB1+DELS*DLIF(21) A2=ALB2+DELS*DLIF(22) POT1=PHSV+DELS*DLIF(23) POT2=PCSV+DELS*DLIF(24) RMASS=RM+DELS*DLIF(25) HL=HLA(IB)+DELS*DLIF(31) CL=CLA(IB)+DELS*DLIF(32) X1=X1A(IB)+DELS*DLIF(33) X2=X2A(IB)+DELS*DLIF(34) y1=y1a(ib) y2=y2a(ib) opsf=opsfa(ib) IF(KSR.EQ.1) GOTO 802 IF(KSR.EQ.2) GOTO 872 IF(LOW(KH).EQ.1) GOTO 802 872 CALL MODLOG(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD, $RMASS,POT1,POT2,G1,G2,A1,A2,N1,N2,FF1,FF2,MOD,XINC,THE,MODE, $SNTHH,CSTHH,SNFIH,CSFIH,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,GLUMP1, $GLUMP2,CSBT1,CSBT2,GMAG1,GMAG2) CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVH,FR1,FR2,HLD, $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POT1,POT2,N1,N2,FF1,FF2,d,hl,cl,x1, $x2,y1,y2,g1,g2,wla(ib),sm1,sm2,tph,tpc,sbrh,sbrc,ifat1,ifat2,t1, $t2,a1,a2,xbol1,xbol2,ybol1,ybol2,phas(ix),rmass,xinc,hot,cool, $snthh,csthh,snfih,csfih,tldh,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) GOTO 801 802 CONTINUE CALL MODLOG(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVL,FR1,FR2,HLD, $RMASS,POT1,POT2,G1,G2,A1,A2,N1L,N2L,FF1,FF2,MOD,XINC,THE,MODE, $SNTHL,CSTHL,SNFIL,CSFIL,GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,GLUMP1, $GLUMP2,CSBT1,CSBT2,GMAG1,GMAG2) CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVL,FR1,FR2,HLD, $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POT1,POT2,N1L,N2L,FF1,FF2,d,hl,cl, $x1,x2,y1,y2,g1,g2,wla(ib),sm1,sm2,tph,tpc,sbrh,sbrc,ifat1,ifat2, $t1,t2,a1,a2,xbol1,xbol2,ybol1,ybol2,phas(ix),rmass,xinc,hot,cool, $snthl,csthl,snfil,csfil,tldl,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) 801 CONTINUE IF(E.NE.0.d0) GOTO 4111 IF(ISYM.EQ.0) GOTO 602 IF(IX.NE.IST) GOTO 602 IF(KSR.EQ.1) GOTO 602 IF(KSR.EQ.2) GOTO 4119 IF(LOW(KH).EQ.1) GOTO 4112 GOTO 4119 4111 IF(IRTE.EQ.1) GOTO 602 IF(KSR.GE.3) GOTO 602 IF(KSR.EQ.2) GOTO 4119 4112 IF(IL.EQ.2) GOTO 4116 SVOL1L=VOL1 SVOL2L=VOL2 SSUM1L=SUMM1 SSUM2L=SUMM2 SSM1L=SM1 SSM2L=SM2 DO 855 IHL=1,MTLL 855 STLDL(IHL)=TLDL(IHL) DO 603 IHL=1,M1L SRVL(IHL)=RV(IHL) SGRXL(IHL)=GRX(IHL) SGRYL(IHL)=GRY(IHL) SGRZL(IHL)=GRZ(IHL) SLMP1L(IHL)=SLUMP1(IHL) SGLM1L(IHL)=GLUMP1(IHL) SGRV1L(IHL)=GRV1(IHL) SXX1L(IHL)=XX1(IHL) SYY1L(IHL)=YY1(IHL) SZZ1L(IHL)=ZZ1(IHL) SGMG1L(IHL)=GMAG1(IHL) SCSB1L(IHL)=CSBT1(IHL) SRF1L(IHL)=RF1(IHL) SFR1L(IHL)=FR1(IHL) 603 CONTINUE DO 606 IHL=1,M2L SRVQL(IHL)=RVQ(IHL) SGRXQL(IHL)=GRXQ(IHL) SGRYQL(IHL)=GRYQ(IHL) SGRZQL(IHL)=GRZQ(IHL) SLMP2L(IHL)=SLUMP2(IHL) SGLM2L(IHL)=GLUMP2(IHL) SGRV2L(IHL)=GRV2(IHL) SXX2L(IHL)=XX2(IHL) SYY2L(IHL)=YY2(IHL) SZZ2L(IHL)=ZZ2(IHL) SGMG2L(IHL)=GMAG2(IHL) SCSB2L(IHL)=CSBT2(IHL) SRF2L(IHL)=RF2(IHL) SFR2L(IHL)=FR2(IHL) 606 CONTINUE GOTO 602 4116 CONTINUE EVOL1L=VOL1 EVOL2L=VOL2 ESUM1L=SUMM1 ESUM2L=SUMM2 ESM1L=SM1 ESM2L=SM2 DO 856 IHL=1,MTLL 856 ETLDL(IHL)=TLDL(IHL) DO 1603 IHL=1,M1L ERVL(IHL)=RV(IHL) EGRXL(IHL)=GRX(IHL) EGRYL(IHL)=GRY(IHL) EGRZL(IHL)=GRZ(IHL) ELMP1L(IHL)=SLUMP1(IHL) EGLM1L(IHL)=GLUMP1(IHL) EGRV1L(IHL)=GRV1(IHL) EXX1L(IHL)=XX1(IHL) EYY1L(IHL)=YY1(IHL) EZZ1L(IHL)=ZZ1(IHL) EGMG1L(IHL)=GMAG1(IHL) ECSB1L(IHL)=CSBT1(IHL) ERF1L(IHL)=RF1(IHL) EFR1L(IHL)=FR1(IHL) 1603 CONTINUE DO 1606 IHL=1,M2L ERVQL(IHL)=RVQ(IHL) EGRXQL(IHL)=GRXQ(IHL) EGRYQL(IHL)=GRYQ(IHL) EGRZQL(IHL)=GRZQ(IHL) ELMP2L(IHL)=SLUMP2(IHL) EGLM2L(IHL)=GLUMP2(IHL) EGRV2L(IHL)=GRV2(IHL) EXX2L(IHL)=XX2(IHL) EYY2L(IHL)=YY2(IHL) EZZ2L(IHL)=ZZ2(IHL) EGMG2L(IHL)=GMAG2(IHL) ECSB2L(IHL)=CSBT2(IHL) ERF2L(IHL)=RF2(IHL) EFR2L(IHL)=FR2(IHL) 1606 CONTINUE GOTO 602 4119 IF(IL.EQ.2) GOTO 4120 SVOL1=VOL1 SVOL2=VOL2 SSUM1=SUMM1 SSUM2=SUMM2 SSM1=SM1 SSM2=SM2 DO 857 IHH=1,MTLH 857 STLDH(IHH)=TLDH(IHH) DO 601 IHH=1,M1H SRV(IHH)=RV(IHH) SGRX(IHH)=GRX(IHH) SGRY(IHH)=GRY(IHH) SGRZ(IHH)=GRZ(IHH) SLMP1(IHH)=SLUMP1(IHH) SGLM1(IHH)=GLUMP1(IHH) SGRV1(IHH)=GRV1(IHH) SXX1(IHH)=XX1(IHH) SYY1(IHH)=YY1(IHH) SZZ1(IHH)=ZZ1(IHH) SGMG1(IHH)=GMAG1(IHH) SCSB1(IHH)=CSBT1(IHH) SRF1(IHH)=RF1(IHH) SFR1(IHH)=FR1(IHH) 601 CONTINUE DO 605 IHH=1,M2H SRVQ(IHH)=RVQ(IHH) SGRXQ(IHH)=GRXQ(IHH) SGRYQ(IHH)=GRYQ(IHH) SGRZQ(IHH)=GRZQ(IHH) SLMP2(IHH)=SLUMP2(IHH) SGLM2(IHH)=GLUMP2(IHH) SGRV2(IHH)=GRV2(IHH) SXX2(IHH)=XX2(IHH) SYY2(IHH)=YY2(IHH) SZZ2(IHH)=ZZ2(IHH) SGMG2(IHH)=GMAG2(IHH) SCSB2(IHH)=CSBT2(IHH) SRF2(IHH)=RF2(IHH) SFR2(IHH)=FR2(IHH) 605 CONTINUE GOTO 602 4120 CONTINUE EVOL1=VOL1 EVOL2=VOL2 ESUM1=SUMM1 ESUM2=SUMM2 ESM1=SM1 ESM2=SM2 DO 858 IHH=1,MTLH 858 ETLDH(IHH)=TLDH(IHH) DO 1601 IHH=1,M1H ERV(IHH)=RV(IHH) EGRX(IHH)=GRX(IHH) EGRY(IHH)=GRY(IHH) EGRZ(IHH)=GRZ(IHH) ELMP1(IHH)=SLUMP1(IHH) EGLM1(IHH)=GLUMP1(IHH) EGRV1(IHH)=GRV1(IHH) EXX1(IHH)=XX1(IHH) EYY1(IHH)=YY1(IHH) EZZ1(IHH)=ZZ1(IHH) EGMG1(IHH)=GMAG1(IHH) ECSB1(IHH)=CSBT1(IHH) ERF1(IHH)=RF1(IHH) EFR1(IHH)=FR1(IHH) 1601 CONTINUE DO 1605 IHH=1,M2H ERVQ(IHH)=RVQ(IHH) EGRXQ(IHH)=GRXQ(IHH) EGRYQ(IHH)=GRYQ(IHH) EGRZQ(IHH)=GRZQ(IHH) ELMP2(IHH)=SLUMP2(IHH) EGLM2(IHH)=GLUMP2(IHH) EGRV2(IHH)=GRV2(IHH) EXX2(IHH)=XX2(IHH) EYY2(IHH)=YY2(IHH) EZZ2(IHH)=ZZ2(IHH) EGMG2(IHH)=GMAG2(IHH) ECSB2(IHH)=CSBT2(IHH) ERF2(IHH)=RF2(IHH) EFR2(IHH)=FR2(IHH) 1605 CONTINUE 602 CONTINUE HTT=HOT IF(MODE.EQ.-1) HTT=0.d0 XR=(HTT+COOL+EL3A(IB))*ELIT+VKM1*VC1+VKM2*VC2 IF(KSR.NE.1) GOTO 710 BL(IX)=XR GOTO 420 710 CONTINUE IF(KSR.NE.2) GOTO 711 BR(IX)=XR II=NMAT+IX OBS(II)=FLUX(IX)-XR if(iss.eq.2) goto 319 GOTO 420 711 CONTINUE XLR=BR(IX) IF(LOW(KH).EQ.1) XLR=BL(IX) IF(IL.NE.2) GOTO 388 XLR=XR GOTO 87 388 XNR=XR 87 khkeep=kh*keep(kh) if(khkeep.ne.11.and.khkeep.ne.14) OBS(II)=(XNR-XLR)/DEL(KH) if(kh.eq.11) dfdap(ix)=(xnr-xlr)/del(kh) if(kh.eq.14) dfdph(ix)=(xnr-xlr)/del(kh) 319 CONTINUE GOTO 420 300 IF(IB.LE.NVC) OBS(II)=(BR(IX)-VGA)/A GOTO 420 308 IF(IB.LE.NVC) OBS(II)=1.d0 GOTO 420 844 obs(ii)=dfdph(ix)/(period+dtime*dpdt) goto 420 845 brac=period+dtime*dpdt obs(ii)=dfdph(ix)*dtime/(brac*period) goto 420 846 dis=dabs(dtime*dpdt/period) if(dis.gt.toldis) goto 848 brac2=2.d0*period+dtime*dpdt dphpd=-2.d0*(dtime/brac2)**2 $+tt*dtime**3*(2.d0*brac2*dpdt-3.d0*dpdt**2*dtime)/brac2**4 goto 849 848 brac=period+dtime*dpdt dphpd=dtime/(brac*dpdt)-(dlog(brac)-dlog(period))/dpdt**2 849 obs(ii)=-dfdph(ix)*dphpd goto 420 847 obs(ii)=dtime*dfdap(ix) goto 420 301 IF(IB.GT.NVC) OBS(II)=1.d0 420 CONTINUE 418 CONTINUE 417 CONTINUE 419 CONTINUE write(6,101) write(6,101) write(6,101) write(6,159) write(6,101) write(6,169) do 298 icv=1,nlvc icvp=icv+1 nbs=knobs(icvp)-knobs(icv) iw=knobs(icv)+1 jstart=nmat+iw jstop=jstart+nbs-1 resq=0.d0 iww=iw-1 do 299 jres=jstart,jstop iww=iww+1 299 resq=resq+wt(iww)*obs(jres)**2 write(6,119) icv,nbs,resq 298 continue write(6,101) write(6,101) GOTO 65 71 JF=0 IF(KO.EQ.0) stop DO 261 J=1,NCOF 261 OBS(J)=HOLD(J) IF(KDISK.EQ.0) GOTO 72 REWIND 9 READ(9,67)(OBS(J),J=1,NCOF) 72 READ(5,20)(KEEP(I),I=1,35),IFDER,IFM,IFR,xlamda IF(KEEP(1).EQ.2) stop DO 232 I=1,35 232 IF(KEP(I).EQ.1) KEEP(I)=1 NOBS=KNOBS(LC1) matrix=31 do 428 ima=1,30 428 matrix=matrix-keep(ima) matrix=matrix+nlc*(5-keep(31)-keep(32)-keep(33)-keep(34)-keep(35)) MAT=MATRIX-1 EM=dfloat(MATRIX-15) KTR=.24d0*EM+2.2d0 IF(EM.LE.1.5d0) KTR=1 IF(EM.GT.12.d0) KTR=5 NCOEFF=MATRIX*NOBS KC=1 NSHIFT(1)=0 DO 59 I=2,36 IF(I.GT.31) KC=NLC KE=0 J=I-1 IF(KEEP(J).GT.KEP(J)) KE=1 59 NSHIFT(I)=NOBS*KE*KC+NSHIFT(J) NOBBS=NOBS DO 30 I=1,36 IF(KEEP(I).EQ.1) GOTO 30 IF(I.GT.30) NOBBS=NOBS*NLC IF(I.EQ.36) NOBBS=NOBS DO 32 J=1,NOBBS JF=JF+1 KX=JF+NSHIFT(I) 32 OBS(JF)=OBS(KX) 30 CONTINUE 65 WRITE(6,102) WRITE(6,20)(KEEP(I),I=1,35),IFDER,IFM,IFR,xlamda NOBS=KNOBS(LC1) WRITE(6,101) IF(IFDER.EQ.0) KTR=5 WRITE(6,82) WRITE(6,101) DO 96 IB=1,NLVC IST=KNOBS(IB)+1 IB1=IB+1 ISP=KNOBS(IB1) DO 96 I=IST,ISP GOTO(5,6,7,8,96),KTR 5 WRITE(6,15)(OBS(J),J=I,NCOEFF,NOBS) GOTO 96 6 WRITE(6,16)(OBS(J),J=I,NCOEFF,NOBS) GOTO 96 7 WRITE(6,17)(OBS(J),J=I,NCOEFF,NOBS) GOTO 96 8 WRITE(6,19)(OBS(J),J=I,NCOEFF,NOBS) 96 CONTINUE 78 CONTINUE IF(KO.LE.1) GOTO 70 IF(IBEF.EQ.1) GOTO 70 DO 62 J=1,NCOEFF 62 HOLD(J)=OBS(J) IF(KDISK.EQ.0) GOTO 73 REWIND 9 WRITE(9,67)(OBS(J),J=1,NCOEFF) 73 CONTINUE 70 WRITE(6,101) DO 97 IB=1,NLVC IST=KNOBS(IB)+1 IB1=IB+1 ISP=KNOBS(IB1) NOIS=NOISE(IB) DO 97 I=IST,ISP IF(IB.GT.NVC) GOTO 444 ROOTWT=dsqrt(WT(I))/(100.d0*SIGMA(IB)) GOTO 445 444 ROOTWT=dsqrt(WT(I))/(100.d0*SIGMA(IB)*dsqrt(FLUX(I))**NOIS) 445 CONTINUE DO 97 LOB=I,NCOEFF,NOBS 97 OBS(LOB)=OBS(LOB)*ROOTWT IF(IFDER.NE.0) WRITE(6,83) IF(IFDER.NE.0) WRITE(6,101) DO 98 I=1,NOBS GOTO(45,46,47,48,98),KTR 45 WRITE(6,15)(OBS(J),J=I,NCOEFF,NOBS) GOTO 98 46 WRITE(6,16)(OBS(J),J=I,NCOEFF,NOBS) GOTO 98 47 WRITE(6,17)(OBS(J),J=I,NCOEFF,NOBS) GOTO 98 48 WRITE(6,19)(OBS(J),J=I,NCOEFF,NOBS) 98 CONTINUE 77 CALL xxsqar(OBS,NOBS,MAT,OUT,sd,xlamda,deter,CN,CNN,cnc,clc, $s,ccl,ll,mm) MSQ=MAT*MAT IF(IFM.EQ.0) GOTO 436 WRITE(6,101) WRITE(6,181) WRITE(6,101) DO 38 JR=1,MAT 38 WRITE(6,37) (CN(JX),JX=JR,MSQ,MAT),CCL(JR) WRITE(6,101) WRITE(6,183) WRITE(6,101) 436 CONTINUE NO1=23 NO2=24 NRM=25 DO 334 IRM=1,24 IF(IRM.LE.23) NO1=NO1-KEEP(IRM) NO2=NO2-KEEP(IRM) 334 NRM=NRM-KEEP(IRM) CORO1=1-KEEP(23) CORO2=1-KEEP(24) CORQ=1-KEEP(25) DO 34 JM=1,MAT DO 33 JQ=1,MAT JT=JM+MAT*(JQ-1) IJM=(MAT+1)*(JM-1)+1 IJQ=(MAT+1)*(JQ-1)+1 33 V(JQ)=CNN(JT)/DSQRT(CNN(IJM)*CNN(IJQ)) co1q=0.d0 co2q=0.d0 if(jm.eq.nrm.and.no1.gt.0) co1q=v(no1)*corq*coro1 if(jm.eq.nrm.and.no2.gt.0) co2q=v(no2)*corq*coro2 34 WRITE(6,37)(V(IM),IM=1,MAT) IF(IFM.EQ.0) GOTO 36 WRITE(6,101) WRITE(6,184) WRITE(6,101) CALL DGMPRD(CN,CNN,CNOUT,MAT,MAT,MAT) DO 116 J8=1,MAT 116 WRITE(6,37)(CNOUT(J7),J7=J8,MSQ,MAT) WRITE(6,101) WRITE(6,185) WRITE(6,101) ANSCH=0.D0 DO 118 J5=1,MAT V(J5)=0.D0 DO 117 J6=1,MAT int=mat*(j6-1) idi=j6+int I9=J5+int 117 V(J5)=OUT(J6)*CN(I9)*dsqrt(cnc(idi))+V(J5) ERR=V(J5)-CCL(J5) 118 ANSCH=ANSCH+DABS(ERR) WRITE(6,137)(V(J4),J4=1,MAT) WRITE(6,101) WRITE(6,138) ANSCH 36 CONTINUE WRITE(6,101) WRITE(6,101) write(6,715) WRITE(6,101) write(6,43) iout=0 do 93 kpar=1,35 imax=1 if(kpar.gt.30) imax=nlc if(keep(kpar).eq.1) goto 93 do 94 kurv=1,imax kcurv=kurv if(kpar.le.30) kcurv=0 iout=iout+1 ipar=kpar if(kpar.gt.30) ipar=30+kurv+(kpar-31)*nlc parout=para(ipar)+out(iout) write(6,615) kpar,kcurv,para(ipar),out(iout),parout,sd(iout) 94 continue 93 continue WRITE(6,101) WRITE(6,101) write(6,716) WRITE(6,101) write(6,43) iout=0 do 53 kpar=1,35 imax=1 if(kpar.gt.30) imax=nlc if(keep(kpar).eq.1) goto 53 do 52 kurv=1,imax kcurv=kurv if(kpar.le.30) kcurv=0 iout=iout+1 ipar=kpar if(kpar.gt.30) ipar=30+kurv+(kpar-31)*nlc parout=para(ipar)+out(iout) write(6,616) kpar,kcurv,para(ipar),out(iout),parout,sd(iout) 52 continue 53 continue WRITE(6,101) 88 RESSQ=0.d0 JST=MAT*NOBS+1 DO 199 JRES=JST,NCOEFF 199 RESSQ=RESSQ+OBS(JRES)**2 WRITE(6,101) WRITE(6,40) WRITE(6,21) RESSQ,S,deter IBEF=1 IF(IFR.EQ.0) GOTO 71 WRITE(6,102) WRITE(6,101) WRITE(6,101) WRITE(6,650) WRITE(6,101) WRITE(6,101) WRITE(6,653) WRITE(6,101) DO1=sd(NO1)*CORO1 DO2=sd(NO2)*CORO2 if(mod.eq.1) do2=do1 if(mod.eq.1) co2q=co1q DQ=sd(NRM)*CORQ COQ=CO1Q F=F1 DP=1.d0-E OME=PHSV DOM=DO1 KOMP=0 925 CONTINUE KOMP=KOMP+1 DO 926 KD=1,4 if(kd.ne.2) goto 928 if(po(komp).ge.omcr(komp)) goto 928 goto 926 928 continue TH=XTHA(KD) FI=XFIA(KD) CALL ROMQ(OME,RM,F,DP,E,TH,FI,R,DRDO,DRDQ,DODQ,KOMP,MODE) DR=dsqrt(DRDQ**2*DQ**2+DRDO**2*DOM**2+2.d0*COQ*DRDQ*DRDO*DQ*DOM) WRITE(6,654)KOMP,ARAD(KD),R,DRDO,DRDQ,DR 926 CONTINUE DO2DQ=DODQ IF(KOMP.EQ.1)DO1DQ=DODQ COQ=CO2Q F=F2 OME=PCSV DOM=DO2 WRITE(6,101) IF(KOMP.EQ.1) GOTO 925 WRITE(6,101) WRITE(6,651) IF(KOMP.EQ.2) WRITE(6,652)DO1DQ,DO2DQ,CO1Q,CO2Q,DO1,DO2,DQ GOTO 71 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