| 1 | #include "erfa.h"
|
|---|
| 2 |
|
|---|
| 3 | int 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-2017, 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 | w = betr*betr + bett*bett;
|
|---|
| 184 | del = - w / (sqrt(1.0 - w) + 1.0);
|
|---|
| 185 | betr = d * betsr + del;
|
|---|
| 186 | bett = d * betst;
|
|---|
| 187 | if (i > 0) {
|
|---|
| 188 | dd = fabs(d - od);
|
|---|
| 189 | ddel = fabs(del - odel);
|
|---|
| 190 | if ((i > 1) && (dd >= odd) && (ddel >= oddel)) break;
|
|---|
| 191 | odd = dd;
|
|---|
| 192 | oddel = ddel;
|
|---|
| 193 | }
|
|---|
| 194 | od = d;
|
|---|
| 195 | odel = del;
|
|---|
| 196 | }
|
|---|
| 197 | if (i >= IMAX) iwarn += 4;
|
|---|
| 198 |
|
|---|
| 199 | /* Replace observed radial velocity with inertial value. */
|
|---|
| 200 | w = (betsr != 0.0) ? d + del / betsr : 1.0;
|
|---|
| 201 | eraSxp(w, usr, ur);
|
|---|
| 202 |
|
|---|
| 203 | /* Replace observed tangential velocity with inertial value. */
|
|---|
| 204 | eraSxp(d, ust, ut);
|
|---|
| 205 |
|
|---|
| 206 | /* Combine the two to obtain the inertial space velocity. */
|
|---|
| 207 | eraPpp(ur, ut, pv[1]);
|
|---|
| 208 |
|
|---|
| 209 | /* Return the status. */
|
|---|
| 210 | return iwarn;
|
|---|
| 211 |
|
|---|
| 212 | }
|
|---|
| 213 | /*----------------------------------------------------------------------
|
|---|
| 214 | **
|
|---|
| 215 | **
|
|---|
| 216 | ** Copyright (C) 2013-2017, NumFOCUS Foundation.
|
|---|
| 217 | ** All rights reserved.
|
|---|
| 218 | **
|
|---|
| 219 | ** This library is derived, with permission, from the International
|
|---|
| 220 | ** Astronomical Union's "Standards of Fundamental Astronomy" library,
|
|---|
| 221 | ** available from http://www.iausofa.org.
|
|---|
| 222 | **
|
|---|
| 223 | ** The ERFA version is intended to retain identical functionality to
|
|---|
| 224 | ** the SOFA library, but made distinct through different function and
|
|---|
| 225 | ** file names, as set out in the SOFA license conditions. The SOFA
|
|---|
| 226 | ** original has a role as a reference standard for the IAU and IERS,
|
|---|
| 227 | ** and consequently redistribution is permitted only in its unaltered
|
|---|
| 228 | ** state. The ERFA version is not subject to this restriction and
|
|---|
| 229 | ** therefore can be included in distributions which do not support the
|
|---|
| 230 | ** concept of "read only" software.
|
|---|
| 231 | **
|
|---|
| 232 | ** Although the intent is to replicate the SOFA API (other than
|
|---|
| 233 | ** replacement of prefix names) and results (with the exception of
|
|---|
| 234 | ** bugs; any that are discovered will be fixed), SOFA is not
|
|---|
| 235 | ** responsible for any errors found in this version of the library.
|
|---|
| 236 | **
|
|---|
| 237 | ** If you wish to acknowledge the SOFA heritage, please acknowledge
|
|---|
| 238 | ** that you are using a library derived from SOFA, rather than SOFA
|
|---|
| 239 | ** itself.
|
|---|
| 240 | **
|
|---|
| 241 | **
|
|---|
| 242 | ** TERMS AND CONDITIONS
|
|---|
| 243 | **
|
|---|
| 244 | ** Redistribution and use in source and binary forms, with or without
|
|---|
| 245 | ** modification, are permitted provided that the following conditions
|
|---|
| 246 | ** are met:
|
|---|
| 247 | **
|
|---|
| 248 | ** 1 Redistributions of source code must retain the above copyright
|
|---|
| 249 | ** notice, this list of conditions and the following disclaimer.
|
|---|
| 250 | **
|
|---|
| 251 | ** 2 Redistributions in binary form must reproduce the above copyright
|
|---|
| 252 | ** notice, this list of conditions and the following disclaimer in
|
|---|
| 253 | ** the documentation and/or other materials provided with the
|
|---|
| 254 | ** distribution.
|
|---|
| 255 | **
|
|---|
| 256 | ** 3 Neither the name of the Standards Of Fundamental Astronomy Board,
|
|---|
| 257 | ** the International Astronomical Union nor the names of its
|
|---|
| 258 | ** contributors may be used to endorse or promote products derived
|
|---|
| 259 | ** from this software without specific prior written permission.
|
|---|
| 260 | **
|
|---|
| 261 | ** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|---|
| 262 | ** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|---|
| 263 | ** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
|
|---|
| 264 | ** FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
|
|---|
| 265 | ** COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
|
|---|
| 266 | ** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
|
|---|
| 267 | ** BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
|
|---|
| 268 | ** LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
|
|---|
| 269 | ** CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|---|
| 270 | ** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
|
|---|
| 271 | ** ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|---|
| 272 | ** POSSIBILITY OF SUCH DAMAGE.
|
|---|
| 273 | **
|
|---|
| 274 | */
|
|---|
