Skip to content

Commit e1dee24

Browse files
msdemleivitcpp
msdemlei
authored and
vitcpp
committed
Now more carefully nulling out null inputs in epoch_prop.
This is so we do not invent values for pm, parallax, or rv. We are actually a bit over-cautious and invalidate both PMs even if only one is missing. This is because PMs mix, and hence there are traces of the invented 0 in the other component of the PM. Similarly, when the parallax is missing or bad, the RV would be heavily contaminated; in these cases, we copy through the original RV, as it may still be useful and certainly will not be grossly wrong. (sorry about a few whitespace diffs introduced by killing trailing whitespace in sql/epochprop.sql and src/epochprop.c)
1 parent 0607488 commit e1dee24

File tree

5 files changed

+73
-72
lines changed

5 files changed

+73
-72
lines changed

‎doc/functions.sgm

+6-3
Original file line numberDiff line numberDiff line change
@@ -867,9 +867,12 @@
867867
<para>
868868
It is an error to have either pos or delta_t NULL. For all
869869
other arguments, NULLs are turned into 0s, except for parallax,
870-
where some very small default is put in. In that case,
871-
both parallax and radial_velocity will be NULL in the output
872-
array.
870+
where some very small default is put in. Whatever is NULL
871+
on the input is NULL on the output. In addition, we null
872+
out a non-NULL input on one component of the proper motion
873+
if the other component is NULL, and we null out the radial
874+
velocity if the parallax is missing, as it would be horribly
875+
off with the propagation algorithm we use here.
873876
</para>
874877

875878
<para>

‎expected/epochprop.out

+31-31
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
1-
SET extra_float_digits = 2;
2-
SELECT
1+
SET extra_float_digits = 1;
2+
SELECT
33
to_char(DEGREES(tp[1]), '999D9999999999'),
44
to_char(DEGREES(tp[2]), '999D9999999999'),
55
to_char(tp[3], '999D999'),
66
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
77
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
88
to_char(tp[6], '999D999')
99
FROM (
10-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
10+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
1111
546.9759,
1212
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
1313
-100) AS tp) AS q;
@@ -16,93 +16,93 @@ FROM (
1616
269.4742714391 | 4.4072939987 | 543.624 | -791.442 | 10235.412 | -110.450
1717
(1 row)
1818

19-
SELECT
19+
SELECT
2020
to_char(DEGREES(tp[1]), '999D9999999999'),
2121
to_char(DEGREES(tp[2]), '999D9999999999'),
2222
to_char(tp[3], '999D999'),
2323
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
2424
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
2525
to_char(tp[6], '999D999')
2626
FROM (
27-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
27+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
2828
0,
2929
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
3030
-100) AS tp) AS q;
31-
to_char | to_char | to_char | to_char | to_char | to_char
32-
-----------------+-----------------+---------+----------+------------+---------
33-
269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 |
31+
to_char | to_char | to_char | to_char | to_char | to_char
32+
-----------------+-----------------+----------+----------+------------+----------
33+
269.4744079540 | 4.4055337210 | .000 | -801.210 | 10361.762 | -110.000
3434
(1 row)
3535

36-
SELECT
36+
SELECT
3737
to_char(DEGREES(tp[1]), '999D9999999999'),
3838
to_char(DEGREES(tp[2]), '999D9999999999'),
3939
to_char(tp[3], '999D999'),
4040
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
4141
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
4242
to_char(tp[6], '999D999')
4343
FROM (
44-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
44+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
4545
NULL,
4646
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
4747
-100) AS tp) AS q;
48-
to_char | to_char | to_char | to_char | to_char | to_char
49-
-----------------+-----------------+---------+----------+------------+---------
50-
269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 |
48+
to_char | to_char | to_char | to_char | to_char | to_char
49+
-----------------+-----------------+---------+----------+------------+----------
50+
269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | -110.000
5151
(1 row)
5252

53-
SELECT
53+
SELECT
5454
to_char(DEGREES(tp[1]), '999D9999999999'),
5555
to_char(DEGREES(tp[2]), '999D9999999999'),
5656
to_char(tp[3], '999D999'),
5757
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
5858
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
5959
to_char(tp[6], '999D999')
6060
FROM (
61-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
61+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
6262
23,
6363
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL,
6464
20) AS tp) AS q;
65-
to_char | to_char | to_char | to_char | to_char | to_char
66-
-----------------+-----------------+----------+----------+------------+----------
67-
269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 | 2.159
65+
to_char | to_char | to_char | to_char | to_char | to_char
66+
-----------------+-----------------+----------+----------+------------+---------
67+
269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 |
6868
(1 row)
6969

