|
| 1 | +/* Handling of (conventional) proper motions. |
| 2 | +
|
| 3 | +This code is largely based on a FORTRAN function written by Lennart Lindegren |
| 4 | +(Lund Obs) in 1995 that implements the procedure described in The Hipparcos |
| 5 | +and Tycho Catalogues (ESA SP-1200), Volume 1, Section 1.5.5, 'Epoch |
| 6 | +Transformation: Rigorous Treatment'; cf. |
| 7 | +<https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf/99adf6e3-6893-4824-8fc2-8d3c9cbba2b5>. |
| 8 | +*/ |
| 9 | + |
| 10 | +#include <math.h> |
| 11 | +#include <pgs_util.h> |
| 12 | + |
| 13 | +#include "point.h" |
| 14 | +#include "epochprop.h" |
| 15 | +#include "vector3d.h" |
| 16 | + |
| 17 | +PG_FUNCTION_INFO_V1(epoch_prop); |
| 18 | + |
| 19 | +/* Astronomical unit in kilometers */ |
| 20 | +#define AU 1.495978701e8 |
| 21 | + |
| 22 | +/* Julian year in seconds */ |
| 23 | +#define J_YEAR (365.25*86400) |
| 24 | + |
| 25 | +/* A_nu as per ESA/SP-1200 */ |
| 26 | +#define A_NU (AU/(J_YEAR)) |
| 27 | + |
| 28 | +/* Following SOFA, we use 1e-7 arcsec as minimal parallax |
| 29 | + ("celestial sphere"); parallax=0 exactly means "infinite distance", which |
| 30 | + leads to all sorts for problems; our parallaxes come in in mas, so: */ |
| 31 | +#define PX_MIN 1e-7*1000 |
| 32 | + |
| 33 | +/* propagate an object at a phase vector over a time difference of delta_t, |
| 34 | +stuffing an updated phase vector in result. |
| 35 | +
|
| 36 | +This does not propagate errors. |
| 37 | +*/ |
| 38 | +static void propagate_phasevec( |
| 39 | + const phasevec *pv, |
| 40 | + const double delta_t, |
| 41 | + phasevec *result) { |
| 42 | + |
| 43 | + double distance_factor, mu0abs, zeta0, parallax; |
| 44 | + |
| 45 | + Vector3D p0, r0, q0; |
| 46 | + Vector3D mu0, pprime, qprime, mu, muprime, u, uprime; |
| 47 | + |
| 48 | + /* for very small or null parallaxes, our algorithm breaks; avoid that |
| 49 | + and, if we did emergency measures, do not talk about parallax and |
| 50 | + radial velocity in the output */ |
| 51 | + if (pv->parallax_valid) { |
| 52 | + parallax = pv->parallax; |
| 53 | + } else { |
| 54 | + parallax = PX_MIN; |
| 55 | + } |
| 56 | + result->parallax_valid = pv->parallax_valid; |
| 57 | + |
| 58 | + /* compute the normal triad as Vector3D-s, eq. (1.2.15)*/ |
| 59 | + spoint_vector3d(&r0, &(pv->pos)); |
| 60 | + |
| 61 | + p0.x = -sin(pv->pos.lng); |
| 62 | + p0.y = cos(pv->pos.lng); |
| 63 | + p0.z = 0; |
| 64 | + |
| 65 | + q0.x = -sin(pv->pos.lat) * cos(pv->pos.lng); |
| 66 | + q0.y = -sin(pv->pos.lat) * sin(pv->pos.lng); |
| 67 | + q0.z = cos(pv->pos.lat); |
| 68 | + |
| 69 | + /* the original proper motion vector */ |
| 70 | + mu0.x = mu0.y = mu0.z = 0; |
| 71 | + vector3d_addwithscalar(&mu0, pv->pm[0], &p0); |
| 72 | + vector3d_addwithscalar(&mu0, pv->pm[1], &q0); |
| 73 | + mu0abs = vector3d_length(&mu0); |
| 74 | + |
| 75 | + /* radial velocity in mas/yr ("change of parallax per year"). eq. (1.5.24) |
| 76 | + We're transforming this to rad/yr so the units work out below */ |
| 77 | + zeta0 = (pv->rv * parallax / A_NU) / 3.6e6 / RADIANS; |
| 78 | + /* distance factor eq. (1.5.25) */ |
| 79 | + distance_factor = 1/sqrt(1 |
| 80 | + + 2 * zeta0 * delta_t |
| 81 | + + (mu0abs * mu0abs + zeta0 * zeta0) * delta_t * delta_t); |
| 82 | + |
| 83 | + /* the propagated proper motion vector, eq. (1.5.28) */ |
| 84 | + muprime.x = muprime.y = muprime.z = 0; |
| 85 | + vector3d_addwithscalar(&muprime, (1 + zeta0 * delta_t), &mu0); |
| 86 | + vector3d_addwithscalar(&muprime, -mu0abs * mu0abs * delta_t, &r0); |
| 87 | + mu.x = mu.y = mu.z = 0; |
| 88 | + vector3d_addwithscalar(&mu, pow(distance_factor, 3), &muprime); |
| 89 | + |
| 90 | + /* parallax, eq. (1.5.27) */ |
| 91 | + result->parallax = distance_factor*parallax; |
| 92 | + /* zeta/rv, eq. (1.5.29); go back from rad to mas, too */ |
| 93 | + result->rv = (zeta0 + (mu0abs * mu0abs + zeta0 * zeta0) * delta_t) |
| 94 | + * distance_factor * distance_factor |
| 95 | + * 3.6e6 * RADIANS |
| 96 | + * A_NU / result->parallax; |
| 97 | + |
| 98 | + /* propagated position, eq. (1.5.26) */ |
| 99 | + uprime.x = uprime.y = uprime.z = 0; |
| 100 | + vector3d_addwithscalar(&uprime, (1 + zeta0 * delta_t), &r0); |
| 101 | + vector3d_addwithscalar(&uprime, delta_t, &mu0); |
| 102 | + u.x = u.y = u.z = 0; |
| 103 | + vector3d_addwithscalar(&u, distance_factor, &uprime); |
| 104 | + vector3d_spoint(&(result->pos), &u); |
| 105 | + |
| 106 | + /* compute a new triad for the propagated position, eq (1.5.15) */ |
| 107 | + pprime.x = -sin(result->pos.lng); |
| 108 | + pprime.y = cos(result->pos.lng); |
| 109 | + pprime.z = 0; |
| 110 | + qprime.x = -sin(result->pos.lat) * cos(result->pos.lng); |
| 111 | + qprime.y = -sin(result->pos.lat) * sin(result->pos.lng); |
| 112 | + qprime.z = cos(result->pos.lat); |
| 113 | + |
| 114 | + /* use it to compute the proper motions, eq. (1.5.32) */ |
| 115 | + result->pm[0] = vector3d_scalar(&pprime, &mu); |
| 116 | + result->pm[1] = vector3d_scalar(&qprime, &mu); |
| 117 | +} |
| 118 | + |
| 119 | + |
| 120 | +/* |
| 121 | + Propagate a position with proper motions and optionally parallax |
| 122 | + and radial velocity. |
| 123 | +
|
| 124 | + Arguments: pos0 (spoint), pm_long, pm_lat (in rad/yr) |
| 125 | + par (parallax, mas), rv (in km/s), delta_t (in years) |
| 126 | +
|
| 127 | + This returns a 6-array of lat, long (in rad), parallax (in mas) |
| 128 | + pmlat, pmlong (in rad/yr), rv (in km/s). |
| 129 | +*/ |
| 130 | +Datum |
| 131 | +epoch_prop(PG_FUNCTION_ARGS) { |
| 132 | + double delta_t; |
| 133 | + phasevec input, output; |
| 134 | + ArrayType *result; |
| 135 | + Datum retvals[6]; |
| 136 | + |
| 137 | + if (PG_ARGISNULL(0)) { |
| 138 | + ereport(ERROR, |
| 139 | + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), |
| 140 | + errmsg("NULL position not supported in epoch propagation"))); } |
| 141 | + memcpy(&(input.pos), (void*)PG_GETARG_POINTER(0), sizeof(SPoint)); |
| 142 | + if (PG_ARGISNULL(1)) { |
| 143 | + input.parallax = 0; |
| 144 | + } else { |
| 145 | + input.parallax = PG_GETARG_FLOAT8(1); |
| 146 | + } |
| 147 | + input.parallax_valid = fabs(input.parallax) > PX_MIN; |
| 148 | + |
| 149 | + if (PG_ARGISNULL(2)) { |
| 150 | + input.pm[0] = 0; |
| 151 | + } else { |
| 152 | + input.pm[0] = PG_GETARG_FLOAT8(2); |
| 153 | + } |
| 154 | + |
| 155 | + if (PG_ARGISNULL(3)) { |
| 156 | + input.pm[1] = 0; |
| 157 | + } else { |
| 158 | + input.pm[1] = PG_GETARG_FLOAT8(3); |
| 159 | + } |
| 160 | + |
| 161 | + if (PG_ARGISNULL(4)) { |
| 162 | + input.rv = 0; |
| 163 | + } else { |
| 164 | + input.rv = PG_GETARG_FLOAT8(4); |
| 165 | + } |
| 166 | + |
| 167 | + if (PG_ARGISNULL(5)) { |
| 168 | + ereport(ERROR, |
| 169 | + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), |
| 170 | + errmsg("NULL delta t not supported in epoch propagation"))); } |
| 171 | + delta_t = PG_GETARG_FLOAT8(5); |
| 172 | + |
| 173 | + propagate_phasevec(&input, delta_t, &output); |
| 174 | + |
| 175 | + /* change to internal units: rad, rad/yr, mas, and km/s */ |
| 176 | + retvals[0] = Float8GetDatum(output.pos.lng); |
| 177 | + retvals[1] = Float8GetDatum(output.pos.lat); |
| 178 | + retvals[2] = Float8GetDatum(output.parallax); |
| 179 | + retvals[3] = Float8GetDatum(output.pm[0]); |
| 180 | + retvals[4] = Float8GetDatum(output.pm[1]); |
| 181 | + retvals[5] = Float8GetDatum(output.rv); |
| 182 | + |
| 183 | + { |
| 184 | + bool isnull[6] = {0, 0, 0, 0, 0, 0}; |
| 185 | + int lower_bounds[1] = {1}; |
| 186 | + int dims[1] = {6}; |
| 187 | +#ifdef USE_FLOAT8_BYVAL |
| 188 | + bool embyval = true; |
| 189 | +#else |
| 190 | + bool embyval = false; |
| 191 | +#endif |
| 192 | + |
| 193 | + if (! output.parallax_valid) { |
| 194 | + /* invalidate parallax and rv */ |
| 195 | + isnull[2] = 1; |
| 196 | + isnull[5] = 1; |
| 197 | + } |
| 198 | + |
| 199 | + result = construct_md_array(retvals, isnull, 1, dims, lower_bounds, |
| 200 | + FLOAT8OID, sizeof(float8), embyval, 'd'); |
| 201 | + } |
| 202 | + PG_RETURN_ARRAYTYPE_P(result); |
| 203 | +} |
0 commit comments