C @(#)precess.for	19.1 (ES0-DMD) 02/25/03 14:27:54
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 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding 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
C @(#)precess.for	19.1 (ESO-IPG) 14:27:54 02/25/03
        PROGRAM PRECES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: PRECESS.FOR
C.PURPOSE: Precess RA and DEC coordinates in table created by CREATE/OPTABLE
C.         command.
C.AUTHOR:  Alessandra Gemmo           Padova Department of Astronomy
C.VERSION: 910826 AG Creation
C------------------------------------------------------------------------------
        IMPLICIT NONE
C
C ***
C
        INTEGER    MADRID(1)
        INTEGER    NRINP
        INTEGER    NCINP
        INTEGER    NROUT
        INTEGER    NCOUT
        INTEGER    NCPAR1
        PARAMETER  (NCOUT=15)
        PARAMETER  (NCPAR1=2)
        PARAMETER  (NROUT=300)
        INTEGER    ISTAT,IAC
        INTEGER    KUN,KNUL
        INTEGER    NACINP,NARINP,NSINP
        INTEGER    OUTTYP,OUTCOL(2),COLNUM(2)
        INTEGER    TIDINP
        INTEGER    I1,I2
C
C ***
C
        DOUBLE PRECISION    RARAD,DECRAD
        DOUBLE PRECISION    PRA,PDEC,EQUI0,EQUI1
        DOUBLE PRECISION    PRARAD,PDECRA
        DOUBLE PRECISION    RA,DEC,PI
C
C ***
C
        CHARACTER*80        STRING(10)
        CHARACTER*60        INPFIL
        CHARACTER*16        LABEL1(NCPAR1),OUTLAB
        CHARACTER*16        UNIT1(NCPAR1),OUTUNI
        CHARACTER*16        OUTFOR
        CHARACTER*16        FORMR8(2)
C
C ***
C
        LOGICAL             NULL
C
C ***
C
        INCLUDE             'MID_INCLUDE:ST_DEF.INC'
        COMMON              /VMR/MADRID
        INCLUDE             'MID_INCLUDE:ST_DAT.INC'
C
        DATA                LABEL1/'PRA     ','PDEC    '/
        DATA                UNIT1 /'HOURS   ','DEGREES '/
C
C ***
C
        DATA                FORMR8/'F12.9','F12.9'/
        DATA                PI/3.14159265358979D0/
C
C *** start the code
C
        CALL STSPRO('OPTAB')
C
C *** get the input table
C
        CALL STKRDC('INPUTFI',1,1,60,IAC,INPFIL,KUN,KNUL,ISTAT)
C
C *** open input table
C
        CALL TBTOPN(INPFIL,F_IO_MODE,TIDINP,ISTAT)
C
C *** get information about the input table
C
        CALL TBIGET(TIDINP,NCINP,NRINP,NSINP,NACINP,NARINP,ISTAT)
        IF(NRINP.EQ.0)THEN         !no rows in table
          STRING(1) = '*** FATAL: There are no data in input table'
          CALL STETER(9,STRING(1))
        ENDIF
C
C *** initialize PRA and PDEC columns
C
        DO I1=1,NCPAR1
          OUTTYP = D_R8_FORMAT
          OUTFOR = FORMR8(I1)
          OUTUNI = UNIT1(I1)
          OUTLAB = LABEL1(I1)
          CALL TBCINI(TIDINP,OUTTYP,1,OUTFOR,OUTUNI,
     2                       OUTLAB,OUTCOL(I1),ISTAT)
        ENDDO
C
C *** fill in PRA and PDEC columns
C
        CALL STKRDD('EQUINOX',1,1,IAC,EQUI0,KUN,KNUL,ISTAT)
        CALL STKRDD('EQUINOX',2,1,IAC,EQUI1,KUN,KNUL,ISTAT)
C
C *** read in RA and DEC columns
C
        DO I2=1,NRINP
          CALL TBLSER(TIDINP,'RA',COLNUM(1),ISTAT)
          CALL TBERDD(TIDINP,I2,COLNUM(1),RA,NULL,ISTAT)
C
          CALL TBLSER(TIDINP,'DEC',COLNUM(2),ISTAT)
          CALL TBERDD(TIDINP,I2,COLNUM(2),DEC,NULL,ISTAT)
C
          RARAD  = RA*1.5D1*(PI/1.8D2)
          DECRAD = DEC*(PI/1.8D2)
C
C *** compute precession
C
          CALL PRE(RARAD,DECRAD,PRARAD,PDECRA,EQUI0,EQUI1)
C
          PRA    = PRARAD*(1.8D2/PI)/15.
          PDEC   = PDECRA*(1.8D2/PI)
C
          CALL TBEWRD(TIDINP,I2,OUTCOL(1),PRA,ISTAT)
          CALL TBEWRD(TIDINP,I2,OUTCOL(2),PDEC,ISTAT)
        ENDDO
C
C *** close the table
C
        CALL TBTCLO(TIDINP,ISTAT)
C
C *** over and out
C
        CALL STSEPI
        END
C------------------------------------------------------------------------------
      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
      DOUBLE PRECISION DT0,DT,DTS,DTC,DC1900,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
      DATA DCSAR/4.848136812D-6/,DC1900/1900.0D0/,DC1M2/0.01D0/,
     * DC1/2304.25D0/,DC2/1.396D0/,DC3/0.302D0/,DC4/0.018D0/,
     * DC5/0.791D0/,DC6/2004.683D0/,DC7/-0.853D0/,DC8/-0.426D0/,
     * DC9/-0.042D0/
C
      PI = 3.1415926535897928
C
C ... first determine the precession angles zeta, z and theta
C
      DT0=(EQX0-DC1900)*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
