alignment/SPK: address code review feedback

Fix undefined behavior in TransformTelescopeToCelestial: spkVtel reads
tar->a and tar->b unconditionally before its internal mode switch, so
leaving them uninitialized in TARG mode was UB. Initialize both to 0.0.

Check the spkAstr return code in UpdateAstrometry. Zero m_Ast before the
call so a failure cannot leave stale data behind; on error, seed eral from
the libnova LST so subsequent spkVtel calls see a coherent ERA rather than
garbage, then return early.

Restore pmfit.c closer to the original vendored source.

Author: ckemper67

libs/alignment/SPKMathPlugin.cpp

@@ -220,6 +220,8 @@ bool SPKMathPlugin::TransformTelescopeToCelestial(const TelescopeDirectionVector
 
     spkTAR tar;
     tar.sys = APPT;
+    tar.a   = 0.0;  // not used by spkVtel in TARG mode, but read unconditionally before the mode switch
+    tar.b   = 0.0;
 
     double tara, tare, tarr, tarp, soln[5];
     int status = spkVtel(TARG, &m_Obs, &m_Opt, &m_PM, &m_Ast, &ax3, &tar, &m_PO,
@@ -267,20 +269,32 @@ void SPKMathPlugin::UpdateAstrometry(double JD)
     utc.mi = date.minutes;
     utc.sec = date.seconds;
 
+    // Compute LST first — needed both for the ERA synchronization below and
+    // as a safe fallback if spkAstr fails.
+    double gmst_hrs = ln_get_apparent_sidereal_time(JD);
+    double lst_rad  = HOURS_TO_RAD(range24(gmst_hrs + RAD_TO_HOURS(m_Obs.slon)));
+    m_LST_rad = lst_rad;
+
     spkEOP eop = {0, 0, 0};
     // pressure=0 disables spkVtel's internal refraction model; temperature and humidity
     // are nominal placeholders.  Atmospheric refraction is not modeled in this pipeline
     // (INDI's coordinate flow stops at "apparent" / JNow; see tar.sys = APPT above).
     spkAIR air = {0.0, 10.0, 0.5};
-    
-    // Use official spkAstr
-    spkAstr(1, utc, 0.0, &eop, &m_Obs, &air, &m_Opt, &m_Ast);
 
-    // Synchronize SOFA internal "clock" (ERA) with libnova LST
+    // Zero m_Ast before calling spkAstr so a failure cannot leave stale data behind.
+    memset(&m_Ast, 0, sizeof(m_Ast));
+    int astrStatus = spkAstr(1, utc, 0.0, &eop, &m_Obs, &air, &m_Opt, &m_Ast);
+    if (astrStatus < 0)
+    {
+        ASSDEBUGF("SPK UpdateAstrometry: spkAstr failed (status=%d), transforms will be inaccurate", astrStatus);
+        // Leave m_Ast zeroed; set eral from libnova LST so subsequent spkVtel
+        // calls see a minimally coherent ERA rather than garbage.
+        m_Ast.astrom.eral = iauAnp(lst_rad);
+        return;
+    }
+
+    // Synchronize SOFA internal "clock" (ERA) with libnova LST.
     // This ensures consistency with the simulator's view of HA/Az.
-    double gmst_hrs = ln_get_apparent_sidereal_time(JD);
-    double lst_rad  = HOURS_TO_RAD(range24(gmst_hrs + RAD_TO_HOURS(m_Obs.slon)));
-    m_LST_rad = lst_rad;
     // eral = ERA + Longitude = LST + EO
     m_Ast.astrom.eral = iauAnp(lst_rad + m_Ast.eo);
 }

libs/alignment/spk/pmfit.c

@@ -1,23 +1,7 @@
 #include <stdio.h>
-#include <math.h>
 #include <stdlib.h>
+#include <math.h>
 #include <string.h>
-#include <stddef.h>
-#include "sofa.h"
-#include "spk.h"
-
-int Bfun ( int, double, char, double, double, double* );
-int Simeqn (int, double*, double* );
-
-int Pmfit ( double, char, int, double*, int,
-            double*, double*, double* );
-
-/* Internal helper Bfun prototype */
-int Bfun ( int, double, char, double, double, double* );
-int Simeqn (int, double*, double* );
-
-/*--------------------------------------------------------------------*/
-
 #include <sofa.h>
 
 int Bfun ( int, double, char, double, double, double* );
@@ -253,7 +237,7 @@ int Pmfit ( double phi, char mount, int n, double* obs, int nt,
    }
 
 /* On-sky variance and RMS. */
-   skyvar = ( 2*n > nt ) ? sumr2 / (double) ( 2*n - nt ) : 0.0;
+   skyvar = sumr2 / (double) (n-nt);
 
 /* Standard deviations of the model coefficients. */
    for ( i = 0; i < nt; i++ ) sigmas[i] = sqrt ( skyvar * a[i*(nt+1)] );
@@ -410,6 +394,27 @@ int Bfun ( int nt, double phi, char mount, double rdem, double pdem,
 ** ---------------------------------------------------
 ** The basis functions at the current point in the sky
 ** ---------------------------------------------------
+**
+** Each term is an [x,y] correction vector at the given point in the
+** sky, where the x-axis is east-west for equatorials or left-right
+** for altazimuths and the y-axis north-south or up-down respectively.
+**
+** The equatorial terms are from Appendix 1 of Wallace, P.T., 1975,
+** Programming the Control Computer of the Anglo-Australian 3.9 Metre
+** Telescope.  In Telescope Automation: proceedings of an International
+** Conference held at MIT, Cambridge, MA USA, April 29 - May 1, 1975,
+** Eds. Maureen K. Huguenin and Thomas B. McCord. Published by MIT,
+** 1975, p.281.
+**
+** The altazimuth terms are from Table 3 in Jeffrey G. Mangum,
+** Jacob W. M. Baars, Albert Greve, Robert Lucas, Ralph C. Snel,
+** Patrick Wallace, Mark Holdaway, 2006, Evaluation of the ALMA
+** Prototype Antennas, Publications of the Astronomical Society of the
+** Pacific, 118, 847, 1260.
+**
+** n.b. The formulas are first-order approximations, not rigorous frame
+** rotations, valid for relatively small values of the pointing-model
+** coefficients.
 */
 
    switch ( (int) mount ) {
@@ -608,4 +613,4 @@ int Simeqn (int n, double *a, double *y )
 **  NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 **  CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 **
-**--------------------------------------------------------------------*/
+**--------------------------------------------------------------------*/

libs/alignment/spk/vtel.c

@@ -1,5 +1,3 @@
-#include <math.h>
-#include <string.h>
 #include "sofa.h"
 #include "spk.h"