1 #include "sofa.h"
2 #include "sofam.h"
3
iauStarpv(double ra,double dec,double pmr,double pmd,double px,double rv,double pv[2][3])4 int iauStarpv(double ra, double dec,
5 double pmr, double pmd, double px, double rv,
6 double pv[2][3])
7 /*
8 ** - - - - - - - - - -
9 ** i a u S t a r p v
10 ** - - - - - - - - - -
11 **
12 ** Convert star catalog coordinates to position+velocity vector.
13 **
14 ** This function is part of the International Astronomical Union's
15 ** SOFA (Standards Of Fundamental Astronomy) software collection.
16 **
17 ** Status: support function.
18 **
19 ** Given (Note 1):
20 ** ra double right ascension (radians)
21 ** dec double declination (radians)
22 ** pmr double RA proper motion (radians/year)
23 ** pmd double Dec proper motion (radians/year)
24 ** px double parallax (arcseconds)
25 ** rv double radial velocity (km/s, positive = receding)
26 **
27 ** Returned (Note 2):
28 ** pv double[2][3] pv-vector (au, au/day)
29 **
30 ** Returned (function value):
31 ** int status:
32 ** 0 = no warnings
33 ** 1 = distance overridden (Note 6)
34 ** 2 = excessive speed (Note 7)
35 ** 4 = solution didn't converge (Note 8)
36 ** else = binary logical OR of the above
37 **
38 ** Notes:
39 **
40 ** 1) The star data accepted by this function are "observables" for an
41 ** imaginary observer at the solar-system barycenter. Proper motion
42 ** and radial velocity are, strictly, in terms of barycentric
43 ** coordinate time, TCB. For most practical applications, it is
44 ** permissible to neglect the distinction between TCB and ordinary
45 ** "proper" time on Earth (TT/TAI). The result will, as a rule, be
46 ** limited by the intrinsic accuracy of the proper-motion and
47 ** radial-velocity data; moreover, the pv-vector is likely to be
48 ** merely an intermediate result, so that a change of time unit
49 ** would cancel out overall.
50 **
51 ** In accordance with normal star-catalog conventions, the object's
52 ** right ascension and declination are freed from the effects of
53 ** secular aberration. The frame, which is aligned to the catalog
54 ** equator and equinox, is Lorentzian and centered on the SSB.
55 **
56 ** 2) The resulting position and velocity pv-vector is with respect to
57 ** the same frame and, like the catalog coordinates, is freed from
58 ** the effects of secular aberration. Should the "coordinate
59 ** direction", where the object was located at the catalog epoch, be
60 ** required, it may be obtained by calculating the magnitude of the
61 ** position vector pv[0][0-2] dividing by the speed of light in
62 ** au/day to give the light-time, and then multiplying the space
63 ** velocity pv[1][0-2] by this light-time and adding the result to
64 ** pv[0][0-2].
65 **
66 ** Summarizing, the pv-vector returned is for most stars almost
67 ** identical to the result of applying the standard geometrical
68 ** "space motion" transformation. The differences, which are the
69 ** subject of the Stumpff paper referenced below, are:
70 **
71 ** (i) In stars with significant radial velocity and proper motion,
72 ** the constantly changing light-time distorts the apparent proper
73 ** motion. Note that this is a classical, not a relativistic,
74 ** effect.
75 **
76 ** (ii) The transformation complies with special relativity.
77 **
78 ** 3) Care is needed with units. The star coordinates are in radians
79 ** and the proper motions in radians per Julian year, but the
80 ** parallax is in arcseconds; the radial velocity is in km/s, but
81 ** the pv-vector result is in au and au/day.
82 **
83 ** 4) The RA proper motion is in terms of coordinate angle, not true
84 ** angle. If the catalog uses arcseconds for both RA and Dec proper
85 ** motions, the RA proper motion will need to be divided by cos(Dec)
86 ** before use.
87 **
88 ** 5) Straight-line motion at constant speed, in the inertial frame,
89 ** is assumed.
90 **
91 ** 6) An extremely small (or zero or negative) parallax is interpreted
92 ** to mean that the object is on the "celestial sphere", the radius
93 ** of which is an arbitrary (large) value (see the constant PXMIN).
94 ** When the distance is overridden in this way, the status,
95 ** initially zero, has 1 added to it.
96 **
97 ** 7) If the space velocity is a significant fraction of c (see the
98 ** constant VMAX), it is arbitrarily set to zero. When this action
99 ** occurs, 2 is added to the status.
100 **
101 ** 8) The relativistic adjustment involves an iterative calculation.
102 ** If the process fails to converge within a set number (IMAX) of
103 ** iterations, 4 is added to the status.
104 **
105 ** 9) The inverse transformation is performed by the function
106 ** iauPvstar.
107 **
108 ** Called:
109 ** iauS2pv spherical coordinates to pv-vector
110 ** iauPm modulus of p-vector
111 ** iauZp zero p-vector
112 ** iauPn decompose p-vector into modulus and direction
113 ** iauPdp scalar product of two p-vectors
114 ** iauSxp multiply p-vector by scalar
115 ** iauPmp p-vector minus p-vector
116 ** iauPpp p-vector plus p-vector
117 **
118 ** Reference:
119 **
120 ** Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.
121 **
122 ** This revision: 2021 May 11
123 **
124 ** SOFA release 2021-05-12
125 **
126 ** Copyright (C) 2021 IAU SOFA Board. See notes at end.
127 */
128 {
129 /* Smallest allowed parallax */
130 static const double PXMIN = 1e-7;
131
132 /* Largest allowed speed (fraction of c) */
133 static const double VMAX = 0.5;
134
135 /* Maximum number of iterations for relativistic solution */
136 static const int IMAX = 100;
137
138 int i, iwarn;
139 double w, r, rd, rad, decd, v, x[3], usr[3], ust[3],
140 vsr, vst, betst, betsr, bett, betr,
141 dd, ddel, ur[3], ut[3],
142 d = 0.0, del = 0.0, /* to prevent */
143 odd = 0.0, oddel = 0.0, /* compiler */
144 od = 0.0, odel = 0.0; /* warnings */
145
146
147 /* Distance (au). */
148 if (px >= PXMIN) {
149 w = px;
150 iwarn = 0;
151 } else {
152 w = PXMIN;
153 iwarn = 1;
154 }
155 r = DR2AS / w;
156
157 /* Radial velocity (au/day). */
158 rd = DAYSEC * rv * 1e3 / DAU;
159
160 /* Proper motion (radian/day). */
161 rad = pmr / DJY;
162 decd = pmd / DJY;
163
164 /* To pv-vector (au,au/day). */
165 iauS2pv(ra, dec, r, rad, decd, rd, pv);
166
167 /* If excessive velocity, arbitrarily set it to zero. */
168 v = iauPm(pv[1]);
169 if (v / DC > VMAX) {
170 iauZp(pv[1]);
171 iwarn += 2;
172 }
173
174 /* Isolate the radial component of the velocity (au/day). */
175 iauPn(pv[0], &w, x);
176 vsr = iauPdp(x, pv[1]);
177 iauSxp(vsr, x, usr);
178
179 /* Isolate the transverse component of the velocity (au/day). */
180 iauPmp(pv[1], usr, ust);
181 vst = iauPm(ust);
182
183 /* Special-relativity dimensionless parameters. */
184 betsr = vsr / DC;
185 betst = vst / DC;
186
187 /* Determine the inertial-to-observed relativistic correction terms. */
188 bett = betst;
189 betr = betsr;
190 for (i = 0; i < IMAX; i++) {
191 d = 1.0 + betr;
192 w = betr*betr + bett*bett;
193 del = - w / (sqrt(1.0 - w) + 1.0);
194 betr = d * betsr + del;
195 bett = d * betst;
196 if (i > 0) {
197 dd = fabs(d - od);
198 ddel = fabs(del - odel);
199 if ((i > 1) && (dd >= odd) && (ddel >= oddel)) break;
200 odd = dd;
201 oddel = ddel;
202 }
203 od = d;
204 odel = del;
205 }
206 if (i >= IMAX) iwarn += 4;
207
208 /* Replace observed radial velocity with inertial value. */
209 w = (betsr != 0.0) ? d + del / betsr : 1.0;
210 iauSxp(w, usr, ur);
211
212 /* Replace observed tangential velocity with inertial value. */
213 iauSxp(d, ust, ut);
214
215 /* Combine the two to obtain the inertial space velocity. */
216 iauPpp(ur, ut, pv[1]);
217
218 /* Return the status. */
219 return iwarn;
220
221 /* Finished. */
222
223 /*----------------------------------------------------------------------
224 **
225 ** Copyright (C) 2021
226 ** Standards Of Fundamental Astronomy Board
227 ** of the International Astronomical Union.
228 **
229 ** =====================
230 ** SOFA Software License
231 ** =====================
232 **
233 ** NOTICE TO USER:
234 **
235 ** BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
236 ** CONDITIONS WHICH APPLY TO ITS USE.
237 **
238 ** 1. The Software is owned by the IAU SOFA Board ("SOFA").
239 **
240 ** 2. Permission is granted to anyone to use the SOFA software for any
241 ** purpose, including commercial applications, free of charge and
242 ** without payment of royalties, subject to the conditions and
243 ** restrictions listed below.
244 **
245 ** 3. You (the user) may copy and distribute SOFA source code to others,
246 ** and use and adapt its code and algorithms in your own software,
247 ** on a world-wide, royalty-free basis. That portion of your
248 ** distribution that does not consist of intact and unchanged copies
249 ** of SOFA source code files is a "derived work" that must comply
250 ** with the following requirements:
251 **
252 ** a) Your work shall be marked or carry a statement that it
253 ** (i) uses routines and computations derived by you from
254 ** software provided by SOFA under license to you; and
255 ** (ii) does not itself constitute software provided by and/or
256 ** endorsed by SOFA.
257 **
258 ** b) The source code of your derived work must contain descriptions
259 ** of how the derived work is based upon, contains and/or differs
260 ** from the original SOFA software.
261 **
262 ** c) The names of all routines in your derived work shall not
263 ** include the prefix "iau" or "sofa" or trivial modifications
264 ** thereof such as changes of case.
265 **
266 ** d) The origin of the SOFA components of your derived work must
267 ** not be misrepresented; you must not claim that you wrote the
268 ** original software, nor file a patent application for SOFA
269 ** software or algorithms embedded in the SOFA software.
270 **
271 ** e) These requirements must be reproduced intact in any source
272 ** distribution and shall apply to anyone to whom you have
273 ** granted a further right to modify the source code of your
274 ** derived work.
275 **
276 ** Note that, as originally distributed, the SOFA software is
277 ** intended to be a definitive implementation of the IAU standards,
278 ** and consequently third-party modifications are discouraged. All
279 ** variations, no matter how minor, must be explicitly marked as
280 ** such, as explained above.
281 **
282 ** 4. You shall not cause the SOFA software to be brought into
283 ** disrepute, either by misuse, or use for inappropriate tasks, or
284 ** by inappropriate modification.
285 **
286 ** 5. The SOFA software is provided "as is" and SOFA makes no warranty
287 ** as to its use or performance. SOFA does not and cannot warrant
288 ** the performance or results which the user may obtain by using the
289 ** SOFA software. SOFA makes no warranties, express or implied, as
290 ** to non-infringement of third party rights, merchantability, or
291 ** fitness for any particular purpose. In no event will SOFA be
292 ** liable to the user for any consequential, incidental, or special
293 ** damages, including any lost profits or lost savings, even if a
294 ** SOFA representative has been advised of such damages, or for any
295 ** claim by any third party.
296 **
297 ** 6. The provision of any version of the SOFA software under the terms
298 ** and conditions specified herein does not imply that future
299 ** versions will also be made available under the same terms and
300 ** conditions.
301 *
302 ** In any published work or commercial product which uses the SOFA
303 ** software directly, acknowledgement (see www.iausofa.org) is
304 ** appreciated.
305 **
306 ** Correspondence concerning SOFA software should be addressed as
307 ** follows:
308 **
309 ** By email: sofa@ukho.gov.uk
310 ** By post: IAU SOFA Center
311 ** HM Nautical Almanac Office
312 ** UK Hydrographic Office
313 ** Admiralty Way, Taunton
314 ** Somerset, TA1 2DN
315 ** United Kingdom
316 **
317 **--------------------------------------------------------------------*/
318 }
319