~ubuntu-branches/ubuntu/raring/starlink-pal/raring-proposed

« back to all changes in this revision

Viewing changes to palRdplan.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2012-03-28 11:00:00 UTC
  • Revision ID: package-import@ubuntu.com-20120328110000-9iw40yuy69wuk7po
Tags: upstream-0.1.0
ImportĀ upstreamĀ versionĀ 0.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*
 
2
*+
 
3
*  Name:
 
4
*     palRdplan
 
5
 
 
6
*  Purpose:
 
7
*     Approximate topocentric apparent RA,Dec of a planet
 
8
 
 
9
*  Language:
 
10
*     Starlink ANSI C
 
11
 
 
12
*  Type of Module:
 
13
*     Library routine
 
14
 
 
15
*  Invocation:
 
16
*     void palRdplan( double date, int np, double elong, double phi,
 
17
*                     double * ra, double * dec, double * diam );
 
18
 
 
19
*  Arguments:
 
20
*     date = double (Given)
 
21
*        MJD of observation (JD-2400000.5) in TDB. For all practical
 
22
*        purposes TT can be used instead of TDB, and for many applications
 
23
*        UT will do (except for the Moon).
 
24
*     np = int (Given)
 
25
*        Planet: 1 = Mercury
 
26
*                2 = Venus
 
27
*                3 = Moon
 
28
*                4 = Mars
 
29
*                5 = Jupiter
 
30
*                6 = Saturn
 
31
*                7 = Uranus
 
32
*                8 = Neptune
 
33
*             else = Sun
 
34
*     elong = double (Given)
 
35
*        Observer's east longitude (radians)
 
36
*     phi = double (Given)
 
37
*        Observer's geodetic latitude (radians)
 
38
*     ra = double * (Returned)
 
39
*        RA (topocentric apparent, radians)
 
40
*     dec = double * (Returned)
 
41
*        Dec (topocentric apparent, radians)
 
42
*     diam = double * (Returned)
 
43
*        Angular diameter (equatorial, radians)
 
44
 
 
45
*  Description:
 
46
*     Approximate topocentric apparent RA,Dec of a planet, and its
 
47
*     angular diameter.
 
48
 
 
49
*  Authors:
 
50
*     TIMJ: Tim Jenness (JAC, Hawaii)
 
51
*     {enter_new_authors_here}
 
52
 
 
53
*  Notes:
 
54
*     - Unlike with slaRdplan, Pluto is not supported.
 
55
*     - The longitude and latitude allow correction for geocentric
 
56
*       parallax.  This is a major effect for the Moon, but in the
 
57
*       context of the limited accuracy of the present routine its
 
58
*       effect on planetary positions is small (negligible for the
 
59
*       outer planets).  Geocentric positions can be generated by
 
60
*       appropriate use of the routines palDmoon and iauPlan94.
 
61
 
 
62
*  History:
 
63
*     2012-03-07 (TIMJ):
 
64
*        Initial version, with some documentation from SLA/F.
 
65
*     {enter_further_changes_here}
 
66
 
 
67
*  Copyright:
 
68
*     Copyright (C) 2012 Science and Technology Facilities Council.
 
69
*     All Rights Reserved.
 
70
 
 
71
*  Licence:
 
72
*     This program is free software; you can redistribute it and/or
 
73
*     modify it under the terms of the GNU General Public License as
 
74
*     published by the Free Software Foundation; either version 3 of
 
75
*     the License, or (at your option) any later version.
 
76
*
 
77
*     This program is distributed in the hope that it will be
 
78
*     useful, but WITHOUT ANY WARRANTY; without even the implied
 
79
*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 
80
*     PURPOSE. See the GNU General Public License for more details.
 
81
*
 
82
*     You should have received a copy of the GNU General Public License
 
83
*     along with this program; if not, write to the Free Software
 
84
*     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 
85
*     MA 02110-1301, USA.
 
86
 
 
87
*  Bugs:
 
88
*     {note_any_bugs_here}
 
89
*-
 
90
*/
 
91
 
 
92
#include <math.h>
 
93
 
 
94
#include "pal.h"
 
95
#include "palmac.h"
 
96
#include "sofa.h"
 
