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/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..b97fd58 100644 --- a/orbfit/dms.c +++ b/orbfit/dms.c @@ -1,25 +1,23 @@ /* $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]; + 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++; @@ -30,15 +28,14 @@ dmsdeg(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); - return(0.); + orbfit_fail("Error interpreting dms ->%s<-",string); } args[iarg] = atof(substring); i += strlen(substring); @@ -51,63 +48,8 @@ dmsdeg(char *string) /* slight variation for hours:*/ double -hmsdeg(char *string) +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 0c760bb..07c98b8 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. @@ -34,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); @@ -95,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 */ @@ -110,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); } } @@ -172,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. **/ /** **/ @@ -183,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() +static 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 @@ -205,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; } /**==========================================================================**/ @@ -389,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; @@ -400,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. */ @@ -433,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); } /*--------------------------------------------------------------------------*/ @@ -474,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"); } /*..........................................................................*/ @@ -519,12 +496,12 @@ 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]) { - 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; @@ -535,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. */ @@ -548,7 +522,6 @@ void Interpolate_State(double Time, for ( i=0 ; i<50 ; i++ ) { A[i] = 0.0; - B[i] = 0.0; } /*--------------------------------------------------------------------------*/ @@ -569,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); } /*--------------------------------------------------------------------------*/ @@ -610,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"); } /*..........................................................................*/ @@ -683,7 +654,7 @@ geocenter_ssbary(double jd, int i; if (!init) { - if (Initialize_Ephemeris()) exit(1); + if (Initialize_Ephemeris()) orbfit_rethrow(); init = 1; } @@ -704,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) { @@ -741,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 @@ -881,9 +852,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; } @@ -930,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; @@ -941,28 +909,30 @@ 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", - 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; 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,22 +943,31 @@ 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; } - 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++; } 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 @@ -1004,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/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/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;kdthetax*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..fcff828 100644 --- a/orbfit/nrutil.c +++ b/orbfit/nrutil.c @@ -1,220 +1,111 @@ +#include #include #include #include "nrutil.h" +#include "orbfit.h" -void nrerror(char error_text[]) -/* Numerical Recipes standard error handler */ -{ - void _exit(int ); +/* 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. + * + * 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. */ - fprintf(stderr,"Numerical Recipes run-time error...\n"); - fprintf(stderr,"%s\n",error_text); - fprintf(stderr,"...now exiting to system...\n"); - _exit(1); -} +typedef union nrblock { + struct { + union nrblock *prev, *next; + } link; + max_align_t align; +} nrblock; -float *vector(int nl,int nh) -/* allocate a float vector with subscript range v[nl..nh] */ +static nrblock live = {{&live, &live}}; + +void *nr_malloc(size_t n) { - float *v; + nrblock *b = malloc(sizeof(nrblock) + n); - v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float))-nl; - if (!v) nrerror("allocation failure in vector()"); - return v; + 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; } -int *ivector(int nl,int nh) -/* allocate an int vector with subscript range v[nl..nh] */ +void nr_free(void *p) { - int *v; + nrblock *b; - v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int))-nl; - if (!v) nrerror("allocation failure in ivector()"); - return v; + 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); } -double *dvector(int nl,int nh) -/* allocate a double vector with subscript range v[nl..nh] */ +void nr_free_all(void) { - double *v; - - v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double))-nl; - if (!v) nrerror("allocation failure in dvector()"); - return v; + while (live.link.next != &live) nr_free(live.link.next + 1); } -float **matrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; -/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ +void nrerror(const char *error_text) +/* Numerical Recipes standard error handler */ { - 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; + orbfit_fail("Numerical Recipes run-time error: %s", error_text); } -double **dmatrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; -/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ +int *ivector(int nl, int nh) +/* allocate an int vector with subscript range v[nl..nh] */ { - int i; - double **m; - - /* allocate pointers to rows */ - m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*))-nrl; - 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; - if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); - } - /* return pointer to array of pointers to rows */ - return m; + if (nl < 0 || nh < nl) nrerror("bad range in ivector()"); + return (int *) nr_malloc((size_t) (nh+1)*sizeof(int)); } -int **imatrix(nrl,nrh,ncl,nch) -int nrl,nrh,ncl,nch; -/* allocate an int matrix with subscript range m[nrl..nrh][ncl..nch] */ +double *dvector(int nl, int nh) +/* allocate a double vector with subscript range v[nl..nh] */ { - 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; + if (nl < 0 || nh < nl) nrerror("bad range in dvector()"); + return (double *) nr_malloc((size_t) (nh+1)*sizeof(double)); } -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; +double **dmatrix(int nrl, int nrh, int ncl, int nch) +/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ { - 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; + int i; + double **m; - /* return pointer to array of pointers to rows */ - return m; -} + if (nrl < 0 || nrh < nrl || ncl < 0 || nch < ncl) + nrerror("bad range in dmatrix()"); -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 */ + /* allocate pointers to rows, then the rows */ + m=(double **) nr_malloc((size_t) (nrh+1)*sizeof(double*)); + for(i=0;i=nrl;i--) free((char*) (m[i]+ncl)); - free((char*) (m+nrl)); + (void) nl; (void) nh; + nr_free(v); } -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--) nr_free(m[i]); + nr_free(m); } diff --git a/orbfit/nrutil.h b/orbfit/nrutil.h index bc6ca08..23b24af 100644 --- a/orbfit/nrutil.h +++ b/orbfit/nrutil.h @@ -1,18 +1,28 @@ -float *vector(int, int); -float **matrix(int, int, int, int); +#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); 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 *) +#ifdef __GNUC__ + __attribute__((noreturn)) +#endif + ; + +void lubksb(double **, int, int *, double[]); +int 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..29adde0 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -1,6 +1,9 @@ /* $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 "orbfit_api.h" #include #include #include @@ -43,6 +46,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*/ @@ -76,10 +80,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 */ @@ -129,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 @@ -226,21 +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); - -/* Read these options */ -int -read_options(int* iarg, - int argc, - char *argv[]); /* preliminary fit to observations */ void @@ -363,3 +348,42 @@ 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); + +/* 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); + +/* 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 e052dfe..50446c3 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 @@ -23,10 +21,8 @@ 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*/ -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.) { @@ -123,7 +117,7 @@ kbo2d(PBASIS *pin, } /* incorporate derivative w.r.t simplified time delay */ if (dx!=NULL) dx[5] += vk[0]*invz*invz/SPEED_OF_LIGHT; - if (dx!=NULL) dy[5] += vk[1]*invz*invz/SPEED_OF_LIGHT; + if (dy!=NULL) dy[5] += vk[1]*invz*invz/SPEED_OF_LIGHT; return distance; } @@ -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); @@ -206,60 +199,99 @@ 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 +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; - 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); - 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) { - fprintf(stderr,"Format error in observation file:\n ->%s\n",inbuff); + 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) { + 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 { + 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)) { + orbfit_error("Format error in satellite position line:\n ->%.*s", LINE(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) { + orbfit_error("Format error in observation file:\n ->%.*s", LINE(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) { + orbfit_error("Missing observatory code:\n ->%.*s", LINE(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) { - 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); 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) { - 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); @@ -292,29 +324,29 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) * set up the tangent point and jd0 to correspond to 1st observation. */ /* Blank & commented lines are skipped */ -/* Input file is stdin if fname==NULL */ +/* Input file is stdin if fname==NULL. obsarray holds MAXOBS entries. + * Returns non-zero on error. */ int read_radec(OBSERVATION obsarray[], char *fname, int *nobs) { FILE *fptr; - OBSERVATION *obs; + OBSERVATION *obs = NULL, scanned; int scan_status_flag; - int i; char inbuff[256]; - double jd,ra,dec,elat,elon; + double elat,elon; if (fname==NULL) fptr = stdin; else if ( (fptr=fopen(fname,"r"))==NULL) { - fprintf(stderr,"Error opening observations file %s\n",fname); - exit(1); + orbfit_error("Error opening observations file %s",fname); + return(1); } *nobs=0; while ( fgets_nocomment(inbuff,255,fptr,NULL)!=NULL) { // obs refers to the previous observation, after the first loop. - scan_status_flag = scan_observation(inbuff, &(obsarray[*nobs]), obs); + scan_status_flag = scan_observation(inbuff, &scanned, obs); // scanned line was 2nd line of two line format so don't advance nobs as all we did was reset the observer x/y if (scan_status_flag == -1) { @@ -324,11 +356,18 @@ 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"); - exit(1); + if (fname!=NULL) fclose(fptr); + return(1); } + if (*nobs >= MAXOBS) { + orbfit_error("More than %d observations in %s", 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); @@ -373,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. */ @@ -395,59 +421,62 @@ read_abg(char *fname, double **covar) { FILE *fptr; - int i; + int i, status = 1; char inbuff[256]; 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"); - return(1); + 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"); - return(1); + 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"); - return(1); + 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"); - return(1); + 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"); - return(1); + 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"); - return(1); + orbfit_error("Error reading coordinate system from %s", fname); + goto done; } lat0 *= DTOR; lon0 *= DTOR; + status = 0; - if (fname!=NULL) fclose(fptr); - return(0); +done: + if (fptr!=stdin) fclose(fptr); + return(status); } /* Take a set of observations and make a preliminary fit using the @@ -513,9 +542,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 */ @@ -555,7 +594,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,26 +630,28 @@ 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; + 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 @@ -619,7 +660,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"); } @@ -633,9 +674,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 +699,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 */ @@ -764,9 +803,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++) { @@ -805,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) @@ -861,6 +880,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 @@ -875,14 +904,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; - double tstep=20.*DAY; /* 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 && @@ -928,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; @@ -939,13 +965,12 @@ 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); } /* 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; @@ -1054,8 +1079,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 */ @@ -1085,8 +1109,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 +1117,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); @@ -1103,52 +1125,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 ( ; *iargg = a[5]; p->gdot = a[6]; + free_dvector(a,1,6); + free_dmatrix(alpha,1,6,1,6); return; } diff --git a/orbfit/orbfit_api.h b/orbfit/orbfit_api.h new file mode 100644 index 0000000..99a5706 --- /dev/null +++ b/orbfit/orbfit_api.h @@ -0,0 +1,47 @@ +/* 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. + * + * On failure the results are NaN and orbfit_last_error() describes the + * 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. */ +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. */ +ORBFIT_API double *predict(char *abg_file, double jdate, int obscode); + +/* Position {x, y, z} (AU) of the target at UTC Julian Date jdate. */ +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. */ +ORBFIT_API double *abg_to_aei(char *abg_file); + +/* Message from the last failed call, or "" */ +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/orbfit_error.c b/orbfit/orbfit_error.c new file mode 100644 index 0000000..beb3ac2 --- /dev/null +++ b/orbfit/orbfit_error.c @@ -0,0 +1,67 @@ +/* Error reporting for the orbfit library. Errors are recorded so the + * Python caller can read them with orbfit_last_error(); orbfit_fail() + * also unwinds to the guard set by the current entry point instead of + * terminating the host process. */ +#include +#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 2025d0b..948e651 100644 --- a/orbfit/orbfitmodule.c +++ b/orbfit/orbfitmodule.c @@ -1,94 +1,90 @@ #include "orbfit.h" +#include "orbfit_api.h" #include #include #include -double *fitradec(char *mpc_filename, char *abg_filename) +/* 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) { + int i; + for (i = 0; i < n; i++) result[i] = NAN; +} - FILE *abg_file ; - FILE *res_file ; +static void +write_abg(const char *abg_filename, PBASIS *p, double **covar) +{ + FILE *abg_file; - OBSERVATION obsarray[MAXOBS]; - int nobs; + if ((abg_file = fopen(abg_filename,"w")) == NULL) + orbfit_fail("Error opening a/b/g output file %s", abg_filename); - ORBIT orbit; - PBASIS p; - XVBASIS xv; + 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); - double d, dd; + 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]; + jmp_buf guard; + OBSERVATION *obsarray; + int nobs; + PBASIS p; + double d, dd; double **covar; double chisq; - int i; int dof; - covar = dmatrix(1,6,1,6); + ORBFIT_GUARD(guard, result, 2); + 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"); - return result ; - } + if (read_radec(obsarray, mpc_filename, &nobs)) orbfit_rethrow(); - /* - for (i=0;i 1.) return 1.; + if (x < -1.) return -1.; + return x; +} + void orbitElements(XVBASIS *xv, ORBIT *orb) { - int i,j,k; + int k; double combinedMass; /* mass of Sun + mass of KBO */ double epochTime; @@ -325,8 +334,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*/ @@ -370,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; @@ -430,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.); } @@ -462,15 +486,14 @@ 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; 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 */ @@ -489,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); @@ -502,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/setup.py b/setup.py index 6a74e81..f67e483 100644 --- a/setup.py +++ b/setup.py @@ -28,6 +28,5 @@ package_dir = {"": "src"}, headers = headers, ext_modules=[Extension('mp_ephem.orbfit', sources, - extra_compile_args=['-Wno-unused-variable'], - )], + extra_compile_args=['-fvisibility=hidden'])], ) 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 new file mode 100644 index 0000000..2bb8c77 --- /dev/null +++ b/test/c/orbfit_driver.c @@ -0,0 +1,88 @@ +/* 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 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 +#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, const int *required, int show) +{ + int i, bad = 0; + if (show) { + printf("%-14s", what); + for (i = 0; i < n; i++) printf(" %.10g", v[i]); + 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); + 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[]) +{ + 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; + + 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, fit_req, i == 0); + r = abg_to_aei(argv[2]); + 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, 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/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_errors.py b/test/test_errors.py new file mode 100644 index 0000000..2b53fa9 --- /dev/null +++ b/test/test_errors.py @@ -0,0 +1,86 @@ +""" +Failures inside liborbfit raise BKOrbitError instead of ending the process, +concurrent calls from threads give the same answers as serial ones, and the +library exports only its public entry points. +""" +import os +import tempfile +import threading +import unittest + +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 + 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) diff --git a/test/test_orbit.py b/test/test_orbit.py index 0e580d4..714656e 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): @@ -91,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!" @@ -211,4 +215,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 new file mode 100644 index 0000000..42e4419 --- /dev/null +++ b/test/test_regression.py @@ -0,0 +1,58 @@ +""" +Pin the orbit fit, prediction and uncertainty for o3o08 so that changes to the +C library that move results are noticed. + +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 + +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.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.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.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.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__': + unittest.main()