70-
SELECT
70+
SELECT
7171
to_char(DEGREES(tp[1]), '999D9999999999'),
7272
to_char(DEGREES(tp[2]), '999D9999999999'),
7373
to_char(tp[3], '999D999'),
7474
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
7575
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
7676
to_char(tp[6], '999D999')
7777
FROM (
78-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
78+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
7979
23,
8080
NULL, RADIANS(10362/3.6e6), -110,
8181
120) AS tp) AS q;
82-
to_char | to_char | to_char | to_char | to_char | to_char
83-
-----------------+-----------------+----------+----------+------------+----------
84-
269.4520769500 | 5.0388680565 | 23.007 | -.000 | 10368.061 | -97.120
82+
to_char | to_char | to_char | to_char | to_char | to_char
83+
-----------------+-----------------+----------+---------+---------+----------
84+
269.4520769500 | 4.6933649660 | 23.007 | | | -110.000
8585
(1 row)
8686

8787
SELECT epoch_prop(NULL,
8888
23,
8989
0.01 , RADIANS(10362/3.6e6), -110,
9090
120);
9191
ERROR: NULL position not supported in epoch propagation
92-
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
92+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
9393
23,
9494
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
9595
20) AS tp;
96-
tp
97-
---------------------------------------------
98-
(4.7027479265831289 , 0.082919450934599334)
96+
tp
97+
-------------------------------------------
98+
(4.702747926583129 , 0.08291945093459933)
9999
(1 row)
100100

101-
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
101+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
102102
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6),
103103
20) AS tp;
104-
tp
105-
---------------------------------------------
106-
(4.7027479306195161 , 0.082919398938087627)
104+
tp
105+
-------------------------------------------
106+
(4.702747930619516 , 0.08291939893808763)
107107
(1 row)
108108

‎sql/epochprop.sql

+13-13
Original file line numberDiff line numberDiff line change
@@ -1,66 +1,66 @@
1-
SET extra_float_digits = 2;
1+
SET extra_float_digits = 1;
22

3-
SELECT
3+
SELECT
44
to_char(DEGREES(tp[1]), '999D9999999999'),
55
to_char(DEGREES(tp[2]), '999D9999999999'),
66
to_char(tp[3], '999D999'),
77
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
88
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
99
to_char(tp[6], '999D999')
1010
FROM (
11-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
11+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
1212
546.9759,
1313
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
1414
-100) AS tp) AS q;
1515

16-
SELECT
16+
SELECT
1717
to_char(DEGREES(tp[1]), '999D9999999999'),
1818
to_char(DEGREES(tp[2]), '999D9999999999'),
1919
to_char(tp[3], '999D999'),
2020
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
2121
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
2222
to_char(tp[6], '999D999')
2323
FROM (
24-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
24+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
2525
0,
2626
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
2727
-100) AS tp) AS q;
2828

29-
SELECT
29+
SELECT
3030
to_char(DEGREES(tp[1]), '999D9999999999'),
3131
to_char(DEGREES(tp[2]), '999D9999999999'),
3232
to_char(tp[3], '999D999'),
3333
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
3434
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
3535
to_char(tp[6], '999D999')
3636
FROM (
37-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
37+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
3838
NULL,
3939
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
4040
-100) AS tp) AS q;
4141

42-
SELECT
42+
SELECT
4343
to_char(DEGREES(tp[1]), '999D9999999999'),
4444
to_char(DEGREES(tp[2]), '999D9999999999'),
4545
to_char(tp[3], '999D999'),
4646
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
4747
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
4848
to_char(tp[6], '999D999')
4949
FROM (
50-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
50+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
5151
23,
5252
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL,
5353
20) AS tp) AS q;
5454

55-
SELECT
55+
SELECT
5656
to_char(DEGREES(tp[1]), '999D9999999999'),
5757
to_char(DEGREES(tp[2]), '999D9999999999'),
5858
to_char(tp[3], '999D999'),
5959
to_char(DEGREES(tp[4])*3.6e6, '999D999'),
6060
to_char(DEGREES(tp[5])*3.6e6, '99999D999'),
6161
to_char(tp[6], '999D999')
6262
FROM (
63-
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
63+
SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)),
6464
23,
6565
NULL, RADIANS(10362/3.6e6), -110,
6666
120) AS tp) AS q;
@@ -70,11 +70,11 @@ SELECT epoch_prop(NULL,
7070
0.01 , RADIANS(10362/3.6e6), -110,
7171
120);
7272

