C @(#)compxy.for	19.1 (ESO-DMD) 02/25/03 14:01:41
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
