Skip to content

Add code for epoch propagation #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jun 16, 2023
Merged

Add code for epoch propagation #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2cfc44d
+ Basic proper motion propagation
msdemlei Apr 14, 2022
44edb5c
Adding documentation for epoch propagation.
msdemlei Sep 5, 2022
7c3e5a4
Now catching NULLs in the delta-t argument of epoch_prop, too.
msdemlei Jun 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 7 additions & 6 deletions Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@ MODULE_big = pg_sphere
OBJS = sscan.o sparse.o sbuffer.o vector3d.o point.o \
euler.o circle.o line.o ellipse.o polygon.o \
path.o box.o output.o gq_cache.o gist.o key.o \
gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o
gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o \
epochprop.o

EXTENSION = pg_sphere
RELEASE_SQL = $(EXTENSION)--$(PGSPHERE_VERSION).sql
Expand All @@ -23,13 +24,13 @@ DATA_built = $(RELEASE_SQL) \
DOCS = README.pg_sphere COPYRIGHT.pg_sphere
REGRESS = init tables points euler circle line ellipse poly path box index \
contains_ops contains_ops_compat bounding_box_gist gnomo healpix \
moc mocautocast
moc mocautocast epochprop

REGRESS_9_5 = index_9.5 # experimental for spoint3

TESTS = init_test tables points euler circle line ellipse poly path box index \
contains_ops contains_ops_compat bounding_box_gist gnomo healpix \
moc mocautocast
moc mocautocast epochprop

ifndef CXXFLAGS
# no support for CXXFLAGS in PGXS before v11
Expand All @@ -39,15 +40,15 @@ endif

EXTRA_CLEAN = $(PGS_SQL) pg_sphere.test.sql

CRUSH_TESTS = init_extended circle_extended
CRUSH_TESTS = init_extended circle_extended

# order of sql files is important
PGS_SQL = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \
pgs_line.sql pgs_ellipse.sql pgs_polygon.sql pgs_path.sql \
pgs_box.sql pgs_contains_ops.sql pgs_contains_ops_compat.sql \
pgs_gist.sql gnomo.sql \
healpix.sql pgs_gist_spoint3.sql pgs_moc_type.sql pgs_moc_compat.sql pgs_moc_ops.sql \
pgs_moc_geo_casts.sql
pgs_moc_geo_casts.sql pgs_epochprop.sql
PGS_SQL_9_5 = pgs_9.5.sql # experimental for spoint3

USE_PGXS = 1
Expand Down Expand Up @@ -199,7 +200,7 @@ ifeq ($(has_parallel), n)
sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL
endif

pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in
pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in pgs_epochprop.sql.in
cat $^ > $@
ifeq ($(has_parallel), n)
sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL
Expand Down
130 changes: 129 additions & 1 deletion doc/functions.sgm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -667,5 +667,133 @@
</example>

</sect2>


<sect2 id="funcs.epochprop">
<title>
Epoch propagation
</title>

<sect3 id="funcs.epochprop.full">
<title>6-Parameter Epoch Propagation</title>
<funcsynopsis>
<funcprototype>
<funcdef><type>double precision[6]</type>
<function>epoch_prop</function></funcdef>
<paramdef>spoint <parameter>pos</parameter></paramdef>
<paramdef>double precision <parameter>parallax</parameter></paramdef>
<paramdef>double precision <parameter>pm_long</parameter></paramdef>
<paramdef>double precision <parameter>pm_lat</parameter></paramdef>
<paramdef>double precision <parameter>radial_velocity</parameter></paramdef>
<paramdef>double precision <parameter>delta_t</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<para>
Propagates a spherical phase vector in time (in particular,
applies proper motion to positions)
</para>
<para>
Following both pg_sphere and, where missing, astronomical
conventions makes units somewhat eclectic here; pm_long and pm_lat
need to be in rad/yr, whereas parallax is in mas, and
radial_velocity in km/s. The time difference must be in
(Julian) years.
</para>

<para>
This function returns a 6-array of [long, lat, parallax,
pm_long, pm_lat, radial_velocity] of the corresponding values
delta_t years after the reference epoch for the original position.
As in the function arguments, long and lat are in rad, pm_lon and
pm_lat are in rad/yr, parallax is in mas, and radial_velocity is
in km/s. If you are only interested in the position, consider
the epoch_prop_pos functions below that have a somewhat less
contorted signature.
</para>

<para>
It is an error to have either pos or delta_t NULL. For all
other arguments, NULLs are turned into 0s, except for parallax,
where some very small default is put in. In that case,
both parallax and radial_velocity will be NULL in the output
array.
</para>

