Skip to content

Commit 2cfc44d

Browse files
msdemleimsdemlei
committed
+ Basic proper motion propagation
This bumps version to 1.2.1, too.
1 parent ef881cd commit 2cfc44d

File tree

8 files changed

+474
-5
lines changed

8 files changed

+474
-5
lines changed

‎Makefile

+9-5
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@ MODULE_big = pg_sphere
77
OBJS = sscan.o sparse.o sbuffer.o vector3d.o point.o \
88
euler.o circle.o line.o ellipse.o polygon.o \
99
path.o box.o output.o gq_cache.o gist.o key.o \
10-
gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o
10+
gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o \
11+
epochprop.o
1112

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

2829
REGRESS_9_5 = index_9.5 # experimental for spoint3
2930

3031
TESTS = init_test tables points euler circle line ellipse poly path box index \
3132
contains_ops contains_ops_compat bounding_box_gist gnomo healpix \
32-
moc mocautocast
33+
moc mocautocast epochprop
3334

3435
ifndef CXXFLAGS
3536
# no support for CXXFLAGS in PGXS before v11
@@ -39,15 +40,15 @@ endif
3940

4041
EXTRA_CLEAN = $(PGS_SQL) pg_sphere.test.sql
4142

42-
CRUSH_TESTS = init_extended circle_extended
43+
CRUSH_TESTS = init_extended circle_extended
4344

4445
# order of sql files is important
4546
PGS_SQL = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \
4647
pgs_line.sql pgs_ellipse.sql pgs_polygon.sql pgs_path.sql \
4748
pgs_box.sql pgs_contains_ops.sql pgs_contains_ops_compat.sql \
4849
pgs_gist.sql gnomo.sql \
4950
healpix.sql pgs_gist_spoint3.sql pgs_moc_type.sql pgs_moc_compat.sql pgs_moc_ops.sql \
50-
pgs_moc_geo_casts.sql
51+
pgs_moc_geo_casts.sql pgs_epochprop.sql
5152
PGS_SQL_9_5 = pgs_9.5.sql # experimental for spoint3
5253

5354
USE_PGXS = 1
@@ -197,6 +198,9 @@ pg_sphere--1.1.5beta4gavo--1.2.0.sql: pgs_moc_ops.sql.in
197198
cat $^ > $@
198199
ifeq ($(has_parallel), n)
199200
sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL
201+
202+
pg_sphere--1.2.0-1.2.1.sql: pgs_epochprop.sql.in
203+
cat upgrade_scripts/$@.in $^ > $@
200204
endif
201205

202206
pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in

‎epochprop.c

+199
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,199 @@
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+
delta_t = PG_GETARG_FLOAT8(5);
168+
169+
propagate_phasevec(&input, delta_t, &output);
170+
171+
/* change to internal units: rad, rad/yr, mas, and km/s */
172+
retvals[0] = Float8GetDatum(output.pos.lng);
173+
retvals[1] = Float8GetDatum(output.pos.lat);
174+
retvals[2] = Float8GetDatum(output.parallax);
175+
retvals[3] = Float8GetDatum(output.pm[0]);
176+
retvals[4] = Float8GetDatum(output.pm[1]);
177+
retvals[5] = Float8GetDatum(output.rv);
178+
179+
{
180+
bool isnull[6] = {0, 0, 0, 0, 0, 0};
181+
int lower_bounds[1] = {1};
182+
int dims[1] = {6};
183+
#ifdef USE_FLOAT8_BYVAL
184+
bool embyval = true;
185+
#else
186+
bool embyval = false;
187+
#endif
188+
189+
if (! output.parallax_valid) {
190+
/* invalidate parallax and rv */
191+
isnull[2] = 1;
192+
isnull[5] = 1;
193+
}
194+
195+
result = construct_md_array(retvals, isnull, 1, dims, lower_bounds,
196+
FLOAT8OID, sizeof(float8), embyval, 'd');
197+
}
198+
PG_RETURN_ARRAYTYPE_P(result);
199+
}

