From 3bd2cab0e125155e5a61fe1615a890dfd32b56d0 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:08:30 -0700 Subject: [PATCH 1/7] Add sanitizer CI job, standalone C driver and o3o08 regression test. Co-authored-by: Cursor --- .github/workflows/python-app.yml | 46 ++++++++++++++++++++++++++ orbfit/orbfit_api.h | 22 ++++++++++++ orbfit/orbfitmodule.c | 1 + test/c/orbfit_driver.c | 52 +++++++++++++++++++++++++++++ test/test_regression.py | 57 ++++++++++++++++++++++++++++++++ 5 files changed, 178 insertions(+) create mode 100644 orbfit/orbfit_api.h create mode 100644 test/c/orbfit_driver.c create mode 100644 test/test_regression.py diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 27808fc..85d0f3a 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -46,3 +46,49 @@ jobs: - name: Test with pytest run: | pytest + + c-checks: + # Strict C23 compile, then the C library under AddressSanitizer, + # LeakSanitizer and UBSan, both standalone and as loaded by Python. + runs-on: ubuntu-latest + env: + ORBIT_EPHEMERIS: ${{ github.workspace }}/src/mp_ephem/data/binEphem.405_64 + ORBIT_OBSERVATORIES: ${{ github.workspace }}/src/mp_ephem/data/observatories.dat + SANITIZE: -fsanitize=address,undefined -fno-omit-frame-pointer -fno-sanitize-recover=undefined + + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.13" + - name: Compile with strict warnings as C23 + run: | + gcc --version + for f in orbfit/*.c test/c/*.c; do + gcc -std=gnu2x -Wall -Wextra -Wno-unused-parameter -Werror -Iorbfit -fsyntax-only "$f" + done + - name: Run the C library standalone under sanitizers + env: + ASAN_OPTIONS: detect_leaks=1 + UBSAN_OPTIONS: print_stacktrace=1 + run: | + gcc -std=gnu2x -g -O1 $SANITIZE -Iorbfit orbfit/*.c test/c/orbfit_driver.c -lm -o orbfit_driver + ./orbfit_driver test/data/o3o08.mpc "$RUNNER_TEMP/o3o08.abg" 3 + ./orbfit_driver test/data/hst.mpc "$RUNNER_TEMP/hst.abg" 2 + - name: Build the extension with sanitizers + run: | + python -m pip install --upgrade pip + pip install pytest + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + CFLAGS="-g -O1 $SANITIZE" LDFLAGS="$SANITIZE" pip install --no-cache-dir . + - name: Test with pytest under sanitizers + env: + # Python itself is not built with LeakSanitizer, so only check for memory errors here. + ASAN_OPTIONS: detect_leaks=0 + UBSAN_OPTIONS: print_stacktrace=1 + run: | + export LD_PRELOAD="$(gcc -print-file-name=libasan.so):$(gcc -print-file-name=libubsan.so)" + pytest test diff --git a/orbfit/orbfit_api.h b/orbfit/orbfit_api.h new file mode 100644 index 0000000..4dd1a32 --- /dev/null +++ b/orbfit/orbfit_api.h @@ -0,0 +1,22 @@ +/* Entry points of the orbfit shared library, as called from mp_ephem/bk_orbit.py + * through ctypes. Each returns a pointer to a static result array that is + * overwritten by the next call; the library is not reentrant. */ +#ifndef ORBFIT_API_H +#define ORBFIT_API_H + +/* Fit the observations in mpc_filename and write the a/b/g orbit to + * abg_filename. Returns {distance, distance uncertainty} in AU. */ +double *fitradec(char *mpc_filename, char *abg_filename); + +/* RA, Dec (deg), error ellipse a, b (arcsec), PA (deg), distance (AU), + * ecliptic lon, lat (deg) at UTC Julian Date jdate from site obscode. */ +double *predict(char *abg_file, double jdate, int obscode); + +/* Position {x, y, z} (AU) of the target at UTC Julian Date jdate. */ +double *predict_helio(char *abg_file, double jdate, int obscode); + +/* a, e, i, Node, peri, T, their uncertainties, epoch, distance and distance + * uncertainty for the orbit in abg_file. */ +double *abg_to_aei(char *abg_file); + +#endif diff --git a/orbfit/orbfitmodule.c b/orbfit/orbfitmodule.c index 2025d0b..94bbf62 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -1,4 +1,5 @@ #include "orbfit.h" +#include "orbfit_api.h" #include #include #include diff --git a/test/c/orbfit_driver.c b/test/c/orbfit_driver.c new file mode 100644 index 0000000..37d38b4 --- /dev/null +++ b/test/c/orbfit_driver.c @@ -0,0 +1,52 @@ +/* Exercise the orbfit library entry points without Python, so the C code can + * be run under AddressSanitizer / LeakSanitizer / UBSan. + * + * orbfit_driver [repeat] + * + * ORBIT_EPHEMERIS and ORBIT_OBSERVATORIES must point at the data files. + * Exits non-zero if any returned value is not finite. */ +#include +#include +#include + +#include "orbfit_api.h" + +static int +check(const char *what, const double *v, int n, int show) +{ + int i, bad = 0; + if (show) printf("%-14s", what); + for (i = 0; i < n; i++) { + if (!isfinite(v[i])) bad = 1; + if (show) printf(" %.10g", v[i]); + } + if (show) printf("\n"); + if (bad) fprintf(stderr, "non-finite value returned by %s\n", what); + return bad; +} + +int +main(int argc, char *argv[]) +{ + double *r, jd; + int i, repeat, bad = 0; + + if (argc < 3) { + fprintf(stderr, "usage: %s [repeat]\n", argv[0]); + return 2; + } + repeat = argc > 3 ? atoi(argv[3]) : 1; + + for (i = 0; i < repeat; i++) { + r = fitradec(argv[1], argv[2]); + bad |= check("fitradec", r, 2, i == 0); + r = abg_to_aei(argv[2]); + bad |= check("abg_to_aei", r, 15, i == 0); + jd = r[12]; + r = predict(argv[2], jd + 365.25, 568); + bad |= check("predict", r, 8, i == 0); + r = predict_helio(argv[2], jd + 365.25, 568); + bad |= check("predict_helio", r, 3, i == 0); + } + return bad; +} diff --git a/test/test_regression.py b/test/test_regression.py new file mode 100644 index 0000000..221d534 --- /dev/null +++ b/test/test_regression.py @@ -0,0 +1,57 @@ +""" +Pin the orbit fit, prediction and uncertainty for o3o08 so that changes to the +C library that move results are noticed. + +Tolerances allow for input precision at the 0.05 arcsec level (the size of the +float-to-double fix in RA/Dec parsing) but catch any real change in the fit. +""" +import os +import unittest + +from astropy import units + +import mp_ephem + +__PATH__ = os.path.dirname(__file__) + + +class O3o08Regression(unittest.TestCase): + + @classmethod + def setUpClass(cls): + observations = mp_ephem.EphemerisReader().read(os.path.join(__PATH__, 'data', 'o3o08.mpc')) + cls.orbit = mp_ephem.BKOrbit(observations) + + def test_elements(self): + orbit = self.orbit + self.assertAlmostEqual(orbit.a.to(units.au).value, 39.34188, delta=2e-3) + self.assertAlmostEqual(orbit.e.value, 0.277796, delta=1e-4) + self.assertAlmostEqual(orbit.inc.to(units.degree).value, 8.048165, delta=1e-3) + self.assertAlmostEqual(orbit.Node.to(units.degree).value, 113.85162, delta=2e-3) + self.assertAlmostEqual(orbit.om.to(units.degree).value, 66.2447, delta=0.05) + self.assertAlmostEqual(orbit.T.to(units.day).value, 2447884.58, delta=5) + self.assertAlmostEqual(orbit.epoch.jd, 2456392.05115, delta=1e-6) + + def test_distance_and_uncertainty(self): + orbit = self.orbit + self.assertAlmostEqual(orbit.distance.to(units.au).value, 31.65719, delta=1e-3) + self.assertAlmostEqual(orbit.da.to(units.au).value, 0.02748, delta=0.003) + + def test_predict_at_first_observation(self): + orbit = self.orbit + orbit.predict(orbit.observations[0].date) + self.assertAlmostEqual(orbit.coordinate.ra.degree, 238.9183567, delta=0.1 / 3600) + self.assertAlmostEqual(orbit.coordinate.dec.degree, -13.4157912, delta=0.1 / 3600) + self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 0.0994, delta=0.01) + self.assertAlmostEqual(orbit.ddec.to(units.arcsec).value, 0.0630, delta=0.01) + + def test_predict_two_years_out(self): + orbit = self.orbit + orbit.predict("2016 04 09.55115") + self.assertAlmostEqual(orbit.coordinate.ra.degree, 245.3271649, delta=2.0 / 3600) + self.assertAlmostEqual(orbit.coordinate.dec.degree, -15.1861396, delta=2.0 / 3600) + self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 2.382, delta=0.2) + + +if __name__ == '__main__': + unittest.main() From 00451a8c1c35764debf4dc3fbe8453a8a7315ab8 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:12:53 -0700 Subject: [PATCH 2/7] Build cleanly as C23: ANSI nrutil, include guards, no hidden warnings. Co-authored-by: Cursor --- orbfit/aeiderivs.c | 3 - orbfit/dms.c | 8 +- orbfit/ephem_earth.c | 7 +- orbfit/ephem_read.h | 50 ----------- orbfit/mrqcof_orbit.c | 7 +- orbfit/mrqmin_orbit.c | 6 +- orbfit/nrutil.c | 193 +++++++----------------------------------- orbfit/nrutil.h | 25 +++--- orbfit/orbfit.h | 15 ++++ orbfit/orbfit1.c | 28 ++---- orbfit/orbfit2.c | 5 +- orbfit/orbfitmodule.c | 12 +-- orbfit/transforms.c | 14 ++- setup.py | 4 +- 14 files changed, 77 insertions(+), 300 deletions(-) delete mode 100644 orbfit/ephem_read.h diff --git a/orbfit/aeiderivs.c b/orbfit/aeiderivs.c index 8bc133b..e2b9c21 100644 --- a/orbfit/aeiderivs.c +++ b/orbfit/aeiderivs.c @@ -1,8 +1,5 @@ /* $Id: aeiderivs.c,v 2.0 2001/09/14 19:05:06 garyb Exp $ */ -#ifndef lint -static char vcid[] = "$Id: aeiderivs.c,v 2.0 2001/09/14 19:05:06 garyb Exp $"; -#endif /* lint */ /* Subroutine to determine partial derivative matrix of orbital params * w.r.t. barycentric phase space params. Incredible algebra by Bharat. */ diff --git a/orbfit/dms.c b/orbfit/dms.c index fc5ac81..02c9ef1 100644 --- a/orbfit/dms.c +++ b/orbfit/dms.c @@ -1,17 +1,15 @@ /* $Id: dms.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -#ifndef lint -static char vcid[] = "$Id: dms.c,v 1.1 2006/11/22 20:31:50 observe Exp $"; -#endif /* lint */ /* convert a string to decimal degrees */ #include #include #include #include #include +#include "orbfit.h" double -dmsdeg(char *string) +dmsdeg(const char *string) { float args[3]; char substring[10]; @@ -51,7 +49,7 @@ dmsdeg(char *string) /* slight variation for hours:*/ double -hmsdeg(char *string) +hmsdeg(const char *string) { return(15. * dmsdeg(string)); } diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index 0c760bb..5eb1a08 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -1,8 +1,4 @@ /* $Id: ephem_earth.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -/* $Id: ephem_earth.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -#ifndef lint -static char vcid[] = "$Id: ephem_earth.c,v 1.1 2006/11/22 20:31:50 observe Exp $"; -#endif /* lint */ /*** ephem_earth.c - I've changed the below to include a routine explicitly *** returning the location of earth geocenter relative to SSBARY. Also *** have eliminated the nutation & libration routines. @@ -524,7 +520,7 @@ void Interpolate_State(double Time, double Position[3], double Velocity[3]) { - double A[50] , B[50] , Cp[50] , P_Sum[3] , V_Sum[3] , Up[50] , + double A[50] , Cp[50] , P_Sum[3] , V_Sum[3] , Up[50] , T_break , T_seg , T_sub , Tc; int i , j; long int C , G , N , offset = 0; @@ -548,7 +544,6 @@ void Interpolate_State(double Time, for ( i=0 ; i<50 ; i++ ) { A[i] = 0.0; - B[i] = 0.0; } /*--------------------------------------------------------------------------*/ diff --git a/orbfit/ephem_read.h b/orbfit/ephem_read.h deleted file mode 100644 index a2a2258..0000000 --- a/orbfit/ephem_read.h +++ /dev/null @@ -1,50 +0,0 @@ -/* $Id: ephem_read.h,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -/**==========================================================================**/ -/** **/ -/** SOURCE FILE: ephem_read.h **/ -/** **/ -/** This file contains C prototype declarations for the functions that **/ -/** are defined in the source file ephem_read.c; the purpose of each of **/ -/** them is explained there. **/ -/** **/ -/** Programmer: David Hoffman/EG5 **/ -/** NASA, Johnson Space Center **/ -/** Houston, TX 77058 **/ -/** e-mail: david.a.hoffman1@jsc.nasa.gov **/ -/** **/ -/**==========================================================================**/ - -#ifndef TYPES_DEFINED -#include "ephem_types.h" -#endif - -/*----------------------------------------------------------------------------*/ -/* Initialize_Ephemeris */ -/*----------------------------------------------------------------------------*/ - - int Initialize_Ephemeris( char *fileName ); - -/*----------------------------------------------------------------------------*/ -/* Interpolate_Libration */ -/*----------------------------------------------------------------------------*/ - - void Interpolate_Libration( double Time, int Target, double Libration[3] ); - -/*----------------------------------------------------------------------------*/ -/* Interpolate_Nutation */ -/*----------------------------------------------------------------------------*/ - - void Interpolate_Nutation( double Time , int Target , double Nutation[2] ); - -/*----------------------------------------------------------------------------*/ -/* Interpolate_Position */ -/*----------------------------------------------------------------------------*/ - - void Interpolate_Position( double Time , int Target , double Position[3] ); - -/*----------------------------------------------------------------------------*/ -/* Interpolate_State */ -/*----------------------------------------------------------------------------*/ - - void Interpolate_State( double Time , int Target , stateType *Planet ); - diff --git a/orbfit/mrqcof_orbit.c b/orbfit/mrqcof_orbit.c index fc0753e..81acd2c 100644 --- a/orbfit/mrqcof_orbit.c +++ b/orbfit/mrqcof_orbit.c @@ -1,8 +1,4 @@ /* $Id: mrqcof_orbit.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -#ifndef lint -static char vcid[] = "$Id: mrqcof_orbit.c,v 1.1 2006/11/22 20:31:50 observe Exp $"; -#endif /* lint */ - /* Altered version of Numerical Recipes "mrqcof" to deal with our situation * in which we are simultaneously fitting the x and y data. Also have * made use of the existing "OBSERVATION" structure as input data. @@ -28,7 +24,6 @@ void mrqcof_orbit(OBSERVATION obsarray[], { int i,j,k,l,m,mfit=0; double xmod,ymod,wtx,wty,sig2x,sig2y,dx,dy,*dyda,*dxda; - double distance; double fb1; PBASIS params; OBSERVATION *oo; @@ -52,7 +47,7 @@ void mrqcof_orbit(OBSERVATION obsarray[], for (i=1;i<=ndata;i++) { oo = &obsarray[i-1]; - distance = kbo2d(¶ms,oo,&xmod,dxda,&ymod,dyda); + kbo2d(¶ms,oo,&xmod,dxda,&ymod,dyda); sig2x=1.0/(oo->dthetax*oo->dthetax); sig2y=1.0/(oo->dthetay*oo->dthetay); dx=oo->thetax-xmod; diff --git a/orbfit/mrqmin_orbit.c b/orbfit/mrqmin_orbit.c index e9e0e49..3cbbc92 100644 --- a/orbfit/mrqmin_orbit.c +++ b/orbfit/mrqmin_orbit.c @@ -16,8 +16,6 @@ void mrqmin_orbit(OBSERVATION obsarray[], int ndata, double a[], int ia[], int ma, double **covar, double **alpha, double *chisq, double *alamda, double energy_wt) { - void covsrt(double **covar, int ma, int ia[], int mfit); - void gaussj(double **a, int n, double **b, int m); void mrqcof_orbit(OBSERVATION obsarray[], int ndata, double a[], int ia[], int ma, double **alpha, double beta[], double *chisq, @@ -64,8 +62,8 @@ void mrqmin_orbit(OBSERVATION obsarray[], int ndata, double a[], int ia[], } for (j=0,l=1;l<=ma;l++) if (ia[l]) atry[l]=a[l]+da[++j]; - mrqcof_orbit(obsarray,ndata,atry,ia,ma,covar,da,chisq, - energy_wt); + mrqcof_orbit(obsarray,ndata,atry,ia,ma,covar,da,chisq, + energy_wt); if (*chisq < ochisq) { *alamda *= 0.1; ochisq=(*chisq); diff --git a/orbfit/nrutil.c b/orbfit/nrutil.c index f8a4568..de6c488 100644 --- a/orbfit/nrutil.c +++ b/orbfit/nrutil.c @@ -2,219 +2,84 @@ #include #include "nrutil.h" -void nrerror(char error_text[]) +/* Numerical Recipes style vectors and matrices indexed v[nl..nh]. Element 0 + * up to nl-1 is allocated but unused, which avoids forming a pointer before + * the start of the allocation; only nl >= 0 is supported. */ + +void nrerror(const char *error_text) /* Numerical Recipes standard error handler */ { - void _exit(int ); - fprintf(stderr,"Numerical Recipes run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); - _exit(1); + exit(1); } -float *vector(int nl,int nh) -/* allocate a float vector with subscript range v[nl..nh] */ -{ - float *v; - - v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float))-nl; - if (!v) nrerror("allocation failure in vector()"); - return v; -} - -int *ivector(int nl,int nh) +int *ivector(int nl, int nh) /* allocate an int vector with subscript range v[nl..nh] */ { int *v; - v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int))-nl; + if (nl < 0 || nh < nl) nrerror("bad range in ivector()"); + v=(int *)malloc((size_t) (nh+1)*sizeof(int)); if (!v) nrerror("allocation failure in ivector()"); return v; } -double *dvector(int nl,int nh) +double *dvector(int nl, int nh) /* allocate a double vector with subscript range v[nl..nh] */ { double *v; - v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double))-nl; + if (nl < 0 || nh < nl) nrerror("bad range in dvector()"); + v=(double *)malloc((size_t) (nh+1)*sizeof(double)); if (!v) nrerror("allocation failure in dvector()"); return v; } -float **matrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; -/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - int i; - float **m; - - /* allocate pointers to rows */ - m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*))-nrl; - if (!m) nrerror("allocation failure 1 in matrix()"); - - /* allocate rows and set pointers to them */ - for(i=nrl;i<=nrh;i++) { - m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float))-ncl; - if (!m[i]) nrerror("allocation failure 2 in matrix()"); - } - /* return pointer to array of pointers to rows */ - return m; -} - -double **dmatrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; +double **dmatrix(int nrl, int nrh, int ncl, int nch) /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ { int i; double **m; + if (nrl < 0 || nrh < nrl || ncl < 0 || nch < ncl) + nrerror("bad range in dmatrix()"); + /* allocate pointers to rows */ - m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*))-nrl; + m=(double **) calloc((size_t) (nrh+1), sizeof(double*)); if (!m) nrerror("allocation failure 1 in dmatrix()"); /* allocate rows and set pointers to them */ for(i=nrl;i<=nrh;i++) { - m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double))-ncl; + m[i]=(double *) malloc((size_t) (nch+1)*sizeof(double)); if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); } /* return pointer to array of pointers to rows */ return m; } -int **imatrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; -/* allocate an int matrix with subscript range m[nrl..nrh][ncl..nch] */ -{ - int i,**m; - - /* allocate pointers to rows */ - m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*))-nrl; - if (!m) nrerror("allocation failure 1 in imatrix()"); - - /* allocate rows and set pointers to them */ - for(i=nrl;i<=nrh;i++) { - m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int))-ncl; - if (!m[i]) nrerror("allocation failure 2 in imatrix()"); - } - /* return pointer to array of pointers to rows */ - return m; -} - -float **submatrix(a, oldrl, oldrh, oldcl, oldch, newrl, newcl) -/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ -float **a; -int oldrl,oldrh,oldcl; -int oldch; -int newrl,newcl; -{ - int i,j; - float **m; - - /* allocate array of pointers to rows */ - m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*))-newrl; - if (!m) nrerror("allocation failure in submatrix()"); - - /* set pointers to rows */ - for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl; - - /* return pointer to array of pointers to rows */ - return m; -} - -float **convert_matrix(a,nrl,nrh,ncl,nch) -float *a; -int nrl,nrh,ncl,nch; -/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix -declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 -and ncol=nch-ncl+1. The routine should be called with the address -&a[0][0] as the first argument. */ -{ - int i,j,nrow,ncol; - float **m; - - nrow=nrh-nrl+1; - ncol=nch-ncl+1; - /* allocate pointers to rows */ - if ((m=(float **) malloc((unsigned) (nrow)*sizeof(float*))) == NULL) - nrerror("allocation failure in convert_matrix()"); - m -= nrl; - - /* set pointers to rows */ - for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl; - /* return pointer to array of pointers to rows */ - return m; -} - -void free_vector(v,nl,nh) -float *v; -int nl,nh; -/* free a float vector allocated with vector() */ -{ - free((char*) (v+nl)); -} - -void free_ivector(v,nl,nh) -int *v,nl,nh; +void free_ivector(int *v, int nl, int nh) /* free an int vector allocated with ivector() */ { - free((char*) (v+nl)); + (void) nl; (void) nh; + free(v); } -void free_dvector(v,nl,nh) -double *v; -int nl,nh; +void free_dvector(double *v, int nl, int nh) /* free a double vector allocated with dvector() */ { - free((char*) (v+nl)); + (void) nl; (void) nh; + free(v); } -void free_matrix(m,nrl,nrh,ncl,nch) -float **m; -int nrl,nrh,ncl,nch; -/* free a float matrix allocated by matrix() */ -{ - int i; - - for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); - free((char*) (m+nrl)); -} - -void free_dmatrix(m,nrl,nrh,ncl,nch) -double **m; -int nrl,nrh,ncl,nch; +void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch) /* free a double matrix allocated by dmatrix() */ { int i; - for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); - free((char*) (m+nrl)); -} - -void free_imatrix(m,nrl,nrh,ncl,nch) -int **m; -int nrl,nrh,ncl,nch; -/* free an int matrix allocated by imatrix() */ -{ - int i; - - for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); - free((char*) (m+nrl)); -} - -void free_submatrix(b,nrl,nrh,ncl,nch) -/* free a submatrix allocated by submatrix() */ -float **b; -int nrl,nrh,ncl,nch; -{ - free((char*) (b+nrl)); -} - -void free_convert_matrix(b,nrl,nrh,ncl,nch) -/* free a matrix allocated by convert_matrix() */ -float **b; -int nrl,nrh,ncl,nch; -{ - free((char*) (b+nrl)); + (void) ncl; (void) nch; + if (m == NULL) return; + for(i=nrh;i>=nrl;i--) free(m[i]); + free(m); } diff --git a/orbfit/nrutil.h b/orbfit/nrutil.h index bc6ca08..0b88655 100644 --- a/orbfit/nrutil.h +++ b/orbfit/nrutil.h @@ -1,18 +1,19 @@ -float *vector(int, int); -float **matrix(int, int, int, int); +#ifndef NRUTIL_H +#define NRUTIL_H + double *dvector(int, int); double **dmatrix(int, int, int, int); int *ivector(int, int); -int **imatrix(int, int, int, int); -float **submatrix(float **, int, int, int, int, int, int); -void free_vector(float *, int, int); void free_dvector(double *, int, int); void free_ivector(int *, int, int); -void free_matrix(float **, int, int, int, int); void free_dmatrix(double **, int, int, int, int); -void free_imatrix(int **, int, int, int, int); -void free_submatrix(float **, int, int, int, int); -void nrerror(char *); -void lubksb(double **, int, int*, double[]); -float **convert_matrix(float *, int, int, int, int); -void free_convert_matrix(float **, int, int, int, int); \ No newline at end of file +void nrerror(const char *); + +void lubksb(double **, int, int *, double[]); +void ludcmp(double **a, int n, int *indx, double *d); +void gaussj(double **a, int n, double **b, int m); +void covsrt(double **covar, int ma, int ia[], int mfit); +float ran1(long *idum); +float gasdev(long *idum); + +#endif diff --git a/orbfit/orbfit.h b/orbfit/orbfit.h index 259be0a..caab7bd 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -1,5 +1,7 @@ /* $Id: orbfit.h,v 1.1 2006/11/22 20:31:50 observe Exp $ */ /* Definitions for the orbit-fitting software */ +#ifndef ORBFIT_H +#define ORBFIT_H #include "nrutil.h" #include #include @@ -363,3 +365,16 @@ is_visible(OBSERVATION *obs); void fake_observation(PBASIS *p, OBSERVATION *obs); + +/* Sexagesimal string to decimal degrees (hmsdeg: hours to degrees) */ +double dmsdeg(const char *string); +double hmsdeg(const char *string); + +/* Invert 1-indexed matrix in[1..dim][1..dim] into out; in is destroyed. + * Returns non-zero if the matrix is singular. */ +int invert_matrix(double **in, double **out, int dim); + +/* Geocentric ICRS x/y/z (km) to SSBary x/y/z (AU) at UTC JD jd */ +void geo_to_ssbary(double jd, double *x, double *y, double *z); + +#endif diff --git a/orbfit/orbfit1.c b/orbfit/orbfit1.c index e052dfe..2497dcc 100644 --- a/orbfit/orbfit1.c +++ b/orbfit/orbfit1.c @@ -1,8 +1,5 @@ /* $Id: orbfit1.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -#ifndef lint -static char vcid[] = "$Id: orbfit1.c,v 1.1 2006/11/22 20:31:50 observe Exp $"; -#endif /* lint */ /* orbfit1.c - program to fit orbital parameters to observed KBO data. This is basically a template program which makes some simplifying assumptions, with the goal of setting up for a more exact fit and for @@ -15,6 +12,7 @@ encounter in the future. #include "orbfit.h" #include +#include #include "ephem_types.h" #include @@ -25,8 +23,6 @@ double jd0; /* Zeropoint of time scale */ double mpc_dtheta=DEFAULT_DTHETA; /*default astrometric error*/ -int invert_matrix(double **, double **, int); -void geo_to_ssbary(double , double *, double *, double *); void set_mpc_dtheta(double d) { mpc_dtheta = d; @@ -40,9 +36,7 @@ kbo3d_helio(PBASIS *pin, { double vk[3],dxk[7],dyk[7],dzk[7]; double xe,ye,ze; - double invz; double distance; - int i; /* Get the Earth position if not already calculated*/ if (obs->xe < -9.) { @@ -187,7 +181,6 @@ char * fgets_nocomment(char *inbuff, int length, FILE *fpin, FILE *fpout) { - int needmore=0; char test[10]; while (1) { if (fgets(inbuff,length,fpin)==NULL) return(NULL); @@ -212,8 +205,6 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) char two_line_flag[80]; // if set then this is the 2nd line of a two line MPC entry char codestring[16]; double jd; - extern double dmsdeg(char *string); - extern double hmsdeg(char *string); /* get date to see which format this is */ /* For an MPC format the first field will have non-numeric characters*/ jd = strtod(inbuff, &endptr); @@ -299,9 +290,8 @@ read_radec(OBSERVATION obsarray[], char *fname, int *nobs) FILE *fptr; OBSERVATION *obs; int scan_status_flag; - int i; char inbuff[256]; - double jd,ra,dec,elat,elon; + double elat,elon; if (fname==NULL) fptr = stdin; @@ -555,7 +545,7 @@ predict_posn(PBASIS *pin, OBSERVATION *obs, double **sigxy) /*this holds xy error matrix*/ { - int i,j,t; + int i,j; double *dx,*dy, distance; dx=dvector(1,6); @@ -591,9 +581,6 @@ int invert_matrix(double **in, double **out, int dim) { - extern void ludcmp(double **a, int n, int *indx, double *d); - extern void ludcmp(double **a, int n, int *indx, double *d); - int *indx,i,j; double *tvec,det; @@ -633,9 +620,7 @@ void mpc3d(double t, /* time is in years here */ double *x, double *y, double *z) { - double xTelEq,yTelEq,zTelEq; double xec, yec, zec; - double xgeo[3]; geo_to_ssbary(t/DAY+jd0, x, y, z); @@ -660,7 +645,7 @@ earth3d(double t, /* time is in years here */ int obscode, double *x, double *y, double *z) { - double x1,y1,z1,xTelEq,yTelEq,zTelEq; + double x1,y1,z1; double xec, yec, zec; /* get observatory posn wrt barycenter */ @@ -878,7 +863,6 @@ kbo3d(PBASIS *pin, static int init=0; static int tdir; double t1,dt,dtv; - double tstep=20.*DAY; /* decide whether we need to reset integrator to t=0*/ if (!init) { @@ -1085,8 +1069,6 @@ fake_observation(PBASIS *p, OBSERVATION *obs) { static long idum=-1; - float gasdev(long *idum); - float distance; /* seed the random number generator*/ if (idum<0) { @@ -1095,7 +1077,7 @@ fake_observation(PBASIS *p, idum = -tp.tv_usec; } - distance = kbo2d(p,obs,&(obs->thetax),NULL,&(obs->thetay),NULL); + kbo2d(p,obs,&(obs->thetax),NULL,&(obs->thetay),NULL); /* Add measurement noise to the positions */ obs->thetax += obs->dthetax*gasdev(&idum); diff --git a/orbfit/orbfit2.c b/orbfit/orbfit2.c index 92f96d4..ee0fbf4 100644 --- a/orbfit/orbfit2.c +++ b/orbfit/orbfit2.c @@ -1,8 +1,5 @@ /* $Id: orbfit2.c,v 2.0 2001/09/14 19:05:06 garyb Exp $ */ -#ifndef lint -static char vcid[] = "$Id: orbfit2.c,v 2.0 2001/09/14 19:05:06 garyb Exp $"; -#endif /* lint */ /***** Routines that execute the fit of orbit to observation data. ****/ /*** assembled 6/14/00 gmb from subroutines scattered about elsewhere. ** will use the mrqmin_orbit and mrqcof_orbit subroutines as well. @@ -51,7 +48,7 @@ fit_observations(OBSERVATION obsarray[], int *dof, FILE *logfile) { - int *ia,i,j; + int *ia,i; double *a; double gbind2; double energy_fit; /*weight given to binding-energy constraint*/ diff --git a/orbfit/orbfitmodule.c b/orbfit/orbfitmodule.c index 94bbf62..43b8b9b 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -8,7 +8,6 @@ double *fitradec(char *mpc_filename, char *abg_filename) { FILE *abg_file ; - FILE *res_file ; OBSERVATION obsarray[MAXOBS]; int nobs; @@ -21,7 +20,6 @@ double *fitradec(char *mpc_filename, char *abg_filename) static double result[2]; double **covar; double chisq; - int i; int dof; covar = dmatrix(1,6,1,6); @@ -145,19 +143,13 @@ double *predict(char *abg_file, double jdate, int obscode) PBASIS p; OBSERVATION futobs; - struct date_time dt; - char inbuff[256], rastring[20], decstring[20]; - char outbuff[256]; - char *f_string; double **covar,**sigxy,a,b,PA,**derivs; double lat,lon,**covecl; double ra,dec, **coveq; - double yr,mo,day,hr,mn,ss; double xx,yy,xy,bovasqrd,det; double distance; static double result[8]; - int i,nfields; - int iarg=1; + int i; sigxy = dmatrix(1,2,1,2); derivs = dmatrix(1,2,1,2); @@ -240,8 +232,6 @@ double *abg_to_aei(char *abg_file) double d, dd; double **covar_abg, **covar_xyz, **derivs, **covar_aei; - int i,j; - covar_abg = dmatrix(1,6,1,6); covar_xyz = dmatrix(1,6,1,6); covar_aei = dmatrix(1,6,1,6); diff --git a/orbfit/transforms.c b/orbfit/transforms.c index 21bcc0c..ef3cf7e 100644 --- a/orbfit/transforms.c +++ b/orbfit/transforms.c @@ -1,8 +1,5 @@ /* $Id: transforms.c,v 1.1 2006/11/22 20:31:50 observe Exp $ */ -#ifndef lint -static char vcid[] = "$Id: transforms.c,v 1.1 2006/11/22 20:31:50 observe Exp $"; -#endif /* lint */ /*********** Coordinate transformation routines *************/ /* All angles assumed to be in radians upon input */ @@ -16,7 +13,7 @@ eq_to_ec( double ra_eq, double *lon_ec, double **partials) { - double sd,cd,cr,se,ce,y,x; + double sd,cd,se,ce,y,x; se = sin(ECL); ce = cos(ECL); @@ -75,7 +72,7 @@ ec_to_eq( double lat_ec, double *dec_eq, double **partials) { - double sd,cd,cr,se,ce,y,x; + double sd,cd,se,ce,y,x; se = sin(-ECL); ce = cos(ECL); @@ -303,7 +300,7 @@ void orbitElements(XVBASIS *xv, ORBIT *orb) { - int i,j,k; + int k; double combinedMass; /* mass of Sun + mass of KBO */ double epochTime; @@ -325,8 +322,7 @@ orbitElements(XVBASIS *xv, double hMagnitude, ascendingNodeMagnitude; /* magnitude of angular momentum */ double ascEccDotProduct, argumentOfPerifocus; double xBar, yBar; - double cosE, sinE, E1, E2, eccentricAnomaly; - /* E1 and E2 are used to decide the quadrant of Eccentric Anomaly */ + double cosE, sinE, eccentricAnomaly; double meanAnomaly, meanMotion, timeOfPerifocalPassage; combinedMass = GM * 1.00134 ; /* Alter GM to account for total SS mass*/ @@ -462,7 +458,7 @@ elements_to_xv(ORBIT *o, double jd, XVBASIS *xv) { - double eccentricAnomaly, r0[3], v0[3], r1[3], v1[3], r2[3], v2[3]; + double r0[3], v0[3], r1[3], v1[3], r2[3], v2[3]; double meanAnomaly; double c, s, t, dt; diff --git a/setup.py b/setup.py index 6a74e81..c96ed34 100644 --- a/setup.py +++ b/setup.py @@ -27,7 +27,5 @@ packages=find_packages(where='src'), package_dir = {"": "src"}, headers = headers, - ext_modules=[Extension('mp_ephem.orbfit', sources, - extra_compile_args=['-Wno-unused-variable'], - )], + ext_modules=[Extension('mp_ephem.orbfit', sources)], ) From 9a53a4454e87368aae8e9886ce1886ed67b55dee Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:17:49 -0700 Subject: [PATCH 3/7] Parse RA/Dec in double precision and write .abg files at full precision. dmsdeg accumulated in float, giving up to 0.05 arcsec of RA error; results for short arcs move by a fraction of their quoted uncertainty. The 80-column MPC reader also dropped the last RA/Dec digit and lost the signs of satellite positions on two-line records. Co-authored-by: Cursor --- orbfit/dms.c | 65 +++------------------------------- orbfit/ephem_earth.c | 17 +++++---- orbfit/orbfit.h | 13 ++++--- orbfit/orbfit1.c | 78 ++++++++++++++++++++++++++++++----------- orbfit/orbfitmodule.c | 4 +-- test/c/orbfit_driver.c | 29 +++++++++------ test/data/o3o08.abg | 16 ++++----- test/test_orbit.py | 13 +++---- test/test_regression.py | 35 +++++++++--------- 9 files changed, 133 insertions(+), 137 deletions(-) diff --git a/orbfit/dms.c b/orbfit/dms.c index 02c9ef1..d3a1b9d 100644 --- a/orbfit/dms.c +++ b/orbfit/dms.c @@ -11,13 +11,13 @@ double dmsdeg(const char *string) { - float args[3]; - char substring[10]; + double args[3]; + char substring[32]; int i,iarg,j,sgn; for (iarg=0;iarg<3;args[iarg++]=0.); /*strip off leading whitespace and look for - sign*/ - for (i=0; isspace( (int) string[i]); i++); + for (i=0; isspace( (unsigned char) string[i]); i++); if (string[i]=='-') { sgn = -1; i++; @@ -28,11 +28,11 @@ dmsdeg(const char *string) for (iarg=0; iarg<3; iarg++) { substring[0] = 0; - j = sscanf(string+i, "%[ :]",substring); + j = sscanf(string+i, "%31[ :]",substring); if (j==EOF || *(string+i)==0) break; /*failing to get EOF at end*/ i += strlen(substring); substring[0] = 0; - j = sscanf(string+i,"%[-+.0-9]",substring); + j = sscanf(string+i,"%31[-+.0-9]",substring); if (j==EOF || *(string+i)==0) break; else if (j!=1) { fprintf(stderr,"Error interpeting dms ->%s<-\n",string); @@ -54,58 +54,3 @@ hmsdeg(const char *string) return(15. * dmsdeg(string)); } -void -degdms(double degr, - char *outbuff) -{ - int ideg,imin,sign=1; - double sec; - - if (degr<0.) { - sign=-1; - degr *= -1.; - } - - ideg = floor(degr); - degr = (degr-ideg)*60.; - imin = floor(degr); - sec = (degr-imin)*60.; - - if (sign>0) - if (ideg<100) - sprintf(outbuff," +%02d:%02d:%06.3f",ideg,imin,sec); - else - sprintf(outbuff,"+%03d:%02d:%06.3f",ideg,imin,sec); - else - if (ideg<100) - sprintf(outbuff," -%02d:%02d:%06.3f",ideg,imin,sec); - else - sprintf(outbuff,"-%03d:%02d:%06.3f",ideg,imin,sec); - return; -} -void -deghms(double degr, - char *outbuff) -{ - int ideg,imin,sign=1; - double sec; - - if (degr<0) { - sign=-1; - degr *= -1.; - } - - degr /= 15.; - ideg = floor(degr); - degr = (degr-ideg)*60.; - imin = floor(degr); - sec = (degr-imin)*60.; - - if (sign>0) - sprintf(outbuff," %02d:%02d:%07.4f",ideg,imin,sec); - else - sprintf(outbuff,"-%02d:%02d:%07.4f",ideg,imin,sec); - - return; -} - diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index 5eb1a08..974044c 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -876,9 +876,8 @@ observatory_geocenter(double jd, } -/* Copy columns [start, end) of a line, or fewer if the line is short */ -static void -column(const char *line, int start, int end, char *out) +void +copy_columns(const char *line, int start, int end, char *out) { int n = strcspn(line, "\r\n"); if (nn) end=n; - strncpy(out, line+start, end-start); + memcpy(out, line+start, end-start); out[end-start]=0; } @@ -936,7 +935,7 @@ read_observatories(char *fname) int nlon, ncos, nsin; char *eol; - column(inbuff, 0, 3, code); + copy_columns(inbuff, 0, 3, code); if ((obscode = obscode_from_string(code))==OBSCODE_INVALID) continue; if (strchr(inbuff, ':')!=NULL && strchr(inbuff, ':') < inbuff+30) { fprintf(stderr,"Observatories file %s is not in the MPC ObsCodes format:\n->%s\n", @@ -955,9 +954,9 @@ read_observatories(char *fname) sss->code = obscode; sss->warned = 0; - column(inbuff, 4, 13, lonstring); - column(inbuff, 13, 21, cosstring); - column(inbuff, 21, 30, sinstring); + copy_columns(inbuff, 4, 13, lonstring); + copy_columns(inbuff, 13, 21, cosstring); + copy_columns(inbuff, 21, 30, sinstring); nlon = sscanf(lonstring, "%lf", &lon); ncos = sscanf(cosstring, "%lf", &(sss->rhocos)); nsin = sscanf(sinstring, "%lf", &(sss->rhosin)); @@ -973,7 +972,7 @@ read_observatories(char *fname) exit(1); } - column(inbuff, 30, 30+79, sss->name); + copy_columns(inbuff, 30, 30+79, sss->name); for (eol=sss->name+strlen(sss->name); eol>sss->name && isspace((unsigned char) eol[-1]); eol--) ; *eol = 0; nsites++; diff --git a/orbfit/orbfit.h b/orbfit/orbfit.h index caab7bd..1596348 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -45,6 +45,7 @@ #define ARCSEC (PI/180./3600.) #define DAY (1./365.25) /*Julian day, 86400 s*/ #define SPEED_OF_LIGHT 63241.06515 /*in AU/YR*/ +#define AU_KM 149597870.7 /*IAU 2012 astronomical unit in km*/ #define ECL (23.43928*PI/180.) /*Obliquity of ecliptic at J2000*/ #define TSTEP (20.*DAY) /*Time step for orbit integrator*/ @@ -78,10 +79,10 @@ struct date_time { int y; int mo; - float d; - float h; - float mn; - float s; + double d; + double h; + double mn; + double s; }; /* Data of an individual observation */ @@ -377,4 +378,8 @@ int invert_matrix(double **in, double **out, int dim); /* Geocentric ICRS x/y/z (km) to SSBary x/y/z (AU) at UTC JD jd */ void geo_to_ssbary(double jd, double *x, double *y, double *z); +/* Copy 0-indexed columns [start, end) of line into out (size end-start+1), + * stopping at the end of the line. */ +void copy_columns(const char *line, int start, int end, char *out); + #endif diff --git a/orbfit/orbfit1.c b/orbfit/orbfit1.c index 2497dcc..107eada 100644 --- a/orbfit/orbfit1.c +++ b/orbfit/orbfit1.c @@ -199,41 +199,79 @@ fgets_nocomment(char *inbuff, int length, * part of OBSERVATION structure, ra & dec in the x & y parts. Format * error returns non-zero. */ +/* Read a signed coordinate from the second line of an MPC satellite + * observation: sign in column sign_col, value in the next 11 columns. */ +static int +mpc_signed_value(const char *line, int sign_col, double *value) +{ + char field[16]; + copy_columns(line, sign_col+1, sign_col+12, field); + if (sscanf(field, "%lf", value) != 1) return 1; + if (line[sign_col] == '-') *value = -*value; + else if (line[sign_col] != '+' && line[sign_col] != ' ') return 1; + return 0; +} + int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) { char rastring[80],decstring[80],*endptr; - char two_line_flag[80]; // if set then this is the 2nd line of a two line MPC entry char codestring[16]; double jd; /* get date to see which format this is */ /* For an MPC format the first field will have non-numeric characters*/ jd = strtod(inbuff, &endptr); - if (jd==0. || *endptr==0 || !isspace(*endptr)) { - /* See if this is perhaps in MPC format */ + if (jd==0. || *endptr==0 || !isspace((unsigned char) *endptr)) { + /* 80-column MPC format */ struct date_time dd; - if (sscanf(inbuff+15,"%d %d %f",&(dd.y),&(dd.mo),&(dd.d))!=3) { + char field[32]; + + if (strcspn(inbuff, "\r\n") > 32 && inbuff[14] == 's') { + /* Second line of a satellite observation: geocentric position of + * the observer for the previous line, in km (column 33 = 1) or + * AU (column 33 = 2). */ + double x, y, z, scale; + if (previous == NULL) { + fprintf(stderr,"Satellite position line without an observation:\n ->%s\n",inbuff); + return(1); + } + if (inbuff[32] == '1') scale = 1.; + else if (inbuff[32] == '2') scale = AU_KM; + else { + fprintf(stderr,"Unknown units in satellite position line:\n ->%s\n",inbuff); + return(1); + } + if (mpc_signed_value(inbuff, 34, &x) || mpc_signed_value(inbuff, 46, &y) + || mpc_signed_value(inbuff, 58, &z)) { + fprintf(stderr,"Format error in satellite position line:\n ->%s\n",inbuff); + return(1); + } + previous->xe = x*scale; + previous->ye = y*scale; + previous->ze = z*scale; + return(-1); + } + + copy_columns(inbuff, 15, 32, field); + if (sscanf(field,"%d %d %lf",&(dd.y),&(dd.mo),&(dd.d))!=3) { fprintf(stderr,"Format error in observation file:\n ->%s\n",inbuff); return(1); } dd.h = dd.mn = dd.s = 0.; - jd = date_to_jd(dd); - sscanf(inbuff+32,"%s", two_line_flag); - - if (strncmp(two_line_flag, "1\0", 2)==0) { - // this is the 2nd line which will contain the x/y/z location of the observations - sscanf(inbuff+32,"%s - %lf - %lf - %lf", two_line_flag, &(previous->xe), &(previous->ye), &(previous->ze)); - return(-1); + if ((jd = date_to_jd(dd)) == 0.) return(1); + copy_columns(inbuff, 32, 44, rastring); + copy_columns(inbuff, 44, 56, decstring); + copy_columns(inbuff, 77, 80, field); + if (sscanf(field,"%15s",codestring) != 1) { + fprintf(stderr,"Missing observatory code:\n ->%s\n",inbuff); + return(1); } - strncpy(rastring,inbuff+32,11); - strncpy(decstring,inbuff+44,11); - sscanf(inbuff+77,"%3s",codestring); obs->obscode = obscode_from_string(codestring); obs->dthetay = mpc_dtheta; } else if (jd<10000.) { /* See if perhaps this was y/m/d instead of JD: */ struct date_time dd; - if (sscanf(inbuff, "%d %d %f %s %s %lf %15s", + if (sscanf(inbuff, "%d %d %lf %79s %79s %lf %15s", &(dd.y), &(dd.mo), &(dd.d), rastring, decstring, &(obs->dthetay), codestring) != 7) { @@ -242,11 +280,11 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) } obs->obscode = obscode_from_string(codestring); dd.h = dd.mn = dd.s = 0.; - jd = date_to_jd(dd); + if ((jd = date_to_jd(dd)) == 0.) return(1); } else { - if (sscanf(inbuff, "%lf %s %s %lf %lf %lf %lf", + if (sscanf(inbuff, "%lf %79s %79s %lf %lf %lf %lf", &jd, rastring, decstring, &(obs->dthetay), &(obs->xe), &(obs->ye), &(obs->ze)) != 7) { - if (sscanf(inbuff, "%lf %s %s %lf %15s", + if (sscanf(inbuff, "%lf %79s %79s %lf %15s", &jd, rastring, decstring, &(obs->dthetay), codestring) != 5) { @@ -288,7 +326,7 @@ int read_radec(OBSERVATION obsarray[], char *fname, int *nobs) { FILE *fptr; - OBSERVATION *obs; + OBSERVATION *obs = NULL; int scan_status_flag; char inbuff[256]; double elat,elon; @@ -606,7 +644,7 @@ print_matrix(FILE *fptr, double **matrix, int xdim, int ydim) int i,j; for (i=1; i<=ydim; i++) { for (j=1; j<=xdim; j++) { - fprintf(fptr,"%11.4e ",matrix[i][j]); + fprintf(fptr,"%24.16e ",matrix[i][j]); } fprintf(fptr,"\n"); } diff --git a/orbfit/orbfitmodule.c b/orbfit/orbfitmodule.c index 43b8b9b..50674b8 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -44,7 +44,7 @@ double *fitradec(char *mpc_filename, char *abg_filename) /* fprintf(stderr, "# Chi-squared of fit: %.2f DOF: %d\n",chisq,dof); */ fprintf(abg_file, "# Exact a, adot, b, bdot, g, gdot:\n"); - fprintf(abg_file, "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f\n",p.a,p.adot,p.b, + fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n",p.a,p.adot,p.b, p.bdot, p.g, p.gdot); pbasis_to_bary(&p, &xv, NULL); @@ -65,7 +65,7 @@ double *fitradec(char *mpc_filename, char *abg_filename) /* Print out information on the coordinate system */ fprintf(abg_file, "# lat0 lon0 xBary yBary zBary JD0\n"); - fprintf(abg_file, "%12.7f %12.7f %10.7f %10.7f %10.7f %.6f\n", + fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n", lat0/DTOR,lon0/DTOR,xBary,yBary,zBary,jd0); fclose(abg_file); diff --git a/test/c/orbfit_driver.c b/test/c/orbfit_driver.c index 37d38b4..093fcd0 100644 --- a/test/c/orbfit_driver.c +++ b/test/c/orbfit_driver.c @@ -4,23 +4,26 @@ * orbfit_driver [repeat] * * ORBIT_EPHEMERIS and ORBIT_OBSERVATORIES must point at the data files. - * Exits non-zero if any returned value is not finite. */ + * Exits non-zero if a returned position or element is not finite. + * Uncertainties may be NaN when the arc leaves the orbit undetermined. */ #include #include #include #include "orbfit_api.h" +/* required[] lists the indices of v that must be finite, ending with -1 */ static int -check(const char *what, const double *v, int n, int show) +check(const char *what, const double *v, int n, const int *required, int show) { int i, bad = 0; - if (show) printf("%-14s", what); - for (i = 0; i < n; i++) { - if (!isfinite(v[i])) bad = 1; - if (show) printf(" %.10g", v[i]); + if (show) { + printf("%-14s", what); + for (i = 0; i < n; i++) printf(" %.10g", v[i]); + printf("\n"); } - if (show) printf("\n"); + for (i = 0; required[i] >= 0; i++) + if (!isfinite(v[required[i]])) bad = 1; if (bad) fprintf(stderr, "non-finite value returned by %s\n", what); return bad; } @@ -28,6 +31,10 @@ check(const char *what, const double *v, int n, int show) int main(int argc, char *argv[]) { + static const int fit_req[] = {0, -1}; + static const int aei_req[] = {0, 1, 2, 3, 4, 5, 12, 13, -1}; + static const int predict_req[] = {0, 1, 5, 6, 7, -1}; + static const int helio_req[] = {0, 1, 2, -1}; double *r, jd; int i, repeat, bad = 0; @@ -39,14 +46,14 @@ main(int argc, char *argv[]) for (i = 0; i < repeat; i++) { r = fitradec(argv[1], argv[2]); - bad |= check("fitradec", r, 2, i == 0); + bad |= check("fitradec", r, 2, fit_req, i == 0); r = abg_to_aei(argv[2]); - bad |= check("abg_to_aei", r, 15, i == 0); + bad |= check("abg_to_aei", r, 15, aei_req, i == 0); jd = r[12]; r = predict(argv[2], jd + 365.25, 568); - bad |= check("predict", r, 8, i == 0); + bad |= check("predict", r, 8, predict_req, i == 0); r = predict_helio(argv[2], jd + 365.25, 568); - bad |= check("predict_helio", r, 3, i == 0); + bad |= check("predict_helio", r, 3, helio_req, i == 0); } return bad; } diff --git a/test/data/o3o08.abg b/test/data/o3o08.abg index da5aeb0..e5826ad 100644 --- a/test/data/o3o08.abg +++ b/test/data/o3o08.abg @@ -1,11 +1,11 @@ # Exact a, adot, b, bdot, g, gdot: - 0.00001882 0.03844611 -0.00000100 -0.00319661 0.03237758 0.00873920 +1.8685477512898321e-05 0.038446089532549368 -1.0260138542225362e-06 -0.0031966413924042919 0.032377591206302167 0.0087421759049221386 # Covariance matrix: - 2.3263e-13 2.9633e-13 -9.4643e-16 -4.9622e-14 3.8138e-13 1.4037e-11 - 2.9633e-13 1.5949e-12 1.9597e-15 -1.9623e-13 1.1534e-12 5.4228e-11 --9.4643e-16 1.9597e-15 9.3316e-14 -9.2562e-14 -1.3215e-15 1.1610e-13 --4.9622e-14 -1.9623e-13 -9.2562e-14 2.1178e-13 -1.5176e-13 -7.6622e-12 - 3.8138e-13 1.1534e-12 -1.3215e-15 -1.5176e-13 1.0685e-12 4.2575e-11 - 1.4037e-11 5.4228e-11 1.1610e-13 -7.6622e-12 4.2575e-11 2.1086e-09 + 2.3263472266901926e-13 2.9632725418606829e-13 -9.4640004148286025e-16 -4.9622609738329566e-14 3.8138296432087947e-13 1.4036773723832638e-11 + 2.9632725418606819e-13 1.5948623846921276e-12 1.9598479056736960e-15 -1.9623042546551733e-13 1.1533904221971157e-12 5.4227553634461145e-11 + -9.4640004148286735e-16 1.9598479056736695e-15 9.3316029791426523e-14 -9.2562061998268248e-14 -1.3214117225642389e-15 1.1610977393780694e-13 + -4.9622609738329541e-14 -1.9623042546551728e-13 -9.2562061998268261e-14 2.1178007773452764e-13 -1.5175824335113173e-13 -7.6623716265248177e-12 + 3.8138296432087942e-13 1.1533904221971157e-12 -1.3214117225642180e-15 -1.5175824335113178e-13 1.0685047142269171e-12 4.2574794601369175e-11 + 1.4036773723832635e-11 5.4227553634461139e-11 1.1610977393780795e-13 -7.6623716265248193e-12 4.2574794601369168e-11 2.1085814475869564e-09 # lat0 lon0 xBary yBary zBary JD0 - 6.8056487 -120.3801450 0.6426369 0.0913121 -0.7649654 2456392.051150 +6.8056498737182647 -120.3801384024556 0.64263700278996438 0.091312060333539641 -0.76496534283853379 2456392.0511500002 diff --git a/test/test_orbit.py b/test/test_orbit.py index 0e580d4..527d2c6 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -22,7 +22,8 @@ def test_parser(self): def test_fitradec(self): orbit = mp_ephem.BKOrbit(None, self.mpc_filename) orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.a.to('au').value, 91.82, 0) + # A 36 minute arc: a is barely constrained (da is about 25 au), so this only pins the fit. + self.assertAlmostEqual(orbit.a.to('au').value, 56.22, 0) # JPL SBDB gives i=3.32 deg for 2020 KQ11 self.assertAlmostEqual(orbit.inc.to('degree').value, 3.3, delta=0.5) @@ -38,7 +39,7 @@ def test_parser(self): def test_fitradec(self): orbit = mp_ephem.BKOrbit(None, self.mpc_filename) orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.a.to('au').value, 43.68817, 3) + self.assertAlmostEqual(orbit.a.to('au').value, 43.76895, 3) class OrbitFit(unittest.TestCase): @@ -65,12 +66,12 @@ def test_orbit(self): :return: """ - self.assertAlmostEqual(self.orbit.a.to(units.AU).value, 39.3419, 3) - self.assertAlmostEqual(self.orbit.e.value, 0.2778, 3) + self.assertAlmostEqual(self.orbit.a.to(units.AU).value, 39.3437, 3) + self.assertAlmostEqual(self.orbit.e.value, 0.2779, 3) self.assertAlmostEqual(self.orbit.inc.to(units.degree).value, 8.05, 2) self.assertAlmostEqual(self.orbit.Node.to(units.degree).value, 113.85, 2) self.assertAlmostEqual(self.orbit.om.to(units.degree).value, 66.24, 2) - self.assertAlmostEqual(self.orbit.T.to(units.day).value, 2447884.5762, 3) + self.assertAlmostEqual(self.orbit.T.to(units.day).value, 2447883.8480, 3) self.assertAlmostEqual(self.orbit.epoch.jd, 2456392.05115, 4) def test_summarize(self): @@ -211,4 +212,4 @@ def test_parser(self): def test_fitradec(self): orbit = mp_ephem.BKOrbit(None, self.mpc_filename) orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.a.to('au').value, 44.48, 2) # Example value, adjust as needed \ No newline at end of file + self.assertAlmostEqual(orbit.a.to('au').value, 44.31, 2) # Example value, adjust as needed \ No newline at end of file diff --git a/test/test_regression.py b/test/test_regression.py index 221d534..42e4419 100644 --- a/test/test_regression.py +++ b/test/test_regression.py @@ -2,8 +2,9 @@ Pin the orbit fit, prediction and uncertainty for o3o08 so that changes to the C library that move results are noticed. -Tolerances allow for input precision at the 0.05 arcsec level (the size of the -float-to-double fix in RA/Dec parsing) but catch any real change in the fit. +Values are from the double-precision RA/Dec parser and full-precision .abg +files. The float parser it replaced shifted a by 0.0018 au and the 2016 RA +by 0.22 arcsec, which these tolerances are tight enough to catch. """ import os import unittest @@ -24,33 +25,33 @@ def setUpClass(cls): def test_elements(self): orbit = self.orbit - self.assertAlmostEqual(orbit.a.to(units.au).value, 39.34188, delta=2e-3) - self.assertAlmostEqual(orbit.e.value, 0.277796, delta=1e-4) - self.assertAlmostEqual(orbit.inc.to(units.degree).value, 8.048165, delta=1e-3) - self.assertAlmostEqual(orbit.Node.to(units.degree).value, 113.85162, delta=2e-3) - self.assertAlmostEqual(orbit.om.to(units.degree).value, 66.2447, delta=0.05) - self.assertAlmostEqual(orbit.T.to(units.day).value, 2447884.58, delta=5) + self.assertAlmostEqual(orbit.a.to(units.au).value, 39.343702, delta=2e-4) + self.assertAlmostEqual(orbit.e.value, 0.2778718, delta=1e-5) + self.assertAlmostEqual(orbit.inc.to(units.degree).value, 8.048180, delta=1e-4) + self.assertAlmostEqual(orbit.Node.to(units.degree).value, 113.851459, delta=2e-4) + self.assertAlmostEqual(orbit.om.to(units.degree).value, 66.23552, delta=2e-3) + self.assertAlmostEqual(orbit.T.to(units.day).value, 2447883.848, delta=0.2) self.assertAlmostEqual(orbit.epoch.jd, 2456392.05115, delta=1e-6) def test_distance_and_uncertainty(self): orbit = self.orbit - self.assertAlmostEqual(orbit.distance.to(units.au).value, 31.65719, delta=1e-3) - self.assertAlmostEqual(orbit.da.to(units.au).value, 0.02748, delta=0.003) + self.assertAlmostEqual(orbit.distance.to(units.au).value, 31.657179, delta=1e-4) + self.assertAlmostEqual(orbit.da.to(units.au).value, 0.027496, delta=3e-4) def test_predict_at_first_observation(self): orbit = self.orbit orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.coordinate.ra.degree, 238.9183567, delta=0.1 / 3600) - self.assertAlmostEqual(orbit.coordinate.dec.degree, -13.4157912, delta=0.1 / 3600) - self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 0.0994, delta=0.01) - self.assertAlmostEqual(orbit.ddec.to(units.arcsec).value, 0.0630, delta=0.01) + self.assertAlmostEqual(orbit.coordinate.ra.degree, 238.9183555, delta=0.02 / 3600) + self.assertAlmostEqual(orbit.coordinate.dec.degree, -13.4157912, delta=0.02 / 3600) + self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 0.0994, delta=0.002) + self.assertAlmostEqual(orbit.ddec.to(units.arcsec).value, 0.0630, delta=0.002) def test_predict_two_years_out(self): orbit = self.orbit orbit.predict("2016 04 09.55115") - self.assertAlmostEqual(orbit.coordinate.ra.degree, 245.3271649, delta=2.0 / 3600) - self.assertAlmostEqual(orbit.coordinate.dec.degree, -15.1861396, delta=2.0 / 3600) - self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 2.382, delta=0.2) + self.assertAlmostEqual(orbit.coordinate.ra.degree, 245.3271034, delta=0.05 / 3600) + self.assertAlmostEqual(orbit.coordinate.dec.degree, -15.1861295, delta=0.05 / 3600) + self.assertAlmostEqual(orbit.dra.to(units.arcsec).value, 2.382, delta=0.02) if __name__ == '__main__': From 24413075208799b29fd040d72b1ef97a5078a8aa Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:21:36 -0700 Subject: [PATCH 4/7] Fix memory leaks and undefined behaviour in the C library. Free all work matrices in predict, predict_helio, mrqfit and on error paths; return NaN rather than stale results on failure; check fopen; bound the observation count; initialise pivots in ludcmp/gaussj; fix the kbo2d/kbo3d NULL checks and ec_to_proj partials; keep orbitElements finite for circular and ecliptic orbits. A degenerate linear fit now falls back to the energy-constrained fit instead of exiting. Co-authored-by: Cursor --- orbfit/gaussj.c | 2 + orbfit/ludcmp.c | 10 ++- orbfit/nrutil.h | 2 +- orbfit/orbfit1.c | 84 +++++++++++++------- orbfit/orbfit2.c | 18 +++-- orbfit/orbfitmodule.c | 173 +++++++++++++++++------------------------- orbfit/transforms.c | 76 +++++++++++++------ 7 files changed, 198 insertions(+), 167 deletions(-) diff --git a/orbfit/gaussj.c b/orbfit/gaussj.c index c938ae0..8fc8db5 100644 --- a/orbfit/gaussj.c +++ b/orbfit/gaussj.c @@ -15,6 +15,7 @@ void gaussj(double **a, int n, double **b, int m) for (j=1;j<=n;j++) ipiv[j]=0; for (i=1;i<=n;i++) { big=0.0; + irow=icol=0; for (j=1;j<=n;j++) if (ipiv[j] != 1) for (k=1;k<=n;k++) { @@ -26,6 +27,7 @@ void gaussj(double **a, int n, double **b, int m) } } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1"); } + if (icol == 0) nrerror("gaussj: matrix is not finite"); ++(ipiv[icol]); if (irow != icol) { for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l]) diff --git a/orbfit/ludcmp.c b/orbfit/ludcmp.c index dfed7b1..7984697 100644 --- a/orbfit/ludcmp.c +++ b/orbfit/ludcmp.c @@ -2,7 +2,8 @@ #include "nrutil.h" #define TINY 1.0e-20; -void ludcmp(double **a, int n, int *indx, double *d) +/* Returns non-zero if a row of a is zero or not finite. */ +int ludcmp(double **a, int n, int *indx, double *d) { int i,imax,j,k; double big,dum,sum,temp; @@ -14,7 +15,10 @@ void ludcmp(double **a, int n, int *indx, double *d) big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; - if (big == 0.0) nrerror("Singular matrix in routine LUDCMP"); + if (big == 0.0 || !isfinite(big)) { + free_dvector(vv,1,n); + return 1; + } vv[i]=1.0/big; } for (j=1;j<=n;j++) { @@ -24,6 +28,7 @@ void ludcmp(double **a, int n, int *indx, double *d) a[i][j]=sum; } big=0.0; + imax=j; for (i=j;i<=n;i++) { sum=a[i][j]; for (k=1;k= MAXOBS) { + fprintf(stderr,"More than %d observations in %s\n", MAXOBS, + fname==NULL ? "input" : fname); + if (fname!=NULL) fclose(fptr); + return(1); + } obs = &(obsarray[*nobs]); + *obs = scanned; (*nobs)++; eq_to_ec(obs->thetax,obs->thetay,&elat,&elon,NULL); @@ -423,7 +432,7 @@ read_abg(char *fname, double **covar) { FILE *fptr; - int i; + int i, status = 1; char inbuff[256]; if (fname==NULL) @@ -436,46 +445,48 @@ read_abg(char *fname, /* Skipping comments, read the a/b/g specs*/ if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { fprintf(stderr,"Data missing from a/b/g/ data file.\n"); - return(1); + goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &(p->a),&(p->adot),&(p->b),&(p->bdot), &(p->g),&(p->gdot)) != 6) { fprintf(stderr,"Error reading a/b/g data\n"); - return(1); + goto done; } for (i=1; i<=6; i++) { if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { fprintf(stderr,"Data missing from a/b/g/ covariance.\n"); - return(1); + goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &covar[i][1],&covar[i][2],&covar[i][3], &covar[i][4],&covar[i][5],&covar[i][6]) != 6) { fprintf(stderr,"Error reading a/b/g covariance\n"); - return(1); + goto done; } } /* Now read the coordinate system info */ if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { fprintf(stderr,"Data missing from a/b/g/ data file.\n"); - return(1); + goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &lat0, &lon0, &xBary, &yBary, &zBary, &jd0) != 6) { fprintf(stderr,"Error reading coord system info\n"); - return(1); + goto done; } lat0 *= DTOR; lon0 *= DTOR; + status = 0; +done: if (fname!=NULL) fclose(fptr); - return(0); + return(status); } /* Take a set of observations and make a preliminary fit using the @@ -541,9 +552,19 @@ prelim_fit(OBSERVATION obsarray[], if (invert_matrix(alpha,covar,5)) { - /* some failure in the inversion...*/ - fprintf(stderr,"Error inverting the alpha matrix\n"); - exit(1); + /* Linear fit is degenerate: flag it with a negative gamma variance so + * fit_observations falls back to the gdot=0, energy-constrained fit. */ + for (i=1; i<=6; i++) + for (j=1; j<=6; j++) covar[i][j]=0.; + covar[5][5] = -1.; + pout->a = pout->adot = pout->b = pout->bdot = pout->gdot = 0.; + pout->g = 0.03; + free_dvector(dx,1,6); + free_dvector(dy,1,6); + free_dvector(soln,1,6); + free_dvector(beta,1,6); + free_dmatrix(alpha,1,6,1,6); + return; } /* Now multiply matrices to get the solution vector */ @@ -619,23 +640,28 @@ int invert_matrix(double **in, double **out, int dim) { - int *indx,i,j; + int *indx,i,j,status=0; double *tvec,det; tvec = dvector(1,dim); indx = ivector(1,dim); - ludcmp(in,dim,indx,&det); - - for (j=1; j<=dim; j++) { - for (i=1; i<=dim; i++) tvec[i]=0.; - tvec[j] = 1.0; - lubksb(in,dim,indx,tvec); - for (i=1; i<=dim; i++) out[i][j]=tvec[i]; + if (ludcmp(in,dim,indx,&det)) { + status = 1; + } else { + for (j=1; j<=dim; j++) { + for (i=1; i<=dim; i++) tvec[i]=0.; + tvec[j] = 1.0; + lubksb(in,dim,indx,tvec); + for (i=1; i<=dim; i++) { + out[i][j]=tvec[i]; + if (!isfinite(tvec[i])) status = 1; + } + } } - free_ivector(indx,1,6); - free_dvector(tvec,1,6); - return(0); + free_ivector(indx,1,dim); + free_dvector(tvec,1,dim); + return(status); } void @@ -967,7 +993,7 @@ kbo3d(PBASIS *pin, vout[i] = v[i] + a[i]*(t-tv); } /* x and y derivatives w.r.t. parameters - inertial approx ok here*/ - if (dx== NULL || dy==NULL || dx==NULL) return; + if (dx== NULL || dy==NULL || dz==NULL) return; for (i=1; i<=6; i++) dx[i]=dy[i]=dz[i]=0.; dx[1] = dy[3] = z0; diff --git a/orbfit/orbfit2.c b/orbfit/orbfit2.c index ee0fbf4..cef2a50 100644 --- a/orbfit/orbfit2.c +++ b/orbfit/orbfit2.c @@ -38,7 +38,8 @@ mrqfit(OBSERVATION *obsarray, * covar should be 6x6 dmatrix upon entry, will hold covariance on exit. * Info on the fit is sent to logfile if it is non-NULL. */ -/* Return value is 5 if energy constraint is used, 6 otherwise */ +/* Return value is the number of fitted parameters: 4 if gdot is fixed, + * 5 if the energy constraint is used, 6 otherwise; -1 if nobs < 2. */ int fit_observations(OBSERVATION obsarray[], int nobs, @@ -49,18 +50,20 @@ fit_observations(OBSERVATION obsarray[], FILE *logfile) { int *ia,i; - double *a; double gbind2; double energy_fit; /*weight given to binding-energy constraint*/ int fitparms=6; - a = dvector(1,6); + + if (nobs<2) { + fprintf(stderr,"ERROR: not enough observations nobs=%d\n",nobs); + return -1; + } + ia = ivector(1,6); energy_fit = 0.; /* Don't even try the no-energy fit unless there are 3 observations*/ - if (nobs<2) { - fprintf(stderr,"ERROR: not enough observations nobs=%d\n",nobs); - } else if (nobs==2) { + if (nobs==2) { fitparms=4; } else { prelim_fit(obsarray,nobs,p,covar); @@ -140,7 +143,6 @@ fit_observations(OBSERVATION obsarray[], mrqfit(obsarray, nobs, p, ia, covar, chisq, energy_fit); } - free_dvector(a,1,6); free_ivector(ia,1,6); if (fitparms==4) { @@ -213,6 +215,8 @@ mrqfit(OBSERVATION *obsarray, p->g = a[5]; p->gdot = a[6]; + free_dvector(a,1,6); + free_dmatrix(alpha,1,6,1,6); return; } diff --git a/orbfit/orbfitmodule.c b/orbfit/orbfitmodule.c index 50674b8..09c72b4 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -4,90 +4,77 @@ #include #include -double *fitradec(char *mpc_filename, char *abg_filename) +static void +fill_nan(double *result, int n) { + int i; + for (i = 0; i < n; i++) result[i] = NAN; +} - FILE *abg_file ; - - OBSERVATION obsarray[MAXOBS]; +double *fitradec(char *mpc_filename, char *abg_filename) +{ + static double result[2]; + FILE *abg_file = NULL; + OBSERVATION *obsarray = NULL; int nobs; - - ORBIT orbit; PBASIS p; - XVBASIS xv; - double d, dd; - static double result[2]; - double **covar; + double **covar = NULL; double chisq; int dof; - covar = dmatrix(1,6,1,6); + fill_nan(result, 2); + obsarray = malloc(MAXOBS*sizeof(OBSERVATION)); + if (obsarray == NULL) { + fprintf(stderr, "Out of memory for observations\n"); + goto cleanup; + } + covar = dmatrix(1,6,1,6); if (read_radec(obsarray, mpc_filename, &nobs)) { fprintf(stderr, "Error reading input observations\n"); - return result ; + goto cleanup; } - /* - for (i=0;i 1.) return 1.; + if (x < -1.) return -1.; + return x; +} + void orbitElements(XVBASIS *xv, ORBIT *orb) @@ -366,38 +378,54 @@ orbitElements(XVBASIS *xv, hMagnitude = sqrt( angularMomentum[1]*angularMomentum[1] + angularMomentum[2]*angularMomentum[2] + angularMomentum[3]*angularMomentum[3] ); - inclination = acos(angularMomentum[3]/hMagnitude); /* in radians here */ + inclination = acos(clamp_unit(angularMomentum[3]/hMagnitude)); /* in radians here */ ascendingNodeMagnitude = sqrt(ascendingNode[1]*ascendingNode[1] + ascendingNode[2]*ascendingNode[2] + ascendingNode[3]*ascendingNode[3]); - longitudeOfAscendingNode = acos(ascendingNode[1]/ascendingNodeMagnitude); + if (ascendingNodeMagnitude <= ELEMENTS_TINY*hMagnitude) { + /* Orbit in the ecliptic: the node is undefined, so measure from +x */ + ascendingNode[1] = 1.; + ascendingNode[2] = ascendingNode[3] = 0.; + ascendingNodeMagnitude = 1.; + } + longitudeOfAscendingNode = acos(clamp_unit(ascendingNode[1]/ascendingNodeMagnitude)); /* Capital Omega in radians here */ if (ascendingNode[2] < 0) longitudeOfAscendingNode = 2*PI - longitudeOfAscendingNode; - /* ???could use atan2 here?? */ ascEccDotProduct = ascendingNode[1]*eccentricityVector[1] + ascendingNode[2]*eccentricityVector[2] + ascendingNode[3]*eccentricityVector[3]; - argumentOfPerifocus = acos(ascEccDotProduct/ - (ascendingNodeMagnitude*eccentricity)); - /* Small omega in radians here */ - if (eccentricityVector[3] < 0) argumentOfPerifocus = - 2*PI - argumentOfPerifocus; - xBar = (semiLatusRectum - rMagnitude)/eccentricity; - yBar = radVelDotProduct*sqrt(semiLatusRectum/combinedMass)/eccentricity; - - /* From here, we assume that the motion is elliptical */ - - cosE = (xBar/semimajor) + eccentricity; - sinE = yBar/(semimajor*sqrt(1-eccentricity*eccentricity)); - /* where semimajor*sqrt(1-eccentricity*eccentricity) is semiminor */ - eccentricAnomaly = atan2(sinE,cosE); - - meanAnomaly = eccentricAnomaly - eccentricity*sinE; /* radians */ - meanMotion = sqrt(combinedMass/(pow(semimajor,3.))); - timeOfPerifocalPassage = epochTime - meanAnomaly/meanMotion/DAY; - /* This comes from M=n(t-T) where t is epoch time and T is time of perifocal - passage, in days */ + if (eccentricity <= ELEMENTS_TINY) { + /* Circular orbit: perihelion is undefined, put it at the node */ + argumentOfPerifocus = 0.; + } else { + argumentOfPerifocus = acos(clamp_unit(ascEccDotProduct/ + (ascendingNodeMagnitude*eccentricity))); + /* Small omega in radians here */ + if (eccentricityVector[3] < 0) argumentOfPerifocus = + 2*PI - argumentOfPerifocus; + } + + if (eccentricity >= 1. || semimajor <= 0.) { + /* Unbound: no time of perihelion from the elliptical relations */ + timeOfPerifocalPassage = NAN; + } else if (eccentricity <= ELEMENTS_TINY) { + timeOfPerifocalPassage = epochTime; + } else { + xBar = (semiLatusRectum - rMagnitude)/eccentricity; + yBar = radVelDotProduct*sqrt(semiLatusRectum/combinedMass)/eccentricity; + + cosE = (xBar/semimajor) + eccentricity; + sinE = yBar/(semimajor*sqrt(1-eccentricity*eccentricity)); + /* where semimajor*sqrt(1-eccentricity*eccentricity) is semiminor */ + eccentricAnomaly = atan2(sinE,cosE); + + meanAnomaly = eccentricAnomaly - eccentricity*sinE; /* radians */ + meanMotion = sqrt(combinedMass/(pow(semimajor,3.))); + timeOfPerifocalPassage = epochTime - meanAnomaly/meanMotion/DAY; + /* This comes from M=n(t-T) where t is epoch time and T is time of perifocal + passage, in days */ + } orb->a=semimajor; orb->e=eccentricity; From 66f2ba46cf399f0aa6c9425df3e0ec548be03f8a Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:29:33 -0700 Subject: [PATCH 5/7] Report library errors to Python instead of exiting; serialize calls The C library no longer calls exit(). Errors deep in a calculation (out-of-range ephemeris dates, Numerical Recipes failures, bad observatory files, unbound orbits in elements_to_xv) call orbfit_fail(), which longjmps back to a guard in each entry point. The guard frees all Numerical Recipes allocations (now tracked in a list), resets the orbit integrator, and returns NaN results. Input errors are recorded with orbfit_error() and propagated through return codes. The message is available from orbfit_last_error(). The ephemeris reader checks its reads, can be retried after a failed initialization, and seeks to records by absolute position. bk_orbit.py loads the library once, raises BKOrbitError with the library's message, and holds a module lock around every call because the library keeps its state in static storage. Remove the unused read_options() command-line parser. Co-authored-by: Cursor --- orbfit/dms.c | 3 +- orbfit/ephem_earth.c | 317 +++++++++++++++++++-------------------- orbfit/nrutil.c | 84 +++++++---- orbfit/nrutil.h | 11 +- orbfit/orbfit.h | 28 +++- orbfit/orbfit1.c | 113 +++++--------- orbfit/orbfit2.c | 2 +- orbfit/orbfit_api.h | 8 +- orbfit/orbfit_error.c | 67 +++++++++ orbfit/orbfitmodule.c | 123 ++++++++------- orbfit/transforms.c | 13 +- src/mp_ephem/bk_orbit.py | 133 ++++++++++------ test/c/orbfit_driver.c | 31 +++- test/test_errors.py | 73 +++++++++ 14 files changed, 609 insertions(+), 397 deletions(-) create mode 100644 orbfit/orbfit_error.c create mode 100644 test/test_errors.py diff --git a/orbfit/dms.c b/orbfit/dms.c index d3a1b9d..b97fd58 100644 --- a/orbfit/dms.c +++ b/orbfit/dms.c @@ -35,8 +35,7 @@ dmsdeg(const char *string) j = sscanf(string+i,"%31[-+.0-9]",substring); if (j==EOF || *(string+i)==0) break; else if (j!=1) { - fprintf(stderr,"Error interpeting dms ->%s<-\n",string); - return(0.); + orbfit_fail("Error interpreting dms ->%s<-",string); } args[iarg] = atof(substring); i += strlen(substring); diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index 974044c..6f048ec 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -91,6 +91,10 @@ set_ephem_file(char *fname) { static recOneType R1; static FILE *Ephemeris_File; static double Coeff_Array[ARRAY_SIZE] , T_beg , T_end , T_span; + static double T_first; /* start of the first data record */ + static double Record_span; /* days covered by each record */ + static long N_records; /* data records in the file */ + static int Ephemeris_ready = 0; static int Debug = FALSE; /* Generates detailed output if true */ @@ -106,57 +110,45 @@ set_ephem_file(char *fname) { /** **/ /**==========================================================================**/ -void Read_Coefficients( double Time ) +static void Read_Coefficients( double Time ) { - double T_delta = 0.0; - long Offset = 0 ; /*** ??? change to long 8/9/99 ***/ - - /*--------------------------------------------------------------------------*/ - /* Find ephemeris data that record contains input time. Note that one, and */ - /* only one, of the following conditional statements will be true (if both */ - /* were false, this function would not have been called). */ - /*--------------------------------------------------------------------------*/ - - if ( Time < T_beg ) /* Compute backwards location offset */ - { - T_delta = T_beg - Time; - Offset = (int) -ceil(T_delta/T_span); /***Needed negative sign here???*/ - } - - if ( Time > T_end ) /* Compute forewards location offset */ - { - T_delta = Time - T_end; - Offset = (int) ceil(T_delta/T_span); - } - - /*--------------------------------------------------------------------------*/ - /* Retrieve ephemeris data from new record. */ - /*--------------------------------------------------------------------------*/ - - fseek(Ephemeris_File,(Offset-1)*ARRAY_SIZE*sizeof(double),SEEK_CUR); - fread(&Coeff_Array,sizeof(double),ARRAY_SIZE,Ephemeris_File); - + double buffer[ARRAY_SIZE]; + long record; + long record_bytes = (long) (ARRAY_SIZE*sizeof(double)); + + if (!Ephemeris_ready) + orbfit_fail("Ephemeris file has not been initialized"); + + /* Records are of equal length and follow the two header records. */ + if (!(Time >= T_first && Time <= T_first + N_records*Record_span)) + orbfit_fail("JD %.2f is out of range of ephemeris file", Time); + record = (long) floor((Time - T_first)/Record_span); + if (record >= N_records) record = N_records - 1; + + if (fseek(Ephemeris_File, (2+record)*record_bytes, SEEK_SET) != 0 + || fread(buffer,sizeof(double),ARRAY_SIZE,Ephemeris_File) != ARRAY_SIZE) + orbfit_fail("Error reading record %ld of ephemeris file", record); + if (Time < buffer[0] || Time > buffer[1]) + orbfit_fail("JD %.2f is not in ephemeris record %ld (%.2f - %.2f)", + Time, record, buffer[0], buffer[1]); + + memcpy(Coeff_Array, buffer, sizeof buffer); T_beg = Coeff_Array[0]; T_end = Coeff_Array[1]; T_span = T_end - T_beg; - if (Time < T_beg || Time > T_end) { - fprintf(stderr,"JD %.2f is out of range of ephemeris file\n",Time); - exit(1); - } /*--------------------------------------------------------------------------*/ /* Debug print (optional) */ /*--------------------------------------------------------------------------*/ if ( Debug ) { - printf("\n In: Read_Coefficients \n"); - printf("\n ARRAY_SIZE = %4d",ARRAY_SIZE); - printf("\n Offset = %3ld",Offset); - printf("\n T_delta = %7.3f",T_delta); - printf("\n T_Beg = %7.3f",T_beg); - printf("\n T_End = %7.3f",T_end); - printf("\n T_Span = %7.3f\n\n",T_span); + fprintf(stderr,"\n In: Read_Coefficients \n"); + fprintf(stderr,"\n ARRAY_SIZE = %4d",ARRAY_SIZE); + fprintf(stderr,"\n Record = %3ld",record); + fprintf(stderr,"\n T_Beg = %7.3f",T_beg); + fprintf(stderr,"\n T_End = %7.3f",T_end); + fprintf(stderr,"\n T_Span = %7.3f\n\n",T_span); } } @@ -168,9 +160,7 @@ void Read_Coefficients( double Time ) /** ephemeris data. It opens the ephemeris data file, reads the header **/ /** data, loads the first coefficient record into a global array, then **/ /** returns a status code that indicates whether or not all of this was **/ -/** done successfully. **/ -/** **/ -/** Input: A character string giving the name of an ephemeris data file. **/ +/** done successfully. A failed call can be retried. **/ /** **/ /** Returns: An integer status code. **/ /** **/ @@ -179,18 +169,13 @@ void Read_Coefficients( double Time ) ** for the file itself, as environment-specified file or as default in ** this directory **/ -int Initialize_Ephemeris() +int Initialize_Ephemeris(void) { int headerID; + long file_bytes, record_bytes = (long) (ARRAY_SIZE*sizeof(double)); char fileName[FNAMESIZE]; - /*** gmb: don't duplicate this call */ - static int init=0; - if (init) return SUCCESS; - init = 1; - /*--------------------------------------------------------------------------*/ - /* Open ephemeris file. */ - /*--------------------------------------------------------------------------*/ + if (Ephemeris_ready) return SUCCESS; /** use previously specified filename, ** or environment-specified file, or the default filename, in that order @@ -201,68 +186,69 @@ int Initialize_Ephemeris() strncpy(fileName, getenv(EPHEM_ENVIRON), FNAMESIZE-1); else strncpy(fileName, DEFAULT_EPHEM_FILE, FNAMESIZE-1); - - fileName[FNAMESIZE-1]=0; - - Ephemeris_File = fopen(fileName,"r"); - - /*--------------------------------------------------------------------------*/ - /* Read header & first coefficient array, then return status code. */ - /*--------------------------------------------------------------------------*/ - - if ( Ephemeris_File == NULL ) /*........................No need to continue */ + Ephemeris_File = fopen(fileName,"rb"); + if ( Ephemeris_File == NULL ) { - printf("\n Unable to open ephemeris file: %s.\n",fileName); + orbfit_error("Unable to open ephemeris file: %s", fileName); return FAILURE; } - else - { /*.................Read first three header records from ephemeris file */ - - fread(&H1,sizeof(double),ARRAY_SIZE,Ephemeris_File); - fread(&H2,sizeof(double),ARRAY_SIZE,Ephemeris_File); - fread(&Coeff_Array,sizeof(double),ARRAY_SIZE,Ephemeris_File); - - /*...............................Store header data in global variables */ - - R1 = H1.data; - - /*..........................................Set current time variables */ - T_beg = Coeff_Array[0]; - T_end = Coeff_Array[1]; - T_span = T_end - T_beg; + /* Read the two header records and the first coefficient record */ + if (fread(&H1,sizeof(double),ARRAY_SIZE,Ephemeris_File) != ARRAY_SIZE + || fread(&H2,sizeof(double),ARRAY_SIZE,Ephemeris_File) != ARRAY_SIZE + || fread(&Coeff_Array,sizeof(double),ARRAY_SIZE,Ephemeris_File) != ARRAY_SIZE) + { + orbfit_error("Error reading header of ephemeris file: %s", fileName); + goto fail; + } - /*..............................Convert header ephemeris ID to integer */ + R1 = H1.data; + T_beg = Coeff_Array[0]; + T_end = Coeff_Array[1]; + T_span = T_end - T_beg; + headerID = (int) R1.DENUM; - headerID = (int) R1.DENUM; - - /*..............................................Debug Print (optional) */ - - if ( Debug ) - { - printf("\n In: Initialize_Ephemeris \n"); - printf("\n ARRAY_SIZE = %4d",ARRAY_SIZE); - printf("\n headerID = %3d",headerID); - printf("\n T_Beg = %7.3f",T_beg); - printf("\n T_End = %7.3f",T_end); - printf("\n T_Span = %7.3f\n\n",T_span); - } - - /*..................................................Return status code */ - - if ( headerID == EPHEMERIS ) - { - return SUCCESS; - } - else - { - printf("\n Opened wrong file: %s",fileName); - printf(" for ephemeris: %d.\n",EPHEMERIS); - return FAILURE; - } + if ( Debug ) + { + fprintf(stderr,"\n In: Initialize_Ephemeris \n"); + fprintf(stderr,"\n ARRAY_SIZE = %4d",ARRAY_SIZE); + fprintf(stderr,"\n headerID = %3d",headerID); + fprintf(stderr,"\n T_Beg = %7.3f",T_beg); + fprintf(stderr,"\n T_End = %7.3f",T_end); + fprintf(stderr,"\n T_Span = %7.3f\n\n",T_span); + } + + if ( headerID != EPHEMERIS ) + { + orbfit_error("Opened wrong file: %s for ephemeris: %d", + fileName, EPHEMERIS); + goto fail; } + if (!(T_span > 0.)) + { + orbfit_error("Bad record time span in ephemeris file: %s", fileName); + goto fail; + } + if (fseek(Ephemeris_File, 0L, SEEK_END) != 0 + || (file_bytes = ftell(Ephemeris_File)) < 0 + || (N_records = file_bytes/record_bytes - 2) < 1) + { + orbfit_error("Error sizing ephemeris file: %s", fileName); + goto fail; + } + + T_first = T_beg; + Record_span = T_span; + Ephemeris_ready = 1; + return SUCCESS; + +fail: + fclose(Ephemeris_File); + Ephemeris_File = NULL; + T_beg = T_end = T_span = 0.; + return FAILURE; } /**==========================================================================**/ @@ -396,11 +382,8 @@ void Interpolate_Position( double Time , int Target , double Position[3] ) /* This function doesn't "do" nutations or librations. */ /*--------------------------------------------------------------------------*/ - if ( Target >= 11 ) /* Also protects against weird input errors */ - { - printf("\n This function does not compute nutations or librations.\n"); - return; - } + if ( Target < 0 || Target >= 11 ) /* Also protects against weird input */ + orbfit_fail("Ephemeris target %d is not a body position", Target); /*--------------------------------------------------------------------------*/ /* Initialize local coefficient array. */ @@ -429,11 +412,11 @@ void Interpolate_Position( double Time , int Target , double Position[3] ) if ( Debug ) { - printf("\n In: Interpolate_Position\n"); - printf("\n Target = %2d",Target); - printf("\n C = %4ld (before)",C); - printf("\n N = %4ld",N); - printf("\n G = %4ld\n",G); + fprintf(stderr,"\n In: Interpolate_Position\n"); + fprintf(stderr,"\n Target = %2d",Target); + fprintf(stderr,"\n C = %4ld (before)",C); + fprintf(stderr,"\n N = %4ld",N); + fprintf(stderr,"\n G = %4ld\n",G); } /*--------------------------------------------------------------------------*/ @@ -470,26 +453,24 @@ void Interpolate_Position( double Time , int Target , double Position[3] ) for (i=C ; i<(C+3*N) ; i++) A[i-C] = Coeff_Array[i]; } else /* Something has gone terribly wrong */ - { - printf("\n Number of granules must be >= 1: check header data.\n"); - } + orbfit_fail("Number of ephemeris granules must be >= 1: check header data"); /*...................................................Debug print (optional) */ if ( Debug ) { - printf("\n C = %4ld (after)",C); - printf("\n offset = %4ld",offset); - printf("\n Time = %12.7f",Time); - printf("\n T_sub = %12.7f",T_sub); - printf("\n T_seg = %12.7f",T_seg); - printf("\n Tc = %12.7f\n",Tc); - printf("\n Array Coefficients:\n"); + fprintf(stderr,"\n C = %4ld (after)",C); + fprintf(stderr,"\n offset = %4ld",offset); + fprintf(stderr,"\n Time = %12.7f",Time); + fprintf(stderr,"\n T_sub = %12.7f",T_sub); + fprintf(stderr,"\n T_seg = %12.7f",T_seg); + fprintf(stderr,"\n Tc = %12.7f\n",Tc); + fprintf(stderr,"\n Array Coefficients:\n"); for ( i=0 ; i<3*N ; i++ ) { - printf("\n A[%2d] = % 22.15e",i,A[i]); + fprintf(stderr,"\n A[%2d] = % 22.15e",i,A[i]); } - printf("\n\n"); + fprintf(stderr,"\n\n"); } /*..........................................................................*/ @@ -531,11 +512,8 @@ void Interpolate_State(double Time, /* This function doesn't "do" nutations or librations. */ /*--------------------------------------------------------------------------*/ - if ( Target >= 11 ) /* Also protects against weird input errors */ - { - printf("\n This function does not compute nutations or librations.\n"); - return; - } + if ( Target < 0 || Target >= 11 ) /* Also protects against weird input */ + orbfit_fail("Ephemeris target %d is not a body position", Target); /*--------------------------------------------------------------------------*/ /* Initialize local coefficient array. */ @@ -564,11 +542,11 @@ void Interpolate_State(double Time, if ( Debug ) { - printf("\n In: Interpolate_State\n"); - printf("\n Target = %2d",Target); - printf("\n C = %4ld (before)",C); - printf("\n N = %4ld",N); - printf("\n G = %4ld\n",G); + fprintf(stderr,"\n In: Interpolate_State\n"); + fprintf(stderr,"\n Target = %2d",Target); + fprintf(stderr,"\n C = %4ld (before)",C); + fprintf(stderr,"\n N = %4ld",N); + fprintf(stderr,"\n G = %4ld\n",G); } /*--------------------------------------------------------------------------*/ @@ -605,26 +583,24 @@ void Interpolate_State(double Time, for (i=C ; i<(C+3*N) ; i++) A[i-C] = Coeff_Array[i]; } else /* Something has gone terribly wrong */ - { - printf("\n Number of granules must be >= 1: check header data.\n"); - } + orbfit_fail("Number of ephemeris granules must be >= 1: check header data"); /*...................................................Debug print (optional) */ if ( Debug ) { - printf("\n C = %4ld (after)",C); - printf("\n offset = %4ld",offset); - printf("\n Time = %12.7f",Time); - printf("\n T_sub = %12.7f",T_sub); - printf("\n T_seg = %12.7f",T_seg); - printf("\n Tc = %12.7f\n",Tc); - printf("\n Array Coefficients:\n"); + fprintf(stderr,"\n C = %4ld (after)",C); + fprintf(stderr,"\n offset = %4ld",offset); + fprintf(stderr,"\n Time = %12.7f",Time); + fprintf(stderr,"\n T_sub = %12.7f",T_sub); + fprintf(stderr,"\n T_seg = %12.7f",T_seg); + fprintf(stderr,"\n Tc = %12.7f\n",Tc); + fprintf(stderr,"\n Array Coefficients:\n"); for ( i=0 ; i<3*N ; i++ ) { - printf("\n A[%2d] = % 22.15e",i,A[i]); + fprintf(stderr,"\n A[%2d] = % 22.15e",i,A[i]); } - printf("\n\n"); + fprintf(stderr,"\n\n"); } /*..........................................................................*/ @@ -678,7 +654,7 @@ geocenter_ssbary(double jd, int i; if (!init) { - if (Initialize_Ephemeris()) exit(1); + if (Initialize_Ephemeris()) orbfit_rethrow(); init = 1; } @@ -924,10 +900,8 @@ read_observatories(char *fname) fileName[FNAMESIZE-1]=0; /* fprintf(stderr,"Loading from %s\n", fileName); */ - if ((sitefile = fopen(fileName,"r"))==NULL) { - fprintf(stderr,"Error opening observatories file %s\n",fileName); - exit(1); - } + if ((sitefile = fopen(fileName,"r"))==NULL) + orbfit_fail("Error opening observatories file %s",fileName); while (fgets_nocomment(inbuff, BUFFSIZE-1, sitefile, NULL)!=NULL) { SITE *sss; @@ -938,17 +912,19 @@ read_observatories(char *fname) copy_columns(inbuff, 0, 3, code); if ((obscode = obscode_from_string(code))==OBSCODE_INVALID) continue; if (strchr(inbuff, ':')!=NULL && strchr(inbuff, ':') < inbuff+30) { - fprintf(stderr,"Observatories file %s is not in the MPC ObsCodes format:\n->%s\n", - fileName, inbuff); - exit(1); + orbfit_error("Observatories file %s is not in the MPC ObsCodes format:\n->%s", + fileName, inbuff); + goto fail; } if (nsites>=maxsites) { maxsites = maxsites ? 2*maxsites : 1024; - if ((sitelist = realloc(sitelist, maxsites*sizeof(SITE)))==NULL) { - fprintf(stderr,"Out of memory reading observatories file %s\n", fileName); - exit(1); + SITE *grown = realloc(sitelist, maxsites*sizeof(SITE)); + if (grown==NULL) { + orbfit_error("Out of memory reading observatories file %s", fileName); + goto fail; } + sitelist = grown; } sss = &(sitelist[nsites]); sss->code = obscode; @@ -967,9 +943,9 @@ read_observatories(char *fname) sss->space = 1; sss->lon = sss->rhocos = sss->rhosin = 0.; } else { - fprintf(stderr,"Bad line in observatories file %s:\n->%s\n", - fileName, inbuff); - exit(1); + orbfit_error("Bad line in observatories file %s:\n->%s", + fileName, inbuff); + goto fail; } copy_columns(inbuff, 30, 30+79, sss->name); @@ -978,11 +954,20 @@ read_observatories(char *fname) nsites++; } fclose(sitefile); - + if (nsites<1) { - fprintf(stderr,"Error: no observatory sites found\n"); - exit(1); + free(sitelist); + sitelist = NULL; + orbfit_fail("No observatory sites found in %s", fileName); } + return; + +fail: + fclose(sitefile); + nsites = 0; + free(sitelist); + sitelist = NULL; + orbfit_rethrow(); } /* Get arbitrary planetary barycenter position. Get the @@ -998,7 +983,7 @@ bodycenter_ssbary(double jd, int i; if (!init) { - if (Initialize_Ephemeris()) exit(1); + if (Initialize_Ephemeris()) orbfit_rethrow(); init = 1; } diff --git a/orbfit/nrutil.c b/orbfit/nrutil.c index de6c488..fcff828 100644 --- a/orbfit/nrutil.c +++ b/orbfit/nrutil.c @@ -1,40 +1,71 @@ +#include #include #include #include "nrutil.h" +#include "orbfit.h" /* Numerical Recipes style vectors and matrices indexed v[nl..nh]. Element 0 * up to nl-1 is allocated but unused, which avoids forming a pointer before - * the start of the allocation; only nl >= 0 is supported. */ + * the start of the allocation; only nl >= 0 is supported. + * + * Every block is kept on a list of live allocations so nr_free_all() can + * release them after orbfit_fail() unwinds out of a calculation. */ + +typedef union nrblock { + struct { + union nrblock *prev, *next; + } link; + max_align_t align; +} nrblock; + +static nrblock live = {{&live, &live}}; + +void *nr_malloc(size_t n) +{ + nrblock *b = malloc(sizeof(nrblock) + n); + + if (!b) nrerror("allocation failure"); + b->link.prev = &live; + b->link.next = live.link.next; + live.link.next->link.prev = b; + live.link.next = b; + return b + 1; +} + +void nr_free(void *p) +{ + nrblock *b; + + if (p == NULL) return; + b = (nrblock *) p - 1; + b->link.prev->link.next = b->link.next; + b->link.next->link.prev = b->link.prev; + free(b); +} + +void nr_free_all(void) +{ + while (live.link.next != &live) nr_free(live.link.next + 1); +} void nrerror(const char *error_text) /* Numerical Recipes standard error handler */ { - fprintf(stderr,"Numerical Recipes run-time error...\n"); - fprintf(stderr,"%s\n",error_text); - fprintf(stderr,"...now exiting to system...\n"); - exit(1); + orbfit_fail("Numerical Recipes run-time error: %s", error_text); } int *ivector(int nl, int nh) /* allocate an int vector with subscript range v[nl..nh] */ { - int *v; - if (nl < 0 || nh < nl) nrerror("bad range in ivector()"); - v=(int *)malloc((size_t) (nh+1)*sizeof(int)); - if (!v) nrerror("allocation failure in ivector()"); - return v; + return (int *) nr_malloc((size_t) (nh+1)*sizeof(int)); } double *dvector(int nl, int nh) /* allocate a double vector with subscript range v[nl..nh] */ { - double *v; - if (nl < 0 || nh < nl) nrerror("bad range in dvector()"); - v=(double *)malloc((size_t) (nh+1)*sizeof(double)); - if (!v) nrerror("allocation failure in dvector()"); - return v; + return (double *) nr_malloc((size_t) (nh+1)*sizeof(double)); } double **dmatrix(int nrl, int nrh, int ncl, int nch) @@ -46,16 +77,11 @@ double **dmatrix(int nrl, int nrh, int ncl, int nch) if (nrl < 0 || nrh < nrl || ncl < 0 || nch < ncl) nrerror("bad range in dmatrix()"); - /* allocate pointers to rows */ - m=(double **) calloc((size_t) (nrh+1), sizeof(double*)); - if (!m) nrerror("allocation failure 1 in dmatrix()"); - - /* allocate rows and set pointers to them */ - for(i=nrl;i<=nrh;i++) { - m[i]=(double *) malloc((size_t) (nch+1)*sizeof(double)); - if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); - } - /* return pointer to array of pointers to rows */ + /* allocate pointers to rows, then the rows */ + m=(double **) nr_malloc((size_t) (nrh+1)*sizeof(double*)); + for(i=0;i=nrl;i--) free(m[i]); - free(m); + for(i=nrh;i>=nrl;i--) nr_free(m[i]); + nr_free(m); } diff --git a/orbfit/nrutil.h b/orbfit/nrutil.h index e145e20..23b24af 100644 --- a/orbfit/nrutil.h +++ b/orbfit/nrutil.h @@ -1,5 +1,10 @@ #ifndef NRUTIL_H #define NRUTIL_H +#include + +void *nr_malloc(size_t n); +void nr_free(void *p); +void nr_free_all(void); double *dvector(int, int); double **dmatrix(int, int, int, int); @@ -7,7 +12,11 @@ int *ivector(int, int); void free_dvector(double *, int, int); void free_ivector(int *, int, int); void free_dmatrix(double **, int, int, int, int); -void nrerror(const char *); +void nrerror(const char *) +#ifdef __GNUC__ + __attribute__((noreturn)) +#endif + ; void lubksb(double **, int, int *, double[]); int ludcmp(double **a, int n, int *indx, double *d); diff --git a/orbfit/orbfit.h b/orbfit/orbfit.h index 1596348..dcb03c9 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -239,12 +239,6 @@ set_ephem_file(char *fname); void set_observatory_file(char *fname); -/* Read these options */ -int -read_options(int* iarg, - int argc, - char *argv[]); - /* preliminary fit to observations */ void prelim_fit(OBSERVATION obsarray[], @@ -382,4 +376,26 @@ void geo_to_ssbary(double jd, double *x, double *y, double *z); * stopping at the end of the line. */ void copy_columns(const char *line, int start, int end, char *out); +/* Error handling (orbfit_error.c). orbfit_error() records and prints a + * message. orbfit_fail() does the same and then longjmps to the guard + * installed by orbfit_begin(), or exits if there is none; orbfit_rethrow() + * jumps with the message already recorded. An entry point calls setjmp() + * on its guard, then orbfit_begin(), and orbfit_end() before returning. */ +#include +#ifdef __GNUC__ +#define ORBFIT_PRINTF __attribute__((format(printf, 1, 2))) +#define ORBFIT_NORETURN __attribute__((noreturn)) +#else +#define ORBFIT_PRINTF +#define ORBFIT_NORETURN +#endif +void orbfit_error(const char *fmt, ...) ORBFIT_PRINTF; +void orbfit_fail(const char *fmt, ...) ORBFIT_PRINTF ORBFIT_NORETURN; +void orbfit_rethrow(void) ORBFIT_NORETURN; +void orbfit_begin(jmp_buf *guard); +void orbfit_end(void); + +/* Forget integrator state left by an interrupted calculation */ +void kbo3d_reset(void); + #endif diff --git a/orbfit/orbfit1.c b/orbfit/orbfit1.c index 7ed0a6a..1485179 100644 --- a/orbfit/orbfit1.c +++ b/orbfit/orbfit1.c @@ -199,6 +199,9 @@ fgets_nocomment(char *inbuff, int length, * part of OBSERVATION structure, ra & dec in the x & y parts. Format * error returns non-zero. */ +/* Length and text of a line without its newline, for "%.*s" */ +#define LINE(s) (int) strcspn((s), "\r\n"), (s) + /* Read a signed coordinate from the second line of an MPC satellite * observation: sign in column sign_col, value in the next 11 columns. */ static int @@ -231,18 +234,18 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) * AU (column 33 = 2). */ double x, y, z, scale; if (previous == NULL) { - fprintf(stderr,"Satellite position line without an observation:\n ->%s\n",inbuff); + orbfit_error("Satellite position line without an observation:\n ->%.*s", LINE(inbuff)); return(1); } if (inbuff[32] == '1') scale = 1.; else if (inbuff[32] == '2') scale = AU_KM; else { - fprintf(stderr,"Unknown units in satellite position line:\n ->%s\n",inbuff); + orbfit_error("Unknown units in satellite position line:\n ->%.*s", LINE(inbuff)); return(1); } if (mpc_signed_value(inbuff, 34, &x) || mpc_signed_value(inbuff, 46, &y) || mpc_signed_value(inbuff, 58, &z)) { - fprintf(stderr,"Format error in satellite position line:\n ->%s\n",inbuff); + orbfit_error("Format error in satellite position line:\n ->%.*s", LINE(inbuff)); return(1); } previous->xe = x*scale; @@ -253,7 +256,7 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) copy_columns(inbuff, 15, 32, field); if (sscanf(field,"%d %d %lf",&(dd.y),&(dd.mo),&(dd.d))!=3) { - fprintf(stderr,"Format error in observation file:\n ->%s\n",inbuff); + orbfit_error("Format error in observation file:\n ->%.*s", LINE(inbuff)); return(1); } dd.h = dd.mn = dd.s = 0.; @@ -262,7 +265,7 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) copy_columns(inbuff, 44, 56, decstring); copy_columns(inbuff, 77, 80, field); if (sscanf(field,"%15s",codestring) != 1) { - fprintf(stderr,"Missing observatory code:\n ->%s\n",inbuff); + orbfit_error("Missing observatory code:\n ->%.*s", LINE(inbuff)); return(1); } obs->obscode = obscode_from_string(codestring); @@ -275,7 +278,7 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) &(dd.y), &(dd.mo), &(dd.d), rastring, decstring, &(obs->dthetay), codestring) != 7) { - fprintf(stderr, "Format error in observation file:\n ->%s\n", inbuff); + orbfit_error("Format error in observation file:\n ->%.*s", LINE(inbuff)); return (1); } obs->obscode = obscode_from_string(codestring); @@ -288,7 +291,7 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) &jd, rastring, decstring, &(obs->dthetay), codestring) != 5) { - fprintf(stderr, "Format error in observation file:\n ->%s\n", inbuff); + orbfit_error("Format error in observation file:\n ->%.*s", LINE(inbuff)); return (1); } obs->obscode = obscode_from_string(codestring); @@ -335,7 +338,7 @@ read_radec(OBSERVATION obsarray[], char *fname, int *nobs) if (fname==NULL) fptr = stdin; else if ( (fptr=fopen(fname,"r"))==NULL) { - fprintf(stderr,"Error opening observations file %s\n",fname); + orbfit_error("Error opening observations file %s",fname); return(1); } @@ -353,14 +356,13 @@ read_radec(OBSERVATION obsarray[], char *fname, int *nobs) // all other non-zero status values indicate an error. if ( scan_status_flag == 1) { - fprintf(stderr,"Quitting on format error\n"); if (fname!=NULL) fclose(fptr); return(1); } if (*nobs >= MAXOBS) { - fprintf(stderr,"More than %d observations in %s\n", MAXOBS, - fname==NULL ? "input" : fname); + orbfit_error("More than %d observations in %s", MAXOBS, + fname==NULL ? "input" : fname); if (fname!=NULL) fclose(fptr); return(1); } @@ -438,46 +440,47 @@ read_abg(char *fname, if (fname==NULL) fptr = stdin; else if ( (fptr=fopen(fname,"r"))==NULL) { - fprintf(stderr,"Error opening a/b/g orbit file %s\n",fname); + orbfit_error("Error opening a/b/g orbit file %s",fname); return(1); } + if (fname==NULL) fname = "input"; /* Skipping comments, read the a/b/g specs*/ if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { - fprintf(stderr,"Data missing from a/b/g/ data file.\n"); + orbfit_error("Data missing from a/b/g data file %s", fname); goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &(p->a),&(p->adot),&(p->b),&(p->bdot), &(p->g),&(p->gdot)) != 6) { - fprintf(stderr,"Error reading a/b/g data\n"); + orbfit_error("Error reading a/b/g data from %s", fname); goto done; } for (i=1; i<=6; i++) { if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { - fprintf(stderr,"Data missing from a/b/g/ covariance.\n"); + orbfit_error("Data missing from a/b/g covariance in %s", fname); goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &covar[i][1],&covar[i][2],&covar[i][3], &covar[i][4],&covar[i][5],&covar[i][6]) != 6) { - fprintf(stderr,"Error reading a/b/g covariance\n"); + orbfit_error("Error reading a/b/g covariance from %s", fname); goto done; } } /* Now read the coordinate system info */ if (fgets_nocomment(inbuff,255,fptr,NULL)==NULL) { - fprintf(stderr,"Data missing from a/b/g/ data file.\n"); + orbfit_error("Coordinate system missing from a/b/g data file %s", fname); goto done; } if (sscanf(inbuff, "%lf %lf %lf %lf %lf %lf", &lat0, &lon0, &xBary, &yBary, &zBary, &jd0) != 6) { - fprintf(stderr,"Error reading coord system info\n"); + orbfit_error("Error reading coordinate system from %s", fname); goto done; } lat0 *= DTOR; @@ -485,7 +488,7 @@ read_abg(char *fname, status = 0; done: - if (fname!=NULL) fclose(fptr); + if (fptr!=stdin) fclose(fptr); return(status); } @@ -813,9 +816,8 @@ matrix_multiply(double **m1, double **m2, double sum; if (x1!=y2) { - fprintf(stderr,"Trying to multiply mismatched matrices, " - " %d x %d times %d x %d\n", y1,x1,y2,x2); - exit(1); + orbfit_fail("Trying to multiply mismatched matrices, " + " %d x %d times %d x %d", y1,x1,y2,x2); } for (i=1; i<=y1; i++) for (k=1; k<=x2; k++) { @@ -910,6 +912,16 @@ accel_nbody(double *xx, return; } +/* kbo3d keeps the integrator state between calls with the same orbit */ +static int kbo3d_init = 0; + +/* Discard that state, e.g. after a failure part way through a step */ +void +kbo3d_reset(void) +{ + kbo3d_init = 0; +} + /* Give the KBO's 3-d position, along with derivatives. */ /* Here is a leapfrog-integrator version.*/ void @@ -924,13 +936,12 @@ kbo3d(PBASIS *pin, int i, restart=1; static double x[3],v[3],a[3], tx, tv, z0; static PBASIS psave; - static int init=0; static int tdir; double t1,dt,dtv; /* decide whether we need to reset integrator to t=0*/ - if (!init) { - init = 1; + if (!kbo3d_init) { + kbo3d_init = 1; restart = 1; } else if (pin->a==psave.a && pin->adot==psave.adot && pin->b==psave.b && pin->bdot==psave.bdot && @@ -1102,8 +1113,7 @@ zenith_angle(OBSERVATION *obs) &xobs, &yobs, &zobs); r = sqrt(xobs*xobs+yobs*yobs+zobs*zobs); if (r<=0.) { - fprintf(stderr,"Non-positive geocentric radius in zenith_angle()\n"); - exit(1); + orbfit_fail("Non-positive geocentric radius in zenith_angle()"); } xobs /=r; yobs/=r; zobs/=r; /* Rotate this ICRS vector into ecliptic, then projected coords */ @@ -1149,52 +1159,3 @@ fake_observation(PBASIS *p, return; } - -/* Parse command-line looking for the standard options. - * iarg is first arg to examine, and is returned as the - * first non-option argument. - * Returns 1 on parse error. - * currently understands m, j, o, and v options. - */ -#define BUFFSIZE 512 -int -read_options(int* iarg, - int argc, - char *argv[]) -{ - double d; - for ( ; *iarg +#include +#include +#include + +#include "orbfit.h" +#include "orbfit_api.h" + +static jmp_buf *active_guard = NULL; +static char last_error[1024] = ""; + +static void +record(const char *fmt, va_list ap) +{ + vsnprintf(last_error, sizeof last_error, fmt, ap); + fprintf(stderr, "%s\n", last_error); +} + +void +orbfit_error(const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + record(fmt, ap); + va_end(ap); +} + +void +orbfit_rethrow(void) +{ + if (last_error[0] == 0) orbfit_error("orbfit: unspecified error"); + if (active_guard != NULL) longjmp(*active_guard, 1); + exit(1); +} + +void +orbfit_fail(const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + record(fmt, ap); + va_end(ap); + orbfit_rethrow(); +} + +void +orbfit_begin(jmp_buf *guard) +{ + last_error[0] = 0; + active_guard = guard; +} + +void +orbfit_end(void) +{ + active_guard = NULL; +} + +const char * +orbfit_last_error(void) +{ + return last_error; +} diff --git a/orbfit/orbfitmodule.c b/orbfit/orbfitmodule.c index 09c72b4..948e651 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -4,6 +4,22 @@ #include #include +/* Each entry point installs a guard so that orbfit_fail() anywhere below + * returns here: the results are left as NaN, every Numerical Recipes + * allocation still live is released and the integrator state is dropped. */ +#define ORBFIT_GUARD(guard, result, n) \ + do { \ + fill_nan(result, n); \ + if (setjmp(guard)) { \ + nr_free_all(); \ + kbo3d_reset(); \ + fill_nan(result, n); \ + orbfit_end(); \ + return result; \ + } \ + orbfit_begin(&guard); \ + } while (0) + static void fill_nan(double *result, int n) { @@ -11,59 +27,54 @@ fill_nan(double *result, int n) for (i = 0; i < n; i++) result[i] = NAN; } +static void +write_abg(const char *abg_filename, PBASIS *p, double **covar) +{ + FILE *abg_file; + + if ((abg_file = fopen(abg_filename,"w")) == NULL) + orbfit_fail("Error opening a/b/g output file %s", abg_filename); + + fprintf(abg_file, "# Exact a, adot, b, bdot, g, gdot:\n"); + fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n",p->a,p->adot,p->b, + p->bdot, p->g, p->gdot); + + fprintf(abg_file, "# Covariance matrix: \n"); + print_matrix(abg_file, covar, 6, 6); + + /* Print out information on the coordinate system */ + fprintf(abg_file, "# lat0 lon0 xBary yBary zBary JD0\n"); + fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n", + lat0/DTOR,lon0/DTOR,xBary,yBary,zBary,jd0); + + if (fclose(abg_file) != 0) + orbfit_fail("Error writing a/b/g output file %s", abg_filename); +} + double *fitradec(char *mpc_filename, char *abg_filename) { static double result[2]; - FILE *abg_file = NULL; - OBSERVATION *obsarray = NULL; + jmp_buf guard; + OBSERVATION *obsarray; int nobs; PBASIS p; double d, dd; - double **covar = NULL; + double **covar; double chisq; int dof; - fill_nan(result, 2); + ORBFIT_GUARD(guard, result, 2); - obsarray = malloc(MAXOBS*sizeof(OBSERVATION)); - if (obsarray == NULL) { - fprintf(stderr, "Out of memory for observations\n"); - goto cleanup; - } + obsarray = nr_malloc(MAXOBS*sizeof(OBSERVATION)); covar = dmatrix(1,6,1,6); - if (read_radec(obsarray, mpc_filename, &nobs)) { - fprintf(stderr, "Error reading input observations\n"); - goto cleanup; - } + if (read_radec(obsarray, mpc_filename, &nobs)) orbfit_rethrow(); /* Call subroutine to do the actual fitting: */ if (fit_observations(obsarray, nobs, &p, covar, &chisq, &dof, NULL) < 0) - goto cleanup; - - if ((abg_file = fopen(abg_filename,"w")) == NULL) { - fprintf(stderr, "Error opening a/b/g output file %s\n", abg_filename); - goto cleanup; - } - - fprintf(abg_file, "# Exact a, adot, b, bdot, g, gdot:\n"); - fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n",p.a,p.adot,p.b, - p.bdot, p.g, p.gdot); - - fprintf(abg_file, "# Covariance matrix: \n"); - print_matrix(abg_file, covar, 6, 6); + orbfit_rethrow(); - /* Print out information on the coordinate system */ - fprintf(abg_file, "# lat0 lon0 xBary yBary zBary JD0\n"); - fprintf(abg_file, "%.17g %.17g %.17g %.17g %.17g %.17g\n", - lat0/DTOR,lon0/DTOR,xBary,yBary,zBary,jd0); - - if (fclose(abg_file) != 0) { - abg_file = NULL; - fprintf(stderr, "Error writing a/b/g output file %s\n", abg_filename); - goto cleanup; - } - abg_file = NULL; + write_abg(abg_filename, &p, covar); /* Barycentric distance and its uncertainty */ d = sqrt(xBary*xBary + yBary*yBary + pow(zBary-1/p.g,2.)); @@ -71,10 +82,9 @@ double *fitradec(char *mpc_filename, char *abg_filename) result[0] = d; result[1] = dd; -cleanup: - if (abg_file != NULL) fclose(abg_file); free_dmatrix(covar,1,6,1,6); - free(obsarray); + nr_free(obsarray); + orbfit_end(); return result; } @@ -86,18 +96,17 @@ double *fitradec(char *mpc_filename, char *abg_filename) double *predict_helio(char *abg_file, double jdate, int obscode) { static double result[3]; + jmp_buf guard; PBASIS p; OBSERVATION futobs; double xk[3]; double **covar; - fill_nan(result, 3); + ORBFIT_GUARD(guard, result, 3); + covar = dmatrix(1,6,1,6); - if (read_abg(abg_file, &p, covar)) { - fprintf(stderr, "Error input alpha/beta/gamma file %s\n", abg_file); - goto cleanup; - } + if (read_abg(abg_file, &p, covar)) orbfit_rethrow(); /* get observatory code */ futobs.obscode = obscode; @@ -111,8 +120,8 @@ double *predict_helio(char *abg_file, double jdate, int obscode) { result[1] = xk[1]; result[2] = xk[2]; -cleanup: free_dmatrix(covar,1,6,1,6); + orbfit_end(); return result; } @@ -124,6 +133,7 @@ double *predict_helio(char *abg_file, double jdate, int obscode) { double *predict(char *abg_file, double jdate, int obscode) { static double result[8]; + jmp_buf guard; PBASIS p; OBSERVATION futobs; double **covar,**sigxy,a,b,PA,**derivs; @@ -132,17 +142,15 @@ double *predict(char *abg_file, double jdate, int obscode) double xx,yy,xy,bovasqrd,det; double distance; - fill_nan(result, 8); + ORBFIT_GUARD(guard, result, 8); + sigxy = dmatrix(1,2,1,2); derivs = dmatrix(1,2,1,2); covar = dmatrix(1,6,1,6); covecl = dmatrix(1,2,1,2); coveq = dmatrix(1,2,1,2); - if (read_abg(abg_file,&p,covar) ) { - fprintf(stderr, "Error input alpha/beta/gamma file %s\n",abg_file); - goto cleanup; - } + if (read_abg(abg_file,&p,covar)) orbfit_rethrow(); /* get observatory code */ futobs.obscode=obscode; @@ -193,12 +201,12 @@ double *predict(char *abg_file, double jdate, int obscode) result[6] = lon; result[7] = lat; -cleanup: free_dmatrix(sigxy,1,2,1,2); free_dmatrix(derivs,1,2,1,2); free_dmatrix(covar,1,6,1,6); free_dmatrix(covecl,1,2,1,2); free_dmatrix(coveq,1,2,1,2); + orbfit_end(); return result; } @@ -206,22 +214,21 @@ double *predict(char *abg_file, double jdate, int obscode) double *abg_to_aei(char *abg_file) { static double result[15]; + jmp_buf guard; PBASIS p; XVBASIS xv; ORBIT orbit; double d, dd; double **covar_abg, **covar_xyz, **derivs, **covar_aei; - fill_nan(result, 15); + ORBFIT_GUARD(guard, result, 15); + covar_abg = dmatrix(1,6,1,6); covar_xyz = dmatrix(1,6,1,6); covar_aei = dmatrix(1,6,1,6); derivs = dmatrix(1,6,1,6); - if (read_abg(abg_file,&p,covar_abg)) { - fprintf(stderr, "Error in input alpha/beta/gamma file %s\n", abg_file); - goto cleanup; - } + if (read_abg(abg_file,&p,covar_abg)) orbfit_rethrow(); /* Transform the orbit basis and get the deriv. matrix */ pbasis_to_bary(&p, &xv, derivs); @@ -257,10 +264,10 @@ double *abg_to_aei(char *abg_file) result[13] = d; result[14] = dd; -cleanup: free_dmatrix(covar_abg,1,6,1,6); free_dmatrix(covar_xyz,1,6,1,6); free_dmatrix(covar_aei,1,6,1,6); free_dmatrix(derivs,1,6,1,6); + orbfit_end(); return result; } diff --git a/orbfit/transforms.c b/orbfit/transforms.c index fbd30a9..a6616ca 100644 --- a/orbfit/transforms.c +++ b/orbfit/transforms.c @@ -454,7 +454,7 @@ date_to_jd(struct date_time date) short dint; if((date.y <= 1900) | (date.y >= 2100)) { - printf("Date out of range. 1900 - 2100 only.\n"); + orbfit_error("Date year %d out of range. 1900 - 2100 only.", date.y); return(0.); } @@ -493,8 +493,7 @@ elements_to_xv(ORBIT *o, double mu = GM*SSMASS; /*use SSMASS, work in AU/yrs*/ if (o->e >= 1. || o->a <=0.) { - fprintf(stderr,"elements_to_xv only for closed orbits now\n"); - exit(1); + orbfit_fail("elements_to_xv only for closed orbits now"); } /* get the eccentric Anomaly from the mean anomaly */ @@ -513,9 +512,8 @@ elements_to_xv(ORBIT *o, f = x1 - o->e * sin(x1) - meanAnomaly; fmid= x2 - o->e * sin(x2) - meanAnomaly; if (f*fmid > 0.0) { - fprintf(stderr,"Error, eccentricAnomaly root not bracketed\n"); - fprintf(stderr,"f, fmid %f %f\n",f,fmid); - exit(1); + orbfit_fail("Error, eccentricAnomaly root not bracketed: f, fmid %f %f", + f, fmid); } rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2); @@ -526,8 +524,7 @@ elements_to_xv(ORBIT *o, if (fabs(dx) < TOLERANCE || fmid == 0.0) break; } if (j>=JMAX) { - fprintf(stderr,"eccentricAnomaly took too long\n"); - exit(1); + orbfit_fail("eccentricAnomaly took too long"); } meanAnomaly = rtb; #undef JMAX diff --git a/src/mp_ephem/bk_orbit.py b/src/mp_ephem/bk_orbit.py index ac942d3..cc67a11 100644 --- a/src/mp_ephem/bk_orbit.py +++ b/src/mp_ephem/bk_orbit.py @@ -7,6 +7,7 @@ import os import random import tempfile +import threading from astropy import units from astropy.coordinates import SkyCoord from astropy.time import Time @@ -28,9 +29,50 @@ def bytes(str, encoding='utf-8'): class BKOrbitError(Exception): - def __init__(self): - super(BKOrbitError, self).__init__( - "Insufficient observations for an orbit.") + def __init__(self, message="Insufficient observations for an orbit."): + super(BKOrbitError, self).__init__(message) + + +# liborbfit is not reentrant: the fit's coordinate system, the orbit +# integrator and the arrays the entry points return are all static storage +# in the library, and ctypes releases the GIL during a call. Every call into +# the library, and the copy of its result, holds this lock. +_ORBFIT_LOCK = threading.Lock() +_ORBFIT_RESULT_SIZES = {'fitradec': 2, 'abg_to_aei': 15, 'predict': 8, 'predict_helio': 3} +_orbfit_lib = None + + +def _load_orbfit(): + """Load liborbfit once and declare the signatures of its entry points.""" + global _orbfit_lib + with _ORBFIT_LOCK: + if _orbfit_lib is not None: + return _orbfit_lib + path = os.path.dirname(__file__) + lib = ctypes.CDLL(glob.glob(os.path.join(path, 'orbfit*.so'))[0]) + for name, size in _ORBFIT_RESULT_SIZES.items(): + getattr(lib, name).restype = ctypes.POINTER(ctypes.c_double * size) + lib.fitradec.argtypes = [ctypes.c_char_p, ctypes.c_char_p] + lib.abg_to_aei.argtypes = [ctypes.c_char_p] + lib.predict.argtypes = [ctypes.c_char_p, ctypes.c_double, ctypes.c_int] + lib.predict_helio.argtypes = [ctypes.c_char_p, ctypes.c_double, ctypes.c_int] + lib.orbfit_last_error.restype = ctypes.c_char_p + lib.orbfit_last_error.argtypes = [] + _orbfit_lib = lib + return lib + + +def _orbfit_call(lib, name, *args): + """Call a liborbfit entry point and return a copy of its result array. + + :raises BKOrbitError: if the library reports an error (the result is then all NaN). + """ + with _ORBFIT_LOCK: + result = list(getattr(lib, name)(*args).contents) + message = lib.orbfit_last_error() + if message: + raise BKOrbitError("{}: {}".format(name, message.decode('utf-8', errors='replace'))) + return result class BKOrbit(object): @@ -47,8 +89,6 @@ def __init__(self, observations, ast_filename=None, abg_file=None): :rtype : Orbfit """ __PATH__ = os.path.dirname(__file__) - # find the orbit.so library - __ORBFIT_LIB__ = glob.glob(os.path.join(__PATH__, 'orbfit*.so'))[0] # Choose a format of JPL binary ephemeris file is_64bits = ctypes.sizeof(ctypes.c_voidp) == 8 @@ -60,8 +100,7 @@ def __init__(self, observations, ast_filename=None, abg_file=None): os.path.join(__PATH__, 'data', 'observatories.dat')) - liborbfit = os.path.join(__ORBFIT_LIB__) - self.orbfit = ctypes.CDLL(liborbfit) + self.orbfit = _load_orbfit() if ast_filename is None: assert isinstance(observations, tuple) or isinstance(observations, list) or isinstance(observations, @@ -120,8 +159,6 @@ def mag(self): def _fit_radec(self, randomize=False): # noqa: C901 # call fit_radec with mpc file and abgfile - self.orbfit.fitradec.restype = ctypes.POINTER(ctypes.c_double * 2) - self.orbfit.fitradec.argtypes = [ctypes.c_char_p, ctypes.c_char_p] build_abg = True if self.abg_filename is not None and os.access(self.abg_filename, os.R_OK): build_abg = False @@ -173,30 +210,36 @@ def _fit_radec(self, randomize=False): # noqa: C901 else: _abg_file = open(self.abg_filename, 'w+') _abg_file_name = _abg_file.name - _ = self.orbfit.fitradec(ctypes.c_char_p(bytes(_mpc_file.name, 'utf-8')), - ctypes.c_char_p(bytes(_abg_file_name, 'utf-8'))) + try: + _orbfit_call(self.orbfit, 'fitradec', + bytes(_mpc_file.name, 'utf-8'), bytes(_abg_file_name, 'utf-8')) + except BKOrbitError: + _abg_file.close() + raise _abg_file.seek(0) # call abg_to_aei to get elliptical elements and their chi^2 uncertainty. - self.orbfit.abg_to_aei.restype = ctypes.POINTER(ctypes.c_double * 15) - self.orbfit.abg_to_aei.argtypes = [ctypes.c_char_p] - result = self.orbfit.abg_to_aei(ctypes.c_char_p(bytes(_abg_file.name, 'utf-8'))) - self._a = result.contents[0] * units.AU - self._da = result.contents[6] * units.AU - self._e = result.contents[1] * units.dimensionless_unscaled - self._de = result.contents[7] * units.dimensionless_unscaled - self._inc = result.contents[2] * units.degree - self._dinc = result.contents[8] * units.degree - self._Node = result.contents[3] * units.degree - self._dNode = result.contents[9] * units.degree - self._om = result.contents[4] * units.degree - self._dom = result.contents[10] * units.degree - self._T = result.contents[5] * units.day - self._dT = result.contents[11] * units.day - self._epoch = Time(result.contents[12] * units.day, scale='utc', format='jd') - self._distance = result.contents[13] * units.au - self._distance_uncertainty = result.contents[14] * units.au + try: + result = _orbfit_call(self.orbfit, 'abg_to_aei', bytes(_abg_file.name, 'utf-8')) + except BKOrbitError: + _abg_file.close() + raise + self._a = result[0] * units.AU + self._da = result[6] * units.AU + self._e = result[1] * units.dimensionless_unscaled + self._de = result[7] * units.dimensionless_unscaled + self._inc = result[2] * units.degree + self._dinc = result[8] * units.degree + self._Node = result[3] * units.degree + self._dNode = result[9] * units.degree + self._om = result[4] * units.degree + self._dom = result[10] * units.degree + self._T = result[5] * units.day + self._dT = result[11] * units.day + self._epoch = Time(result[12] * units.day, scale='utc', format='jd') + self._distance = result[13] * units.au + self._distance_uncertainty = result[14] * units.au _abg_file.seek(0) self.abg = _abg_file.read() _abg_file.close() @@ -548,12 +591,9 @@ def predict_helio(self, date, obs_code=568, abg_file=None, minimum_delta=None): abg_file.write(bytes(self.abg, 'utf-8')) abg_file.seek(0) - self.orbfit.predict_helio.restype = ctypes.POINTER(ctypes.c_double * 3) - self.orbfit.predict_helio.argtypes = [ctypes.c_char_p, ctypes.c_double, ctypes.c_int] - predict = self.orbfit.predict_helio(ctypes.c_char_p(bytes(abg_file.name, 'utf-8')), - jd, - ctypes.c_int(obscode_to_int(obs_code))) - self.helio = numpy.array((predict.contents[0], predict.contents[1], predict.contents[2])) * units.au + predict = _orbfit_call(self.orbfit, 'predict_helio', bytes(abg_file.name, 'utf-8'), + jd, obscode_to_int(obs_code)) + self.helio = numpy.array(predict) * units.au def predict(self, date, obs_code=568, abg_file=None, minimum_delta=None): """ @@ -593,20 +633,17 @@ def predict(self, date, obs_code=568, abg_file=None, minimum_delta=None): abg_file.write(bytes(self.abg, 'utf-8')) abg_file.seek(0) - self.orbfit.predict.restype = ctypes.POINTER(ctypes.c_double * 8) - self.orbfit.predict.argtypes = [ctypes.c_char_p, ctypes.c_double, ctypes.c_int] - predict = self.orbfit.predict(ctypes.c_char_p(bytes(abg_file.name, 'utf-8')), - jd, - ctypes.c_int(obscode_to_int(obs_code))) + predict = _orbfit_call(self.orbfit, 'predict', bytes(abg_file.name, 'utf-8'), + jd, obscode_to_int(obs_code)) self._coordinate = None - self._ra = predict.contents[0] * units.degree - self._dec = predict.contents[1] * units.degree - self._lon = predict.contents[6] * units.degree - self._lat = predict.contents[7] * units.degree - self._dra = predict.contents[2] * units.arcsec - self._ddec = predict.contents[3] * units.arcsec - self._pa = predict.contents[4] * units.degree - self._distance = predict.contents[5] * units.AU + self._ra = predict[0] * units.degree + self._dec = predict[1] * units.degree + self._lon = predict[6] * units.degree + self._lat = predict[7] * units.degree + self._dra = predict[2] * units.arcsec + self._ddec = predict[3] * units.arcsec + self._pa = predict[4] * units.degree + self._distance = predict[5] * units.AU self._date = str(_date) self._time = _date self._set_coordinates() diff --git a/test/c/orbfit_driver.c b/test/c/orbfit_driver.c index 093fcd0..2bb8c77 100644 --- a/test/c/orbfit_driver.c +++ b/test/c/orbfit_driver.c @@ -4,7 +4,8 @@ * orbfit_driver [repeat] * * ORBIT_EPHEMERIS and ORBIT_OBSERVATORIES must point at the data files. - * Exits non-zero if a returned position or element is not finite. + * Exits non-zero if a returned position or element is not finite, or if a + * call that should fail does not report an error. * Uncertainties may be NaN when the arc leaves the orbit undetermined. */ #include #include @@ -25,9 +26,24 @@ check(const char *what, const double *v, int n, const int *required, int show) for (i = 0; required[i] >= 0; i++) if (!isfinite(v[required[i]])) bad = 1; if (bad) fprintf(stderr, "non-finite value returned by %s\n", what); + if (orbfit_last_error()[0]) { + fprintf(stderr, "%s reported an error: %s\n", what, orbfit_last_error()); + bad = 1; + } return bad; } +/* A failed call must return NaN and leave a message */ +static int +check_failed(const char *what, const double *v) +{ + if (!isnan(v[0]) || orbfit_last_error()[0] == 0) { + fprintf(stderr, "%s did not report a failure\n", what); + return 1; + } + return 0; +} + int main(int argc, char *argv[]) { @@ -54,6 +70,19 @@ main(int argc, char *argv[]) bad |= check("predict", r, 8, predict_req, i == 0); r = predict_helio(argv[2], jd + 365.25, 568); bad |= check("predict_helio", r, 3, helio_req, i == 0); + + /* Failures part way through a calculation must not exit or leak, and + * the next call must work. */ + bad |= check_failed("predict out of ephemeris range", + predict(argv[2], 1.0e7, 568)); + bad |= check_failed("predict_helio out of ephemeris range", + predict_helio(argv[2], 1.0e7, 568)); + bad |= check_failed("abg_to_aei missing file", + abg_to_aei("/nonexistent/orbfit.abg")); + bad |= check_failed("fitradec missing file", + fitradec("/nonexistent/orbfit.mpc", argv[2])); + r = predict(argv[2], jd + 365.25, 568); + bad |= check("predict after failure", r, 8, predict_req, 0); } return bad; } diff --git a/test/test_errors.py b/test/test_errors.py new file mode 100644 index 0000000..e1b3f94 --- /dev/null +++ b/test/test_errors.py @@ -0,0 +1,73 @@ +""" +Failures inside liborbfit raise BKOrbitError instead of ending the process, +and concurrent calls from threads give the same answers as serial ones. +""" +import os +import tempfile +import threading +import unittest + +from astropy.time import Time + +import mp_ephem + +__PATH__ = os.path.dirname(__file__) + + +class LibraryErrors(unittest.TestCase): + + @classmethod + def setUpClass(cls): + observations = mp_ephem.EphemerisReader().read(os.path.join(__PATH__, 'data', 'o3o08.mpc')) + cls.orbit = mp_ephem.BKOrbit(observations) + + def test_date_outside_ephemeris(self): + with self.assertRaisesRegex(mp_ephem.BKOrbitError, 'out of range of ephemeris'): + self.orbit.predict(Time('2300-01-01', scale='utc')) + with self.assertRaisesRegex(mp_ephem.BKOrbitError, 'out of range of ephemeris'): + self.orbit.predict_helio(Time('2300-01-01', scale='utc')) + + def test_recovers_after_failure(self): + date = Time('2016-04-09', scale='utc') + self.orbit.predict(date) + expected = self.orbit.coordinate.ra.degree + with self.assertRaises(mp_ephem.BKOrbitError): + self.orbit.predict(Time('2300-01-01', scale='utc')) + self.orbit.predict(date) + self.assertEqual(self.orbit.coordinate.ra.degree, expected) + + def test_bad_abg_file(self): + with tempfile.NamedTemporaryFile(mode='w', suffix='.abg') as abg: + abg.write('# not an orbit\n1 2 3\n') + abg.flush() + with self.assertRaisesRegex(mp_ephem.BKOrbitError, 'a/b/g'): + mp_ephem.BKOrbit(None, os.path.join(__PATH__, 'data', 'o3o08.mpc'), abg_file=abg.name) + + def test_threads_match_serial(self): + dates = [Time(2457000.5 + 40 * i, format='jd', scale='utc') for i in range(8)] + serial = [] + for date in dates: + self.orbit.predict(date) + serial.append((self.orbit.ra.degree, self.orbit.dec.degree)) + + results = {} + errors = [] + + def work(index): + try: + orbit = mp_ephem.BKOrbit(None, os.path.join(__PATH__, 'data', 'o3o08.mpc')) + for _ in range(3): + orbit.predict(dates[index]) + results[index] = (orbit.ra.degree, orbit.dec.degree) + except Exception as ex: # pragma: no cover - reported below + errors.append(ex) + + threads = [threading.Thread(target=work, args=(i,)) for i in range(len(dates))] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + self.assertEqual(errors, []) + for index, (ra, dec) in enumerate(serial): + self.assertAlmostEqual(results[index][0], ra, delta=1e-9) + self.assertAlmostEqual(results[index][1], dec, delta=1e-9) From 7746ce0d9b2057552cb407f3bf9b72e1c92fb8ed Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:31:57 -0700 Subject: [PATCH 6/7] Export only the public entry points of the shared library Build with -fvisibility=hidden and mark the entry points used through ctypes (and utc_jd_to_tt, used by the tests) with ORBFIT_API in orbfit_api.h, so the generic Numerical Recipes names and other internals cannot collide with other libraries loaded in the Python process. Make file-local helpers and state static, remove the unused accel() and write_obs(), and send the ephemeris debug output to stderr. Co-authored-by: Cursor --- orbfit/ephem_earth.c | 14 +++++++------- orbfit/orbfit.h | 14 +------------- orbfit/orbfit1.c | 38 ++------------------------------------ orbfit/orbfit2.c | 4 ++-- orbfit/orbfit_api.h | 31 +++++++++++++++++++++++++------ orbfit/transforms.c | 4 ++-- setup.py | 3 ++- test/test_errors.py | 15 ++++++++++++++- 8 files changed, 55 insertions(+), 68 deletions(-) diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index 6f048ec..07c98b8 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -30,14 +30,14 @@ typedef struct { } SITE; static SITE *sitelist=NULL; -char observatory_file[FNAMESIZE]=""; +static char observatory_file[FNAMESIZE]=""; void set_observatory_file(char *fname) { strncpy(observatory_file, fname, FNAMESIZE-1); observatory_file[FNAMESIZE-1]=0; } -char ephem_file[FNAMESIZE]=""; +static char ephem_file[FNAMESIZE]=""; void set_ephem_file(char *fname) { strncpy(ephem_file, fname, FNAMESIZE-1); @@ -169,7 +169,7 @@ static void Read_Coefficients( double Time ) ** for the file itself, as environment-specified file or as default in ** this directory **/ -int Initialize_Ephemeris(void) +static int Initialize_Ephemeris(void) { int headerID; long file_bytes, record_bytes = (long) (ARRAY_SIZE*sizeof(double)); @@ -371,7 +371,7 @@ double utc_jd_to_tt(double jd_utc) /** **/ /**==========================================================================**/ -void Interpolate_Position( double Time , int Target , double Position[3] ) +static void Interpolate_Position( double Time , int Target , double Position[3] ) { double A[50] , Cp[50] , sum[3] , T_break , T_seg , T_sub , Tc; int i , j; @@ -496,7 +496,7 @@ void Interpolate_Position( double Time , int Target , double Position[3] ) return; } -void Interpolate_State(double Time, +static void Interpolate_State(double Time, int Target, double Position[3], double Velocity[3]) @@ -675,7 +675,7 @@ geocenter_ssbary(double jd, #define FLATTEN 0.003352813 /* flattening of earth, 1/298.257 */ #define EQUAT_RAD 6378137. /* equatorial radius of earth, meters */ -double +static double lst(double jd, double longit) { @@ -712,7 +712,7 @@ lst(double jd, return(sid_g); } -void +static void topo(double lmst, double rhocos, double rhosin, double *x_geo, double *y_geo, double *z_geo) /* computes the geocentric equatorial vector (AU) of a site from its diff --git a/orbfit/orbfit.h b/orbfit/orbfit.h index dcb03c9..29adde0 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -3,6 +3,7 @@ #ifndef ORBFIT_H #define ORBFIT_H #include "nrutil.h" +#include "orbfit_api.h" #include #include #include @@ -132,10 +133,6 @@ void earth3d(double t, /* time is in years here */ int obscode, double *x, double *y, double *z); -/* UTC Julian Date to TT Julian Date. DE405 is evaluated at TT. */ -double -utc_jd_to_tt(double jd_utc); - /*ICRS vector from SSBary to observatory location. * jd is a UTC Julian Date. */ void @@ -229,15 +226,6 @@ int read_radec(OBSERVATION obsarray[], char *fname, int *nobs); -/* reset astrometric error assigned to MPC obs.*/ -void -set_mpc_dtheta(double d); - -/* set filename for ephemeris or observatory data */ -void -set_ephem_file(char *fname); -void -set_observatory_file(char *fname); /* preliminary fit to observations */ void diff --git a/orbfit/orbfit1.c b/orbfit/orbfit1.c index 1485179..50446c3 100644 --- a/orbfit/orbfit1.c +++ b/orbfit/orbfit1.c @@ -21,7 +21,7 @@ double lat0, lon0; /* ecliptic lat & lon of tangent point */ double xBary, yBary, zBary; /*Posn of barycenter in our system*/ double jd0; /* Zeropoint of time scale */ -double mpc_dtheta=DEFAULT_DTHETA; /*default astrometric error*/ +static double mpc_dtheta=DEFAULT_DTHETA; /*default astrometric error*/ void set_mpc_dtheta(double d) { @@ -412,19 +412,6 @@ read_radec(OBSERVATION obsarray[], char *fname, int *nobs) return(0); } -/* Write an observation (in arcseconds) to already-opened file - * descriptor. - */ -void -write_obs(FILE *fptr, - OBSERVATION *obs) -{ - fprintf(fptr,"%10.7f %9.2f %5.2f %9.2f %5.2f %3d\n", - obs->obstime, obs->thetax/ARCSEC, obs->dthetax/ARCSEC, - obs->thetay/ARCSEC, obs->dthetay/ARCSEC, obs->obscode); - return; -} - /* Read an orbit in alpha, beta, gamma format & covar matrix from file. Also reads the coordinate frame specifications. */ @@ -856,31 +843,12 @@ covar_map(double **covar_in, return; } -/* calculate acceleration given barycentric coords x */ -/* This version just does the barycentric approximation. */ -void -accel(double *xx, - double *a, - double t) -{ - double x,y,z,acc,r2; - x = xx[0] - xBary; - y = xx[1] - yBary; - z = xx[2] - zBary; - r2 = x*x+y*y+z*z; - acc = -GM*SSMASS*pow( r2 , -1.5); - a[0] = x*acc; - a[1] = y*acc; - a[2] = z*acc; - return; -} - /* calculate acceleration given barycentric coords x */ /* Here is a more sophisticated version that includes perturbations*/ /**** ??? Note that we could cache the planet positions if this turns *** out to be slower than desired. ***/ -void +static void accel_nbody(double *xx, double *a, double t) @@ -987,7 +955,6 @@ kbo3d(PBASIS *pin, /* jump velocity to middle of time step */ accel_nbody(x,a,tx); - /*accel(x,a,tx);*/ dtv = 0.5*dt + (tx-tv); for (i=0; i<3; i++) v[i] += dtv*a[i]; tv += dtv; @@ -998,7 +965,6 @@ kbo3d(PBASIS *pin, /* Now take the last step from tx to t */ accel_nbody(x,a,tx); - /*accel(x,a,tx);*/ for (i=0; i<3; i++) { xout[i] = x[i] + v[i]*(t-tx) + a[i]*(t-tx)*(0.5*(t+tx)-tv); vout[i] = v[i] + a[i]*(t-tv); diff --git a/orbfit/orbfit2.c b/orbfit/orbfit2.c index 0413bfc..dea70ee 100644 --- a/orbfit/orbfit2.c +++ b/orbfit/orbfit2.c @@ -24,7 +24,7 @@ mrqmin_orbit(OBSERVATION obsarray[], int ndata, double a[], int ia[], int ma, double **covar, double **alpha, double *chisq, double *alamda, double energy_wt); -void +static void mrqfit(OBSERVATION *obsarray, int nobs, PBASIS *p, @@ -158,7 +158,7 @@ fit_observations(OBSERVATION obsarray[], /* Subroutine which executes the mrqmin optimization */ -void +static void mrqfit(OBSERVATION *obsarray, int nobs, PBASIS *p, diff --git a/orbfit/orbfit_api.h b/orbfit/orbfit_api.h index e75892e..99a5706 100644 --- a/orbfit/orbfit_api.h +++ b/orbfit/orbfit_api.h @@ -3,26 +3,45 @@ * overwritten by the next call; the library is not reentrant. * * On failure the results are NaN and orbfit_last_error() describes the - * problem; it is the empty string after a successful call. */ + * problem; it is the empty string after a successful call. + * + * The library is built with hidden symbol visibility; only the functions + * marked ORBFIT_API here are exported. */ #ifndef ORBFIT_API_H #define ORBFIT_API_H +#if defined(__GNUC__) +#define ORBFIT_API __attribute__((visibility("default"))) +#else +#define ORBFIT_API +#endif + /* Fit the observations in mpc_filename and write the a/b/g orbit to * abg_filename. Returns {distance, distance uncertainty} in AU. */ -double *fitradec(char *mpc_filename, char *abg_filename); +ORBFIT_API double *fitradec(char *mpc_filename, char *abg_filename); /* RA, Dec (deg), error ellipse a, b (arcsec), PA (deg), distance (AU), * ecliptic lon, lat (deg) at UTC Julian Date jdate from site obscode. */ -double *predict(char *abg_file, double jdate, int obscode); +ORBFIT_API double *predict(char *abg_file, double jdate, int obscode); /* Position {x, y, z} (AU) of the target at UTC Julian Date jdate. */ -double *predict_helio(char *abg_file, double jdate, int obscode); +ORBFIT_API double *predict_helio(char *abg_file, double jdate, int obscode); /* a, e, i, Node, peri, T, their uncertainties, epoch, distance and distance * uncertainty for the orbit in abg_file. */ -double *abg_to_aei(char *abg_file); +ORBFIT_API double *abg_to_aei(char *abg_file); /* Message from the last failed call, or "" */ -const char *orbfit_last_error(void); +ORBFIT_API const char *orbfit_last_error(void); + +/* Override the ephemeris and observatory files (otherwise taken from the + * ORBIT_EPHEMERIS and ORBIT_OBSERVATORIES environment variables) and the + * default astrometric uncertainty (arcsec) of MPC-format observations. */ +ORBFIT_API void set_ephem_file(char *fname); +ORBFIT_API void set_observatory_file(char *fname); +ORBFIT_API void set_mpc_dtheta(double d); + +/* TT Julian Date for a UTC Julian Date */ +ORBFIT_API double utc_jd_to_tt(double jd_utc); #endif diff --git a/orbfit/transforms.c b/orbfit/transforms.c index a6616ca..c23c61c 100644 --- a/orbfit/transforms.c +++ b/orbfit/transforms.c @@ -128,9 +128,9 @@ xyz_ec_to_eq(double x_ec, double y_ec, double z_ec, *** of projection as static variables to update as needed ****/ -double old_lat0=-999., old_lon0, clat0, slat0, clon0, slon0; +static double old_lat0=-999., old_lon0, clat0, slat0, clon0, slon0; -void +static void check_latlon0(double lat0, double lon0) { diff --git a/setup.py b/setup.py index c96ed34..f67e483 100644 --- a/setup.py +++ b/setup.py @@ -27,5 +27,6 @@ packages=find_packages(where='src'), package_dir = {"": "src"}, headers = headers, - ext_modules=[Extension('mp_ephem.orbfit', sources)], + ext_modules=[Extension('mp_ephem.orbfit', sources, + extra_compile_args=['-fvisibility=hidden'])], ) diff --git a/test/test_errors.py b/test/test_errors.py index e1b3f94..2b53fa9 100644 --- a/test/test_errors.py +++ b/test/test_errors.py @@ -1,6 +1,7 @@ """ Failures inside liborbfit raise BKOrbitError instead of ending the process, -and concurrent calls from threads give the same answers as serial ones. +concurrent calls from threads give the same answers as serial ones, and the +library exports only its public entry points. """ import os import tempfile @@ -10,10 +11,22 @@ from astropy.time import Time import mp_ephem +from mp_ephem import bk_orbit __PATH__ = os.path.dirname(__file__) +class LibraryExports(unittest.TestCase): + + def test_only_api_exported(self): + lib = bk_orbit._load_orbfit() + for name in ('fitradec', 'predict', 'predict_helio', 'abg_to_aei', 'orbfit_last_error', + 'set_ephem_file', 'set_observatory_file', 'set_mpc_dtheta', 'utc_jd_to_tt'): + self.assertTrue(hasattr(lib, name), name) + for name in ('gaussj', 'dvector', 'nrerror', 'kbo3d', 'read_abg', 'orbfit_fail'): + self.assertFalse(hasattr(lib, name), name) + + class LibraryErrors(unittest.TestCase): @classmethod From d4b4939f98d3fa072e3db82ab54a741651ca9dd6 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 20:41:53 -0700 Subject: [PATCH 7/7] Compare the loaded and fitted o3o08 orbits to 1e-10 instead of exactly The .abg file is now written at full precision, so a fit on Linux x86-64 differs from the committed macOS-generated file by up to ~4e-13 relative instead of being hidden by 8-decimal rounding. Co-authored-by: Cursor --- test/test_orbit.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_orbit.py b/test/test_orbit.py index 527d2c6..714656e 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -92,15 +92,18 @@ def test_data(self): def test_abg_load(self): """ - Test that loading an abg file returns the same results as calling fit_radec + Test that loading an abg file returns the same results as calling fit_radec. + The abg file is written at full precision, so the fit on another + platform can differ from it in the last few bits. :return: """ orbit1 = mp_ephem.BKOrbit(self.observations, abg_file=self.abg_filename) orbit2 = mp_ephem.BKOrbit(self.observations) for attr in ['a', 'e', 'Node', 'inc', 'om', 'T', 'distance']: - self.assertEqual(getattr(orbit1, attr), - getattr(orbit2, attr)) + value1 = getattr(orbit1, attr).value + value2 = getattr(orbit2, attr).value + self.assertAlmostEqual(value1, value2, delta=1e-10 * abs(value2), msg=attr) def test_OSSOSParser(self): mpc_line = " O13BL3UV C2013 08 02.50855 01 00 04.549+04 59 01.53 24.3 r 568 O 1645236p27 L3UV Y 106.35 4301.85 0.20 0 24.31 0.15 % hurrah!"