C @(#)precess.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:14 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 17.1.1.1 (ESO-IPG) 17:55:14 01/25/02 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