‎epochprop.h

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
/* Definitions for handling proper motions */
2+
3+
#include <postgres.h>
4+
5+
6+
Datum epoch_prop(PG_FUNCTION_ARGS);
7+
8+
9+
/* a cartesian point; this is like geo_decl's point, but you can't
10+
have both geo_decls and pg_sphere right now (both define a type Point,
11+
not to mention they have different ideas on EPSILON */
12+
typedef struct s_cpoint {
13+
double x, y;
14+
} CPoint;
15+
16+
typedef struct s_phasevec {
17+
SPoint pos; /* Position as an SPoint */
18+
double pm[2]; /* Proper motion long/lat in rad/year,
19+
PM in longitude has cos(lat) applied */
20+
double parallax; /* in rad */
21+
double rv; /* radial velocity in km/s */
22+
int parallax_valid; /* 1 if the parallax really is a NULL */
23+
} phasevec;
24+
25+
26+

‎expected/epochprop.out

+107
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
SELECT
2+
to_char(DEGREES(tp[1]), '999D9999999999'),
3+
to_char(DEGREES(tp[2]), '999D9999999999'),
4+
to_char(tp[3], '999D999'),
5+
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
6+
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
7+
to_char(tp[6], '999D999')
8+
FROM (
9+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
10+
546.9759,
11+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
12+
-100) AS tp) AS q;
13+
to_char | to_char | to_char | to_char | to_char | to_char
14+
-----------------+-----------------+----------+----------+------------+----------
15+
269.4742714391 | 4.4072939987 | 543.624 | -791.442 | 10235.412 | -110.450
16+
(1 row)
17+
18+
SELECT
19+
to_char(DEGREES(tp[1]), '999D9999999999'),
20+
to_char(DEGREES(tp[2]), '999D9999999999'),
21+
to_char(tp[3], '999D999'),
22+
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
23+
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
24+
to_char(tp[6], '999D999')
25+
FROM (
26+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
27+
0,
28+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
29+
-100) AS tp) AS q;
30+
to_char | to_char | to_char | to_char | to_char | to_char
31+
-----------------+-----------------+---------+----------+------------+---------
32+
269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 |
33+
(1 row)
34+
35+
SELECT
36+
to_char(DEGREES(tp[1]), '999D9999999999'),
37+
to_char(DEGREES(tp[2]), '999D9999999999'),
38+
to_char(tp[3], '999D999'),
39+
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
40+
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
41+
to_char(tp[6], '999D999')
42+
FROM (
43+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
44+
NULL,
45+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
46+
-100) AS tp) AS q;
47+
to_char | to_char | to_char | to_char | to_char | to_char
48+
-----------------+-----------------+---------+----------+------------+---------
49+
269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 |
50+
(1 row)
51+
52+
SELECT
53+
to_char(DEGREES(tp[1]), '999D9999999999'),
54+
to_char(DEGREES(tp[2]), '999D9999999999'),
55+
to_char(tp[3], '999D999'),
56+
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
57+
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
58+
to_char(tp[6], '999D999')
59+
FROM (
60+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
61+
23,
62+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL,
63+
20) AS tp) AS q;
64+
to_char | to_char | to_char | to_char | to_char | to_char
65+
-----------------+-----------------+----------+----------+------------+----------
66+
269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 | 2.159
67+
(1 row)
68+
69+
SELECT
70+
to_char(DEGREES(tp[1]), '999D9999999999'),
71+
to_char(DEGREES(tp[2]), '999D9999999999'),
72+
to_char(tp[3], '999D999'),
73+
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
74+
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
75+
to_char(tp[6], '999D999')
76+
FROM (
77+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
78+
23,
79+
NULL, RADIANS(10362/3.6e6), -110,
80+
120) AS tp) AS q;
81+
to_char | to_char | to_char | to_char | to_char | to_char
82+
-----------------+-----------------+----------+----------+------------+----------
83+
269.4520769500 | 5.0388680565 | 23.007 | -.000 | 10368.061 | -97.120
84+
(1 row)
85+
86+
SELECT epoch_prop(NULL,
87+
23,
88+
0.01 , RADIANS(10362/3.6e6), -110,
89+
120);
90+
ERROR: NULL position not supported in epoch propagation
91+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
92+
23,
93+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
94+
20) AS tp;
95+
tp
96+
-----------------------------------------
97+
(4.70274792658313 , 0.0829194509345993)
98+
(1 row)
99+
100+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
101+
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6),
102+
20) AS tp;
103+
tp
104+
-----------------------------------------
105+
(4.70274793061952 , 0.0829193989380876)
106+
(1 row)
107+

0 commit comments

Comments
 (0)