source: branches/FACT++_lidctrl_usb/erfa/src/starpv.c@ 19875

Last change on this file since 19875 was 18711, checked in by tbretz, 10 years ago
Updated to ERFA 1.3.0 (no relevant code change except the leap second at the beginning of 2017)
File size: 10.3 KB
Line 
1#include "erfa.h"
2
3int eraStarpv(double ra, double dec,
4 double pmr, double pmd, double px, double rv,
5 double pv[2][3])
6/*
7** - - - - - - - - - -
8** e r a S t a r p v
9** - - - - - - - - - -
10**
11** Convert star catalog coordinates to position+velocity vector.
12**
13** Given (Note 1):
14** ra double right ascension (radians)
15** dec double declination (radians)
16** pmr double RA proper motion (radians/year)
17** pmd double Dec proper motion (radians/year)
18** px double parallax (arcseconds)
19** rv double radial velocity (km/s, positive = receding)
20**
21** Returned (Note 2):
22** pv double[2][3] pv-vector (AU, AU/day)
23**
24** Returned (function value):
25** int status:
26** 0 = no warnings
27** 1 = distance overridden (Note 6)
28** 2 = excessive speed (Note 7)
29** 4 = solution didn't converge (Note 8)
30** else = binary logical OR of the above
31**
32** Notes:
33**
34** 1) The star data accepted by this function are "observables" for an
35** imaginary observer at the solar-system barycenter. Proper motion
36** and radial velocity are, strictly, in terms of barycentric
37** coordinate time, TCB. For most practical applications, it is
38** permissible to neglect the distinction between TCB and ordinary
39** "proper" time on Earth (TT/TAI). The result will, as a rule, be
40** limited by the intrinsic accuracy of the proper-motion and
41** radial-velocity data; moreover, the pv-vector is likely to be
42** merely an intermediate result, so that a change of time unit
43** would cancel out overall.
44**
45** In accordance with normal star-catalog conventions, the object's
46** right ascension and declination are freed from the effects of
47** secular aberration. The frame, which is aligned to the catalog
48** equator and equinox, is Lorentzian and centered on the SSB.
49**
50** 2) The resulting position and velocity pv-vector is with respect to
51** the same frame and, like the catalog coordinates, is freed from
52** the effects of secular aberration. Should the "coordinate
53** direction", where the object was located at the catalog epoch, be
54** required, it may be obtained by calculating the magnitude of the
55** position vector pv[0][0-2] dividing by the speed of light in
56** AU/day to give the light-time, and then multiplying the space
57** velocity pv[1][0-2] by this light-time and adding the result to
58** pv[0][0-2].
59**
60** Summarizing, the pv-vector returned is for most stars almost
61** identical to the result of applying the standard geometrical
62** "space motion" transformation. The differences, which are the
63** subject of the Stumpff paper referenced below, are:
64**
65** (i) In stars with significant radial velocity and proper motion,
66** the constantly changing light-time distorts the apparent proper
67** motion. Note that this is a classical, not a relativistic,
68** effect.
69**
70** (ii) The transformation complies with special relativity.
71**
72** 3) Care is needed with units. The star coordinates are in radians
73** and the proper motions in radians per Julian year, but the
74** parallax is in arcseconds; the radial velocity is in km/s, but
75** the pv-vector result is in AU and AU/day.
76**
77** 4) The RA proper motion is in terms of coordinate angle, not true
78** angle. If the catalog uses arcseconds for both RA and Dec proper
79** motions, the RA proper motion will need to be divided by cos(Dec)
80** before use.
81**
82** 5) Straight-line motion at constant speed, in the inertial frame,
83** is assumed.
84**
85** 6) An extremely small (or zero or negative) parallax is interpreted
86** to mean that the object is on the "celestial sphere", the radius
87** of which is an arbitrary (large) value (see the constant PXMIN).
88** When the distance is overridden in this way, the status,
89** initially zero, has 1 added to it.
90**
91** 7) If the space velocity is a significant fraction of c (see the
92** constant VMAX), it is arbitrarily set to zero. When this action
93** occurs, 2 is added to the status.
94**
95** 8) The relativistic adjustment involves an iterative calculation.
96** If the process fails to converge within a set number (IMAX) of
97** iterations, 4 is added to the status.
98**
99** 9) The inverse transformation is performed by the function
100** eraPvstar.
101**
102** Called:
103** eraS2pv spherical coordinates to pv-vector
104** eraPm modulus of p-vector
105** eraZp zero p-vector
106** eraPn decompose p-vector into modulus and direction
107** eraPdp scalar product of two p-vectors
108** eraSxp multiply p-vector by scalar
109** eraPmp p-vector minus p-vector
110** eraPpp p-vector plus p-vector
111**
112** Reference:
113**
114** Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.
115**
116** Copyright (C) 2013-2016, NumFOCUS Foundation.
117** Derived, with permission, from the SOFA library. See notes at end of file.
118*/
119{
120/* Smallest allowed parallax */
121 static const double PXMIN = 1e-7;
122
123/* Largest allowed speed (fraction of c) */
124 static const double VMAX = 0.5;
125
126/* Maximum number of iterations for relativistic solution */
127 static const int IMAX = 100;
128
129 int i, iwarn;
130 double w, r, rd, rad, decd, v, x[3], usr[3], ust[3],
131 vsr, vst, betst, betsr, bett, betr,
132 dd, ddel, ur[3], ut[3],
133 d = 0.0, del = 0.0, /* to prevent */
134 odd = 0.0, oddel = 0.0, /* compiler */
135 od = 0.0, odel = 0.0; /* warnings */
136
137
138/* Distance (AU). */
139 if (px >= PXMIN) {
140 w = px;
141 iwarn = 0;
142 } else {
143 w = PXMIN;
144 iwarn = 1;
145 }
146 r = ERFA_DR2AS / w;
147
148/* Radial velocity (AU/day). */
149 rd = ERFA_DAYSEC * rv * 1e3 / ERFA_DAU;
150
151/* Proper motion (radian/day). */
152 rad = pmr / ERFA_DJY;
153 decd = pmd / ERFA_DJY;
154
155/* To pv-vector (AU,AU/day). */
156 eraS2pv(ra, dec, r, rad, decd, rd, pv);
157
158/* If excessive velocity, arbitrarily set it to zero. */
159 v = eraPm(pv[1]);
160 if (v / ERFA_DC > VMAX) {
161 eraZp(pv[1]);
162 iwarn += 2;
163 }
164
165/* Isolate the radial component of the velocity (AU/day). */
166 eraPn(pv[0], &w, x);
167 vsr = eraPdp(x, pv[1]);
168 eraSxp(vsr, x, usr);
169
170/* Isolate the transverse component of the velocity (AU/day). */
171 eraPmp(pv[1], usr, ust);
172 vst = eraPm(ust);
173
174/* Special-relativity dimensionless parameters. */
175 betsr = vsr / ERFA_DC;
176 betst = vst / ERFA_DC;
177
178/* Determine the inertial-to-observed relativistic correction terms. */
179 bett = betst;
180 betr = betsr;
181 for (i = 0; i < IMAX; i++) {
182 d = 1.0 + betr;
183 del = sqrt(1.0 - betr*betr - bett*bett) - 1.0;
184 betr = d * betsr + del;
185 bett = d * betst;
186 if (i > 0) {
187 dd = fabs(d - od);
188 ddel = fabs(del - odel);
189 if ((i > 1) && (dd >= odd) && (ddel >= oddel)) break;
190 odd = dd;
191 oddel = ddel;
192 }
193 od = d;
194 odel = del;
195 }
196 if (i >= IMAX) iwarn += 4;
197
198/* Replace observed radial velocity with inertial value. */
199 w = (betsr != 0.0) ? d + del / betsr : 1.0;
200 eraSxp(w, usr, ur);
201
202/* Replace observed tangential velocity with inertial value. */
203 eraSxp(d, ust, ut);
204
205/* Combine the two to obtain the inertial space velocity. */
206 eraPpp(ur, ut, pv[1]);
207
208/* Return the status. */
209 return iwarn;
210
211}
212/*----------------------------------------------------------------------
213**
214**
215** Copyright (C) 2013-2016, NumFOCUS Foundation.
216** All rights reserved.
217**
218** This library is derived, with permission, from the International
219** Astronomical Union's "Standards of Fundamental Astronomy" library,
220** available from http://www.iausofa.org.
221**
222** The ERFA version is intended to retain identical functionality to
223** the SOFA library, but made distinct through different function and
224** file names, as set out in the SOFA license conditions. The SOFA
225** original has a role as a reference standard for the IAU and IERS,
226** and consequently redistribution is permitted only in its unaltered
227** state. The ERFA version is not subject to this restriction and
228** therefore can be included in distributions which do not support the
229** concept of "read only" software.
230**
231** Although the intent is to replicate the SOFA API (other than
232** replacement of prefix names) and results (with the exception of
233** bugs; any that are discovered will be fixed), SOFA is not
234** responsible for any errors found in this version of the library.
235**
236** If you wish to acknowledge the SOFA heritage, please acknowledge
237** that you are using a library derived from SOFA, rather than SOFA
238** itself.
239**
240**
241** TERMS AND CONDITIONS
242**
243** Redistribution and use in source and binary forms, with or without
244** modification, are permitted provided that the following conditions
245** are met:
246**
247** 1 Redistributions of source code must retain the above copyright
248** notice, this list of conditions and the following disclaimer.
249**
250** 2 Redistributions in binary form must reproduce the above copyright
251** notice, this list of conditions and the following disclaimer in
252** the documentation and/or other materials provided with the
253** distribution.
254**
255** 3 Neither the name of the Standards Of Fundamental Astronomy Board,
256** the International Astronomical Union nor the names of its
257** contributors may be used to endorse or promote products derived
258** from this software without specific prior written permission.
259**
260** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
261** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
262** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
263** FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
264** COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
265** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
266** BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
267** LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
268** CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
269** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
270** ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
271** POSSIBILITY OF SUCH DAMAGE.
272**
273*/
Note: See TracBrowser for help on using the repository browser.