C From interf!uunet!sbast1.ess.sunysb.edu!DPETERSON Sat Sep 12 08:44:15 1992 C Return-Path: C Received: from interf.UUCP by pleiades (4.1/SMI-4.1) C id AA19778; Sat, 12 Sep 92 08:44:13 EDT C From: interf!uunet!sbast1.ess.sunysb.edu!DPETERSON C Received: from uunet.UUCP by (4.0/SMI-4.0) C id AA07524; Sat, 12 Sep 92 08:48:48 EDT C Received: from sbast1.ess.sunysb.edu by relay1.UU.NET with SMTP C (5.61/UUNET-internet-primary) id AA21377; Sat, 12 Sep 92 08:49:43 -0400 C Date: Sat, 12 Sep 1992 8:49:55 EDT C Message-Id: <920912084955.34a00c27@sbast1.ess.sunysb.edu> C To: atlas!tarmstr@interf.com C X-Vmsmail-To: SMTP%"atlas!tarmstr@interf.com" C Status: R C From Deane Peterson (SUNY Stony Brook) 12 Sept. 1992 PROGRAM BINPAC C PROGRAM TO DERIVE A BINARY ORBIT COMMON /ANOMAL/ TRUEA,ECCA,MEANA,CYCLES COMMON /ELEMNT/ PERIOD,TPERIA,AAU,PARALX,INCL,ECCEN,ARGPER,LGNODE, 1 M2M1,VSYS(6),KAMP(2) REAL INCL,LGNODE,M2M1,KAMP,MEANA COMMON /EQNORM/ XNORM(1444),DRIVE(38),CHISQ,NORD,MAPDUM(38) COMMON /TRIG/ COSV,SINV,COSW,SINW,COSVW,SINVW,COSI,SINI C RADIAL VELOCITY DATA DIMENSION VOBS(400,3),SIGV(400,3),TVOBS(400),IDVOBS(400),SIGRAD(3) DIMENSION SIG3V(3,6) C DATA EXCHANGED WITH SPECTROSCOPIC CALCULATION DIMENSION VR(2),DVDEL(6,2) C STORAGE FOR LEAST SQUARES C NOTE THAT ICNVRG NOW CONTAINS ROOM FOR A SWITCH (ICNVRG(16)): C WHEN = 1 THE NEXT RECORD CONTAINS SIGMAS FOR THE VELOCITIES AS C A FUNCTION OF OBSERVER C THE SIGMA ASSIGNED TO EACH VELOCITY IS TAKEN IN THE FOLLOWING PRIORITY: C 1. THE SIGMA APPEARING AFTER THE INDIVIDUAL OBSERVATIONS, IF NE.0 C 2. THE SIGMA FROM THE SET READ BY ICNVRG(16), IF AVAILABLE C 3. THE SIGMA FROM THE THE VELOCITY HEADER RECORD C 4. IF AFTER ALL THAT THE SIGMA REMAINS ZERO, THE OBSERVATION C IS IGNORED IN THE CALCULATION C C AN EPHEMERIS OPTION HAS BEEN ADDED TO THE END SEQUENCE C Least Squares Matrix and Vectors to Map Variables In and Out. DIMENSION ICNVRG(16),COREC(15),PARTAL(15),CORMTX(225),SIGMA(15), 1SCALE(15),COEFDR(15),CSIGMA(15),MAPCVG(15) C BOOKKEEPING DIMENSION IFDATA(3),NDATA(3),ELEM(15),ELNAME(15),ELINTL(15), 1HNAME(15) EQUIVALENCE (ELEM(1),PERIOD) C STORAGE FOR ERROR PROPAGATION DIMENSION INDX(15),PFPTHT(15),INDXMS(3),PFPTM1(3),PFPTM2(3) C STORAGE FOR VISUAL DATA DIMENSION TVISOB(300),RHOOBS(300),SIGRHO(300),THTOBS(300), 1SIGTH(300) C STORAGE FOR OCCULTATION DATA DIMENSION TLUOCC(20),RHOASP(20),ASPECT(20),SIGRHA(20), 1 SIGASP(20),COROAS(20),SIGRA(20) C STORAGE FOR RESIDUALS DIMENSION RESVR(400,3),RESRHO(300),RESTH(300),RESRHA(20), 1RVR(3),SVR(3) C STORAGE FOR THE EPHEMERIS CALCULATION DIMENSION DELDAY(100) CHARACTER*39 FILNAM C OBSERVER ID'S CHARACTER*4 VOBSR(300),LOBSR(20),VO,LO,BLANK,AKEY*1 REAL*8 TJD,DPERIA,DIME DATA TJD/2440000.0D0/,SCALE/4*1.,6.28,1.,2*6.28,7*1./ DATA TWOPI/6.283185/,T2000/11544.533/ DATA RADIAN /57.29578/,CONVRG/1.E-4/ DATA ELNAME/' P ',' TP ',' AU ',' PI ','INCL','ECC ',' W ', 1'OMGA','M2M1','V0 1','V0 2','V0 3','V0 4','V0 5','V0 6'/ DATA INDX/1,3,9,5,6,10*0/,INDXMS/1,3,9/,BLANK/' '/ DATA NDIM/15/ C FUNCTION CONVERTS BESSELIAN YEARS TO JD-2440000 YRTOJD(X)=365.2422*(X-1950.0)-6717.577 C Let's go 10 PRINT *, 'FILENAME' READ 15,FILNAM 15 FORMAT (A39) IF (FILNAM.EQ.'EXIT'.OR.FILNAM.EQ.'END') STOP IF (FILNAM.EQ.'exit'.OR.FILNAM.EQ.'end') STOP OPEN (UNIT=7,FILE=FILNAM,STATUS='OLD') C READ INITIAL ELEMENTS READ (7,20,END=1000) PERIOD,DPERIA,(ELEM(I),I=3,9),PRECES C P TP AU PI I ECC W,OM M2M1 PRECES 20 FORMAT (F10.4,F12.3,F8.4,F5.3,F6.1,F5.3,2F6.1,F5.3,F7.4) READ (7,21,END=1000) (ELEM(I),I=10,15) C VSYS 21 FORMAT (6F6.1) C CONVERT ANGLES TO RADIAN AND SUBTRACT TJD FROM DATES TPERIA=DPERIA-TJD INCL=INCL/RADIAN ARGPER=ARGPER/RADIAN LGNODE=LGNODE/RADIAN C PRECESSION CHANGED FROM PER ANNUM TO PER DAY PRECES=PRECES/(RADIAN*365.2422) C SET UP SCALING PARAMETERS SCALE(1)=PERIOD SCALE(2)=PERIOD SCALE(3)=AAU SCALE(4)=PARALX DO 25 I=1,6 25 SCALE(9+I)=VSYS(I) C READ SWITCHES INDICATING WHICH ELEMENTS TO CONVERGE READ (7,30) ICNVRG 30 FORMAT (16I2) C WRITE (*,'(1X,16I2)') ICNVRG NAUTHR=MAX(1,ICNVRG(16)) NORD=0 C SOLVIT REQUIRES A "FILLED" MATRIX. THIS SETS UP THE INDICES THAT C MAP THE ELEMENTS INTO AND OUT OF THAT MATRIX. DO 40 I=1,NDIM MAPDUM(I)=I ELINTL(I)=ELEM(I) IF (ICNVRG(I).EQ.0) GO TO 40 NORD=NORD+1 MAPCVG(NORD)=I ICNVRG(I)=NORD HNAME(NORD)=ELNAME(I) 40 CONTINUE DO 50 I=1,3 NDATA(I)=0 50 IFDATA(I)=0 C IF ICNVRG(16) GT 0, GET 3 SIGMAS FOR EACH VSYS IF (ICNVRG(16).GT.0) READ (7,55) ((SIG3V(I,J),I=1,3),J=1,NAUTHR) 55 FORMAT (18F4.1) C IF (ICNVRG(16).GT.0) WRITE (*,'(1X,18F4.1)') ((SIG3V(I,J),I=1,3), C 1J=1,NAUTHR) C READ HEADER CARD 60 READ (7,70) ID,ND,SIGRAD 70 FORMAT (2I5,3F5.3) IF (ID.LE.0.OR.ID.GT.3) GO TO 200 C ANOTHER KIND OF DATA IFDATA(ID)=1 NDATA(ID)=ND C ID=1, RADIAL VELOCITIES C ID=2, "VISUAL" OBSERVATIONS C ID=3, LUNAR OCCULTATIONS GO TO (80,110,140),ID C READ IN RADIAL VELOCITY DATA 80 DO 100 N=1,ND READ (7,90) IDVOBS(N),DIME,(VOBS(N,I),I=1,3),(SIGV(N,I),I=1,3) 90 FORMAT (I3,F12.3,3F7.2,3F5.1) TVOBS(N)=DIME-TJD C IF SMALL, GOT BESSELIAN DATE IF (DIME.LE.2500.) TVOBS(N)=YRTOJD(SNGL(DIME)) DO 100 I=1,3 IF (SIGV(N,I).NE.0.) GO TO 95 IF (ICNVRG(16).GT.0) SIGV(N,I)=SIG3V(I,IDVOBS(N)) IF (SIGV(N,I).EQ.0.) SIGV(N,I)=SIGRAD(I) 95 IF (VOBS(N,I).EQ.0) SIGV(N,I)=0. 100 CONTINUE GO TO 60 C READ IN VISUAL DATA 110 DO 130 N=1,ND READ (7,120) VOBSR(N),DIME,RHOOBS(N),SIGRHO(N),THTOBS(N), 1 SIGTH(N) 120 FORMAT (A4,F11.3,2F6.3,2F6.1) TVISOB(N)=DIME-TJD IF (DIME.LE.2500.) TVISOB(N)=YRTOJD(SNGL(DIME)) THTOBS(N)=THTOBS(N)/RADIAN SIGTH(N)=SIGTH(N)/RADIAN 130 CONTINUE GO TO 60 C READ IN LUNAR OCCULTATION DATA 140 DO 160 N=1,ND READ (7,150) LOBSR(N),DIME,RHOASP(N),ASPECT(N),SIGRHA(N), 1 SIGASP(N),COROAS(N) 150 FORMAT (A4,F11.3,F8.5,F7.2,F8.5,F6.2,F7.4) TLUOCC(N)=DIME-TJD IF (DIME.LE.2500.) TLUOCC(N)=YRTOJD(SNGL(DIME)) ASPECT(N)=ASPECT(N)/RADIAN SIGASP(N)=SIGASP(N)/RADIAN 160 CONTINUE C SET UP THE NLLS RUN C NUMBER OF ITERATIONS 200 CLOSE (UNIT=7) IF (IFDATA(1)+IFDATA(2)+IFDATA(3).EQ.0) GO TO 800 WRITE (*,'('' ITERATIONS =''\)') ITERAT=1 READ *,NITER IF (NITER.GT.0) ITERAT=NITER C THE ITERATIONS DO 440 ITER=1,ITERAT NUMEQC=0 C CLEAR ACCUMULATORS CALL ACUMLS(PARTAL,RESID,SIG,1) C DO OVER VARIOUS TYPES OF DATA DO 400 ITYP=1,3 ID=IFDATA(ITYP) IF (ID.EQ.0) GO TO 400 ND=NDATA(ITYP) GO TO (220,250,300),ITYP C DO RADIAL VELOCITY DATA 220 DO 230 I=1,NDIM 230 PARTAL(I)=0. KAMP(2)=10879.*AAU*SIN(INCL)/PERIOD/ 1 SQRT(1.-ECCEN**2)/(1.+M2M1) KAMP(1)=KAMP(2)*M2M1 DO 240 N=1,ND CALL KEPLER (TVOBS(N)) CALL SPECTR(VR,DVDEL,IDVOBS(N)) DO 240 ICOMP=1,3 RESVR(N,ICOMP)=0. IF (SIGV(N,ICOMP).EQ.0.) GO TO 240 DO 231 I=1,NDIM 231 COEFDR(I)=0. IF (ICOMP.EQ.3) GO TO 232 C VELOCITIES OF INDIVIDUAL COMPONENTS C DO INDIVIDUAL COMPONENTS DV=VR(ICOMP)-VSYS(IDVOBS(N)) RESID=VOBS(N,ICOMP)-VR(ICOMP) C PERIOD COEFDR(1)=-DV/PERIOD+DVDEL(6,ICOMP) C PERIASTRON EPOCH COEFDR(2)=DVDEL(5,ICOMP) C A(AU) COEFDR(3)=DV/AAU C INCLINATION COEFDR(5)=DV/TAN(INCL) C ECCENTRICITY COEFDR(6)=DV*ECCEN/(1.-ECCEN**2)+DVDEL(3,ICOMP) C ARGUMENT OF PERIASTRON COEFDR(7)=DVDEL(4,ICOMP) C MASS RATIO COEFDR(9)=-DV/(1.+M2M1) IF (ICOMP.EQ.1) COEFDR(9)=-COEFDR(9)/M2M1 C SYSTEMIC VELOCITY COEFDR(9+IDVOBS(N))=DVDEL(1,ICOMP) GO TO 233 C DO RELATIVE VELOCITY 232 DV=VR(2)-VR(1) RESID=VOBS(N,3)-DV OPM2M1=1.+M2M1 C PERIOD COEFDR(1)=-DV/PERIOD+OPM2M1*DVDEL(6,2) C PERIASTRON PASSAGE COEFDR(2)=OPM2M1*DVDEL(5,2) C A(AU) COEFDR(3)=DV/AAU C INCLINATION COEFDR(5)=DV/TAN(INCL) C ECCENTRICITY COEFDR(6)=DV*ECCEN/(1.-ECCEN**2)+OPM2M1*DVDEL(3,2) C ARGUMENT OF PERIASTRON COEFDR(7)=OPM2M1*DVDEL(4,2) 233 NUMEQC=NUMEQC+1 RESVR(N,ICOMP)=RESID DO 235 I=1,NDIM IN=ICNVRG(I) IF (IN.EQ.0) GO TO 235 C PRESCALE THE LEAST SQUARES MATRIX PARTAL(IN)=COEFDR(I)*SCALE(I) 235 CONTINUE CALL ACUMLS(PARTAL,RESID,SIGV(N,ICOMP),2) 240 CONTINUE GO TO 400 C VISUAL BINARY DATA 250 DO 290 N=1,ND CALL KEPLER(TVISOB(N)) C FIRST DO SEPARATIONS IF (SIGRHO(N).EQ.0.) GO TO 270 NUMEQC=NUMEQC+1 CALL RHO(RHOC,COEFDR) DO 260 I=1,NDIM IC=ICNVRG(I) IF (IC.EQ.0) GO TO 260 PARTAL(IC)=COEFDR(I)*SCALE(I) 260 CONTINUE RESID=RHOOBS(N)-RHOC RESRHO(N)=RESID CALL ACUMLS(PARTAL,RESID,SIGRHO(N),2) C NEXT ANGLE 270 IF(SIGTH(N).EQ.0.) GO TO 290 NUMEQC=NUMEQC+1 CALL THETA(THC,COEFDR) C CORRECT THETA TO T=2000 RESID=THTOBS(N)-THC+PRECES*(TVISOB(N)-T2000) C CHECK IF SUBTRACTING THROUGH 2PI IF (ABS(ABS(RESID)-TWOPI).GT.1.5708) GO TO 275 RES=RESID IF (RES.LT.0.) RESID=RESID+TWOPI IF (RES.GT.0.) RESID=RESID-TWOPI 275 DO 280 I=1,NDIM IC=ICNVRG(I) IF (IC.EQ.0) GO TO 280 PARTAL(IC)=COEFDR(I)*SCALE(I) 280 CONTINUE RESTH(N)=RESID CALL ACUMLS(PARTAL,RESID,SIGTH(N),2) 290 CONTINUE GO TO 400 C LUNAR OCCULTATION DATA 300 DO 310 I=1,NDIM 310 PARTAL(I)=0. DO 330 N=1,ND CALL KEPLER(TLUOCC(N)) C CORRECT ASPECT TO T=2000 ASPCOR=ASPECT(N)+PRECES*(TVISOB(N)-T2000) CALL LUNOCC(ASPCOR,RHOAC,RHOSIN,COEFDR) NUMEQC=NUMEQC+1 DO 320 I=1,NDIM IC=ICNVRG(I) IF (IC.EQ.0) GO TO 320 PARTAL(IC)=COEFDR(I)*SCALE(I) 320 CONTINUE RESID=RHOASP(N)-RHOAC RESRHA(N)=RESID SIGRA(N)=SQRT(SIGRHA(N)**2+(SIGASP(N)*RHOSIN)**2+2.* 1 SIGRHA(N)*SIGASP(N)*RHOSIN*COROAS(N)) SIGRA(N)=AMAX1(SIGRA(N),.001) CALL ACUMLS(PARTAL,RESID,SIGRA(N),2) 330 CONTINUE GO TO 400 C ECLIPSE DATA, WHEN I HAVE THE TIME 400 CONTINUE IF (NITER.EQ.0) GO TO 620 C SOLVE FOR CORRECTIONS CALL SOLLSQ(COREC) DEGFRE=NUMEQC-NORD IFERR=0 WRITE (6,410) ITER,CHISQ,DEGFRE 410 FORMAT ('0ITER=',I3,5X,'CHISQ=',1PE11.4,5X, 1 'DEGREES OF FREEDOM=', 2 0PF5.0/' ELEMENT',8X,'CORRECTION') DO 430 I=1,NDIM IC=ICNVRG(I) IF (IC.EQ.0) GO TO 430 IF (ABS(COREC(IC)).GT.CONVRG) IFERR=1 C UNDO THE SCALING COREC(IC)=COREC(IC)*SCALE(I) ELEM(I)=ELEM(I)+COREC(IC) 415 WRITE (6,420) ELNAME(I),ELEM(I),COREC(IC) 420 FORMAT (1X,A4,1PE12.4,E13.3) 430 CONTINUE WRITE (*,'('' PAUSE: Hit any key''\)') READ (*,'(A1)') AKEY IF (IFERR.LE.0) GO TO 460 440 CONTINUE WRITE (6,450) 450 FORMAT (' DID NOT CONVERGE') 460 WRITE (6,470) 470 FORMAT ('0',8X,'FINAL VALUES',12X,'ERRORS'/' ELEM INITIAL FIN 1AL',7X,'SIGMA',6X,'REL'/) CALL INVRT(CSIGMA,CORMTX) UNITER=SQRT(CHISQ/DEGFRE) DO 480 I=1,NDIM 480 SIGMA(I)=0. C MAP THE ELEMENTS BACK TO THE ORIGINAL ORDER DO 490 I=1,NORD M=MAPCVG(I) 490 SIGMA(M)=CSIGMA(I) C WRITE OUT THE CONVERGED ELEMENTS DO 510 I=1,NDIM C UNDO THE SCALING SIGMA(I)=SIGMA(I)*SCALE(I)*UNITER REL=0. IF (ELEM(I).NE.0.) REL=ABS(SIGMA(I)/ELEM(I)) WRITE (6,500) ELNAME(I),ELINTL(I),ELEM(I),SIGMA(I),REL 500 FORMAT (1X,A4,1PE11.4,E13.6,E10.3,0PF10.7) 510 CONTINUE DPERIA=ELEM(2)+TJD WRITE (6,520) DPERIA,SIGMA(2) 520 FORMAT (' FINAL T(PERIASTRON)= ',F12.3,'+/-',F7.3) WRITE (*,'('' PAUSE: Hit any key''\)') READ (*,'(A1)') AKEY C WRITE OUT THE CORRELATION MATRIX WRITE (6,530) (HNAME(I),I=1,NORD) 530 FORMAT ('0',10X,'CORRELATION MATRIX'/' ELEM ',15(2X,A4)) DO 550 N=1,NORD IO=NORD*(N-1) WRITE (6,540) HNAME(N),(CORMTX(IO+I),I=1,NORD) 540 FORMAT (1X,A4,2X,15F6.3) 550 CONTINUE C SOLVE FOR AUXILLARY QUANTITIES C FIRST MASSES XM1=AAU**3/(PERIOD/365.24)**2/(1.+M2M1) XM2=XM1*M2M1 PFPTM1(1)=-2.*SIGMA(1)/PERIOD PFPTM2(1)=PFPTM1(1) PFPTM1(2)=3.*SIGMA(3)/AAU PFPTM2(2)=PFPTM1(2) PFPTM1(3)=-SIGMA(9)/(1.+M2M1) PFPTM2(3)=SIGMA(9)/M2M1/(1.+M2M1) SSIGM1=SQRT(CCRLAT(INDXMS,PFPTM1,3,INDXMS,PFPTM1,3,CORMTX,ICNVRG, 1NORD)) SSIGM2=SQRT(CCRLAT(INDXMS,PFPTM2,3,INDXMS,PFPTM2,3,CORMTX,ICNVRG, 1NORD)) RHOM12=CCRLAT(INDXMS,PFPTM1,3,INDXMS,PFPTM2,3,CORMTX,ICNVRG,NORD)/ 1SSIGM1/SSIGM2 C NEXT, VELOCITY AMPLITUDES KAMP(2)=10879.*AAU*SIN(INCL)/PERIOD/SQRT(1.-ECCEN**2)/(1.+M2M1) KAMP(1)=KAMP(2)*M2M1 PFPTHT(1)=-SIGMA(1)/PERIOD PFPTHT(2)=SIGMA(3)/AAU PFPTHT(3)=-1./(1.+M2M1)*SIGMA(9) PFPTHT(4)=SIGMA(5)*COSI/SINI PFPTHT(5)=SIGMA(6)*ECCEN/(1.-ECCEN**2) SIGK2=PRPGAT(INDX,PFPTHT,5,CORMTX,ICNVRG,NORD)*ABS(KAMP(2)) PFPTHT(3)=-1./M2M1*PFPTHT(3) SIGK1=PRPGAT(INDX,PFPTHT,5,CORMTX,ICNVRG,NORD)*ABS(KAMP(1)) WRITE (6,560) XM1,SSIGM1*XM1,XM2,SSIGM2*XM2,RHOM12,KAMP(1),SIGK1, 1KAMP(2),SIGK2 560 FORMAT ('0M1=',F6.3,'+/-',F5.3,' M2=',F6.3,'+/-',F5.3,' RHO(M1,M 12)=',F7.4/' K1=',F6.2,'+/-',F5.2,' K2=',F6.2,'+/-',F5.2) C NOW, THE NICETIES WRITE (*,'('' DO YOU WANT TO SAVE THIS? [N]'')') READ (*,'(A1)') AKEY IF (AKEY.NE.'Y'.AND.AKEY.NE.'y') GO TO 620 PRINT *, 'FILENAME: ' READ 15,FILNAM OPEN (UNIT=7,FILE=FILNAM,STATUS='NEW') C SAVE CONVERGED PARAMETERS WRITE (7,580) (ELEM(I),SIGMA(I),I=1,15) 580 FORMAT (1PE14.6,E11.3,E14.6,E11.3,E14.6,E11.3) WRITE (7,590) ICNVRG,NORD 590 FORMAT (15I2,I5) DO 610 N=1,NORD NO=NORD*(N-1) WRITE (7,600) (CORMTX(NO+I),I=1,NORD) 600 FORMAT (15F5.3) 610 CONTINUE CLOSE (UNIT=7) C PRINT OUT RESIDUALS 620 WRITE ( *,'('' DO YOU WANT RESIDUALS? [N] '')') READ (*,'(A1)') AKEY IF (AKEY.NE.'Y'.AND.AKEY.NE.'y') GO TO 800 C RADIAL VELOCITIES IF (NDATA(1).EQ.0) GO TO 670 WRITE (*,'('' WILL YOU WANT THIS SAVED TO A FILE [N] ''\)') READ (*,'(A1)') AKEY IF (AKEY.EQ.'Y'.OR.AKEY.EQ.'y') THEN IFPLOT=1 WRITE (*,'('' NAME FOR FILE: '')') READ (*,'(A39)') FILNAM OPEN (UNIT=7,FILE=FILNAM,STATUS='UNKNOWN') ELSE IFPLOT=0 END IF WRITE (6,630) 630 FORMAT (15X,'R A D I A L V E L O C I T I E S'/ 11X,'OBS PHASE V1 DV1 SIG V2 DV2 SIG VTOT',4X, 2'DVTOT SIG') NPAUSE=0 DO 660 N=1,NDATA(1) NPAUSE=NPAUSE+1 P1=(TVOBS(N)-TPERIA)/PERIOD DO 640 I=1,3 SVR(I)=SIGV(N,I) RVR(I)=RESVR(N,I) IF (VOBS(N,I).NE.0.) GO TO 640 SVR(I)=0. RVR(I)=0. 640 CONTINUE PRINT 650, IDVOBS(N),P1,(VOBS(N,I),RVR(I),SVR(I),I=1,3) 650 FORMAT (I4,F9.2,2(2F7.1,F4.1),2F8.2,F5.2) IF (IFPLOT.EQ.1) THEN WRITE (7,651) IDVOBS(N),P1,(VOBS(N,I),RVR(I),SVR(I),I=1,3) 651 FORMAT (I5,F10.3,2(2F7.1,F4.1),2F8.2,F5.2) END IF IF (NPAUSE.GT.22) THEN WRITE (*,'('' PAUSE: Hit any key ''\)') READ (*,'(A1)') AKEY NPAUSE=0 END IF 660 CONTINUE IF (IFPLOT.EQ.1) CLOSE (UNIT=7) C OPTICAL AND OCCULTATION DATA 670 NDMAX=MAX(NDATA(2),NDATA(3)) IF (NDMAX.LE.0) GO TO 800 IF (NPAUSE.GE.20) THEN WRITE (*,'('' PAUSE: Hit any key ''\)') READ (*,'(A1)') AKEY NPAUSE=1 ELSE NPAUSE=NPAUSE+2 END IF PRINT 680 680 FORMAT (20X,'O P T I C A L',17X,'L U N A R'/ 1' OBS PHASE RHO DRHO SIG THETA DTH SIG', 2' OBS PHASE ASPCT RHOTH DRTH SIG') DO 760 N=1,NDMAX NPAUSE=NPAUSE+1 IF (NDATA(2).LT.N) GO TO 710 VO=VOBSR(N) P2=(TVISOB(N)-TPERIA)/PERIOD IDR=RESRHO(N)*1000. ISR=SIGRHO(N)*1000. IRO=RHOOBS(N)*1000. IF (RHOOBS(N).NE.0.) GO TO 700 IDR=0 ISR=0 IRO=0 700 DT=RESTH(N)*RADIAN ST=SIGTH(N)*RADIAN TH=THTOBS(N)*RADIAN IF (THTOBS(N).NE.0.) GO TO 720 DT=0. ST=0. TH=0. GO TO 720 710 IDR=0 ISR=0 IRO=0 DT=0. ST=0. TH=0. P2=0. VO=BLANK C LUNAR OCCULTATIONS 720 IF (NDATA(3).LT.N) GO TO 730 LO=LOBSR(N) P3=(TLUOCC(N)-TPERIA)/PERIOD DRA=RESRHA(N)*1000. SRA=SIGRA(N)*1000. RA=RHOASP(N)*1000. THA=ASPECT(N)*RADIAN GO TO 740 730 P3=0. DRA=0. SRA=0. RA=0. THA=0. LO=BLANK 740 WRITE (6,750) VO,P2,IRO,IDR,ISR,TH,DT,ST,LO,P3,THA,RA,DRA,SRA 750 FORMAT (1X,A4,F6.2,2I5,I4,2F6.1,F5.1,1X,A4,F6.2,2F6.1,2F5.1) IF (NPAUSE.GT.22) THEN WRITE (*,'('' PAUSE: Hit any key ''\)') READ (*,'(A1)') AKEY NPAUSE=0 END IF 760 CONTINUE C PUT ECLIPSE DATA HERE C CODE FOR EPHEMERIS 800 WRITE (*,'('' DO YOU WANT TO RUN AN EPHEMERIS? [N]''\)') READ (*,'(A1)') AKEY IF (AKEY.NE.'Y'.AND.AKEY.NE.'y') GO TO 10 PRINT *,'DATES HAND ENTERED (1), OR FROM A FILE (0)?' READ *,INO IF (INO.EQ.1) GO TO 809 C GET A FILE PRINT *,'FILE NAME:' READ 15,FILNAM IF (FILNAM.EQ.'EXIT'.OR.FILNAM.EQ.'END') STOP IF (FILNAM.EQ.'exit'.OR.FILNAM.EQ.'end') STOP OPEN (UNIT=7,FILE=FILNAM,STATUS='OLD') N=0 801 N=N+1 READ (7,'(F12.3)') DIME IF (DIME.GE.0.D0) THEN IF (DABS(DIME).GT.10.*PERIOD) THEN DELDAY(N)=DIME-DPERIA ELSE DELDAY(N)=DIME END IF NEPH=N GO TO 801 ELSE CLOSE (UNIT=7) GO TO 820 END IF 809 PRINT *,'GIVE JD OF FIRST EPOCH (0 = TP, .LT.0 FOR EXIT)' READ *,DIME IF (DIME.LT.0) GO TO 1000 DTIME=DIME-DPERIA IF (DIME.EQ.0) DTIME=0. NEPH=1 DELDAY(1)=DTIME KAMP(2)=10879.*AAU*SIN(INCL)/PERIOD/SQRT(1.-ECCEN**2)/(1.+M2M1) KAMP(1)=KAMP(2)*M2M1 810 PRINT *,'GIVE OFFSET (IN DAYS) TO NEXT POINT (- = DONE)' READ *, DEL IF (DEL.LT.0) GO TO 820 NEPH=NEPH+1 DELDAY(NEPH)=DEL+DTIME GO TO 810 820 WRITE (*,'('' WILL YOU WANT THIS SAVED FOR PLOTTING? (Y/N) ''\)') READ (*,'(A1)') AKEY IF (AKEY.EQ.'Y'.OR.AKEY.EQ.'y') THEN IFPLOT=1 WRITE (*,'('' NAME FOR PLOT FILE: ''\)') READ (*,'(A39)') FILNAM OPEN (UNIT=7,FILE=FILNAM,STATUS='UNKNOWN') ELSE IFPLOT=0 END IF STVSYS=VSYS(1) VSYS(1)=0. ID=1 PRINT *,' T-TP JUL DAY RHO THETA V1(G=0) V2-V1 X 1 Y' PRINT *,' DAYS -2440000. MAS DEGREES KM/SEC KM/SEC RSIN 1T RCOST' NPAUSE=0 DO 840 N=1,NEPH NPAUSE=NPAUSE+1 CALL KEPLER(DELDAY(N)+TPERIA) CALL SPECTR(VR,DVDEL,ID) CALL RHO(RHOC,COEFDR) CALL THETA(THC,COEFDR) RHOC=RHOC*1000. Y=RHOC*COS(THC) X=RHOC*SIN(THC) THC=THC*RADIAN DV=VR(2)-VR(1) PRINT 830,DELDAY(N),DELDAY(N)+TPERIA,RHOC,THC,VR(1),DV,X,Y 830 FORMAT (F10.3,F11.3,F7.1,F8.1,F9.2,F8.2,2F7.1) IF (IFPLOT.GT.0) THEN WRITE (7,'(F10.3,2F9.2,2F7.1)') DELDAY(N),VR(1),VR(2),X,Y END IF IF (NPAUSE.GT.21) THEN WRITE (*,'('' PAUSE: Hit any key ''\)') READ (*,'(A1)') AKEY NPAUSE=0 END IF 840 CONTINUE IF (IFPLOT.GT.0) CLOSE (UNIT=7) VSYS(1)=STVSYS GO TO 10 1000 CLOSE (UNIT=7) STOP END SUBROUTINE KEPLER(TIME) C THIS SOLVES KEPLERS EQUATION AND PROVIDES THE TRUE, ECCENTRIC, AND C MEAN ANOMALIES. C PERIOD IN DAYS, SEPARATIONS IN AU AND ANGLES IN RADIANS COMMON /ANOMAL/ TRUEA,ECCA,MEANA,CYCLES COMMON /ELEMNT/ PERIOD,TPERIA,AAU,PARALX,INCL,ECCEN,ARGPER,LGNODE, 1 M2M1,VSYS(6),KAMP(2) REAL INCL,LGNODE,M2M1,KAMP,MEANA COMMON /TRIG/ COSV,SINV,COSW,SINW,COSVW,SINVW,COSI,SINI CYCLES=(TIME-TPERIA)/PERIOD ICYC=CYCLES XMEANA=6.283185*(CYCLES-ICYC) MEANA=6.283185*CYCLES C NOW SOLVE KEPLERS EQUATION ECCA=XMEANA DO 10 I=1,50 SINE=SIN(ECCA) COSE=COS(ECCA) XNUM=ECCA-ECCEN*SINE-XMEANA ECCA=ECCA-XNUM/(1.-ECCEN*COSE) IF (ABS(XNUM).LT.1.E-5) GO TO 30 10 CONTINUE WRITE (6,20) XMEANA,ECCEN,ECCA,XNUM 20 FORMAT ('1KEPLER DIDNT CONVERGE',4F9.5) STOP 30 TRUEA=2.*ATAN(SQRT((1.+ECCEN)/(1.-ECCEN))*TAN(ECCA/2.)) COSW=COS(ARGPER) SINW=SIN(ARGPER) COSVW=COS(TRUEA+ARGPER) SINVW=SIN(TRUEA+ARGPER) SINV=SIN(TRUEA) COSV=COS(TRUEA) SINI=SIN(INCL) COSI=COS(INCL) RETURN END SUBROUTINE SPECTR(VR,DVDEL,ID) C THIS GENERATES THE SPECTROSCOPIC ORBIT AND DERIVATIVES FOR BOTH C PRIMARY AND SECONDARY C C IN ORDER THE DERIVATIVES, DVDEL, ARE C DV/DVSYS,DV/DKAMP,DV/DECCEN,DV/DARGPER,DV/DTPERIA, AND DV/DPERIOD C C NOTE WE TREAT THE AMPLITUDE AS AN INDEPENDENT VARIABLE C C FOR CONSISTENCY WITH THE VISUAL ORBIT, THE ARGUMENT OF PERIASTRON C IS THAT FOR THE SECONDARY. WE EXPLICITLY ALLOW FOR THE 180 DEGREE C DIFFERENCE IN THE CALCULATION FOR THE PRIMARY COMMON /ANOMAL/ TRUEA,ECCA,MEANA,CYCLES COMMON /ELEMNT/ PERIOD,TPERIA,AAU,PARALX,INCL,ECCEN,ARGPER,LGNODE, 1 M2M1,VSYS(6),KAMP(2) REAL INCL,LGNODE,M2M1,KAMP,MEANA COMMON /TRIG/ COSV,SINV,COSW,SINW,COSVW,SINVW,COSI,SINI DIMENSION VR(2),DVDEL(6,2),SIGN(2) DATA SIGN/-1.,1./ OME2=1.-ECCEN**2 OME232=SQRT(OME2)**3 TERM=ECCEN*COSW+COSVW TERM2=SINVW*(1.+ECCEN*COSV)**2 DO 10 I=1,2 SK=SIGN(I)*KAMP(I) VR(I)=VSYS(ID)+SK*TERM C DV/DVSYS DVDEL(1,I)=1. C DV/DKAMP DVDEL(2,I)=SIGN(I)*TERM C DV/DECCEN DVDEL(3,I)=SK*(COSW-SINVW*SINV/OME2*(2.+ECCEN*COSV)) C DV/DARGPER DVDEL(4,I)=-SK*(ECCEN*SINW+SINVW) C DV/DTPERIA DVDEL(5,I)=SK*TERM2*6.283185/PERIOD/OME232 C DV/DPERIOD DVDEL(6,I)=SK*TERM2*MEANA/PERIOD/OME232 10 CONTINUE RETURN END SUBROUTINE THETA(THC,DTHDEL) C THIS GENERATES THE ANGLE VARIABLE OF THE VISUAL ORBIT AS WELL AS C THE DERIVATIVES W.R.T. PERIOD(1),TPERIA(2),INCL(5),ECCEN(6), C ARGPER(7), AND LGNODE(8) COMMON /ANOMAL/ TRUEA,ECCA,MEANA,CYCLES COMMON /ELEMNT/ PERIOD,TPERIA,AAU,PARALX,INCL,ECCEN,ARGPER,LGNODE, 1 M2M1,VSYS(6),KAMP(2) REAL INCL,LGNODE,M2M1,KAMP,MEANA COMMON /TRIG/ COSV,SINV,COSW,SINW,COSVW,SINVW,COSI,SINI DIMENSION DTHDEL(15) DATA TWOPI/6.283185/,NDIM/15/ TH=ATAN2(SINVW*COSI,COSVW)+LGNODE THC=TH IF (TH.LT.0.) THC=TH+TWOPI IF (TH.GT.TWOPI) THC=TH-TWOPI DTHDV=COSI/(1.-(SINVW*SINI)**2) TERM1=1.+ECCEN*COSV TERM2=TERM1**2/SQRT((1.-ECCEN**2)**3) DO 10 I=1,NDIM 10 DTHDEL(I)=0. C DTH/D(PERIOD) DTHDEL(1)=-DTHDV*MEANA/PERIOD*TERM2 C DTH/D(TPERIA) DTHDEL(2)=-DTHDV*TWOPI/PERIOD*TERM2 C DTH/D(INCL) DTHDEL(5)=-DTHDV/COSI*SINI*COSVW*SINVW C DTH/D(ECCEN) DTHDEL(6)=DTHDV*SINV/(1.-ECCEN**2)*(1.+TERM1) C DTH/D(ARGPER) DTHDEL(7)=DTHDV C DTH/D(LGNODE) DTHDEL(8)=1. RETURN END SUBROUTINE RHO(RHOC,DRDEL) C THIS GENERATES THE SEPARATIONS FOR THE VISUAL ORBIT AS WELL AS THE C DERIVATIVES W.R.T. PERIOD(1),TPERIA(2),AAU(3),PI(4),INCL(5), C ECCEN(6), AND ARGPER(7) COMMON /ANOMAL/ TRUEA,ECCA,MEANA,CYCLES COMMON /ELEMNT/ PERIOD,TPERIA,AAU,PARALX,INCL,ECCEN,ARGPER,LGNODE, 1 M2M1,VSYS(6),KAMP(2) REAL INCL,LGNODE,M2M1,KAMP,MEANA COMMON /TRIG/ COSV,SINV,COSW,SINW,COSVW,SINVW,COSI,SINI DIMENSION DRDEL(15) DATA NDIM/15/,TWOPI/6.283185/ DO 10 I=1,NDIM 10 DRDEL(I)=0. AP=AAU*PARALX DRDA=(1.-ECCEN**2)/(1.+ECCEN*COSV) RAD=AP*DRDA DRDE=-AP*COSV DRDP=-AP*ECCEN*MEANA/PERIOD/SQRT(1.-ECCEN**2)*SINV DRDT=-AP*ECCEN*TWOPI/PERIOD/SQRT(1.-ECCEN**2)*SINV DRHODR=SQRT(1.-(SINVW*SINI)**2) DRHODV=-RAD*SINVW*COSVW*SINI**2/DRHODR DVDP=-MEANA/PERIOD*(1.+ECCEN*COSV)**2/SQRT((1.-ECCEN**2)**3) DVDT=-TWOPI/PERIOD*(1.+ECCEN*COSV)**2/SQRT((1.-ECCEN**2)**3) RHOC=RAD*DRHODR C DRHO/D(PERIOD) DRDEL(1)=DRHODV*DVDP+DRHODR*DRDP C DRHO/D(TPERIA) DRDEL(2)=DRHODV*DVDT+DRHODR*DRDT C DRHO/D(AAU) DRHODA=DRHODR*DRDA DRDEL(3)=PARALX*DRHODA C DRHO/D(PI) DRDEL(4)=AAU*DRHODA C DRHO/D(INCL) DRDEL(5)=DRHODV*SINVW*COSI/COSVW/SINI C DRHO/D(ECCEN) DRDEL(6)=DRHODV*SINV/(1.-ECCEN**2)*(2.+ECCEN*COSV)+DRHODR*DRDE C DRHO/D(ARGPER) DRDEL(7)=DRHODV RETURN END SUBROUTINE LUNOCC(ASPECT,RHOAC,RHOSIN,COEFDR) C THIS CALCULATES THE EXPECTED PROJECTION OF THE SEPARATION C ALONG THE ASPECT DIRECTION. C C THE PARTIAL DERIVATIVES WITH RESPECT TO: PERIOD(1),TPERIA(2), C AAU(3),PI(4),INCL(5),ECCEN(6),ARGPER(7) AND LGNODE(8) ARE C ALSO OBTAINED. THIS IS ACCOMPLISHED BY USING THE RHO AND C THETA SUBROUTINES. C C ASSUMES KEPLER HAS ALREADY BEEN CALLED DIMENSION DTHDEL(15),DRDEL(15),COEFDR(15) DATA NDIM/15/ CALL THETA(THC,DTHDEL) CALL RHO(RHOC,DRDEL) COSATH=COS(ASPECT-THC) SINATH=SIN(ASPECT-THC) RHOSIN=RHOC*SINATH RHOAC=RHOC*COSATH DO 10 N=1,NDIM 10 COEFDR(N)=DRDEL(N)*COSATH+DTHDEL(N)*RHOSIN RETURN END C Error Propagation Package FUNCTION CCRLAT(IDX1,PFPT1,N1,IDX2,PFPT2,N2,CORMTX,ICNVRG,NORD) C THIS EVALUATES THE CORRELATION BETWEEN TWO QUANTITIES DERIVED FROM C VALUES DERIVED IN THE LEAST SQUARES ADJUSTMENT. C C IF INDX1=INDX2 AND PFPT1=PFPT2, THEN YOU GET SIGMA**2 C C OTHERWISE, YOU GET SIGMA(QTY1,QTY2) IN DIMENSIONLESS UNITS C C TO GET THE CORRELATION, SIMPLY DIVIDE THE RESULT BY C SIGMA(QTY1)/QTY1 AND SIGMA(QTY2)/QTY2 C C ASSUMES THAT THE PFPTN'S ARE PARTIALS OF LN(F(THETA)) W.R.T. THETA C TIMES SIGMA(THETA) DIMENSION IDX1(1),IDX2(1),PFPT1(1),PFPT2(1),CORMTX(1),ICNVRG(1) CCRLAT=0. IF (N1.LE.0) RETURN IF (N2.LE.0) RETURN DO 10 I=1,N1 IX1=IDX1(I) IC1=ICNVRG(IX1) DO 10 J=1,N2 IX2=IDX2(J) IC2=ICNVRG(IX2) IJC=(IC2-1)*NORD+IC1 10 CCRLAT=CCRLAT+PFPT1(I)*PFPT2(J)*CORMTX(IJC) RETURN END FUNCTION PRPGAT(INDX,PFPTHT,NUMBER,CORMTX,ICNVRG,NORD) C PROPAGATES ERRORS OF MULTIPLICATIVE QUANTITIES C ASSUMES THE PARTIAL OF LN(F) W.R.T. THETA HAS SIGMA(THETA) IN IT DIMENSION CORMTX(1),INDX(1),PFPTHT(1),ICNVRG(1) PRPGAT=0. DO 20 I=1,NUMBER IX=INDX(I) IC=ICNVRG(IX) PRPGAT=PRPGAT+PFPTHT(I)**2 IF (I.EQ.NUMBER) GO TO 20 I1=I+1 DO 10 J=I1,NUMBER JX=INDX(J) JC=ICNVRG(JX) IJC=(JC-1)*NORD+IC 10 PRPGAT=PRPGAT+2.*PFPTHT(I)*PFPTHT(J)*CORMTX(IJC) 20 CONTINUE PRPGAT=SQRT(PRPGAT) RETURN END C Non-linear Least Squares Package SUBROUTINE ACUMLS (PARTAL,RESID,SIG,MODE) C THIS SETS UP THE NORMAL EQUATIONS COMMON /EQNORM/ XNORM(1444),DRIVE(38),CHISQ,NDIM,MAPCVG(38) DIMENSION PARTAL(38) GO TO (10,30),MODE C CLEAR ACCUMULATORS 10 CHISQ=0. DO 20 I=1,NDIM DRIVE(I)=0. DO 20 K=1,NDIM J=I+NDIM*(K-1) XNORM(J)=0. 20 CONTINUE RETURN C ACCUMULATE 30 WT=1./SIG**2 CHISQ=CHISQ+RESID**2*WT DO 40 I=1,NDIM IM=MAPCVG(I) DRIVE(I)=DRIVE(I)+WT*RESID*PARTAL(IM) DO 40 K=1,NDIM KM=MAPCVG(K) J=I+NDIM*(K-1) XNORM(J)=XNORM(J)+WT*PARTAL(IM)*PARTAL(KM) 40 CONTINUE RETURN END SUBROUTINE SOLLSQ(COREC) C SOLVES THE LEAST SQUARES MATRIX FOR CORRECTIONS COMMON /EQNORM/ XNORM(1444),DRIVE(38),CHISQ,NDIM,MAPCVG(38) DIMENSION COREC(38) COMMON /DUMMY/ DMTX(1444),IDUM(38) DO 10 I=1,NDIM COREC(I)=DRIVE(I) DO 10 J=1,NDIM K=I+NDIM*(J-1) 10 DMTX(K)=XNORM(K) CALL SOLVIT(DMTX,NDIM,COREC,IDUM) RETURN END SUBROUTINE INVRT(SIGMA,CORMTX) C INVERTS NORMAL EQUATIONS, DERIVES SIGMAS, AND NORMALIZES CORMTX COMMON /EQNORM/ XNORM(1444),DRIVE(38),CHISQ,NDIM,MAPCVG(38) COMMON /DUMMY/ DMTX(1444),IDUM(38) DIMENSION SIGMA(38),CORMTX(1444) DO 20 I=1,NDIM L=1+NDIM*(I-1) DO 10 J=1,NDIM M=J+NDIM*(I-1) CORMTX(M)=0. DO 10 K=1,NDIM N=J+NDIM*(K-1) DMTX(N)=XNORM(N) 10 CONTINUE N=I+NDIM*(I-1) CORMTX(N)=1. CALL SOLVIT(DMTX,NDIM,CORMTX(L),IDUM) 20 CONTINUE DO 30 I=1,NDIM J=I+NDIM*(I-1) SIGMA(I)=SQRT(CORMTX(J)) 30 CONTINUE DO 40 I=1,NDIM DO 40 J=1,NDIM K=I+NDIM*(J-1) CORMTX(K)=CORMTX(K)/SIGMA(I)/SIGMA(J) 40 CONTINUE RETURN END SUBROUTINE SOLVIT(A,N,B,IPIVOT) C SOLVES LINEAR EQUATIONS C A IS A COMPLETELY FILLED N BY N ARRAY WHICH IS DESTROYED. C B IS THE RIGHT SIDE VECTOR OF LENGTH N AND RETURNS AS THE SOLUTION C IPIVOT IS A SCRATCH AREA OF LENGTH N. DIMENSION A(300),B(30),IPIVOT(30) EQUIVALENCE(AMAX,SWAP,PIVOT,T) DO 20 J=1,N 20 IPIVOT(J)=0 DO 550 I=1,N AMAX=0. DO 105 J=1,N IF(IPIVOT(J).EQ.1)GO TO 105 JK=J-N DO 100 K=1,N JK=JK+N IF(IPIVOT(K).EQ.1)GO TO 100 AA=ABS(A(JK)) IF(AMAX.GE.AA)GO TO 100 IROW=J ICOLUM=K AMAX=AA 100 CONTINUE 105 CONTINUE IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1 IF(IROW.EQ.ICOLUM)GO TO 260 IRL=IROW-N ICL=ICOLUM-N DO 200 L=1,N IRL=IRL+N SWAP=A(IRL) ICL=ICL+N A(IRL)=A(ICL) 200 A(ICL)=SWAP SWAP=B(IROW) B(IROW)=B(ICOLUM) B(ICOLUM)=SWAP 260 ICIC=ICOLUM*N+ICOLUM-N PIVOT=A(ICIC) A(ICIC)=1. ICL=ICOLUM-N DO 350 L=1,N ICL=ICL+N 350 A(ICL)=A(ICL)/PIVOT B(ICOLUM)=B(ICOLUM)/PIVOT L1IC=ICOLUM*N-N DO 550 L1=1,N L1IC=L1IC+1 IF(L1.EQ.ICOLUM)GO TO 550 T=A(L1IC) A(L1IC)=0. L1L=L1-N ICL=ICOLUM-N DO 450 L=1,N L1L=L1L+N ICL=ICL+N 450 A(L1L)=A(L1L)-A(ICL)*T B(L1)=B(L1)-B(ICOLUM)*T 550 CONTINUE RETURN END