Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 9626) +++ /trunk/psastro/src/Makefile.am (revision 9627) @@ -13,4 +13,5 @@ libpsastro_la_SOURCES = \ psastroArguments.c \ +psastroCleanup.c \ psastroParseCamera.c \ psastroDataLoad.c \ Index: /trunk/psastro/src/psastro.c =================================================================== --- /trunk/psastro/src/psastro.c (revision 9626) +++ /trunk/psastro/src/psastro.c (revision 9627) @@ -2,5 +2,5 @@ static void usage (void) { - psErrorStackPrint(stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)"); + psErrorStackPrint(stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)\n"); exit (2); } @@ -20,5 +20,5 @@ // load identify the data sources if (!psastroParseCamera (config)) { - psErrorStackPrint(stderr, "error setting up the camera"); + psErrorStackPrint(stderr, "error setting up the camera\n"); exit (1); } @@ -27,5 +27,5 @@ // select subset of stars for astrometry if (!psastroDataLoad (config)) { - psErrorStackPrint(stderr, "error loading input data"); + psErrorStackPrint(stderr, "error loading input data\n"); exit (1); } @@ -34,5 +34,5 @@ // apply the initial guess if (!psastroAstromGuess (config)) { - psErrorStackPrint(stderr, "failed to determine initial astrometry guess"); + psErrorStackPrint(stderr, "failed to determine initial astrometry guess\n"); exit (1); } @@ -41,11 +41,11 @@ psArray *refs = psastroLoadRefstars (config); if (!refs) { - psErrorStackPrint(stderr, "failed to load reference data"); + psErrorStackPrint(stderr, "failed to load reference data\n"); exit (1); } // choose reference stars corresponding to the selected chips - if (psastroChooseRefstars (config, refs)) { - psErrorStackPrint(stderr, "failed to select reference data for chips"); + if (!psastroChooseRefstars (config, refs)) { + psErrorStackPrint(stderr, "failed to select reference data for chips\n"); exit (1); } @@ -54,10 +54,10 @@ if (mosastro == NULL) { if (!psastroChipAstrom (config, refs)) { - psErrorStackPrint(stderr, "failed to perform single chip astrometry"); + psErrorStackPrint(stderr, "failed to perform single chip astrometry\n"); exit (1); } } else { if (!psastroMosaicAstrom (config, refs)) { - psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry"); + psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry\n"); exit (1); } @@ -69,5 +69,5 @@ // write out the results if (!psastroDataSave (config)) { - psErrorStackPrint(stderr, "failed to write out data"); + psErrorStackPrint(stderr, "failed to write out data\n"); exit (1); } @@ -75,16 +75,5 @@ psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); - psFree (refs); - psFree (config); - - psTimerStop (); - psMemCheckCorruption (true); - pmModelGroupCleanup (); - psTimeFinalize (); - pmConceptsDone (); - pmConfigDone (); - // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); - - exit (0); + psastroCleanup (config, refs); + exit (EXIT_SUCCESS); } Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 9626) +++ /trunk/psastro/src/psastro.h (revision 9627) @@ -31,4 +31,5 @@ pmConfig *psastroArguments (int argc, char **argv); +void psastroCleanup (pmConfig *config, psArray *refs); bool psastroParseCamera (pmConfig *config); bool psastroDataLoad (pmConfig *config); @@ -44,7 +45,10 @@ psPlaneDistort *psPlaneDistortIdentity (); -psArray *psastroLoadReferences (pmConfig *config); +psArray *psastroLoadRefstars (pmConfig *config); bool psastroChipAstrom (pmConfig *config, psArray *refs); bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refset, psArray *rawset, psMetadata *recipe, psMetadata *updates); +bool psastroChooseRefstars (pmConfig *config, psArray *refs); +bool psastroRefstarSubset (pmReadout *readout); +pmLumFunc *psastroLuminosityFunction (psArray *stars); // utility functions: @@ -56,5 +60,4 @@ psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicAstrom (pmConfig *config, psArray *refs); -bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs); bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 9626) +++ /trunk/psastro/src/psastroArguments.c (revision 9627) @@ -27,4 +27,8 @@ // load config data from default locations pmConfig *config = pmConfigRead(&argc, argv); + if (config == NULL) { + psError(PSASTRO_ERR_CONFIG, false, "Can't read site configuration"); + return NULL; + } // save these recipe options until we have loaded the options Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 9626) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 9627) @@ -80,5 +80,5 @@ p_psDeproject (raw->sky, raw->TP, fpa->projection); - if ((i < 10) && (psTraceGetLevel("psastro.guess") > 5)) { + if ((i < 10) && (psTraceGetLevel("psastro") > 5)) { fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", raw->chip->x, raw->chip->y, @@ -100,4 +100,8 @@ tp->x, tp->y, raw->sky->r, raw->sky->d); + + psFree (fp); + psFree (tp); + psFree (ch); } Index: /trunk/psastro/src/psastroChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroChipAstrom.c (revision 9626) +++ /trunk/psastro/src/psastroChipAstrom.c (revision 9627) @@ -54,4 +54,5 @@ pmAstromWriteWCS (chip->toFPA, fpa->toSky, updates, plateScale); psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_DATA_METADATA, "psastro header stats", updates); + psFree (updates); } } Index: /trunk/psastro/src/psastroChooseRefstars.c =================================================================== --- /trunk/psastro/src/psastroChooseRefstars.c (revision 9626) +++ /trunk/psastro/src/psastroChooseRefstars.c (revision 9627) @@ -79,8 +79,7 @@ psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); + psFree (refstars); psastroRefstarSubset (readout); - - psFree (refstars); } } Index: /trunk/psastro/src/psastroErrorCodes.dat =================================================================== --- /trunk/psastro/src/psastroErrorCodes.dat (revision 9626) +++ /trunk/psastro/src/psastroErrorCodes.dat (revision 9627) @@ -7,7 +7,5 @@ ARGUMENTS Incorrect arguments CONFIG Problem in configure files -FITS Problem in FITS I/O -FITS_WCS Error interpreting FITS WCS information -PHOTOM Problem in photometry -PSF Problem in PSF -APERTURE Problem with aperture photometry +IO Problem in FITS I/O +WCS Error interpreting FITS WCS information +DATA Problem in data values Index: /trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- /trunk/psastro/src/psastroLoadRefstars.c (revision 9626) +++ /trunk/psastro/src/psastroLoadRefstars.c (revision 9627) @@ -20,4 +20,5 @@ float DECmin = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MIN"); float DECmax = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MAX"); + float MAGmax = psMetadataLookupF32(NULL, recipe, "MAG_MAX"); char *CATDIR = psMetadataLookupStr(NULL, recipe, "DVO.CATDIR"); @@ -38,17 +39,19 @@ close (fd); + psTimerStart ("psastro"); + // use fork to add timeout capability // use cfitsio |filename format to avoid tempFile # if ELIXIR_MODE if (CATDIR) { - sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT elixir -D CATMODE mef -region %f %f %f %f -o %s", CATDIR, RAmin, DECmin, RAmax, DECmax, tempFile); + sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT elixir -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", CATDIR, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); } else { - sprintf (tempLine, "getstar -D CATFORMAT elixir -D CATMODE mef -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile); + sprintf (tempLine, "getstar -D CATFORMAT elixir -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); } # else if (CATDIR) { - sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT panstarrs -D CATMODE mef -region %f %f %f %f -o %s", CATDIR, RAmin, DECmin, RAmax, DECmax, tempFile); + sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT panstarrs -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", CATDIR, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); } else { - sprintf (tempLine, "getstar -D CATFORMAT panstarrs -D CATMODE mef -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile); + sprintf (tempLine, "getstar -D CATFORMAT panstarrs -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); } # endif @@ -61,4 +64,6 @@ } + psLogMsg ("psastro", 3, "ran getstar : %f sec\n", psTimerMark ("psastro")); + // the output from getstar is a file with the Average table psFits *fits = psFitsOpen (tempFile, "r"); @@ -70,10 +75,14 @@ # endif + psTimerStart ("psastro"); psMetadata *header = psFitsReadHeader (NULL, fits); psArray *table = psFitsReadTable (fits); psFitsClose (fits); + unlink (tempFile); + psLogMsg ("psastro", 3, "read getstar output table : %f sec\n", psTimerMark ("psastro")); // convert the Average table to the pmAstromObj entries + psTimerStart ("psastro"); psArray *refs = psArrayAlloc (table->n); for (int i = 0; i < table->n; i++) { @@ -98,7 +107,9 @@ psFree (header); psFree (table); + psLogMsg ("psastro", 3, "converted table to pmAstromObj : %f sec\n", psTimerMark ("psastro")); psTrace ("psastro", 3, "loaded %ld reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n", refs->n, RAmin, DECmin, RAmax, DECmax); + return refs; } Index: /trunk/psastro/src/psastroLuminosityFunction.c =================================================================== --- /trunk/psastro/src/psastroLuminosityFunction.c (revision 9626) +++ /trunk/psastro/src/psastroLuminosityFunction.c (revision 9627) @@ -15,18 +15,18 @@ psMemSetDeallocator(func, (psFreeFunc) pmLumFuncFree); - func->mMin = 0; - func->mMax = 0; - func->slope = 0; - func->offset = 0; + func->mMin = mMin; + func->mMax = mMax; + func->slope = slope; + func->offset = offset; return func; } -pmLumFunc *psastroRefstarSubset (psArray *stars) { +pmLumFunc *psastroLuminosityFunction (psArray *stars) { // determine the max and min magnitude for the array of stars pmAstromObj *star = stars->data[0]; - mMin = star->Mag; - mMax = star->Mag; + double mMin = star->Mag; + double mMax = star->Mag; for (int i = 0; i < stars->n; i++) { star = stars->data[i]; @@ -35,5 +35,5 @@ } - int nBin = 1 + (nMax - nMin) / dMag; + int nBin = 1 + (mMax - mMin) / dMag; if (nBin <= 1) { psError(PSASTRO_ERR_DATA, true, "magnitude range of 0.0\n"); @@ -45,8 +45,20 @@ // bin[i] = mMin + i*dMag : mMin + (i+1)*dMag psVector *nMags = psVectorAlloc (nBin, PS_TYPE_F32); + nMags->n = nBin; + psVectorInit (nMags, 0); + for (int i = 0; i < stars->n; i++) { star = stars->data[i]; int bin = (star->Mag - mMin) / dMag; nMags->data.F32[bin] += 1.0; + } + + // find the peak and position + int iPeak = 0; + int nPeak = 0; + for (int i = 0; i < nMags->n; i++) { + if (nMags->data.F32[i] < nPeak) continue; + iPeak = i; + nPeak = nMags->data.F32[i]; } @@ -56,13 +68,21 @@ // create 2 vectors represnting the dlogN/dlogS line // exclude bins with no stars + // exclude points after the peak with N < 0.8*peak int n = 0; - for (int i = 0; i < stars->n; i++) { + for (int i = 0; i < nMags->n; i++) { if (nMags->data.F32[i] < 1) continue; + if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue; lnMag->data.F32[n] = log10 (nMags->data.F32[i]); Mag->data.F32[n] = mMin + (i + 0.5)*dMag; n++; } + lnMag->n = n; + Mag->n = n; + psLogMsg ("psastro", 4, "fitting %d points to luminosity function\n", n); psVector *mask = psVectorAlloc (nBin, PS_TYPE_MASK); + mask->n = Mag->n; + psVectorInit (mask, 0); + psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); @@ -70,9 +90,9 @@ stats->clipIter = 3; - line = psVectorClipFitPolynomial1D(line, stats, mask, 1, lnMag, NULL, Mag); + line = psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag); // find min and max unmasked magnitudes - mMinValid = NAN; - mMaxValid = NAN; + double mMinValid = NAN; + double mMaxValid = NAN; for (int i = 0; i < Mag->n; i++) { if (mask->data.U8[i]) continue; @@ -89,4 +109,12 @@ pmLumFunc *lumFunc = pmLumFuncAlloc (mMinValid, mMaxValid, line->coeff[0], line->coeff[1]); + + psFree (lnMag); + psFree (nMags); + psFree (Mag); + psFree (mask); + psFree (line); + psFree (stats); + return (lumFunc); } Index: /trunk/psastro/src/psastroRefstarSubset.c =================================================================== --- /trunk/psastro/src/psastroRefstarSubset.c (revision 9626) +++ /trunk/psastro/src/psastroRefstarSubset.c (revision 9627) @@ -5,9 +5,15 @@ // select the raw objects for this readout psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); - if (rawstars == NULL) { continue; } + if (rawstars == NULL) { + psError(PSASTRO_ERR_DATA, false, "missing rawstars in psastroRefstarSubset\n"); + return false; + } // select the raw objects for this readout psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); - if (refstars == NULL) { continue; } + if (refstars == NULL) { + psError(PSASTRO_ERR_DATA, false, "missing refstars in psastroRefstarSubset\n"); + return false; + } // calculate luminosity functions for rawstars and refstars @@ -15,14 +21,27 @@ // is needed... pmLumFunc *rawfunc = psastroLuminosityFunction (rawstars); + if (rawfunc == NULL) { + psError(PSASTRO_ERR_DATA, false, "error measuring rawstar luminosity function\n"); + return NULL; + } pmLumFunc *reffunc = psastroLuminosityFunction (refstars); + if (reffunc == NULL) { + psError(PSASTRO_ERR_DATA, false, "error measuring refstar luminosity function\n"); + return NULL; + } + + // XXXX test + // psFree (rawfunc); + // psFree (reffunc); + // return true; // what is the offset between the two lines at the average magnitude? - mRef = 0.5*(reffunc->mMin + reffunc->mMax); - logRho = mRef * reffunc->slope + reffunc->offset; - mRaw = (logRho - rawfunc->offset) / rawfunc->slope; + double mRef = 0.5*(reffunc->mMin + reffunc->mMax); + double logRho = mRef * reffunc->slope + reffunc->offset; + double mRaw = (logRho - rawfunc->offset) / rawfunc->slope; psLogMsg ("psastro", 4, "mRef: %f, logRho: %f, mRaw: %f\n", mRef, logRho, mRaw); - mRefMax = rawfunc->mMax - mRaw + mRef; + double mRefMax = rawfunc->mMax - mRaw + mRef; psLogMsg ("psastro", 4, "clipping stars fainter than %f\n", mRefMax); @@ -31,12 +50,16 @@ pmAstromObj *ref = refstars->data[i]; if (ref->Mag > mRefMax) continue; - psArrayAdd (subset, ref); + psArrayAdd (subset, 100, ref); } - psLogMsg ("psastro", 4, "keeping %d of %d reference stars\n", subset->n, refstars->n); + psLogMsg ("psastro", 4, "keeping %ld of %ld reference stars\n", subset->n, refstars->n); psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS"); psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", subset); + psFree (rawfunc); + psFree (reffunc); + psFree (subset); + return true; }