aboutsummaryrefslogtreecommitdiff
path: root/math/slalib/mapqk.f
blob: 0e1997ff944b780a0577ac454577c9e019967c7c (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
SUBROUTINE slMAPQ (RM, DM, PR, PD, PX, RV, AMPRMS, RA, DA)
*+
*     - - - - - -
*      M A P Q
*     - - - - - -
*
*  Quick mean to apparent place:  transform a star RA,Dec from
*  mean place to geocentric apparent place, given the
*  star-independent parameters.
*
*  Use of this routine is appropriate when efficiency is important
*  and where many star positions, all referred to the same equator
*  and equinox, are to be transformed for one epoch.  The
*  star-independent parameters can be obtained by calling the
*  slMAPA routine.
*
*  If the parallax and proper motions are zero the slMAPZ
*  routine can be used instead.
*
*  The reference frames and timescales used are post IAU 1976.
*
*  Given:
*     RM,DM    d      mean RA,Dec (rad)
*     PR,PD    d      proper motions:  RA,Dec changes per Julian year
*     PX       d      parallax (arcsec)
*     RV       d      radial velocity (km/sec, +ve if receding)
*
*     AMPRMS   d(21)  star-independent mean-to-apparent parameters:
*
*       (1)      time interval for proper motion (Julian years)
*       (2-4)    barycentric position of the Earth (AU)
*       (5-7)    heliocentric direction of the Earth (unit vector)
*       (8)      (grav rad Sun)*2/(Sun-Earth distance)
*       (9-11)   barycentric Earth velocity in units of c
*       (12)     sqrt(1-v**2) where v=modulus(ABV)
*       (13-21)  precession/nutation (3,3) matrix
*
*  Returned:
*     RA,DA    d      apparent RA,Dec (rad)
*
*  References:
*     1984 Astronomical Almanac, pp B39-B41.
*     (also Lederle & Schwan, Astron. Astrophys. 134,
*      1-6, 1984)
*
*  Notes:
*
*  1)  The vectors AMPRMS(2-4) and AMPRMS(5-7) are referred to
*      the mean equinox and equator of epoch EQ.
*
*  2)  Strictly speaking, the routine is not valid for solar-system
*      sources, though the error will usually be extremely small.
*      However, to prevent gross errors in the case where the
*      position of the Sun is specified, the gravitational
*      deflection term is restrained within about 920 arcsec of the
*      centre of the Sun's disc.  The term has a maximum value of
*      about 1.85 arcsec at this radius, and decreases to zero as
*      the centre of the disc is approached.
*
*  Called:
*     slDS2C       spherical to Cartesian
*     slDVDV        dot product
*     slDMXV        matrix x vector
*     slDC2S       Cartesian to spherical
*     slDA2P      normalize angle 0-2Pi
*
*  P.T.Wallace   Starlink   15 January 2000
*
*  Copyright (C) 2000 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION RM,DM,PR,PD,PX,RV,AMPRMS(21),RA,DA

*  Arc seconds to radians
      DOUBLE PRECISION AS2R
      PARAMETER (AS2R=0.484813681109535994D-5)

*  Km/s to AU/year
      DOUBLE PRECISION VF
      PARAMETER (VF=0.21094502D0)

      INTEGER I

      DOUBLE PRECISION PMT,GR2E,AB1,EB(3),EHN(3),ABV(3),
     :                 Q(3),PXR,W,EM(3),P(3),PN(3),PDE,PDEP1,
     :                 P1(3),P1DV,P2(3),P3(3)

      DOUBLE PRECISION slDVDV,slDA2P



*  Unpack scalar and vector parameters
      PMT = AMPRMS(1)
      GR2E = AMPRMS(8)
      AB1 = AMPRMS(12)
      DO I=1,3
         EB(I) = AMPRMS(I+1)
         EHN(I) = AMPRMS(I+4)
         ABV(I) = AMPRMS(I+8)
      END DO

*  Spherical to x,y,z
      CALL slDS2C(RM,DM,Q)

*  Space motion (radians per year)
      PXR = PX*AS2R
      W = VF*RV*PXR
      EM(1) = -PR*Q(2)-PD*COS(RM)*SIN(DM)+W*Q(1)
      EM(2) =  PR*Q(1)-PD*SIN(RM)*SIN(DM)+W*Q(2)
      EM(3) =          PD*COS(DM)        +W*Q(3)

*  Geocentric direction of star (normalized)
      DO I=1,3
         P(I) = Q(I)+PMT*EM(I)-PXR*EB(I)
      END DO
      CALL slDVN(P,PN,W)

*  Light deflection (restrained within the Sun's disc)
      PDE = slDVDV(PN,EHN)
      PDEP1 = PDE+1D0
      W = GR2E/MAX(PDEP1,1D-5)
      DO I=1,3
         P1(I) = PN(I)+W*(EHN(I)-PDE*PN(I))
      END DO

*  Aberration (normalization omitted)
      P1DV = slDVDV(P1,ABV)
      W = 1D0+P1DV/(AB1+1D0)
      DO I=1,3
         P2(I) = AB1*P1(I)+W*ABV(I)
      END DO

*  Precession and nutation
      CALL slDMXV(AMPRMS(13),P2,P3)

*  Geocentric apparent RA,Dec
      CALL slDC2S(P3,RA,DA)
      RA = slDA2P(RA)

      END