From 70b821c7fac01d45777c6c7b5ffc3e3960ceb34b Mon Sep 17 00:00:00 2001 From: Markus Demleitner Date: Wed, 15 May 2024 15:35:20 +0200 Subject: [PATCH] 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) --- doc/functions.sgm | 9 ++++-- expected/epochprop.out | 62 +++++++++++++++++++++--------------------- sql/epochprop.sql | 26 +++++++++--------- src/epochprop.c | 37 +++++++++++++++---------- src/epochprop.h | 11 +------- 5 files changed, 73 insertions(+), 72 deletions(-) diff --git a/doc/functions.sgm b/doc/functions.sgm index 232b755..047b7f1 100644 --- a/doc/functions.sgm +++ b/doc/functions.sgm @@ -867,9 +867,12 @@ 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. + where some very small default is put in. Whatever is NULL + on the input is NULL on the output. In addition, we null + out a non-NULL input on one component of the proper motion + if the other component is NULL, and we null out the radial + velocity if the parallax is missing, as it would be horribly + off with the propagation algorithm we use here. diff --git a/expected/epochprop.out b/expected/epochprop.out index 3a52832..1111f1a 100644 --- a/expected/epochprop.out +++ b/expected/epochprop.out @@ -1,5 +1,5 @@ -SET extra_float_digits = 2; -SELECT +SET extra_float_digits = 1; +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -7,7 +7,7 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 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; @@ -16,7 +16,7 @@ FROM ( 269.4742714391 | 4.4072939987 | 543.624 | -791.442 | 10235.412 | -110.450 (1 row) -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -24,16 +24,16 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 0, 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.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+----------+------------+---------- + 269.4744079540 | 4.4055337210 | .000 | -801.210 | 10361.762 | -110.000 (1 row) -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -41,16 +41,16 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), NULL, 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.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+---------+----------+------------+---------- + 269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | -110.000 (1 row) -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -58,16 +58,16 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 23, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL, 20) AS tp) AS q; - to_char | to_char | to_char | to_char | to_char | to_char ------------------+-----------------+----------+----------+------------+---------- - 269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 | 2.159 + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+----------+------------+--------- + 269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 | (1 row) -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -75,13 +75,13 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 23, NULL, RADIANS(10362/3.6e6), -110, 120) AS tp) AS q; - to_char | to_char | to_char | to_char | to_char | to_char ------------------+-----------------+----------+----------+------------+---------- - 269.4520769500 | 5.0388680565 | 23.007 | -.000 | 10368.061 | -97.120 + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+---------+---------+---------- + 269.4520769500 | 4.6933649660 | 23.007 | | | -110.000 (1 row) SELECT epoch_prop(NULL, @@ -89,20 +89,20 @@ SELECT epoch_prop(NULL, 0.01 , RADIANS(10362/3.6e6), -110, 120); ERROR: NULL position not supported in epoch propagation -SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), 23, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, 20) AS tp; - tp ---------------------------------------------- - (4.7027479265831289 , 0.082919450934599334) + tp +------------------------------------------- + (4.702747926583129 , 0.08291945093459933) (1 row) -SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), +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.7027479306195161 , 0.082919398938087627) + tp +------------------------------------------- + (4.702747930619516 , 0.08291939893808763) (1 row) diff --git a/sql/epochprop.sql b/sql/epochprop.sql index d8ae2b7..a6c69dd 100644 --- a/sql/epochprop.sql +++ b/sql/epochprop.sql @@ -1,6 +1,6 @@ -SET extra_float_digits = 2; +SET extra_float_digits = 1; -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -8,12 +8,12 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 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; -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -21,12 +21,12 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 0, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, -100) AS tp) AS q; -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -34,12 +34,12 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), NULL, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, -100) AS tp) AS q; -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -47,12 +47,12 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 23, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL, 20) AS tp) AS q; -SELECT +SELECT to_char(DEGREES(tp[1]), '999D9999999999'), to_char(DEGREES(tp[2]), '999D9999999999'), to_char(tp[3], '999D999'), @@ -60,7 +60,7 @@ SELECT to_char(DEGREES(tp[5])*3.6e6, '99999D999'), to_char(tp[6], '999D999') FROM ( - SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), 23, NULL, RADIANS(10362/3.6e6), -110, 120) AS tp) AS q; @@ -70,11 +70,11 @@ SELECT epoch_prop(NULL, 0.01 , RADIANS(10362/3.6e6), -110, 120); -SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), 23, RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, 20) AS tp; -SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), 20) AS tp; diff --git a/src/epochprop.c b/src/epochprop.c index ab9abc8..12ad621 100644 --- a/src/epochprop.c +++ b/src/epochprop.c @@ -133,6 +133,7 @@ epoch_prop(PG_FUNCTION_ARGS) { phasevec input, output; ArrayType *result; Datum retvals[6]; + bool output_null[6] = {0, 0, 0, 0, 0, 0}; if (PG_ARGISNULL(0)) { ereport(ERROR, @@ -141,25 +142,29 @@ epoch_prop(PG_FUNCTION_ARGS) { memcpy(&(input.pos), (void*)PG_GETARG_POINTER(0), sizeof(SPoint)); if (PG_ARGISNULL(1)) { input.parallax = 0; + output_null[2] = 1; + /* The way we do our computation, with a bad parallax the RV + will be horribly off, too, so null this out, too; if avaialble, + we will fiddle in the original RV below again. */ + output_null[5] = 1; } 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)) { + if (PG_ARGISNULL(2) || PG_ARGISNULL(3)) { + input.pm[0] = 0; input.pm[1] = 0; + output_null[3] = 1; + output_null[4] = 1; } else { + input.pm[0] = PG_GETARG_FLOAT8(2); input.pm[1] = PG_GETARG_FLOAT8(3); } if (PG_ARGISNULL(4)) { input.rv = 0; + output_null[5] = 1; } else { input.rv = PG_GETARG_FLOAT8(4); } @@ -172,6 +177,15 @@ epoch_prop(PG_FUNCTION_ARGS) { propagate_phasevec(&input, delta_t, &output); + /* If we have an invalid parallax but a good RV, preserve the original, + untransformed RV on output. See + https://github.com/ivoa-std/udf-catalogue/pull/20#issuecomment-2115053757 + for the rationale. */ + if (!PG_ARGISNULL(4) && !input.parallax_valid) { + output_null[5] = 0; + output.rv = input.rv; + } + /* change to internal units: rad, rad/yr, mas, and km/s */ retvals[0] = Float8GetDatum(output.pos.lng); retvals[1] = Float8GetDatum(output.pos.lat); @@ -181,7 +195,6 @@ epoch_prop(PG_FUNCTION_ARGS) { 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 @@ -190,13 +203,7 @@ epoch_prop(PG_FUNCTION_ARGS) { 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, + result = construct_md_array(retvals, output_null, 1, dims, lower_bounds, FLOAT8OID, sizeof(float8), embyval, 'd'); } PG_RETURN_ARRAYTYPE_P(result); diff --git a/src/epochprop.h b/src/epochprop.h index a93e4c3..3b61a02 100644 --- a/src/epochprop.h +++ b/src/epochprop.h @@ -6,15 +6,6 @@ extern Datum epoch_prop(PG_FUNCTION_ARGS); -/* a cartesian point; this is like geo_decl's point, but you can't -have both geo_decls and pg_sphere right now (both define a type Point, -not to mention they have different ideas on EPSILON */ -typedef struct s_cpoint -{ - double x, - y; -} CPoint; - typedef struct s_phasevec { SPoint pos; /* Position as an SPoint */ @@ -22,5 +13,5 @@ typedef struct s_phasevec * longitude has cos(lat) applied */ double parallax; /* in rad */ double rv; /* radial velocity in km/s */ - int parallax_valid; /* 1 if the parallax really is a NULL */ + int parallax_valid; /* 1 if we accept the parallax as physical */ } phasevec;