C @(#)compxy.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:56 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM COMPXY C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C original FORTRAN code: D. Gillet ESO - Garching C modified and adapted to MIDAS environment: D. Baade ST-ECF, Garching C program BARYCORR version 1.00 860203 C original FORTRAN code: D. Gillet ESO - Garching C modified and adapted to MIDAS environment: D. Baade ST-ECF, Garching C program PRECESS version 1.0 860201 C original FORTRAN code: D. Gillet ESO - Garching C modified and adapted to MIDAS environment: D. Baade ST-ECF, Garching C C this is the merged version of old COMPST.FOR, COMPUT.FOR, BARYCORR.FOR C K.Banse ESO - IPG C program COMPXY version 2.00 890519 C 2.10 900213 2.15 900711 2.20 901008 C 2.25 910429 2.3 911127 C 3.35 920602 C 3.50 920902 PJG 3.60 931112 C C.COMMENT C The program was changed so that to input longitude now conform C to IAU standard namely east-longitude is positive!!! C Internally east-longitudes are converted back to west! PJG 920902 C C.KEYWORDS C a,b) Julian date, conversion from universal time to local sidereal time C c) barycentric and heliocentric corrections of times and radial velocities C d) precession of equatorial coordinates C C.PURPOSE C a) calculate Julian date, convert from UT to local ST C b) calculate Julian date, convert from local ST to UT C c) correct times for the light travel time in the solar system, i.e. C reduce to either barycentric or heliocentric times C and C correct radial velocities for the earth's orbital and diurnal velocities, C i.e. reduce to either barycentric or heliocentric radial velocities C d) precess equatorial coordinates from one epoch to another C C.ALGORITHM C C Reference : the Astronomical Almanac 1983, p B7 C in 1983: 1 mean solar day = 1.00273790931 mean sidereal days C ST without equation of equinoxes correction => accuracy +/- 1 sec C C.INPUT/OUTPUT C the following keywords are used: C ACTION/C/1/2 a,b) US for UT -> ST, SU for ST -> UT C c) BA for barycorr stuff C d) PR for precession stuff C C for COMPUTE/UT /ST C INPUTR/R/1/3 date: year,month,day C INPUTR/R/4/3 universal or sidereal time: hour,min,sec C INPUTR/R/7/3 EAST longitude of observatory: degree,min,sec !! NOTE C OUTPUTR/R/1/3 sidereal or universal time : hour,min,sec C OUTPUTD/D/1/1 Julian date C C for COMPUTE/BARYCORR C INPUTR/R/1/3 date: year,month,day C INPUTR/R/4/3 universal time: hour,min,sec C INPUTR/R/7/3 EAST longitude of observatory: degree,min,sec !! NOTE C INPUTR/R/10/3 latitude of observatory: degree,min,sec C INPUTR/R/13/3 right ascension: hour,min,sec C INPUTR/R/16/3 declination: degree,min,sec C OUTPUTD/D/1/1 barycentric correction to time (days) C OUTPUTD/D/2/1 heliocentric correction to time (days) C OUTPUTR/R/1/1 barycentric correction to radial velocity (km/s) C OUTPUTR/R/2/1 heliocentric correction to radial velocity (km/s) C C for COMPUTE/PRECESSION C INPUTR/R/1/3 original epoch: year,month,day C INPUTR/R/4/3 new epoch: year,month,day C INPUTR/R/7/3 original right ascension: hour,min,sec C INPUTR/R/10/3 original declination: degree,min,sec C OUTPUTR/R/1/3 precessed right ascension: hour,min,sec C OUTPUTR/R/4/3 precessed declination: degree,min,sec C OUTPUTR/R/7/1 new epoch: real years C OUTPUTD/D/1/2 old epoch, new epoch C C------------------------------------------------------------------------------ C IMPLICIT NONE C INTEGER IAV,STAT,KUN(1),KNUL,N INTEGER MADRID C DOUBLE PRECISION UTR,STR,T0,DL,THETA0,PE,ST0HG,STG,GAST,R1 DOUBLE PRECISION JD,JD0H,JD00,ZERO DOUBLE PRECISION DCORB(3),DCORH(3),DVELB(3),DVELH(3) DOUBLE PRECISION ALP,BCT,BEOV,BERV,DEL,EDV DOUBLE PRECISION HAR,HCT,HEOV,HERV,PHI,PI DOUBLE PRECISION EQX0,EQX1 DOUBLE PRECISION A0R,A1R,D0R,D1R DOUBLE PRECISION DSMALL,DTEMP(3) C REAL DATE0(3),DATE1(3),DATE00(3),A0(3),A1(3),D0(3),D1(3) REAL DATE(3),UT(3),OLONG(3),ST(3) REAL OLAT(3),ALPHA(3),DELTA(3) REAL RBUF(20) C CHARACTER ACTIO*2,SIGNS*3,INPSGN*3 C COMMON /VMR/MADRID(1) C EQUIVALENCE (RBUF(1),DATE(1)),(RBUF(7),OLONG(1)) EQUIVALENCE (RBUF(10),OLAT(1)),(RBUF(13),ALPHA(1)) EQUIVALENCE (RBUF(16),DELTA(1)) C DATA PI /3.1415926535897928D0/ DATA DSMALL /1.D-38/ C C ... set up MIDAS environment C CALL STSPRO('COMPXY') ZERO = 0.D0 DO 100, N=1,3 DATE0(N) = 0.0 DATE1(N) = 0.0 DATE00(N) = 0.0 A0(N) = 0.0 A1(N) = 0.0 D0(N) = 0.0 D1(N) = 0.0 DATE(N) = 0.0 UT(N) = 0.0 ST(N) = 0.0 OLONG(N) = 0.0 OLAT(N) = 0.0 ALPHA(N) = 0.0 DELTA(N) = 0.0 100 CONTINUE C C ... read keywords for input C CALL STKRDC('ACTION',1,1,2,IAV,ACTIO,KUN,KNUL,STAT) CALL UPCAS(ACTIO,ACTIO) CALL STKRDR('INPUTR',1,20,IAV,RBUF,KUN,KNUL,STAT) CALL STKRDC('INPUTC',1,1,3,IAV,INPSGN,KUN,KNUL,STAT) C C compute/barycorr is done in section 1000, compute/precession in section 2000 C IF (ACTIO.EQ.'BA') GOTO 1000 IF (ACTIO.EQ.'PR') GOTO 2000 C C ***************************************************** C C here for COMPUTE/UT or COMPUTE/ST C C ***************************************************** C IF (ACTIO.EQ.'US') THEN DO 200, N=1,3 UT(N) = RBUF(N+3) 200 CONTINUE UTR = UT(1)+UT(2)/60.D0+UT(3)/3600.D0 !convert UT to real hours ELSE C DO 300, N=1,3 ST(N) = RBUF(N+3) 300 CONTINUE STR = ST(1)+ST(2)/60.D0+ST(3)/3600.D0 !convert ST to real hours ENDIF C C ... likewise convert longitude of observatory to real hours, C take care of sign. NOTE: east longitude is assumed for input !! C SIGNS = '+++' IF ((OLONG(1).LT.0.0) .OR. (OLONG(2).LT.0.0) .OR. + (OLONG(3).LT.0.0) .OR. (INPSGN(1:1).EQ.'-')) THEN SIGNS(1:1) = '-' OLONG(1) = ABS(OLONG(1)) OLONG(2) = ABS(OLONG(2)) OLONG(3) = ABS(OLONG(3)) ENDIF DTEMP(1) = DBLE(OLONG(1)) DTEMP(2) = DBLE(OLONG(2)) DTEMP(3) = DBLE(OLONG(3)) DL = DTEMP(1)+DTEMP(2)/60.D0+DTEMP(3)/3600.D0 IF (SIGNS(1:1).EQ.'-') DL = -DL ! negative longitude DL = -DL*24.D0/360.D0 !convert back to west longitude C C ... compute ST at 0 hours UT (ST0HG) C ... REFERENCE : MEEUS J.,1980,ASTRONOMICAL FORMULAE FOR CALCULATORS C CALL JULDAT(DATE,ZERO,JD0H) DTEMP(1) = 2415020.D0 T0=(JD0H-DTEMP(1))/36525.D0 DTEMP(1) = 0.276919398D0 DTEMP(2) = 100.0021359D0 DTEMP(3) = 0.000001075D0 THETA0 = DTEMP(1) + DTEMP(2)*T0 + DTEMP(3)*T0*T0 PE=DINT(THETA0) THETA0=THETA0-PE ST0HG=THETA0*24.D0 C IF (ACTIO.EQ.'US') THEN C C now do the conversion UT -> ST (MEAN SIDEREAL TIME) STG=ST0HG + (UTR*1.00273790931D0) IF (STG.LT.DL) STG=STG+24.D0 STR=STG-DL 400 IF (STR.GE.24.D0) THEN STR=STR-24.D0 GOTO 400 ENDIF ST(1)=IDINT(STR) R1=(STR-ST(1))*60.D0 ST(2)=IDINT(R1) ST(3)=(R1-ST(2))*60.D0 ELSE C C now do the conversion ST -> UT (MEAN SIDEREAL TIME) GAST=STR+DL IF(GAST.LT.ST0HG)GAST=GAST+24.D0 UTR=(GAST-ST0HG)/1.00273790931D0 IF (UTR.GE.24.D0) UTR=UTR-24.D0 UT(1)=IDINT(UTR) R1=(UTR-UT(1))*60.D0 UT(2)=IDINT(R1) UT(3)=(R1-UT(2))*60.D0 ENDIF C C ... compute Julian date C CALL JULDAT(DATE,UTR,JD) C C ... write output to keyword C CALL STKWRD('OUTPUTD',JD,1,1,KUN,STAT) IF (ACTIO.EQ.'US') THEN CALL STKWRR('OUTPUTR',ST,1,3,KUN,STAT) ELSE CALL STKWRR('OUTPUTR',UT,1,3,KUN,STAT) ENDIF GOTO 9000 C C ***************************************************** C C here we do COMPUTE/BARYCORR C C ***************************************************** C 1000 SIGNS = '+++' DO 1100 N=1,3 UT(N) = RBUF(N+3) 1100 CONTINUE C C ... convert UT to real hours, calculate Julian date C UTR = UT(1)+UT(2)/60.D0+UT(3)/3600.D0 CALL JULDAT(DATE,UTR,JD) C C ... likewise convert longitude and latitude of observatory to real hours C ... and degrees, respectively; take care of signs C ... NOTE: east longitude is assumed for input !! C IF ((OLONG(1).LT.0.0) .OR. (OLONG(2).LT.0.0) .OR. + (OLONG(3).LT.0.0) .OR. (INPSGN(1:1).EQ.'-')) THEN SIGNS(1:1) = '-' OLONG(1) = ABS(OLONG(1)) OLONG(2) = ABS(OLONG(2)) OLONG(3) = ABS(OLONG(3)) ENDIF DL = OLONG(1)+OLONG(2)/60.D0+OLONG(3)/3600.D0 IF (SIGNS(1:1).EQ.'-') DL = -DL ! negative longitude DL = -DL*24.D0/360.D0 ! convert back to west longitude C IF ((OLAT(1).LT.0.0) .OR. (OLAT(2).LT.0.0) .OR. + (OLAT(3).LT.0.0) .OR. (INPSGN(2:2).EQ.'-')) THEN SIGNS(2:2) = '-' OLAT(1) = ABS(OLAT(1)) OLAT(2) = ABS(OLAT(2)) OLAT(3) = ABS(OLAT(3)) ENDIF PHI = OLAT(1)+OLAT(2)/60.D0+OLAT(3)/3600.D0 IF (SIGNS(2:2).EQ.'-') PHI = -PHI ! negative latitude PHI = PHI*PI/180.D0 C C ... convert right ascension and declination to real radians C ALP = (ALPHA(1)*3600D0+ALPHA(2)*60D0+ALPHA(3))*PI/(12.D0*3600.D0) C IF ((DELTA(1).LT.0.0) .OR. (DELTA(2).LT.0.0) .OR. + (DELTA(3).LT.0.0) .OR. (INPSGN(3:3).EQ.'-')) THEN SIGNS(3:3) = '-' DELTA(1) = ABS(DELTA(1)) DELTA(2) = ABS(DELTA(2)) DELTA(3) = ABS(DELTA(3)) ENDIF DEL = (DELTA(1)*3600.D0 + DELTA(2)*60.D0 + DELTA(3)) + * PI/(3600.D0*180.D0) IF (SIGNS(3:3).EQ.'-') DEL = -DEL ! negative declination C C ... calculate earth's orbital velocity in rectangular coordinates X,Y,Z C ... for both heliocentric and barycentric frames (DVELH, DVELB) C ... Note that setting the second argument of BARVEL to zero as done below C ... means that the input coordinates will not be corrected for precession. C CALL BARVEL(JD,0.0D0,DVELH,DVELB) C C ... with the rectangular velocity components known, the respective projections C ... HEOV and BEOV on a given line of sight (ALP,DEL) can be determined: C C ... REFERENCE: THE ASTRONOMICAL ALMANAC 1982 PAGE:B17 C BEOV=DVELB(1)*DCOS(ALP)*DCOS(DEL)+ 1 DVELB(2)*DSIN(ALP)*DCOS(DEL)+ 2 DVELB(3)*DSIN(DEL) HEOV=DVELH(1)*DCOS(ALP)*DCOS(DEL)+ 1 DVELH(2)*DSIN(ALP)*DCOS(DEL)+ 2 DVELH(3)*DSIN(DEL) C C ... For determination also of the contribution due to the diurnal rotation of C ... the earth (EDV), the hour angle (HAR) is needed at which the observation C ... was made which requires conversion of UT to sidereal time (ST). C C ... Therefore, first compute ST at 0 hours UT (ST0HG) C C ... REFERENCE : MEEUS J.,1980,ASTRONOMICAL FORMULAE FOR CALCULATORS C CALL JULDAT(DATE,ZERO,JD0H) T0=(JD0H-2415020.D0)/36525.D0 THETA0=0.276919398D0+100.0021359D0*T0+0.000001075D0*T0*T0 PE=DINT(THETA0) THETA0=THETA0-PE ST0HG=THETA0*24.D0 C C ... now do the conversion UT -> ST (MEAN SIDEREAL TIME) C C ... REFERENCE : THE ASTRONOMICAL ALMANAC 1983, P B7 C ... IN 1983: 1 MEAN SOLAR DAY = 1.00273790931 MEAN SIDEREAL DAYS C ... ST WITHOUT EQUATION OF EQUINOXES CORRECTION => ACCURACY +/- 1 SEC C STG=ST0HG+UTR*1.00273790931D0 IF (STG.LT.DL) STG=STG+24.D0 STR=STG-DL IF (STR.GE.24.D0) STR=STR-24.D0 STR = STR*PI/12.D0 ! ST in radians C HAR=STR-ALP ! hour angle of observation EDV=-0.4654D0*DSIN(HAR)*DCOS(DEL)*DCOS(PHI) C C ... the total correction (in km/s) is the sum of orbital and diurnal components C HERV=HEOV+EDV BERV=BEOV+EDV C C ... Calculation of the barycentric and heliocentric correction times C ... (BCT and HCT) requires knowledge of the earth's position in its C ... orbit. Subroutine BARCOR returns the rectangular barycentric (DCORB) C ... and heliocentric (DCORH) coordinates. C CALL BARCOR(DCORH,DCORB) C C ... from this, the correction times (in days) can be determined: C ... (REFERENCE: THE ASTRONOMICAL ALMANAC 1982 PAGE:B16) C BCT=+0.0057756D0*(DCORB(1)*DCOS(ALP)*DCOS(DEL)+ 1 DCORB(2)*DSIN(ALP)*DCOS(DEL)+ 2 DCORB(3)* DSIN(DEL)) HCT=+0.0057756D0*(DCORH(1)*DCOS(ALP)*DCOS(DEL)+ 1 DCORH(2)*DSIN(ALP)*DCOS(DEL)+ 2 DCORH(3)* DSIN(DEL)) C C ... write results to keywords C CALL STKWRD('OUTPUTD',BCT,1,1,KUN,STAT) ! barycentric correction time CALL STKWRD('OUTPUTD',HCT,2,1,KUN,STAT) ! heliocentric correction time RBUF(1) = BERV ! barocentric RV correction RBUF(2) = HERV ! heliocentric RV correction C ... (note that EDV is already contained in both BERV and HERV) RBUF(3) = EDV ! diurnal RV correction CALL STKWRR('OUTPUTR',RBUF,1,3,KUN,STAT) GOTO 9000 C C ***************************************************** C C here we handle COMPUTE/PRECESSION C C ***************************************************** C 2000 SIGNS = '++' DO 2100 N=1,3 DATE0(N) = RBUF(N) DATE1(N) = RBUF(N+3) A0(N) = RBUF(N+6) D0(N) = RBUF(N+9) 2100 CONTINUE C C ... convert equinoxes to real years C CALL JULDAT(DATE0,ZERO,JD) DATE00(1)=DATE0(1) CALL JULDAT(DATE00,ZERO,JD00) EQX0=DATE0(1)+(JD-JD00)/365.25D0 C CALL JULDAT(DATE1,ZERO,JD) DATE00(1)=DATE1(1) CALL JULDAT(DATE00,ZERO,JD00) EQX1=DATE1(1)+(JD-JD00)/365.25D0 C C ... check if equinoxes are in fact different C CALL STKWRD('OUTPUTD',EQX0,1,1,KUN,STAT) CALL STKWRD('OUTPUTD',EQX1,2,1,KUN,STAT) IF (ABS(EQX0-EQX1).LT.DSMALL) THEN CALL STTPUT + ('Equinoxes identical - no precession applied',STAT) CALL STKWRR('OUTPUTR',RBUF(7),1,7,KUN,STAT) GOTO 9000 ENDIF C C ... convert right ascension and declination to real radians C A0R = (A0(1)*3600D0+A0(2)*60D0+A0(3))*PI/(12.*3600.) IF ((D0(1).LT.0.0) .OR. + (D0(2).LT.0.0) .OR. + (D0(3).LT.0.0) .OR. + (INPSGN(1:1).EQ.'-')) THEN !take care of -0 ...! SIGNS(1:1) = '-' D0(1) = ABS(D0(1)) D0(2) = ABS(D0(2)) D0(3) = ABS(D0(3)) ENDIF D0R = (D0(1)*3600.D0+D0(2)*60.D0+D0(3)) + * PI/(3600.0D0*180.0D0) IF (SIGNS(1:1).EQ.'-') D0R = -D0R C C ... subroutine PRE computes the precession C CALL PRE(A0R,D0R,A1R,D1R,EQX0,EQX1) C C ... convert precessed right ascension to hours, minutes, seconds C A1R = A1R*12/PI ! convert R.A. to real hours A1(1) = DINT(A1R) R1 = (A1R-A1(1))*60. A1(2) = DINT(R1) A1(3) = (R1-A1(2))*60 C C ... convert precessed declination to degrees, minutes, seconds C ... take care of sign C IF (D1R.LT.0.D0) THEN SIGNS(1:1) = '-' D1R = -D1R ELSE SIGNS(1:1) = '+' ENDIF D1R = D1R*180./PI ! convert decl. to real degrees D1(1) = DINT(D1R) R1 = (D1R-D1(1))*60. D1(2) = DINT(R1) D1(3) = (R1-D1(2))*60 IF (SIGNS(1:1).EQ.'-') THEN ! negative declination IF (D1(1).GT.DSMALL) THEN D1(1) = -D1(1) ELSE IF (D1(2).GT.DSMALL) THEN D1(2) = -D1(2) ELSE D1(3) = -D1(3) ENDIF ENDIF C C ... write results to keyword C CALL STKWRC('SIGN',1,SIGNS(1:1),1,1,KUN,STAT) !communicate sign CALL STKWRR('OUTPUTR',A1,1,3,KUN,STAT) CALL STKWRR('OUTPUTR',D1,4,3,KUN,STAT) C C ... finish off C 9000 CALL STSEPI END SUBROUTINE PRE(A0,D0,A1,D1,EQX0,EQX1) C C ... calculate general precession for two given epochs C C ... References: STUMPFF P.,ASTRON. ASTROPHYS. SUPPL. 41, P. 1 (1980) C ... EXPLAN. SUPPL. ASTRON. EPHEREMERIS, P. 28 (1961) C IMPLICIT NONE !REAL*8 (A-H,O-Z) C DOUBLE PRECISION A0,D0,D1,EQX0,EQX1 C DOUBLE PRECISION DT0,DT,DTS,DTC,DC1900,DC1M2 DOUBLE PRECISION DT0,DT,DTS,DTC,DC2000,DC1M2 DOUBLE PRECISION DC1,DC2,DC3,DC4,DC5,DC6,DC7,DC8,DC9 DOUBLE PRECISION DZED,DZETA,DCSAR,DTHETA DOUBLE PRECISION PI,SINR,COSR,A1,R C C DATA DCSAR/4.848136812D-6/,DC1900/1900.0D0/, C * DC1M2/0.01D0/, C * DC1/2304.25D0/,DC2/1.396D0/,DC3/0.302D0/,DC4/0.018D0/, C * DC5/0.791D0/,DC6/2004.683D0/,DC7/-0.853D0/,DC8/-0.426D0/, C * DC9/-0.042D0/ DATA DCSAR/4.84813681109535994D-6/, * DC2000/2000.0D0/, * DC1M2/0.01D0/, * DC1/2306.2181D0/, * DC2/1.39656D0/, * DC3/0.30188D0/, * DC4/0.018203D0/, * DC5/0.791D0/, * DC6/2004.3109D0/, * DC7/-0.85330D0/, * DC8/-0.42665D0/, * DC9/-0.041833D0/ C PI = 3.1415926535897928D0 C C ... first determine the precession angles zeta, z and theta C C DT0=(EQX0-DC1900)*DC1M2 DT0=(EQX0-DC2000)*DC1M2 DT=(EQX1-EQX0)*DC1M2 DTS=DT*DT DTC=DTS*DT DZETA=((DC1+DC2*DT0)*DT+DC3*DTS+DC4*DTC)*DCSAR DZED=DZETA + DC5*DTS*DCSAR DTHETA=((DC6+DC7*DT0)*DT+DC8*DTS+DC9*DTC)*DCSAR C C ... with these and A0,D0 compute the precessed coordinates A1,D1 C D1 = ASIN(DCOS(DTHETA)*DSIN(D0)+DSIN(DTHETA)*DCOS(D0)* 1 DCOS(A0+DZETA)) SINR = (DCOS(D0)*DSIN(A0+DZETA))/DCOS(D1) COSR = (DCOS(DTHETA)*DCOS(D0)*DCOS(A0+DZETA)- 1 DSIN(DTHETA)*DSIN(D0))/DCOS(D1) R = DASIN(SINR) IF (COSR.LT.0) R=PI-R ! determine quadrant of R A1 = R+DZED ! this is the precessed R.A. IF (A1.GT.2.0D0*PI) A1=A1-2.0D0*PI IF (A1.LT.0.0D0) A1=A1+2.0D0*PI C RETURN END SUBROUTINE BARVEL(DJE,DEQ,DVELH,DVELB) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C PURPOSE : compute rectangular heliocentric and barycentric components of C the earth's orbital velocity C C REFERENCE : STUMPFF P. ASTRON. ASTOPHYS. SUPPL. 41,1,1980 C MODIFICATION : D. GILLET 1983-9-15 C C OTHER ROUTINES CALLED : PRE (calculates matrix of general precession) C C------------------------------------------------------------------------------- C IMPLICIT NONE !REAL*8 (A-H,O-Z) C DOUBLE PRECISION DJE,DEQ,DVELH(3),DVELB(3),SN(4) DOUBLE PRECISION DT,DTL,DCT0,DCJUL,DTSQ,DLOCAL,DC2PI,CC2PI DOUBLE PRECISION DRD,DRLD,DCSLD,DC1 DOUBLE PRECISION DXBD,DYBD,DZBD,DZHD,DXHD,DYHD DOUBLE PRECISION DYAHD,DZAHD,DYABD,DZABD DOUBLE PRECISION DML,DEPS,PHI,PHID,PSID,DPARAM,PARAM DOUBLE PRECISION CCFDI,CCKM,CCMLD,PLON,POMG,PECC DOUBLE PRECISION PERTL,PERTLD,PERTRD,PERTP,PERTR,PERTPD DOUBLE PRECISION SINA,CCSGD,DC1MME,TL DOUBLE PRECISION CCSEC3,COSA,ESQ DOUBLE PRECISION DCFEL(3,8),DCEPS(3),CCSEL(3,17),DCARGS(2,15) DOUBLE PRECISION CCAMPS(5,15),CCSEC(3,4),DCARGM(2,3) DOUBLE PRECISION CCAMPM(4,3),CCPAMV(4) DOUBLE PRECISION A,B,E,F,G,SINF,COSF,T,TSQ,TWOE,TWOG C DOUBLE PRECISION DPREMA(3,3),DPSI,D1PDRO,DSINLS DOUBLE PRECISION DCOSLS,DSINEP,DCOSEP DOUBLE PRECISION FORBEL(7),SORBEL(17),SINLP(4),COSLP(4) DOUBLE PRECISION SINLM,COSLM,SIGMA C INTEGER IDEQ,K,N C COMMON /BARXYZ/ DPREMA,DPSI,D1PDRO,DSINLS,DCOSLS, + DSINEP,DCOSEP,FORBEL,SORBEL,SINLP, + COSLP,SINLM,COSLM,SIGMA,IDEQ C EQUIVALENCE (SORBEL(1),E),(FORBEL(1),G) C DATA DC2PI/6.2831853071796D0/,CC2PI/6.283185/, *DC1/1.0D0/,DCT0/2415020.0D0/,DCJUL/36525.0D0/ C DATA DCFEL/ 1.7400353D+00, 6.2833195099091D+02, 5.2796D-06, * 6.2565836D+00, 6.2830194572674D+02,-2.6180D-06, * 4.7199666D+00, 8.3997091449254D+03,-1.9780D-05, * 1.9636505D-01, 8.4334662911720D+03,-5.6044D-05, * 4.1547339D+00, 5.2993466764997D+01, 5.8845D-06, * 4.6524223D+00, 2.1354275911213D+01, 5.6797D-06, * 4.2620486D+00, 7.5025342197656D+00, 5.5317D-06, * 1.4740694D+00, 3.8377331909193D+00, 5.6093D-06/ C DATA DCEPS/ 4.093198D-01,-2.271110D-04,-2.860401D-08/ C DATA CCSEL/ 1.675104D-02,-4.179579D-05,-1.260516D-07, * 2.220221D-01, 2.809917D-02, 1.852532D-05, * 1.589963D+00, 3.418075D-02, 1.430200D-05, * 2.994089D+00, 2.590824D-02, 4.155840D-06, * 8.155457D-01, 2.486352D-02, 6.836840D-06, * 1.735614D+00, 1.763719D-02, 6.370440D-06, * 1.968564D+00, 1.524020D-02,-2.517152D-06, * 1.282417D+00, 8.703393D-03, 2.289292D-05, * 2.280820D+00, 1.918010D-02, 4.484520D-06, * 4.833473D-02, 1.641773D-04,-4.654200D-07, * 5.589232D-02,-3.455092D-04,-7.388560D-07, * 4.634443D-02,-2.658234D-05, 7.757000D-08, * 8.997041D-03, 6.329728D-06,-1.939256D-09, * 2.284178D-02,-9.941590D-05, 6.787400D-08, * 4.350267D-02,-6.839749D-05,-2.714956D-07, * 1.348204D-02, 1.091504D-05, 6.903760D-07, * 3.106570D-02,-1.665665D-04,-1.590188D-07/ C DATA DCARGS/ 5.0974222D+00,-7.8604195454652D+02, * 3.9584962D+00,-5.7533848094674D+02, * 1.6338070D+00,-1.1506769618935D+03, * 2.5487111D+00,-3.9302097727326D+02, * 4.9255514D+00,-5.8849265665348D+02, * 1.3363463D+00,-5.5076098609303D+02, * 1.6072053D+00,-5.2237501616674D+02, * 1.3629480D+00,-1.1790629318198D+03, * 5.5657014D+00,-1.0977134971135D+03, * 5.0708205D+00,-1.5774000881978D+02, * 3.9318944D+00, 5.2963464780000D+01, * 4.8989497D+00, 3.9809289073258D+01, * 1.3097446D+00, 7.7540959633708D+01, * 3.5147141D+00, 7.9618578146517D+01, * 3.5413158D+00,-5.4868336758022D+02/ C DATA CCAMPS/ *-2.279594D-5, 1.407414D-5, 8.273188D-6, 1.340565D-5,-2.490817D-7, *-3.494537D-5, 2.860401D-7, 1.289448D-7, 1.627237D-5,-1.823138D-7, * 6.593466D-7, 1.322572D-5, 9.258695D-6,-4.674248D-7,-3.646275D-7, * 1.140767D-5,-2.049792D-5,-4.747930D-6,-2.638763D-6,-1.245408D-7, * 9.516893D-6,-2.748894D-6,-1.319381D-6,-4.549908D-6,-1.864821D-7, * 7.310990D-6,-1.924710D-6,-8.772849D-7,-3.334143D-6,-1.745256D-7, *-2.603449D-6, 7.359472D-6, 3.168357D-6, 1.119056D-6,-1.655307D-7, *-3.228859D-6, 1.308997D-7, 1.013137D-7, 2.403899D-6,-3.736225D-7, * 3.442177D-7, 2.671323D-6, 1.832858D-6,-2.394688D-7,-3.478444D-7, * 8.702406D-6,-8.421214D-6,-1.372341D-6,-1.455234D-6,-4.998479D-8, *-1.488378D-6,-1.251789D-5, 5.226868D-7,-2.049301D-7, 0.0D0, *-8.043059D-6,-2.991300D-6, 1.473654D-7,-3.154542D-7, 0.0D0, * 3.699128D-6,-3.316126D-6, 2.901257D-7, 3.407826D-7, 0.0D0, * 2.550120D-6,-1.241123D-6, 9.901116D-8, 2.210482D-7, 0.0D0, *-6.351059D-7, 2.341650D-6, 1.061492D-6, 2.878231D-7, 0.0D0/ C DATA CCSEC3/-7.757020D-08/ C DATA CCSEC/ 1.289600D-06, 5.550147D-01, 2.076942D+00, * 3.102810D-05, 4.035027D+00, 3.525565D-01, * 9.124190D-06, 9.990265D-01, 2.622706D+00, * 9.793240D-07, 5.508259D+00, 1.559103D+01/ C DATA DCSLD/ 1.990987D-07/, CCSGD/ 1.990969D-07/ C DATA CCKM/3.122140D-05/, CCMLD/2.661699D-06/, CCFDI/2.399485D-07/ C DATA DCARGM/ 5.1679830D+00, 8.3286911095275D+03, * 5.4913150D+00,-7.2140632838100D+03, * 5.9598530D+00, 1.5542754389685D+04/ C DATA CCAMPM/ * 1.097594D-01, 2.896773D-07, 5.450474D-02, 1.438491D-07, * -2.223581D-02, 5.083103D-08, 1.002548D-02,-2.291823D-08, * 1.148966D-02, 5.658888D-08, 8.249439D-03, 4.063015D-08/ C DATA CCPAMV/8.326827D-11,1.843484D-11,1.988712D-12,1.881276D-12/, * DC1MME/0.99999696D0/ C IDEQ=DEQ DT=(DJE-DCT0)/DCJUL T=DT DTSQ=DT*DT TSQ=DTSQ DO 100, K=1,8 DLOCAL=DMOD(DCFEL(1,K)+DT*DCFEL(2,K)+DTSQ*DCFEL(3,K),DC2PI) IF (K.EQ.1) DML=DLOCAL IF (K.NE.1) FORBEL(K-1)=DLOCAL 100 CONTINUE DEPS=DMOD(DCEPS(1)+DT*DCEPS(2)+DTSQ*DCEPS(3), DC2PI) DO 200, K=1,17 SORBEL(K)=DMOD(CCSEL(1,K)+T*CCSEL(2,K)+TSQ*CCSEL(3,K),CC2PI) 200 CONTINUE DO 300, K=1,4 A=DMOD(CCSEC(2,K)+T*CCSEC(3,K),CC2PI) SN(K)=DSIN(A) 300 CONTINUE PERTL = CCSEC(1,1) *SN(1) +CCSEC(1,2)*SN(2) * +(CCSEC(1,3)+T*CCSEC3)*SN(3) +CCSEC(1,4)*SN(4) PERTLD=0.0 PERTR =0.0 PERTRD=0.0 C DO 400, K=1,15 A=DMOD(DCARGS(1,K)+DT*DCARGS(2,K), DC2PI) COSA=DCOS(A) SINA=DSIN(A) PERTL =PERTL+CCAMPS(1,K)*COSA+CCAMPS(2,K)*SINA PERTR =PERTR+CCAMPS(3,K)*COSA+CCAMPS(4,K)*SINA IF (K.GE.11) GO TO 400 PERTLD=PERTLD+(CCAMPS(2,K)*COSA-CCAMPS(1,K)*SINA)*CCAMPS(5,K) PERTRD=PERTRD+(CCAMPS(4,K)*COSA-CCAMPS(3,K)*SINA)*CCAMPS(5,K) 400 CONTINUE ESQ=E*E DPARAM=DC1-ESQ PARAM=DPARAM TWOE=E+E TWOG=G+G PHI=TWOE*((1.0-ESQ*0.125D0)*DSIN(G)+E*0.625D0*DSIN(TWOG) * +ESQ*0.5416667D0*DSIN(G+TWOG) ) F=G+PHI SINF=DSIN(F) COSF=DCOS(F) DPSI=DPARAM/(DC1+E*COSF) PHID=TWOE*CCSGD*((1.0+ESQ*1.5D0)*COSF+E*(1.25D0-SINF*SINF*0.5D0)) PSID=CCSGD*E*SINF/SQRT(PARAM) D1PDRO=(DC1+PERTR) DRD=D1PDRO*(PSID+DPSI*PERTRD) DRLD=D1PDRO*DPSI*(DCSLD+PHID+PERTLD) DTL=DMOD(DML+PHI+PERTL, DC2PI) DSINLS=DSIN(DTL) DCOSLS=DCOS(DTL) DXHD = DRD*DCOSLS-DRLD*DSINLS DYHD = DRD*DSINLS+DRLD*DCOSLS PERTL =0.0 PERTLD=0.0 PERTP =0.0 PERTPD=0.0 C DO 500 K=1,3 A=DMOD(DCARGM(1,K)+DT*DCARGM(2,K), DC2PI) SINA =DSIN(A) COSA =DCOS(A) PERTL =PERTL +CCAMPM(1,K)*SINA PERTLD=PERTLD+CCAMPM(2,K)*COSA PERTP =PERTP +CCAMPM(3,K)*COSA PERTPD=PERTPD-CCAMPM(4,K)*SINA 500 CONTINUE TL=FORBEL(2)+PERTL SINLM=DSIN(TL) COSLM=DCOS(TL) SIGMA=CCKM/(1.0+PERTP) A=SIGMA*(CCMLD+PERTLD) B=SIGMA*PERTPD DXHD=DXHD+A*SINLM+B*COSLM DYHD=DYHD-A*COSLM+B*SINLM DZHD= -SIGMA*CCFDI*DCOS(FORBEL(3)) DXBD=DXHD*DC1MME DYBD=DYHD*DC1MME DZBD=DZHD*DC1MME DO 600 K=1,4 PLON=FORBEL(K+3) POMG=SORBEL(K+1) PECC=SORBEL(K+9) TL=DMOD(PLON+2.0*PECC*DSIN(PLON-POMG), CC2PI) SINLP(K)=DSIN(TL) COSLP(K)=DCOS(TL) DXBD=DXBD+CCPAMV(K)*(SINLP(K)+PECC*DSIN(POMG)) DYBD=DYBD-CCPAMV(K)*(COSLP(K)+PECC*DCOS(POMG)) DZBD=DZBD-CCPAMV(K)*SORBEL(K+13)*DCOS(PLON-SORBEL(K+5)) 600 CONTINUE DCOSEP=DCOS(DEPS) DSINEP=DSIN(DEPS) DYAHD=DCOSEP*DYHD-DSINEP*DZHD DZAHD=DSINEP*DYHD+DCOSEP*DZHD DYABD=DCOSEP*DYBD-DSINEP*DZBD DZABD=DSINEP*DYBD+DCOSEP*DZBD DVELH(1)=DXHD DVELH(2)=DYAHD DVELH(3)=DZAHD DVELB(1)=DXBD DVELB(2)=DYABD DVELB(3)=DZABD DO 800 N=1,3 DVELH(N)=DVELH(N)*1.4959787D8 DVELB(N)=DVELB(N)*1.4959787D8 800 CONTINUE C RETURN END SUBROUTINE BARCOR(DCORH,DCORB) C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C PURPOSE : calculate rectangular barycentric and heliocentric coordinates C of the earth C C REFERENCE : STUMPFF,P.,ASTRON. ASTROPHYS. SUPPL. 41,1,1980 C MODIFICATION : D. GILLET 1983-9-15 C C AUXILIARY ROUTINES CALLED : BARVEL,PRE C C------------------------------------------------------------------------------- C IMPLICIT NONE !REAL*8 (A-H,O-Z) C DOUBLE PRECISION DCORH(3),DCORB(3),CCPAM(4) DOUBLE PRECISION DR,CCIM,FLAT,FLATM DOUBLE PRECISION DXH,DYH,DZH,DXB,DYB,DZB DOUBLE PRECISION DYAH,DZAH,DYAB,DZAB DOUBLE PRECISION A,B DOUBLE PRECISION DC1MME C DOUBLE PRECISION DPREMA(3,3),DPSI,D1PDRO,DSINLS,DCOSLS DOUBLE PRECISION DSINEP,DCOSEP,FORBEL(7),SORBEL(17),SINLP(4) DOUBLE PRECISION COSLP(4),SINLM,COSLM,SIGMA C INTEGER IDEQ,K,N C COMMON /BARXYZ/ DPREMA,DPSI,D1PDRO,DSINLS,DCOSLS, + DSINEP,DCOSEP,FORBEL,SORBEL,SINLP, + COSLP,SINLM,COSLM,SIGMA,IDEQ C DATA CCPAM/4.960906D-3,2.727436D-3,8.392311D-4,1.556861D-3/, * CCIM /8.978749D-2/, DC1MME/0.99999696D0/ C DR=DPSI*D1PDRO FLATM=CCIM*DSIN(FORBEL(3)) A=SIGMA*DCOS(FLATM) DXH=DR*DCOSLS-A*COSLM DYH=DR*DSINLS-A*SINLM DZH= -SIGMA*DSIN(FLATM) DXB=DXH*DC1MME DYB=DYH*DC1MME DZB=DZH*DC1MME DO 10, K=1,4 FLAT=SORBEL(K+13)*DSIN(FORBEL(K+3)-SORBEL(K+5)) A=CCPAM(K)*(1.0-SORBEL(K+9)*DCOS(FORBEL(K+3)-SORBEL(K+1))) B=A*DCOS(FLAT) DXB=DXB-B*COSLP(K) DYB=DYB-B*SINLP(K) DZB=DZB-A*DSIN(FLAT) 10 CONTINUE DYAH=DCOSEP*DYH-DSINEP*DZH DZAH=DSINEP*DYH+DCOSEP*DZH DYAB=DCOSEP*DYB-DSINEP*DZB DZAB=DSINEP*DYB+DCOSEP*DZB IF(IDEQ.NE.0) GO TO 20 C DCORH(1)=DXH DCORH(2)=DYAH DCORH(3)=DZAH DCORB(1)=DXB DCORB(2)=DYAB DCORB(3)=DZAB RETURN C 20 DO 30, N=1,3 DCORH(N)=DXH*DPREMA(N,1)+DYAH*DPREMA(N,2)+DZAH*DPREMA(N,3) DCORB(N)=DXB*DPREMA(N,1)+DYAB*DPREMA(N,2)+DZAB*DPREMA(N,3) 30 CONTINUE C RETURN END SUBROUTINE JULDAT(INDATE,UTR,JD) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C FORTRAN subroutine JULDAT version 1.0 870102 C original coding: D. Gillet ESO - Garching C variables renamed and restructured: D. Baade ST-ECF, Garching C C.KEYWORDS C geocentric Julian date C C.PURPOSE C calculate geocentric Julian date for any civil date (time in UT) C C.ALGORITHM C adapted from MEEUS J.,1980, ASTRONOMICAL FORMULAE FOR CALCULATORS C C.INPUT/OUTPUT C the following are passed from and to the calling program: C INDATE(3) : civil date as year,month,day OR year.fraction C UT : universal time expressed in real hours C JD : real geocentric Julian date C C.REVISIONS C made to accept also REAL dates D. Baade 910408 C C------------------------------------------------------------------------------- C IMPLICIT NONE C DOUBLE PRECISION YP,P,C,A,UT DOUBLE PRECISION UTR,JD C INTEGER STAT,IA,IB,IC,ND,DATE(3) C REAL INDATE(3),FRAC C UT=UTR / 24.0D0 C C CHECK FORMAT OF DATE: may be either year,month,date OR year.fraction,0,0 C (Note that the fraction of the year must NOT include fractions of a day.) C For all other formats exit and terminate also calling command sequence. C IF ((INDATE(1)-INT(INDATE(1))).GT.1.0E-6) THEN IF ((INDATE(2).GT.1.0E-6).OR.(INDATE(3).GT.1.0E-6)) + CALL STETER(1,'Error: Date was entered in wrong format.') C C copy date input buffer copy to other buffer so that calling program C does not notice any changes C C FIRST CASE: format was year.fraction C DATE(1)=INT(INDATE(1)) FRAC=INDATE(1)-DATE(1) DATE(2)=1 DATE(3)=1 ELSE C C SECOND CASE: format was year,month,day C DATE(1)=NINT(INDATE(1)) FRAC=0 DATE(2)=NINT(INDATE(2)) DATE(3)=NINT(INDATE(3)) IF ((DATE(2).EQ.0).AND.(DATE(3).EQ.0)) THEN DATE(2)=1 DATE(3)=1 ENDIF IF ((DATE(2).LT.1).OR.(DATE(2).GT.12)) + CALL STETER(1,'Error: such a month does not exist') IF ((DATE(3).LT.1).OR.(DATE(3).GT.31)) + CALL STETER(1,'Error: such a day does not exist') ENDIF C C from here on, the normal procedure applies which is based on the C format year,month,day: C IF (DATE(2) .GT. 2) THEN YP=DATE(1) P=DATE(2) ELSE YP=DATE(1)-1 P=DATE(2)+12.0 ENDIF C C = DATE(1) + DATE(2)*1.D-2 + DATE(3)*1.D-4 + UT*1.D-6 IF (C .GE. 1582.1015D0) THEN IA=IDINT(YP/100.D0) A=DBLE(IA) IB=2-IA+IDINT(A/4.D0) ELSE IB=0 ENDIF C JD = DINT(365.25D0*YP) + DINT(30.6001D0*(P+1.D0)) + DATE(3) + UT * + DBLE(IB) + 1720994.5D0 C C finally, take into account fraction of year (if any), respect leap C year conventions C IF (FRAC.GT.1.0E-6) THEN ND=365 IF (C.GE.1582.1015D0) THEN IC = MOD(DATE(1),4) IF (IC.EQ.0) THEN ND=366 IC = MOD(DATE(1),100) IF (IC.EQ.0) THEN IC = MOD(DATE(1),400) IF (IC.NE.0) ND=365 ENDIF ENDIF ENDIF C IF (ABS(FRAC*ND-NINT(FRAC*ND)).GT.0.3) THEN CALL STTPUT + ('Warning: Fraction of year MAY not correspond to ',STAT) CALL STTPUT(' integer number of days.',STAT) ENDIF C JD = JD+NINT(FRAC*ND) ENDIF C RETURN END