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)
postgrespro · Jun 10, 2024 · e1dee24 · e1dee24
1 parent 0607488
commit e1dee24
Show file tree

Hide file tree

Showing 5 changed files with 73 additions and 72 deletions.
diff --git a/doc/functions.sgm b/doc/functions.sgm
@@ -867,9 +867,12 @@
             <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.
+              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.
             </para>
 
             <para>

diff --git a/expected/epochprop.out b/expected/epochprop.out
@@ -1,13 +1,13 @@
-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'),
        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)), 
+       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,93 +16,93 @@ 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'),
        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)), 
+       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'),
        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)), 
+       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'),
        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)), 
+       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'),
        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)), 
+       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,
                23,
                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
@@ -1,66 +1,66 @@
-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'),
        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)), 
+       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'),
        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)), 
+       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'),
        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)), 
+       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'),
        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)), 
+       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'),
        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)), 
+       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
@@ -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
@@ -6,21 +6,12 @@
 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 */
 	double		pm[2];			/* Proper motion long/lat in rad/year, PM in
 								 * 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;