<para>
This uses the rigorous method derived in "The Hipparcos and Tycho
Catalogues", ESA Special Publication 1200 (1997), p 94f. It does
not take into account relativistic effects, and it also does not
account for secular aberration.
</para>
<example>
<title>Propagating Barnard's star into the past</title>
<programlisting><![CDATA[
SELECT
to_char(DEGREES(tp[1]), '999D9999999999'),
to_char(DEGREES(tp[2]), '999D9999999999'),
to_char(tp[3], '999D999'),
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
to_char(tp[6], '999D999')
FROM (
SELECT epoch_prop(
spoint(radians(269.45207695), radians(4.693364966)), 546.9759,
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
-100) AS tp) AS q;
to_char | to_char | to_char | to_char | to_char | to_char
-----------------+-----------------+----------+----------+------------+----------
269.4742714391 | 4.4072939987 | 543.624 | -791.442 | 10235.412 | -110.450
]]></programlisting>
</example>

</sect3>
<sect3>
<title>Epoch Propagation of Positions Only</title>
<funcsynopsis>
<funcprototype>
<funcdef><type>spoint</type>
<function>epoch_prop_pos</function></funcdef>
<paramdef>spoint <parameter>pos</parameter></paramdef>
<paramdef>double precision <parameter>parallax</parameter></paramdef>
<paramdef>double precision <parameter>pm_long</parameter></paramdef>
<paramdef>double precision <parameter>pm_lat</parameter></paramdef>
<paramdef>double precision <parameter>radial_velocity</parameter></paramdef>
<paramdef>double precision <parameter>delta_t</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<funcsynopsis>
<funcprototype>
<funcdef><type>spoint</type>
<function>epoch_prop_pos</function></funcdef>
<paramdef>spoint <parameter>pos</parameter></paramdef>
<paramdef>double precision <parameter>pm_long</parameter></paramdef>
<paramdef>double precision <parameter>pm_lat</parameter></paramdef>
<paramdef>double precision <parameter>delta_t</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<para>
These are simplified versions of epoch_prop returning only spoints;
the propagated values for the other coordinates are discarded
(but still internallay computed; these functions do not
run any faster than epoch_prop itself).
</para>

<para>
As with epoch_prop itself, missing values (except for pos and
delta_t) are substituted by 0 (or a very small value in the
case of parallax).
</para>
<example>
<title>Barnard's star, position and proper motion</title>
<programlisting><![CDATA[
SELECT epoch_prop_pos(
spoint(radians(269.45207695), radians(4.693364966)),
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6),
20) AS tp;
tp
-----------------------------------------
(4.70274793061952 , 0.0829193989380876)
]]></programlisting>
</example>
</sect3>
</sect2>
</sect1>
203 changes: 203 additions & 0 deletions epochprop.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
/* Handling of (conventional) proper motions.

This code is largely based on a FORTRAN function written by Lennart Lindegren
(Lund Obs) in 1995 that implements the procedure described in The Hipparcos
and Tycho Catalogues (ESA SP-1200), Volume 1, Section 1.5.5, 'Epoch
Transformation: Rigorous Treatment'; cf.
<https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf/99adf6e3-6893-4824-8fc2-8d3c9cbba2b5>.
*/

#include <math.h>
#include <pgs_util.h>

#include "point.h"
#include "epochprop.h"
#include "vector3d.h"

PG_FUNCTION_INFO_V1(epoch_prop);

/* Astronomical unit in kilometers */
#define AU 1.495978701e8

/* Julian year in seconds */
#define J_YEAR (365.25*86400)

/* A_nu as per ESA/SP-1200 */
#define A_NU (AU/(J_YEAR))

/* Following SOFA, we use 1e-7 arcsec as minimal parallax
("celestial sphere"); parallax=0 exactly means "infinite distance", which
leads to all sorts for problems; our parallaxes come in in mas, so: */
#define PX_MIN 1e-7*1000