73-
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
73+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
7474
23,
7575
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110,
7676
20) AS tp;
7777

78-
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
78+
SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)),
7979
RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6),
8080
20) AS tp;

‎src/epochprop.c

+22-15
Original file line numberDiff line numberDiff line change
@@ -133,6 +133,7 @@ epoch_prop(PG_FUNCTION_ARGS) {
133133
phasevec input, output;
134134
ArrayType *result;
135135
Datum retvals[6];
136+
bool output_null[6] = {0, 0, 0, 0, 0, 0};
136137

137138
if (PG_ARGISNULL(0)) {
138139
ereport(ERROR,
@@ -141,25 +142,29 @@ epoch_prop(PG_FUNCTION_ARGS) {
141142
memcpy(&(input.pos), (void*)PG_GETARG_POINTER(0), sizeof(SPoint));
142143
if (PG_ARGISNULL(1)) {
143144
input.parallax = 0;
145+
output_null[2] = 1;
146+
/* The way we do our computation, with a bad parallax the RV
147+
will be horribly off, too, so null this out, too; if avaialble,
148+
we will fiddle in the original RV below again. */
149+
output_null[5] = 1;
144150
} else {
145151
input.parallax = PG_GETARG_FLOAT8(1);
146152
}
147153
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-
}
154154

155-
if (PG_ARGISNULL(3)) {
155+
if (PG_ARGISNULL(2) || PG_ARGISNULL(3)) {
156+
input.pm[0] = 0;
156157
input.pm[1] = 0;
158+
output_null[3] = 1;
159+
output_null[4] = 1;
157160
} else {
161+
input.pm[0] = PG_GETARG_FLOAT8(2);
158162
input.pm[1] = PG_GETARG_FLOAT8(3);
159163
}
160164

161165
if (PG_ARGISNULL(4)) {
162166
input.rv = 0;
167+
output_null[5] = 1;
163168
} else {
164169
input.rv = PG_GETARG_FLOAT8(4);
165170
}
@@ -172,6 +177,15 @@ epoch_prop(PG_FUNCTION_ARGS) {
172177

173178
propagate_phasevec(&input, delta_t, &output);
174179

180+
/* If we have an invalid parallax but a good RV, preserve the original,
181+
untransformed RV on output. See
182+
https://github.com/ivoa-std/udf-catalogue/pull/20#issuecomment-2115053757
183+
for the rationale. */
184+
if (!PG_ARGISNULL(4) && !input.parallax_valid) {
185+
output_null[5] = 0;
186+
output.rv = input.rv;
187+
}
188+
175189
/* change to internal units: rad, rad/yr, mas, and km/s */
176190
retvals[0] = Float8GetDatum(output.pos.lng);
177191
retvals[1] = Float8GetDatum(output.pos.lat);
@@ -181,7 +195,6 @@ epoch_prop(PG_FUNCTION_ARGS) {
181195
retvals[5] = Float8GetDatum(output.rv);
182196

183197
{
184-
bool isnull[6] = {0, 0, 0, 0, 0, 0};
185198
int lower_bounds[1] = {1};
186199
int dims[1] = {6};
187200
#ifdef USE_FLOAT8_BYVAL
@@ -190,13 +203,7 @@ epoch_prop(PG_FUNCTION_ARGS) {
190203
bool embyval = false;
191204
#endif
192205

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,
206+
result = construct_md_array(retvals, output_null, 1, dims, lower_bounds,
200207
FLOAT8OID, sizeof(float8), embyval, 'd');
201208
}
202209
PG_RETURN_ARRAYTYPE_P(result);

‎src/epochprop.h

+1-10
Original file line numberDiff line numberDiff line change
@@ -6,21 +6,12 @@
66
extern Datum epoch_prop(PG_FUNCTION_ARGS);
77

88

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-
{
14-
double x,
15-
y;
16-
} CPoint;
17-
189
typedef struct s_phasevec
1910
{
2011
SPoint pos; /* Position as an SPoint */
2112
double pm[2]; /* Proper motion long/lat in rad/year, PM in
2213
* longitude has cos(lat) applied */
2314
double parallax; /* in rad */
2415
double rv; /* radial velocity in km/s */
25-
int parallax_valid; /* 1 if the parallax really is a NULL */
16+
int parallax_valid; /* 1 if we accept the parallax as physical */
2617
} phasevec;

0 commit comments

Comments
 (0)