← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/optopus/src/precess.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/optopus/src/precess.for
 
1
C @(#)precess.for       19.1 (ES0-DMD) 02/25/03 14:27:54
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C @(#)precess.for       19.1 (ESO-IPG) 14:27:54 02/25/03
 
30
        PROGRAM PRECES
 
31
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
32
C.IDENTIFICATION: PRECESS.FOR
 
33
C.PURPOSE: Precess RA and DEC coordinates in table created by CREATE/OPTABLE
 
34
C.         command.
 
35
C.AUTHOR:  Alessandra Gemmo           Padova Department of Astronomy
 
36
C.VERSION: 910826 AG Creation
 
37
C------------------------------------------------------------------------------
 
38
        IMPLICIT NONE
 
39
C
 
40
C ***
 
41
C
 
42
        INTEGER    MADRID(1)
 
43
        INTEGER    NRINP
 
44
        INTEGER    NCINP
 
45
        INTEGER    NROUT
 
46
        INTEGER    NCOUT
 
47
        INTEGER    NCPAR1
 
48
        PARAMETER  (NCOUT=15)
 
49
        PARAMETER  (NCPAR1=2)
 
50
        PARAMETER  (NROUT=300)
 
51
        INTEGER    ISTAT,IAC
 
52
        INTEGER    KUN,KNUL
 
53
        INTEGER    NACINP,NARINP,NSINP
 
54
        INTEGER    OUTTYP,OUTCOL(2),COLNUM(2)
 
55
        INTEGER    TIDINP
 
56
        INTEGER    I1,I2
 
57
C
 
58
C ***
 
59
C
 
60
        DOUBLE PRECISION    RARAD,DECRAD
 
61
        DOUBLE PRECISION    PRA,PDEC,EQUI0,EQUI1
 
62
        DOUBLE PRECISION    PRARAD,PDECRA
 
63
        DOUBLE PRECISION    RA,DEC,PI
 
64
C
 
65
C ***
 
66
C
 
67
        CHARACTER*80        STRING(10)
 
68
        CHARACTER*60        INPFIL
 
69
        CHARACTER*16        LABEL1(NCPAR1),OUTLAB
 
70
        CHARACTER*16        UNIT1(NCPAR1),OUTUNI
 
71
        CHARACTER*16        OUTFOR
 
72
        CHARACTER*16        FORMR8(2)
 
73
C
 
74
C ***
 
75
C
 
76
        LOGICAL             NULL
 
77
C
 
78
C ***
 
79
C
 
80
        INCLUDE             'MID_INCLUDE:ST_DEF.INC'
 
81
        COMMON              /VMR/MADRID
 
82
        INCLUDE             'MID_INCLUDE:ST_DAT.INC'
 
83
C
 
84
        DATA                LABEL1/'PRA     ','PDEC    '/
 
85
        DATA                UNIT1 /'HOURS   ','DEGREES '/
 
86
C
 
87
C ***
 
88
C
 
89
        DATA                FORMR8/'F12.9','F12.9'/
 
90
        DATA                PI/3.14159265358979D0/
 
91
C
 
92
C *** start the code
 
93
C
 
94
        CALL STSPRO('OPTAB')
 
95
C
 
96
C *** get the input table
 
97
C
 
98
        CALL STKRDC('INPUTFI',1,1,60,IAC,INPFIL,KUN,KNUL,ISTAT)
 
99
C
 
100
C *** open input table
 
101
C
 
102
        CALL TBTOPN(INPFIL,F_IO_MODE,TIDINP,ISTAT)
 
103
C
 
104
C *** get information about the input table
 
105
C
 
106
        CALL TBIGET(TIDINP,NCINP,NRINP,NSINP,NACINP,NARINP,ISTAT)
 
107
        IF(NRINP.EQ.0)THEN         !no rows in table
 
108
          STRING(1) = '*** FATAL: There are no data in input table'
 
109
          CALL STETER(9,STRING(1))
 
110
        ENDIF
 
111
C
 
112
C *** initialize PRA and PDEC columns
 
113
C
 
114
        DO I1=1,NCPAR1
 
115
          OUTTYP = D_R8_FORMAT
 
116
          OUTFOR = FORMR8(I1)
 
117
          OUTUNI = UNIT1(I1)
 
118
          OUTLAB = LABEL1(I1)
 
119
          CALL TBCINI(TIDINP,OUTTYP,1,OUTFOR,OUTUNI,
 
120
     2                       OUTLAB,OUTCOL(I1),ISTAT)
 
121
        ENDDO
 
122
C
 
123
C *** fill in PRA and PDEC columns
 
124
C
 
125
        CALL STKRDD('EQUINOX',1,1,IAC,EQUI0,KUN,KNUL,ISTAT)
 
126
        CALL STKRDD('EQUINOX',2,1,IAC,EQUI1,KUN,KNUL,ISTAT)
 
127
C
 
128
C *** read in RA and DEC columns
 
129
C
 
130
        DO I2=1,NRINP
 
131
          CALL TBLSER(TIDINP,'RA',COLNUM(1),ISTAT)
 
132
          CALL TBERDD(TIDINP,I2,COLNUM(1),RA,NULL,ISTAT)
 
133
C
 
134
          CALL TBLSER(TIDINP,'DEC',COLNUM(2),ISTAT)
 
135
          CALL TBERDD(TIDINP,I2,COLNUM(2),DEC,NULL,ISTAT)
 
136
C
 
137
          RARAD  = RA*1.5D1*(PI/1.8D2)
 
138
          DECRAD = DEC*(PI/1.8D2)
 
139
C
 
140
C *** compute precession
 
141
C
 
142
          CALL PRE(RARAD,DECRAD,PRARAD,PDECRA,EQUI0,EQUI1)
 
143
C
 
144
          PRA    = PRARAD*(1.8D2/PI)/15.
 
145
          PDEC   = PDECRA*(1.8D2/PI)
 
146
C
 
147
          CALL TBEWRD(TIDINP,I2,OUTCOL(1),PRA,ISTAT)
 
148
          CALL TBEWRD(TIDINP,I2,OUTCOL(2),PDEC,ISTAT)
 
149
        ENDDO
 
150
C
 
151
C *** close the table
 
152
C
 
153
        CALL TBTCLO(TIDINP,ISTAT)
 
154
C
 
155
C *** over and out
 
156
C
 
157
        CALL STSEPI
 
158
        END
 
159
C------------------------------------------------------------------------------
 
160
      SUBROUTINE PRE(A0,D0,A1,D1,EQX0,EQX1)
 
161
C
 
162
C ... calculate general precession for two given epochs
 
163
C
 
164
C ... References: STUMPFF P.,ASTRON. ASTROPHYS. SUPPL. 41, P. 1 (1980)
 
165
C ...             EXPLAN. SUPPL. ASTRON. EPHEREMERIS, P. 28 (1961)
 
166
C
 
167
      IMPLICIT NONE             !REAL*8 (A-H,O-Z)
 
168
C
 
169
      DOUBLE PRECISION A0,D0,D1,EQX0,EQX1
 
170
      DOUBLE PRECISION DT0,DT,DTS,DTC,DC1900,DC1M2
 
171
      DOUBLE PRECISION DC1,DC2,DC3,DC4,DC5,DC6,DC7,DC8,DC9
 
172
      DOUBLE PRECISION DZED,DZETA,DCSAR,DTHETA
 
173
      DOUBLE PRECISION PI,SINR,COSR,A1,R
 
174
C
 
175
      DATA DCSAR/4.848136812D-6/,DC1900/1900.0D0/,DC1M2/0.01D0/,
 
176
     * DC1/2304.25D0/,DC2/1.396D0/,DC3/0.302D0/,DC4/0.018D0/,
 
177
     * DC5/0.791D0/,DC6/2004.683D0/,DC7/-0.853D0/,DC8/-0.426D0/,
 
178
     * DC9/-0.042D0/
 
179
C
 
180
      PI = 3.1415926535897928
 
181
C
 
182
C ... first determine the precession angles zeta, z and theta
 
183
C
 
184
      DT0=(EQX0-DC1900)*DC1M2
 
185
      DT=(EQX1-EQX0)*DC1M2
 
186
      DTS=DT*DT
 
187
      DTC=DTS*DT
 
188
      DZETA=((DC1+DC2*DT0)*DT+DC3*DTS+DC4*DTC)*DCSAR
 
189
      DZED=DZETA + DC5*DTS*DCSAR
 
190
      DTHETA=((DC6+DC7*DT0)*DT+DC8*DTS+DC9*DTC)*DCSAR
 
191
C
 
192
C ... with these and A0,D0 compute the precessed coordinates A1,D1
 
193
C
 
194
      D1 = ASIN(DCOS(DTHETA)*DSIN(D0)+DSIN(DTHETA)*DCOS(D0)*
 
195
     1                                              DCOS(A0+DZETA))
 
196
      SINR = (DCOS(D0)*DSIN(A0+DZETA))/DCOS(D1)
 
197
      COSR = (DCOS(DTHETA)*DCOS(D0)*DCOS(A0+DZETA)-
 
198
     1                         DSIN(DTHETA)*DSIN(D0))/DCOS(D1)
 
199
      R = DASIN(SINR)
 
200
      IF (COSR.LT.0) R=PI-R                ! determine quadrant of R
 
201
      A1 = R+DZED                      ! this is the precessed R.A.
 
202
      IF (A1.GT.2.0D0*PI) A1=A1-2.0D0*PI
 
203
      IF (A1.LT.0.0D0) A1=A1+2.0D0*PI
 
204
C
 
205
      RETURN
 
206
      END