/* propagate an object at a phase vector over a time difference of delta_t,
stuffing an updated phase vector in result.

This does not propagate errors.
*/
static void propagate_phasevec(
const phasevec *pv,
const double delta_t,
phasevec *result) {

double distance_factor, mu0abs, zeta0, parallax;

Vector3D p0, r0, q0;
Vector3D mu0, pprime, qprime, mu, muprime, u, uprime;

/* for very small or null parallaxes, our algorithm breaks; avoid that
and, if we did emergency measures, do not talk about parallax and
radial velocity in the output */
if (pv->parallax_valid) {
parallax = pv->parallax;
} else {
parallax = PX_MIN;
}
result->parallax_valid = pv->parallax_valid;

/* compute the normal triad as Vector3D-s, eq. (1.2.15)*/
spoint_vector3d(&r0, &(pv->pos));

p0.x = -sin(pv->pos.lng);
p0.y = cos(pv->pos.lng);
p0.z = 0;

q0.x = -sin(pv->pos.lat) * cos(pv->pos.lng);
q0.y = -sin(pv->pos.lat) * sin(pv->pos.lng);
q0.z = cos(pv->pos.lat);

/* the original proper motion vector */
mu0.x = mu0.y = mu0.z = 0;
vector3d_addwithscalar(&mu0, pv->pm[0], &p0);
vector3d_addwithscalar(&mu0, pv->pm[1], &q0);
mu0abs = vector3d_length(&mu0);

/* radial velocity in mas/yr ("change of parallax per year"). eq. (1.5.24)
We're transforming this to rad/yr so the units work out below */
zeta0 = (pv->rv * parallax / A_NU) / 3.6e6 / RADIANS;
/* distance factor eq. (1.5.25) */
distance_factor = 1/sqrt(1
+ 2 * zeta0 * delta_t
+ (mu0abs * mu0abs + zeta0 * zeta0) * delta_t * delta_t);

/* the propagated proper motion vector, eq. (1.5.28) */
muprime.x = muprime.y = muprime.z = 0;
vector3d_addwithscalar(&muprime, (1 + zeta0 * delta_t), &mu0);
vector3d_addwithscalar(&muprime, -mu0abs * mu0abs * delta_t, &r0);
mu.x = mu.y = mu.z = 0;
vector3d_addwithscalar(&mu, pow(distance_factor, 3), &muprime);

/* parallax, eq. (1.5.27) */
result->parallax = distance_factor*parallax;
/* zeta/rv, eq. (1.5.29); go back from rad to mas, too */
result->rv = (zeta0 + (mu0abs * mu0abs + zeta0 * zeta0) * delta_t)
* distance_factor * distance_factor
* 3.6e6 * RADIANS
* A_NU / result->parallax;

/* propagated position, eq. (1.5.26) */
uprime.x = uprime.y = uprime.z = 0;
vector3d_addwithscalar(&uprime, (1 + zeta0 * delta_t), &r0);
vector3d_addwithscalar(&uprime, delta_t, &mu0);
u.x = u.y = u.z = 0;
vector3d_addwithscalar(&u, distance_factor, &uprime);
vector3d_spoint(&(result->pos), &u);

/* compute a new triad for the propagated position, eq (1.5.15) */
pprime.x = -sin(result->pos.lng);
pprime.y = cos(result->pos.lng);
pprime.z = 0;
qprime.x = -sin(result->pos.lat) * cos(result->pos.lng);
qprime.y = -sin(result->pos.lat) * sin(result->pos.lng);
qprime.z = cos(result->pos.lat);

/* use it to compute the proper motions, eq. (1.5.32) */
result->pm[0] = vector3d_scalar(&pprime, &mu);
result->pm[1] = vector3d_scalar(&qprime, &mu);
}


/*
Propagate a position with proper motions and optionally parallax
and radial velocity.

Arguments: pos0 (spoint), pm_long, pm_lat (in rad/yr)
par (parallax, mas), rv (in km/s), delta_t (in years)

This returns a 6-array of lat, long (in rad), parallax (in mas)
pmlat, pmlong (in rad/yr), rv (in km/s).
*/
Datum
epoch_prop(PG_FUNCTION_ARGS) {
double delta_t;
phasevec input, output;
ArrayType *result;
Datum retvals[6];

if (PG_ARGISNULL(0)) {
ereport(ERROR,
(errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED),
errmsg("NULL position not supported in epoch propagation"))); }
memcpy(&(input.pos), (void*)PG_GETARG_POINTER(0), sizeof(SPoint));
if (PG_ARGISNULL(1)) {
input.parallax = 0;
} else {
input.parallax = PG_GETARG_FLOAT8(1);
}
input.parallax_valid = fabs(input.parallax) > PX_MIN;

if (PG_ARGISNULL(2)) {
input.pm[0] = 0;
} else {
input.pm[0] = PG_GETARG_FLOAT8(2);
}

if (PG_ARGISNULL(3)) {
input.pm[1] = 0;
} else {
input.pm[1] = PG_GETARG_FLOAT8(3);
}

if (PG_ARGISNULL(4)) {
input.rv = 0;
} else {
input.rv = PG_GETARG_FLOAT8(4);
}

if (PG_ARGISNULL(5)) {
ereport(ERROR,
(errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED),
errmsg("NULL delta t not supported in epoch propagation"))); }
delta_t = PG_GETARG_FLOAT8(5);

propagate_phasevec(&input, delta_t, &output);

/* change to internal units: rad, rad/yr, mas, and km/s */
retvals[0] = Float8GetDatum(output.pos.lng);
retvals[1] = Float8GetDatum(output.pos.lat);
retvals[2] = Float8GetDatum(output.parallax);
retvals[3] = Float8GetDatum(output.pm[0]);
retvals[4] = Float8GetDatum(output.pm[1]);
retvals[5] = Float8GetDatum(output.rv);

{
bool isnull[6] = {0, 0, 0, 0, 0, 0};
int lower_bounds[1] = {1};
int dims[1] = {6};
#ifdef USE_FLOAT8_BYVAL
bool embyval = true;
#else
bool embyval = false;
#endif

if (! output.parallax_valid) {
/* invalidate parallax and rv */
isnull[2] = 1;
isnull[5] = 1;
}

result = construct_md_array(retvals, isnull, 1, dims, lower_bounds,
FLOAT8OID, sizeof(float8), embyval, 'd');
}
PG_RETURN_ARRAYTYPE_P(result);
}
Loading