97
 
 
98
void palRdplan( double date, int np, double elong, double phi,
 
99
                double * ra, double * dec, double * diam ) {
 
100
 
 
101
  /* AU in km */
 
102
  const double AUKM = 1.49597870e8;
 
103
 
 
104
  /* Equatorial radii (km) */
 
105
  const double EQRAU[] = {
 
106
    696000.0, /* Sun */
 
107
      2439.7,
 
108
      6051.9,
 
109
      1738,
 
110
      3397,
 
111
     71492,
 
112
     60268,
 
113
     25559,
 
114
     24764
 
115
  };
 
116
 
 
117
  /* Local variables */
 
118
  int i, j;
 
119
  double stl;
 
120
  double vgm[6];
 
121
  double v[6];
 
122
  double rmat[3][3];
 
123
  double vse[6];
 
124
  double vsg[6];
 
125
  double vsp[6];
 
126
  double vgo[6];
 
127
  double dx,dy,dz,r,tl;
 
128
 
 
129
  /* Classify np */
 
130
  if (np < 0 || np > 8 ) np=0;  /* Sun */
 
131
 
 
132
  /* Approximate local sidereal time */
 
133
  stl = palGmst( date - palDt( palEpj(date)) / 86400.0) + elong;
 
134
 
 
135
  /* Geocentre to Moon (mean of date) */
 
136
  palDmoon( date, v );
 
137
 
 
138
  /* Nutation to true of date */
 
139
  palNut( date, rmat );
 
140
  iauRxp( rmat, v, vgm );
 
141
  iauRxp( rmat, &(v[3]), &(vgm[3]) );
 
142
 
 
143
  /* Moon? */
 
144
  if (np == 3) {
 
145
 
 
146
    /* geocentre to Moon (true of date) */
 
147
    for (i=0; i<6; i++) {
 
148
      v[i] = vgm[i];
 
149
    }
 
150
 
 
151
  } else {
 
152
 
 
153
    /* Not moon: precession/nutation matrix J2000 to date */
 
154
    palPrenut( 2000.0, date, rmat );
 
155
 
 
156
    /* Sun to Earth-Moon Barycentre (J2000) */
 
157
    palPlanet( date, 3, v, &j );
 
158
 
 
159
    /* Precession and nutation to date */
 
160
    iauRxp( rmat, v, vse );
 
161
    iauRxp( rmat, &(v[3]), &(vse[3]) );
 
162
 
 
163
    /* Sun to geocentre (true of date) */
 
164
    for (i=0; i<6; i++) {
 
165
      vsg[i] = vse[i] - 0.012150581 * vgm[i];
 
166
    }
 
167
 
 
168
    /* Sun ? */
 
169
    if (np == 0) {
 
170
 
 
171
      /* Geocentre to Sun */
 
172
      for (i=0; i<6; i++) {
 
173
        v[i] = -vsg[i];
 
174
      }
 
175
 
 
176
    } else {
 
177
 
 
178
      /* Sun to Planet (J2000) */
 
179
      palPlanet( date, np, v, &j );
 
180
 
 
181
      /* Precession and nutation to date */
 
182
      iauRxp( rmat, v, vsp );
 
183
      iauRxp( rmat, &(v[3]), &(vsp[3]) );
 
184
 
 
185
      /* Geocentre to planet */
 
186
      for (i=0; i<6; i++) {
 
187
        v[i] = vsp[i] - vsg[i];
 
188
      }
 
189
 
 
190
    }
 
191
 
 
192
  }
 
193
 
 
194
  /* Refer to origina at the observer */
 
195
  palPvobs( phi, 0.0, stl, vgo );
 
196
  for (i=0; i<6; i++) {
 
197
    v[i] -= vgo[i];
 
198
  }
 
199
 
 
200
  /* Geometric distance (AU) */
 
201
  dx = v[0];
 
202
  dy = v[1];
 
203
  dz = v[2];
 
204
  r = sqrt( dx*dx + dy*dy + dz*dz );
 
205
 
 
206
  /* Light time */
 
207
  tl = PAL__CR * r;
 
208
 
 
209
  /* Correct position for planetary aberration */
 
210
  for (i=0; i<3; i++) {
 
211
    v[i] -= tl * v[i+3];
 
212
  }
 
213
 
 
214
  /* To RA,Dec */
 
215
  iauC2s( v, ra, dec );
 
216
  *ra = iauAnp( *ra );
 
217
 
 
218
  /* Angular diameter (radians) */
 
219
  *diam = 2.0 * asin( EQRAU[np] / (r * AUKM ) );
 
220
 
 
221
}