Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 38601) +++ trunk/Ohana/src/relastro/Makefile (revision 38986) @@ -20,5 +20,7 @@ relastro_client: $(BIN)/relastro_client.$(ARCH) -install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client $(DESTBIN)/testparallax +install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client + +# $(DESTBIN)/testparallax RELASTRO = \ @@ -26,8 +28,9 @@ $(SRC)/FitChip.$(ARCH).o \ $(SRC)/FitMosaic.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ -$(SRC)/FitPar.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ -$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitPosPMfixed.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -99,8 +102,9 @@ RELASTRO_CLIENT = \ $(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ -$(SRC)/FitPar.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ -$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitPosPMfixed.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ Index: trunk/Ohana/src/relastro/doc/bright.sh =================================================================== --- trunk/Ohana/src/relastro/doc/bright.sh (revision 38986) +++ trunk/Ohana/src/relastro/doc/bright.sh (revision 38986) @@ -0,0 +1,45 @@ + +macro go.example + if ($0 != 6) + echo "USAGE: go.example (exptime) (zp) (Mo) (scale) color" + break + end + + local exptime zp Mo scale + + $exptime = $1 + $zp = $2 + $Mo = $3 + $scale = $4 + $color = $5 + + create mag 10 25 0.05 + + set flux = $exptime*ten(0.4*($zp - mag)) + set dflux = sqrt(flux) + set dmag = dflux / flux + set mag:inst = mag - $zp - 2.5*log($exptime) + + set dPosB = 0.335 / (1.0 + exp($scale*(mag - $Mo))) + set dPos = sqrt(0.015^2 + dmag^2 + dPosB^2) + + dev -n 0 + plot mag dPos -c $color + + dev -n 1 + plot mag:inst dPos -c $color + + echo Mo,inst = {$Mo - $zp - 2.5*log($exptime)} +end + +macro go.filters + + dev -n 0 ; lim 9.5 25.5 0.0 0.35; clear; box + dev -n 1 ; lim -20.0 -4.0 0.0 0.35; clear; box + + go.example 43.0 24.56 13.00 1.3 blue + go.example 40.0 24.75 12.00 1.3 green + go.example 45.0 24.61 11.75 1.3 black + go.example 30.0 24.25 11.25 1.8 violet + go.example 30.0 23.32 11.00 2.0 red +end Index: trunk/Ohana/src/relastro/doc/memory.txt =================================================================== --- trunk/Ohana/src/relastro/doc/memory.txt (revision 38986) +++ trunk/Ohana/src/relastro/doc/memory.txt (revision 38986) @@ -0,0 +1,67 @@ + +2015.08.28 + +-update-offsets -hostID 90 -D CATDIR /data/ipp094.0/eugene/pv3.cam.20150607/catdir.master -hostdir /data/ipp090.0/eugene/pv3.3pi.20150612/catdir.090 -region 5.000000 20.000000 5.000000 20.000000 -statmode INNER_MEAN -minerror 0.001000 -v -reset -update -max-density 3000.000000 -basic-image-search -exclude-bogus 2.000000 -D USE_FIXED_PIXCOORDS 1 -photcode SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9,WISE_W1,WISE_W2,WISE_W3,WISE_W4,SYNTH.g,SYNTH.r,SYNTH.i,SYNTH.z,SYNTH.y -dcr-blue-color g i -dcr-red-color z y + + +10:30 : trying to constrain memory usage in update-offsets mode. + + relastro_client.c:75 : 161928 memory blocks allocated (22732118 bytes total) + relastro_client.c:81 : 161933 memory blocks allocated (8754713623 bytes total) + load_images.c:27 : 161935 memory blocks allocated (8759201623 bytes total) + AstromOffsetMapIO.c:41 : 161939 memory blocks allocated (15643735064 bytes total) + AstromOffsetMapIO.c:42 : 161939 memory blocks allocated (15643735064 bytes total) + AstromOffsetMapIO.c:51 : 198076971 memory blocks allocated (19599617648 bytes total) + AstromOffsetMapIO.c:55 : 198076968 memory blocks allocated (12715087087 bytes total) + select_images.c:266 : 198076973 memory blocks allocated (12817856495 bytes total) + UpdateObjectOffsets.c:53 : 198081293 memory blocks allocated (13202313980 bytes total) + +** we have two big memory wastes: + 1) conversion from AstromOffsetMap_Disk_6x6 to AstromOffsetTable temporarily doubles storage used for maps + 2) select_images for a full-sky analysis returns a full copy of the image table. + + +16:10 : updated AstromOffsetTableLoad to read FITS table in chunks of 100k rows + + load_images.c:31 : 161935 memory blocks allocated (8759201623 bytes total) + load_images.c:33 : 198076968 memory blocks allocated (12715087087 bytes total) + +17:00 : the memory usage claimed by ohana_alloc is ok, but the linux kernel is claiming to hang on to a huge amount. + + I've tried to pre-allocate the be arrays which were being + re-alloced, but it does not help much. Now I'm going to test + not actually generating some of the data: + + 100 loops puts the memory at 17722496k + + without the call to AstromOffsetMapAppendToTable, we have 8602060k after 100 loops + + It looks like on each loop of AstromOffsetMapAppendToTable I'm adding ~100M to the footprint though I can only account for + +20:30 : yay, I figured it out. there is a minimal page size. since I + was allocating each row of the correction images, I was + getting hit by the 900,000 allocs per block of + AstromOffsetMapAppendToTable. By changing the image structure + to store the image as a single 1D array I drop that to 300,000 + allocs per block, with the memory only hitting 16856392k after + 220 loops (top now claims 8.5g which must be wrong somehow). + +2015.08.29 + +18:00 : I've updated relastro to avoid using the subset images. Now I'm trying to test with: + + /data/ipp064.0/eugene/relastro.20140923/tests.fake/catdir.sky.i1.v00 + + memory usage notes: + + load_images.c:15 : 161930 memory blocks allocated (66916795 bytes total) + load_images.c:34 : 161936 memory blocks allocated (67567719 bytes total) + load_images.c:50 : 161938 memory blocks allocated (109481535 bytes total) +eugene 21282 8.0 0.2 183916 124864 pts/1 t 18:08 0:23 /data/ippc19.0/home/eugene/src/psconfig/ipp-dev.lin64/bin/relastro + + with -use-all-images + load_images.c:15 : 161930 memory blocks allocated (66916795 bytes total) + load_images.c:34 : 161936 memory blocks allocated (67567719 bytes total) + load_images.c:56 : 161937 memory blocks allocated (68431383 bytes total) +eugene 22249 0.2 0.1 143760 84512 pts/1 t 18:22 0:00 /data/ippc19.0/home/eugene/src/psconfig/ipp-dev.lin64/bin/relastro + Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 38601) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 38986) @@ -28,4 +28,12 @@ typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget; + +typedef enum { + FIT_RESULT_RA, + FIT_RESULT_DEC, + FIT_RESULT_uR, + FIT_RESULT_uD, + FIT_RESULT_PLX, +} FitAstromResultMode; typedef enum { @@ -115,13 +123,55 @@ double Ro, dRo; double Do, dDo; - double uR, duR; double uD, duD; - - double p, dp; + double p, dp; double chisq; int Nfit; -} PMFit; +} FitAstromResult; + +typedef struct { + double **A; + double **B; + int Nterms; +} FitAstromData; + +typedef struct { + double X, dX; + double Y, dY; + double R, dR; + double D, dD; + double T, dT; + double pR; + double pD; + double C_blue; + double C_red; + int measure; +} FitAstromPoint; + +typedef struct { + off_t Nave; + off_t Npm; + off_t Npar; + off_t Nskip; + off_t Noffset; + + double *values; + FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats + int Nfit; + int NfitAlloc; + + FitAstromPoint *points; + FitAstromPoint *sample; + int Npoints; + int NpointsAlloc; + + FitAstromData *fitdataPos; + FitAstromData *fitdataPM; + FitAstromData *fitdataPar; + + Coords coords; + time_t T2000; +} FitStats; typedef struct { @@ -150,12 +200,4 @@ typedef struct { - off_t Nave; - off_t Npm; - off_t Npar; - off_t Nskip; - off_t Noffset; -} FitStats; - -typedef struct { double scale; double **Roff; @@ -180,4 +222,6 @@ int Nicrfobj; } ICRFobj; + +# define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM /* global variables set in parameter file */ @@ -224,4 +268,6 @@ int PARALLEL_OUTPUT; +int VERBOSE_IMAGE; + int VERBOSE; int VERBOSE2; @@ -234,4 +280,6 @@ int USE_ICRF_SHFIT; int USE_ICRF_POLE; + +int USE_ALL_IMAGES; int RESET; @@ -252,4 +300,6 @@ int CHIPORDER; int CHIPMAP; + +int N_BOOTSTRAP_SAMPLES; int MaxDensityUse; @@ -369,5 +419,5 @@ int liststats_pos PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB)); Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile)); -int load_images PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap)); +int load_images PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap, int UseAllImages)); Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap)); @@ -448,5 +498,5 @@ int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop); int UpdateSimple (Catalog *catalog, int Ncatalog); -int UpdateChips (Catalog *catalog, int Ncatalog); +int UpdateChips (Catalog *catalog, int Ncatalog, int Nloop); int UpdateMosaic (Catalog *catalog, int Ncatalog); int UpdateMeasures (Catalog *catalog, int Ncatalog); @@ -457,8 +507,8 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius); -int ParFactor (double *pR, double *pD, double RA, double Dec, double Time, double Tmean); -int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB); -int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts); -int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB); +int ParFactor (double *pR, double *pD, double RA, double Dec, double Time); +int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); +int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); +int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); Mosaic *getMosaicForImage (off_t N); @@ -661,2 +711,23 @@ int client_logger_init (char *dirname); int client_logger_message (char *format,...); + +int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode); +double VectorFractionInterpolate (double *values, float fraction, int Npts); +int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode); +int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints); +int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog); +int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange); +int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack); +int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff); + +void FitAstromResultInit (FitAstromResult *fit); +void FitAstromPointInit (FitAstromPoint *object); +void FitAstromDataFree (FitAstromData *fit); +FitAstromData *FitAstromDataInit (int Nterms); +void FitStatsFree (FitStats *fitStats); +void FitStatsSum (FitStats *src, FitStats *tgt); +void FitStatsReset (FitStats *tgt); +FitStats *FitStatsInit (int Nmax, int Nboot); +int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average); +void AstromErrorSetLoop (int Nloop, int isImageMode); Index: trunk/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 38986) +++ trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 38986) @@ -0,0 +1,148 @@ +# include "relastro.h" + +// I am modifying FitPM with an eye to (a) threaded operations and (b) bootstrap resampling tests. + +FitStats *FitStatsInit (int Nmax, int Nboot) { + + FitStats *fitStats = NULL; + ALLOCATE (fitStats, FitStats, 1); + + // counters to record successful fits or failures + fitStats->Nave = 0; + fitStats->Npm = 0; + fitStats->Npar = 0; + fitStats->Nskip = 0; + fitStats->Noffset = 0; + + fitStats->values = NULL; // pre-allocated array for median & robust sigma + fitStats->fit = NULL; // container to hold the fit results (Nboot > 1 for bootstrap resampling) + fitStats->Nfit = 0; + fitStats->NfitAlloc = Nboot; + if (Nboot > 0) { + ALLOCATE (fitStats->fit, FitAstromResult, Nboot); + ALLOCATE (fitStats->values, double, Nboot); + } + + // containers to hold the measurements for a given star + fitStats->points = NULL; + fitStats->sample = NULL; + fitStats->Npoints = 0; + fitStats->NpointsAlloc = Nmax; + if (Nmax > 0) { + ALLOCATE (fitStats->points, FitAstromPoint, Nmax); + ALLOCATE (fitStats->sample, FitAstromPoint, Nmax); + } + + // pre-allocated fit matrices for the 3 fit options + fitStats->fitdataPos = NULL; + fitStats->fitdataPM = NULL; + fitStats->fitdataPar = NULL; + + if (Nmax > 0) { + fitStats->fitdataPos = FitAstromDataInit (2); + fitStats->fitdataPM = FitAstromDataInit (4); + fitStats->fitdataPar = FitAstromDataInit (5); + } + + /* project coordinates to a plane centered on the object with units of arcsec */ + InitCoords (&fitStats->coords, "DEC--SIN"); + fitStats->coords.cdelt1 = fitStats->coords.cdelt2 = 1.0 / 3600.0; + + // use J2000 as a reference time + fitStats->T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); + return fitStats; +} + +void FitStatsReset (FitStats *tgt) { + tgt->Nave = 0; + tgt->Npm = 0; + tgt->Npar = 0; + tgt->Nskip = 0; + tgt->Noffset = 0; + return; +} + +void FitStatsSum (FitStats *src, FitStats *tgt) { + tgt->Nave += src->Nave ; + tgt->Npm += src->Npm ; + tgt->Npar += src->Npar ; + tgt->Nskip += src->Nskip ; + tgt->Noffset += src->Noffset ; + return; +} + +void FitStatsFree (FitStats *fitStats) { + if (!fitStats) return; + + FREE (fitStats->fit); + FREE (fitStats->values); + FREE (fitStats->points); + FREE (fitStats->sample); + + FitAstromDataFree (fitStats->fitdataPos); + FitAstromDataFree (fitStats->fitdataPM); + FitAstromDataFree (fitStats->fitdataPar); + + free (fitStats); +} + +FitAstromData *FitAstromDataInit (int Nterms) { + + FitAstromData *fit = NULL; + ALLOCATE (fit, FitAstromData, 1); + + /* do I need to do this as 2 2x2 matrix equations? */ + fit->A = array_init (Nterms, Nterms); + fit->B = array_init (Nterms, 1); + fit->Nterms = Nterms; + + return fit; +} + +void FitAstromDataFree (FitAstromData *fit) { + + if (!fit) return; + + array_free (fit->A, fit->Nterms); + array_free (fit->B, fit->Nterms); + free (fit); + return; +} + +void FitAstromPointInit (FitAstromPoint *object) { + object->X = 0.0; + object->Y = 0.0; + object->R = 0.0; + object->D = 0.0; + object->T = 0.0; + object->dX = 0.0; + object->dY = 0.0; + object->dR = 0.0; + object->dD = 0.0; + object->dT = 0.0; + object->pR = 0.0; + object->pD = 0.0; + object->C_blue = 0.0; + object->C_red = 0.0; + object->measure= -1; + return; +} + +void FitAstromResultInit (FitAstromResult *fit) { + + fit->Ro = 0.0; + fit->dRo = 0.0; + fit->Do = 0.0; + fit->dDo = 0.0; + fit->uR = 0.0; + fit->duR = 0.0; + fit->uD = 0.0; + fit->duD = 0.0; + fit->p = 0.0; + fit->dp = 0.0; + + fit->chisq = NAN; + fit->Nfit = 0; + + return; +} Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 38601) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 38986) @@ -81,5 +81,7 @@ } - if (VERBOSE2) fprintf (stderr, "using %d for %s\n", order_use, image[0].name); + if (VERBOSE2) { + fprintf (stderr, "using %d for %s\n", order_use, image[0].name); + } // when fitting the map, first fit a linear model (below? change Npolyterms to -1) @@ -92,5 +94,5 @@ for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); + fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, ref[i].dPos); } @@ -310,9 +312,10 @@ int i, N; - float *x, *y, *dX, *dY; + float *x, *y, *dX, *dY, *dP; ALLOCATE (x, float, Npts); ALLOCATE (y, float, Npts); ALLOCATE (dX, float, Npts); ALLOCATE (dY, float, Npts); + ALLOCATE (dP, float, Npts); N = 0; @@ -323,4 +326,5 @@ dX[N] = ref[i].X - raw[i].X; dY[N] = ref[i].Y - raw[i].Y; + dP[N] = ref[i].dPos; N++; } @@ -329,6 +333,6 @@ // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY) - AstromOffsetMapFit (map, x, y, dX, N, TRUE); - AstromOffsetMapFit (map, x, y, dY, N, FALSE); + AstromOffsetMapFit (map, x, y, dX, dP, N, TRUE); + AstromOffsetMapFit (map, x, y, dY, dP, N, FALSE); AstromOffsetMapRepair (map, TRUE); @@ -339,4 +343,5 @@ free (dX); free (dY); + free (dP); return TRUE; Index: trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM.c (revision 38601) +++ trunk/Ohana/src/relastro/src/FitPM.c (revision 38986) @@ -1,88 +1,69 @@ # include "relastro.h" -/* do we want an init function which does the alloc and a clear function to free? */ -int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB) { +// initial values of *fit are ignored +int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { int i; - double **A, **B; double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; - double chisq, Xf, Yf; - /* do I need to do this as 2 2x2 matrix equations? */ - A = array_init (4, 4); - B = array_init (4, 1); + myAssert (data->Nterms == 4, "invalid fit arrays"); Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; - for (i = 0; i < Npts; i++) { + + for (i = 0; i < Npoints; i++) { /* handle case where dX or dY = 0.0 */ - wx = 1.0 / SQ(dX[i]); - wy = 1.0 / SQ(dY[i]); + wx = 1.0 / SQ(points[i].dX); + wy = 1.0 / SQ(points[i].dY); Wx += wx; Wy += wy; - Tx += T[i]*wx; - Ty += T[i]*wy; + Tx += points[i].T*wx; + Ty += points[i].T*wy; - Tx2 += SQ(T[i])*wx; - Ty2 += SQ(T[i])*wy; + Tx2 += SQ(points[i].T)*wx; + Ty2 += SQ(points[i].T)*wy; - Xs += X[i]*wx; - Ys += Y[i]*wy; + Xs += points[i].X*wx; + Ys += points[i].Y*wy; - XT += X[i]*T[i]*wx; - YT += Y[i]*T[i]*wy; + XT += points[i].X*points[i].T*wx; + YT += points[i].Y*points[i].T*wy; } - A[0][0] = Wx; - A[0][1] = Tx; + data->A[0][0] = Wx; + data->A[0][1] = Tx; - A[1][0] = Tx; - A[1][1] = Tx2; + data->A[1][0] = Tx; + data->A[1][1] = Tx2; - A[2][2] = Wy; - A[2][3] = Ty; + data->A[2][2] = Wy; + data->A[2][3] = Ty; - A[3][2] = Ty; - A[3][3] = Ty2; + data->A[3][2] = Ty; + data->A[3][3] = Ty2; - B[0][0] = Xs; - B[1][0] = XT; - B[2][0] = Ys; - B[3][0] = YT; + data->B[0][0] = Xs; + data->B[1][0] = XT; + data->B[2][0] = Ys; + data->B[3][0] = YT; - dgaussjordan (A, B, 4, 1); + dgaussjordan (data->A, data->B, 4, 1); - fit[0].Ro = B[0][0]; - fit[0].uR = B[1][0]; - fit[0].Do = B[2][0]; - fit[0].uD = B[3][0]; - fit[0].p = 0.0; + fit->Ro = data->B[0][0]; + fit->uR = data->B[1][0]; + fit->Do = data->B[2][0]; + fit->uD = data->B[3][0]; + fit->p = 0.0; - fit[0].dRo = sqrt(A[0][0]); - fit[0].duR = sqrt(A[1][1]); - fit[0].dDo = sqrt(A[2][2]); - fit[0].duD = sqrt(A[3][3]); - fit[0].dp = 0.0; + fit->dRo = sqrt(data->A[0][0]); + fit->duR = sqrt(data->A[1][1]); + fit->dDo = sqrt(data->A[2][2]); + fit->duD = sqrt(data->A[3][3]); + fit->dp = 0.0; - array_free (A, 4); - array_free (B, 4); + fit->Nfit = Npoints; - // add up the chi square for the fit - chisq = 0.0; - for (i = 0; i < Npts; i++) { - Xf = fit[0].Ro + fit[0].uR*T[i]; - Yf = fit[0].Do + fit[0].uD*T[i]; - chisq += SQ(X[i] - Xf) / SQ(dX[i]); - chisq += SQ(Y[i] - Yf) / SQ(dY[i]); - if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq); - } - fit[0].Nfit = Npts; - - // the reduced chisq is divided by (Ndof = 2*Npts - 4) - fit[0].chisq = chisq / (2.0*Npts - 4.0); return (TRUE); } - -// XXX this function should (optionally?) iterate and clip outlier detections Index: trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 38601) +++ trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 38986) @@ -1,113 +1,94 @@ # include "relastro.h" -/* do we want an init function which does the alloc and a clear function to free? */ -int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB) { +// initial values of *fit are ignored +int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { int i; - double **A, **B; double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2; - double chisq, Xf, Yf; - A = array_init (5, 5); - B = array_init (5, 1); + myAssert (data->Nterms == 5, "invalid fit arrays"); PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0; Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; - for (i = 0; i < Npts; i++) { + + for (i = 0; i < Npoints; i++) { /* handle case where dX or dY = 0.0 */ - wx = 1.0 / SQ(dX[i]); - wy = 1.0 / SQ(dY[i]); + wx = 1.0 / SQ(points[i].dX); + wy = 1.0 / SQ(points[i].dY); Wx += wx; Wy += wy; - Tx += T[i]*wx; - Ty += T[i]*wy; + Tx += points[i].T*wx; + Ty += points[i].T*wy; - Tx2 += SQ(T[i])*wx; - Ty2 += SQ(T[i])*wy; + Tx2 += SQ(points[i].T)*wx; + Ty2 += SQ(points[i].T)*wy; - PR += pR[i]*wx; - PD += pD[i]*wy; + PR += points[i].pR*wx; + PD += points[i].pD*wy; - PRT += pR[i]*T[i]*wx; - PDT += pD[i]*T[i]*wy; + PRT += points[i].pR*points[i].T*wx; + PDT += points[i].pD*points[i].T*wy; - PRX += pR[i]*X[i]*wx; - PDY += pD[i]*Y[i]*wy; + PRX += points[i].pR*points[i].X*wx; + PDY += points[i].pD*points[i].Y*wy; - PR2 += SQ(pR[i])*wx; - PD2 += SQ(pD[i])*wy; + PR2 += SQ(points[i].pR)*wx; + PD2 += SQ(points[i].pD)*wy; - Xs += X[i]*wx; - Ys += Y[i]*wy; + Xs += points[i].X*wx; + Ys += points[i].Y*wy; - XT += X[i]*T[i]*wx; - YT += Y[i]*T[i]*wy; + XT += points[i].X*points[i].T*wx; + YT += points[i].Y*points[i].T*wy; } - A[0][0] = Wx; - A[0][1] = Tx; - A[0][4] = PR; + data->A[0][0] = Wx; + data->A[0][1] = Tx; + data->A[0][4] = PR; - A[1][0] = Tx; - A[1][1] = Tx2; - A[1][4] = PRT; + data->A[1][0] = Tx; + data->A[1][1] = Tx2; + data->A[1][4] = PRT; - A[2][2] = Wy; - A[2][3] = Ty; - A[2][4] = PD; + data->A[2][2] = Wy; + data->A[2][3] = Ty; + data->A[2][4] = PD; - A[3][2] = Ty; - A[3][3] = Ty2; - A[3][4] = PDT; + data->A[3][2] = Ty; + data->A[3][3] = Ty2; + data->A[3][4] = PDT; - A[4][0] = PR; - A[4][1] = PRT; - A[4][2] = PD; - A[4][3] = PDT; - A[4][4] = PR2 + PD2; + data->A[4][0] = PR; + data->A[4][1] = PRT; + data->A[4][2] = PD; + data->A[4][3] = PDT; + data->A[4][4] = PR2 + PD2; - B[0][0] = Xs; - B[1][0] = XT; - B[2][0] = Ys; - B[3][0] = YT; - B[4][0] = PRX + PDY; + data->B[0][0] = Xs; + data->B[1][0] = XT; + data->B[2][0] = Ys; + data->B[3][0] = YT; + data->B[4][0] = PRX + PDY; - dgaussjordan (A, B, 5, 1); + dgaussjordan (data->A, data->B, 5, 1); - fit[0].Ro = B[0][0]; - fit[0].uR = B[1][0]; - fit[0].Do = B[2][0]; - fit[0].uD = B[3][0]; - fit[0].p = B[4][0]; + fit->Ro = data->B[0][0]; + fit->uR = data->B[1][0]; + fit->Do = data->B[2][0]; + fit->uD = data->B[3][0]; + fit->p = data->B[4][0]; - fit[0].dRo = sqrt(A[0][0]); - fit[0].duR = sqrt(A[1][1]); - fit[0].dDo = sqrt(A[2][2]); - fit[0].duD = sqrt(A[3][3]); - fit[0].dp = sqrt(A[4][4]); + fit->dRo = sqrt(data->A[0][0]); + fit->duR = sqrt(data->A[1][1]); + fit->dDo = sqrt(data->A[2][2]); + fit->duD = sqrt(data->A[3][3]); + fit->dp = sqrt(data->A[4][4]); - array_free (A, 5); - array_free (B, 5); + fit->Nfit = Npoints; - /* get the chisq from the matrix values */ - - // add up the chi square for the fit - chisq = 0.0; - for (i = 0; i < Npts; i++) { - Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i]; - Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i]; - chisq += SQ(X[i] - Xf) / SQ(dX[i]); - chisq += SQ(Y[i] - Yf) / SQ(dY[i]); - if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq); - - } - fit[0].Nfit = Npts; - - // the reduced chisq is divided by (Ndof = 2*Npts - 5) - fit[0].chisq = chisq / (2.0*Npts - 5.0); return (TRUE); } Index: trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 38986) +++ trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 38986) @@ -0,0 +1,40 @@ +# include "relastro.h" + +int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { + + int i; + + double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys; + + myAssert (data->Nterms == 2, "invalid fit arrays"); + + Wx = Wy = Tx = Ty = Xs = Ys = 0.0; + + for (i = 0; i < Npoints; i++) { + /* handle case where dX or dY = 0.0 */ + wx = 1.0 / SQ(points[i].dX); + wy = 1.0 / SQ(points[i].dY); + + Wx += wx; + Wy += wy; + + Tx += points[i].T*wx; + Ty += points[i].T*wy; + + Xs += points[i].X*wx; + Ys += points[i].Y*wy; + } + + fit->Ro = (Xs - fit->uR*Tx) / Wx; + fit->Do = (Ys - fit->uD*Ty) / Wy; + + fit->p = 0.0; + + // should be ~ 1.0 / Wx + fit->dRo = sqrt(1.0 / Wx); + fit->dDo = sqrt(1.0 / Wy); + + fit->Nfit = Npoints; + + return (TRUE); +} Index: trunk/Ohana/src/relastro/src/FrameCorrection.c =================================================================== --- trunk/Ohana/src/relastro/src/FrameCorrection.c (revision 38601) +++ trunk/Ohana/src/relastro/src/FrameCorrection.c (revision 38986) @@ -479,6 +479,6 @@ } - AstromOffsetMapFit (map, Xfit, Yfit, dXfit, Nkeep, TRUE); - AstromOffsetMapFit (map, Xfit, Yfit, dYfit, Nkeep, FALSE); + AstromOffsetMapFit (map, Xfit, Yfit, dXfit, NULL, Nkeep, TRUE); + AstromOffsetMapFit (map, Xfit, Yfit, dYfit, NULL, Nkeep, FALSE); for (i = 0; i < Npts; i++) { Index: trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c =================================================================== --- trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 38601) +++ trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 38986) @@ -240,10 +240,10 @@ float *buffer = (float *)matrix->buffer; - float **value = (raDirection) ? map->dXv : map->dYv; + float *value = (raDirection) ? map->dXv : map->dYv; int ix, iy; for (ix = 0; ix < map->Nx; ix++) { for (iy = 0; iy < map->Ny; iy++) { - buffer[ix + map->Nx*iy] = value[ix][iy]; + buffer[ix + map->Nx*iy] = value[ix + map->Nx*iy]; } } @@ -265,10 +265,10 @@ float *buffer = (float *) matrix->buffer; - float **value = (raDirection) ? map->dXv : map->dYv; + float *value = (raDirection) ? map->dXv : map->dYv; int ix, iy; for (ix = 0; ix < map->Nx; ix++) { for (iy = 0; iy < map->Ny; iy++) { - value[ix][iy] = buffer[ix + map->Nx*iy]; + value[ix + map->Nx*iy] = buffer[ix + map->Nx*iy]; } } Index: trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- trunk/Ohana/src/relastro/src/GetAstromError.c (revision 38601) +++ trunk/Ohana/src/relastro/src/GetAstromError.c (revision 38986) @@ -1,4 +1,17 @@ # include "relastro.h" # define WEIGHTED_ERRORS 1 + +// XXX hard-wire the trends identified by CZW +static float BrightMo[] = {-15.6, -16.8, -17.0, -16.7, -16.0}; +static float BrightMs[] = {1.3, 1.3, 1.3, 1.8, 2.0}; + +static int Nloop = -1; +static int isImage = FALSE; + +// Nloop is used to modify the per detection errors +void AstromErrorSetLoop (int N, int isImageMode) { + Nloop = N; + isImage = isImageMode; +} float GetAstromErrorTiny (MeasureTiny *measure, int mode) { @@ -41,5 +54,93 @@ dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM)); + // for GPC1 data, we have a bright end model: + if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) { + int Np = ((int) (measure[0].photcode / 100)) % 100; + myAssert (Np >= 0, "oops"); + myAssert (Np <= 4, "oops"); + + float Minst = measure[0].M - measure[0].dt - 25.0; + float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np]))); + dPtotal = hypot(dPtotal, dPbright); + } dPtotal = MAX (dPtotal, MIN_ERROR); + + // early on, we want 2MASS and Tycho to have a very high weight. This will force images + // to match the 2MASS / Tycho / ICRS reference frame. As Nloop gets higher, the weight + // needs to drop to allow the ps1 measurements to drive the solution + int is2MASS = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2011) && (measure[0].photcode <= 2013); + int isTycho = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2020) && (measure[0].photcode <= 2021); + int has2MASS = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS); + + // serious hack: if the object has 2MASS, we set this internal bit and adjust the + // weight to ensure the image is tied down to the 2mass frame + + if (has2MASS) { + switch (Nloop) { + case 0: + dPtotal = dPtotal / 1000.0; + break; + case 1: + dPtotal = dPtotal / 300.0; + break; + case 2: + dPtotal = dPtotal / 100.0; + break; + case 3: + dPtotal = dPtotal / 30.0; + break; + case 4: + dPtotal = dPtotal / 10.0; + break; + case 5: + dPtotal = dPtotal / 10.0; + break; + default: + break; + } + } + if (is2MASS) { + switch (Nloop) { + case 0: + dPtotal = dPtotal / 1000.0; + break; + case 1: + dPtotal = dPtotal / 300.0; + break; + case 2: + dPtotal = dPtotal / 100.0; + break; + case 3: + dPtotal = dPtotal / 30.0; + break; + case 4: + dPtotal = dPtotal / 10.0; + break; + case 5: + dPtotal = dPtotal / 10.0; + break; + default: + break; + } + } + if (isTycho) { + switch (Nloop) { + case 0: + case 1: + dPtotal = dPtotal / 200.0; + break; + case 2: + case 3: + dPtotal = dPtotal / 100.0; + break; + case 4: + case 5: + dPtotal = dPtotal / 50.0; + break; + default: + break; + } + } + return (dPtotal); } @@ -84,4 +185,15 @@ dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM)); + // for GPC1 data, we have a bright end model: + if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) { + int Np = ((int) (measure[0].photcode / 100)) % 100; + myAssert (Np >= 0, "oops"); + myAssert (Np <= 4, "oops"); + + float Minst = measure[0].M - measure[0].dt - 25.0; + float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np]))); + dPtotal = hypot(dPtotal, dPbright); + } + dPtotal = MAX (dPtotal, MIN_ERROR); return (dPtotal); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 38601) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 38986) @@ -1,3 +1,4 @@ # include "relastro.h" +int isGPC1chip (int photcode); # define USE_IMAGE_ID 1 @@ -690,4 +691,12 @@ raw[i].mask |= MARK_NAN_POS_ERROR; } + + // XXX A TEST: can we use only 2MASS measurements to fit the images? + // XXX Do NOT apply this for the real calibration + if (FALSE && !(measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS)) { + fprintf (stderr, "@"); + raw[i].mask |= MARK_BIG_OFFSET; + } + raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is @@ -751,7 +760,17 @@ ref[i].R = catalog[c].average[n].R; ref[i].D = catalog[c].average[n].D; + + // if we are applying the galaxy model, move the reference position... + if (USE_GALAXY_MODEL) { + // apply proper-motion from average position to measure epoch: + float dTime = (measure[0].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years + + // XXX do this in a better way? + ref[i].R += dTime * catalog[c].average[n].uR / 3600.0 / cos(ref[i].D*RAD_DEG); + ref[i].D += dTime * catalog[c].average[n].uD / 3600.0; + } // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position - if (USE_GALAXY_MODEL) { + if (0) { myAssert (!isnan(measure[0].RoffGAL), "oops"); myAssert (!isnan(measure[0].DoffGAL), "oops"); @@ -981,7 +1000,7 @@ if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); free (R); - free(dR); - free(D); - free(dD); + free (dR); + free (D); + free (dD); } @@ -1110,9 +1129,9 @@ if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); free (R); - free(dR); - free(D); - free(dD); - free(d2); - free(index); + free (dR); + free (D); + free (dD); + free (d2); + free (index); } @@ -1215,5 +1234,6 @@ } - if (MinBadQF > 0.0) { + if ((MinBadQF > 0.0) && isGPC1chip(measure[0].photcode)) { + if (!isfinite(measure[0].psfQF)) return FALSE; if (measure[0].psfQF < MinBadQF) return FALSE; } Index: trunk/Ohana/src/relastro/src/ParFactor.c =================================================================== --- trunk/Ohana/src/relastro/src/ParFactor.c (revision 38601) +++ trunk/Ohana/src/relastro/src/ParFactor.c (revision 38986) @@ -2,49 +2,9 @@ # define J2000 2451545. /* Julian date at standard epoch */ -# if (0) -/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ -/* ra and dec are returned as decimal hours and decimal degrees. */ -void lpsun (double jd, double *ra, double *dec) { +/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ +// jdoff is days since J2000 +int sun_ecliptic (double jdoff, double *lambda, double *beta, double *epsilon, double *Radius) { - double n, L, g, lambda,epsilon,alpha,delta,x,y,z; - - n = jd - J2000; - L = 280.460 + 0.9856474 * n; - g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN; - lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN; - epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN; - - // this is the conversion from ecliptic to celestial coords - x = cos(lambda); - y = cos(epsilon)*sin(lambda); - z = sin(epsilon)*sin(lambda); - - *ra = (atan_circ(x,y))*HRS_IN_RADIAN; - *dec = (asin(z))*DEG_IN_RADIAN; -} -# endif - -# if (0) -/* code borrowed from Skycalc : fix this stuff XXX */ -/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ -int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) { - - double n, L, g; - - - n = jd - J2000; - L = 280.460 + 0.9856474 * n; - g = (357.528 + 0.9856003 * n)*RAD_DEG; - *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees - *beta = 0.0; // approx latitude - *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees - return TRUE; -} -# endif - -/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ -int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) { - - double n = jd - J2000; // day number + double n = jdoff; // day number double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly @@ -58,17 +18,14 @@ /* given RA, DEC, Time, calculate the parallax factor */ -// Time is relative to Tmean, Tmean is years relative to J2000 -int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) { +// Time is years since J2000 +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { - double jd, lambda, beta, epsilon, Radius; + double lambda, beta, epsilon, Radius; - /* given a Time relative to Tmean, Tmean in years since J2000, determine the solar - longitude S */ + /* given a Time in years since J2000, determine the solar longitude S */ - // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); - jd = 365.25*(Time + Tmean) + J2000; - // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd); + double jdoff = 365.25*Time; - sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius); + sun_ecliptic (jdoff, &lambda, &beta, &epsilon, &Radius); double lambda_rad = lambda*RAD_DEG; @@ -99,2 +56,42 @@ return TRUE; } + +# if (0) +/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ +/* ra and dec are returned as decimal hours and decimal degrees. */ +void lpsun (double jd, double *ra, double *dec) { + + double n, L, g, lambda,epsilon,alpha,delta,x,y,z; + + n = jd - J2000; + L = 280.460 + 0.9856474 * n; + g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN; + lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN; + epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN; + + // this is the conversion from ecliptic to celestial coords + x = cos(lambda); + y = cos(epsilon)*sin(lambda); + z = sin(epsilon)*sin(lambda); + + *ra = (atan_circ(x,y))*HRS_IN_RADIAN; + *dec = (asin(z))*DEG_IN_RADIAN; +} + +/* code borrowed from Skycalc : fix this stuff XXX */ +/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ +int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) { + + double n, L, g; + + + n = jd - J2000; + L = 280.460 + 0.9856474 * n; + g = (357.528 + 0.9856003 * n)*RAD_DEG; + *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees + *beta = 0.0; // approx latitude + *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees + return TRUE; +} +# endif + Index: trunk/Ohana/src/relastro/src/Shutdown.c =================================================================== --- trunk/Ohana/src/relastro/src/Shutdown.c (revision 38601) +++ trunk/Ohana/src/relastro/src/Shutdown.c (revision 38986) @@ -21,6 +21,4 @@ va_end (argp); - SetProtect (TRUE); - gfits_db_close (db); fprintf (stderr, "ERROR: relastro halted\n"); exit (1); Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 38601) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 38986) @@ -163,12 +163,12 @@ } - if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax); if (dLmax > DPOS_MAX) { - if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax); return (FALSE); } if (dMmax > DPOS_MAX) { - if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax); return (FALSE); } Index: trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateChips.c (revision 38601) +++ trunk/Ohana/src/relastro/src/UpdateChips.c (revision 38986) @@ -29,5 +29,5 @@ // update astrometry of all chips relative to the average positions // if NTHREADS is non-zero, call the threaded version of this function -int UpdateChips (Catalog *catalog, int Ncatalog) { +int UpdateChips (Catalog *catalog, int Ncatalog, int Nloop) { off_t Nskip, Nmosaic, NnewFit, NoldFit; @@ -40,4 +40,6 @@ char *mode; + AstromErrorSetLoop (Nloop, TRUE); + if (NTHREADS) { UpdateChips_threaded (catalog, Ncatalog); @@ -57,4 +59,6 @@ // each chip is fitted independently, so we could do N at once in parallel for (i = 0; i < Nimage; i++) { + + VERBOSE_IMAGE = !strcmp(image[i].name, "o5745g0516o.356887.cm.982631.smf[XY54]"); // XXX looks like everything below is thread safe : we can unroll this into a set of @@ -91,4 +95,10 @@ } + if (VERBOSE_IMAGE) { + dump_stardata_pts (raw, Nraw, "testimage.raw.dat"); + dump_stardata_pts (ref, Nref, "testimage.ref.dat"); + fprintf (stderr, "dumped test image\n"); + } + // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions. assert (Nraw == Nref); @@ -105,5 +115,5 @@ // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n", i, Nraw); if (!FitChip (raw, ref, Nraw, &image[i])) { - if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); + if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused %d of %d\n", image[i].name, i, Nraw, image[i].nFitAstrom, image[i].nstar); if (1) { @@ -124,5 +134,5 @@ if (!checkStarMap (i)) { - if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused: %d\n", image[i].name, i, Nraw, image[i].nFitAstrom); if (1) { @@ -148,4 +158,5 @@ setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC); if (USE_GALAXY_MODEL) { + // XXX DEPRECATE? image[i].flags |= ID_IMAGE_ASTROM_GMM; } @@ -311,5 +322,5 @@ // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n", i, Nraw); if (!FitChip (raw, ref, Nraw, &image[i])) { - if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); + if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused %d of %d\n", image[i].name, i, Nraw, image[i].nFitAstrom, image[i].nstar); // restore status quo ante (replace truMap with tmpMap) @@ -328,5 +339,5 @@ if (!checkStarMap (i)) { - if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused: %d\n", image[i].name, i, Nraw, image[i].nFitAstrom); // restore status quo ante (replace truMap with tmpMap) Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 38601) +++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 38986) @@ -97,5 +97,5 @@ } -# if (0) +# if (1) UpdateObjectOffsets_parallel_table (table, sky); @@ -155,6 +155,11 @@ char *command = NULL; - strextend (&command, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", - group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); + strextend (&command, "relastro_client -update-offsets"); + strextend (&command, "-hostID %d", group->hosts[i][0].hostID); + strextend (&command, "-D CATDIR %s", CATDIR); + strextend (&command, "-hostdir %s", group->hosts[i][0].pathname); + strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + strextend (&command, "-statmode %s", STATMODE); + strextend (&command, "-minerror %f", MIN_ERROR); if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm"); @@ -176,4 +181,5 @@ if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); + if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); @@ -184,4 +190,6 @@ if (PhotFlagPoor) strextend (&command, "+photflagpoor %d", PhotFlagPoor); // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { @@ -303,4 +311,6 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 38601) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 38986) @@ -2,104 +2,9 @@ # define PAR_TOOFEW 5 -int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset); -int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); - -static off_t Nmax; -static double *X, *dX; -static double *Y, *dY; -static double *R, *dR; -static double *D, *dD; -static double *pX; -static double *pY; -static double *T; -static double *dT; -static double *C_blue; -static double *C_red; - -static Coords coords; - -static time_t T2000; - -void initFitStats (FitStats *fitStats) { - fitStats->Nave = 0; - fitStats->Npm = 0; - fitStats->Npar = 0; - fitStats->Nskip = 0; - fitStats->Noffset = 0; - return; -} - -void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) { - tgtFitStats->Nave += srcFitStats->Nave ; - tgtFitStats->Npm += srcFitStats->Npm ; - tgtFitStats->Npar += srcFitStats->Npar ; - tgtFitStats->Nskip += srcFitStats->Nskip ; - tgtFitStats->Noffset += srcFitStats->Noffset ; - return; -} - -void initObjectData (Catalog *catalog, int Ncatalog) { - - off_t i, j; - - Nmax = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); - } - } - - ALLOCATE (R, double, MAX (1, Nmax)); - ALLOCATE (D, double, MAX (1, Nmax)); - ALLOCATE (T, double, MAX (1, Nmax)); - ALLOCATE (X, double, MAX (1, Nmax)); - ALLOCATE (Y, double, MAX (1, Nmax)); - - ALLOCATE (dR, double, MAX (1, Nmax)); - ALLOCATE (dD, double, MAX (1, Nmax)); - ALLOCATE (dT, double, MAX (1, Nmax)); - ALLOCATE (dX, double, MAX (1, Nmax)); - ALLOCATE (dY, double, MAX (1, Nmax)); - - ALLOCATE (pX, double, MAX (1, Nmax)); - ALLOCATE (pY, double, MAX (1, Nmax)); - - ALLOCATE (C_blue, double, MAX (1, Nmax)); - ALLOCATE (C_red, double, MAX (1, Nmax)); - - /* project coordinates to a plane centered on the object with units of arcsec */ - InitCoords (&coords, "DEC--SIN"); - coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; - - // use J2000 as a reference time - T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); -} - -void freeObjectData () { - - free (R); - free (D); - free (T); - free (X); - free (Y); - - free (dR); - free (dD); - free (dT); - free (dX); - free (dY); - - free (pX); - free (pY); - - free (C_blue); - free (C_red); -} +int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog); // This function operates on both Measure and MeasureTiny. In the big stages, this should // be called with just MeasureTiny set and Measure == NULL int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) { - - initObjectData (catalog, Ncatalog); // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? @@ -109,6 +14,14 @@ } - FitStats sumStatsChips; initFitStats (&sumStatsChips); - FitStats sumStatsStack; initFitStats (&sumStatsStack); + int NmeasureMax = CatalogMaxNmeasure (catalog, Ncatalog); + + // allocate summary stats with Nmax = 0, Nboot = 0 + FitStats *sumStatsChips = FitStatsInit (0, 0); + FitStats *sumStatsStack = FitStatsInit (0, 0); + + FitStats *fitStatsChips = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES); + FitStats *fitStatsStack = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES); + + AstromErrorSetLoop (Nloop, FALSE); int i; @@ -117,6 +30,6 @@ if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); - FitStats fitStatsChips; initFitStats (&fitStatsChips); - FitStats fitStatsStack; initFitStats (&fitStatsStack); + FitStatsReset (fitStatsChips); + FitStatsReset (fitStatsStack); off_t j; @@ -129,16 +42,24 @@ SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt]; - UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack); - UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop); - } - if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsChips.Nave, fitStatsChips.Npm, fitStatsChips.Npar, fitStatsChips.Nskip, fitStatsChips.Noffset); - if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsStack.Nave, fitStatsStack.Npm, fitStatsStack.Npar, fitStatsStack.Nskip, fitStatsStack.Noffset); - sumFitStats (&fitStatsChips, &sumStatsChips); - sumFitStats (&fitStatsStack, &sumStatsStack); - } - freeObjectData (); - - if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar), sumStatsChips.Nave, sumStatsChips.Npm, sumStatsChips.Npar, sumStatsChips.Nskip, sumStatsChips.Noffset); - if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar), sumStatsStack.Nave, sumStatsStack.Npm, sumStatsStack.Npar, sumStatsStack.Nskip, sumStatsStack.Noffset); + UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack); + UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m); + } + if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsChips->Nave, fitStatsChips->Npm, fitStatsChips->Npar, fitStatsChips->Nskip, fitStatsChips->Noffset); + if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsStack->Nave, fitStatsStack->Npm, fitStatsStack->Npar, fitStatsStack->Nskip, fitStatsStack->Noffset); + FitStatsSum (fitStatsChips, sumStatsChips); + FitStatsSum (fitStatsStack, sumStatsStack); + } + + // DumpObjectsWith2MASS (catalog, Ncatalog); + + FitStatsFree (fitStatsChips); + FitStatsFree (fitStatsStack); + + if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsChips->Nave + sumStatsChips->Npm + sumStatsChips->Npar), sumStatsChips->Nave, sumStatsChips->Npm, sumStatsChips->Npar, sumStatsChips->Nskip, sumStatsChips->Noffset); + if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsStack->Nave + sumStatsStack->Npm + sumStatsStack->Npar), sumStatsStack->Nave, sumStatsStack->Npm, sumStatsStack->Npar, sumStatsStack->Nskip, sumStatsStack->Noffset); + + FitStatsFree (sumStatsChips); + FitStatsFree (sumStatsStack); + return (TRUE); } @@ -146,14 +67,37 @@ // This function operates on both Measure and MeasureTiny. In the big stages, this should // be called with just MeasureTiny set and Measure == NULL -int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) { - - int setRefColor = areImagesMatched(); +int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog) { + + int i; + for (i = 0; i < Ncatalog; i++) { + off_t j; + for (j = 0; j < catalog[i].Naverage; j++) { + /* calculate the average value of R,D for a single star */ + off_t m = catalog[i].average[j].measureOffset; + + off_t k; + for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { + MeasureTiny *measure = &catalog[i].measureT[m+k]; + if (measure->dbFlags & ID_MEAS_OBJECT_HAS_2MASS) { + fprintf (stderr, "0x%08x 0x%08x : %12.8f %12.8f %5d\n", catalog[i].average[j].objID, catalog[i].average[j].catID, measure->R, measure->D, measure->photcode); + } + } + } + } + return (TRUE); +} + +// This function operates on both Measure and MeasureTiny. In the big stages, this should +// be called with just MeasureTiny set and Measure == NULL +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff) { + + int k; /* calculate the average value of R,D for a single star */ - PMFit fit; memset (&fit, 0, sizeof(fit)); - PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN; - PMFit fitPM; memset (&fitPM, 0, sizeof(fitPM)); fitPM.chisq = NAN; - PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN; + FitAstromResult fitPos, fitPM, fitPar; + FitAstromResultInit (&fitPos); + FitAstromResultInit (&fitPM); + FitAstromResultInit (&fitPar); // if we fail to fit the astrometry for some reason, we need to set/reset these @@ -167,8 +111,4 @@ if (average[0].Nmeasure == 0) return TRUE; - int NcBlue = 0; - int NcRed = 0; - int N = 0; - int mode = FIT_MODE; // start with the globally-defined fit mode @@ -177,123 +117,9 @@ XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST); - // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) - off_t k; - for (k = 0; k < average[0].Nmeasure; k++) { - - if (XVERB) { - char *date = ohana_sec_to_date (measure[k].t); - int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; - fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); - free (date); - } - - // SKIP gpc1 stack data - if (isGPC1stack(measure[k].photcode)) continue; - - // SKIP gpc1 forced-warp data - if (isGPC1warp(measure[k].photcode)) continue; - - // reset the bit to note that a detection was used (or not) - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } - - // does the measurement pass the supplied filtering constraints? - // MeasFilterTestTiny does not test psfQF - // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst - int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); - if (!keepMeasure) { - continue; - } - - double Ri = getMeanR (&measure[k], average, secfilt); - double Di = getMeanD (&measure[k], average, secfilt); - - // if we are correcting for the Galaxy Motion Model, only should apply it here - // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and - // (b) if we think the measure R,D is already at the image epoch position - if (USE_GALAXY_MODEL && applyGalaxyOffset) { - Ri -= measure[k].RoffGAL / 3600.0; - Di -= measure[k].DoffGAL / 3600.0; - } - - // XXX add in dR,dD GAL here - - // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius - if (ExcludeBogus) { - coords.crval1 = average[0].R; - coords.crval2 = average[0].D; - double Xi, Yi; - RD_to_XY (&Xi, &Yi, Ri, Di, &coords); - double radius = hypot(Xi, Yi); - if (radius > ExcludeBogusRadius) { - measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; - if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; } - continue; - } - measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; } - } - - // outlier rejection - if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { - continue; - } - - R[N] = Ri; - D[N] = Di; - - // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000. - // T[] is time in years since J2000 (jd = 2451545) - T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years - - // dX, dY : error in arcsec -- - dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); - dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); - - // allow a given photcode or measurement to be - // ignored if the error is NAN (for photcode, set astromErrSys to NaN) - if (isnan(dX[N])) continue; - if (isnan(dY[N])) continue; - - // add systematic error in quadrature, if desired - // only do this after the fit has converged (or you will never improve the poor images) - // if (INCLUDE_SYS_ERR) { - // float dRsys = FromShortPixels(measure[k].dRsys); - // dX[N] = hypot(dX[N], dRsys); - // dY[N] = hypot(dY[N], dRsys); - // } - - // dX[N] = 0.1; - // dY[N] = 0.1; - - dT[N] = measure[k].dt; - - // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in - // arcseconds. dR, dD are the errors in those directions in degrees. IF we have - // non-circular errors (different values for X and Y), then dR and dD will be - // incorrect: they would need to be rotated to take out the position angle - dR[N] = dX[N] / 3600.0; - dD[N] = dY[N] / 3600.0; - - if (setRefColor) { - float colorBlue = getColorBlue (m+k, i); - if (!isnan(colorBlue)) { - C_blue[NcBlue] = colorBlue; - NcBlue++; - } - float colorRed = getColorRed (m+k, i); - if (!isnan(colorRed)) { - C_red[NcRed] = colorRed; - NcRed++; - } - } - - measure[k].dbFlags |= ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } - - N++; - } // loop over measurements : average[0].Nmeasure - - if (N < 1) { + // select the measurements to be used in this analysis + UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, FALSE); + + // if there are no exposure detections, use the stack position + if (fitStats->Npoints < 1) { if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) { average[0].R = average[0].Rstk; @@ -301,44 +127,16 @@ average[0].dR = average[0].dRstk; average[0].dD = average[0].dDstk; + average[0].flags |= ID_STACK_ASTROM; } return FALSE; } - // if we have too few good detections for the desired fit, or too limited a - // baseline, use a fit with fewer parameters. XXX if we have too few measurements - // for even the average position, consider including the lower-quality detections? - - // find Tmin & Tmax from the list of accepted measurements - double Tmean = 0.0; - double Tmin = T[0]; - double Tmax = T[0]; - for (k = 0; k < N; k++) { - Tmin = MIN(Tmin, T[k]); - Tmax = MAX(Tmax, T[k]); - Tmean += T[k]; - } - double Trange = Tmax - Tmin; - - if (RELASTRO_OP == OP_HIGH_SPEED) { - Tmean = 0.5*(Tmax - Tmin); - } else { - Tmean /= (float) N; - } - - /* we need to do the fit in a locally linear space; choose a ref coordinate */ - coords.crval1 = R[0]; - coords.crval2 = D[0]; + double Tmean, Trange, parRange; + FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange); // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq - // project all of the R,D coordinates to a plane centered on this coordinate. set - // the times to be relative to Tmean (this is required for parallax as well) - for (k = 0; k < N; k++) { - RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); - T[k] -= Tmean; - if (XVERB) { - fprintf (stderr, OFF_T_FMT" %f %f %f %f %f +/- %f %f\n", k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); - } - } + // if we have too few good detections for the desired fit, or too limited a baseline, + // use a fit with fewer parameters. // *** first fit for the proper motion (skip fit if Trange or Npts is too small) *** @@ -346,20 +144,31 @@ if (Trange < PM_DT_MIN) { mode = FIT_AVERAGE; - goto skipPM; - } - if (N <= PM_TOOFEW) { + goto justPosition; + } + if (fitStats->Npoints <= PM_TOOFEW) { mode = FIT_AVERAGE; - goto skipPM; - } - - FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB); - - if (XVERB) fprintf (stderr, "fitted PM: %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq); + goto justPosition; + } + + if (fitStats->NfitAlloc == 1) { + // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: + FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints); + } else { + fitStats->Nfit = 0; + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue; + fitStats->Nfit ++; + } + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + } + FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY); // project Ro, Do back to RA,DEC - XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords); - if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p); + XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords); if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - // XXX : does this make sense at 0,360 boundary? fitPM.p = fitPM.dp = 0.0; @@ -374,84 +183,119 @@ } -skipPM: // fit the parallax + proper-motion model // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only if (mode == FIT_PM_AND_PAR) { - if (Trange < PM_DT_MIN) { - mode = FIT_PM_ONLY; - goto skipPAR; - } - if (N <= PAR_TOOFEW) { - mode = FIT_PM_ONLY; - goto skipPAR; - } - float pXmin = +2.0; - float pXmax = -2.0; - float pYmin = +2.0; - float pYmax = -2.0; - for (k = 0; k < N; k++) { - ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean); - pXmin = MIN (pXmin, pX[k]); - pXmax = MAX (pXmax, pX[k]); - pYmin = MIN (pYmin, pY[k]); - pYmax = MAX (pYmax, pY[k]); - } - float dXRange = pXmax - pXmin; - float dYRange = pYmax - pYmin; - float parRange = hypot (dXRange, dYRange); - if (parRange < PAR_FACTOR_MIN) { mode = FIT_PM_ONLY; - goto skipPAR; - } - - FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB); - if (XVERB) fprintf (stderr, "fitted PM+PAR: %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq); - - XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); + goto justPosition; + } + if (fitStats->Npoints <= PAR_TOOFEW) { + mode = FIT_PM_ONLY; + goto justPosition; + } + + if (fitStats->NfitAlloc == 1) { + // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: + FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints); + } else { + fitStats->Nfit = 0; + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints); + fitStats->Nfit ++; + } + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + } + FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR); + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); + if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + average[0].flags |= ID_STAR_FIT_PAR; fitStats->Npar ++; - if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - // XXX a hard-wired hack... - if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { + if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) { mode = FIT_PM_ONLY; } } -skipPAR: +justPosition: { - // ALWAYS fit the average model - StatType statsR, statsD; - liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) - liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) + // use bootstrap resampling to check the error distribution + // if we only have one point, this is silly... + + if (fitStats->NfitAlloc == 1) { + FitAstromResultSetPM (&fitPos, 1, average); + FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints); + } else { + fitStats->Nfit = 0; + FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average); + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints); + fitStats->Nfit ++; + } + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + } + FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE); // project Ro, Do back to RA,DEC - XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords); - if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do); - - fitAve.dRo = statsR.sigma; - fitAve.dDo = statsD.sigma; - - fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN; - fitAve.Nfit = N; - - fitAve.uR = fitAve.duR = 0.0; - fitAve.uD = fitAve.duD = 0.0; - fitAve.p = fitAve.dp = 0.0; + XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); average[0].flags |= ID_STAR_FIT_AVE; fitStats->Nave ++; } + // update the bit flags of which points were used + for (k = 0; k < fitStats->Npoints; k++) { + int Nm = fitStats->points[k].measure; + myAssert (Nm >= 0, "oops"); + measure[Nm].dbFlags |= ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; } + } + + // we can set the star reference-image color only if we have loaded the image data + int setRefColor = areImagesMatched(); if (setRefColor) { + float *C_blue = NULL; + float *C_red = NULL; + ALLOCATE (C_blue, float, fitStats->Npoints); + ALLOCATE (C_red, float, fitStats->Npoints); + + int NcBlue = 0; + int NcRed = 0; + + for (k = 0; k < fitStats->Npoints; k++) { + int Nm = fitStats->points[k].measure; + float colorBlue = getColorBlue (measOff + Nm, cat); + if (!isnan(colorBlue)) { + C_blue[NcBlue] = colorBlue; + NcBlue++; + } + float colorRed = getColorRed (measOff + Nm, cat); + if (!isnan(colorRed)) { + C_red[NcRed] = colorRed; + NcRed++; + } + } + + // need to reassign here if isfinite() float colorMedian; - dsort (C_blue, NcBlue); + fsort (C_blue, NcBlue); colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN; average[0].refColorBlue = colorMedian; - dsort (C_red, NcRed); + fsort (C_red, NcRed); colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN; average[0].refColorRed = colorMedian; + + free (C_blue); + free (C_red); } @@ -459,9 +303,11 @@ // XXXX for now, just use the mode as the result: int result = mode; + FitAstromResult fit; + FitAstromResultInit (&fit); switch (result) { case FIT_AVERAGE: average[0].flags |= ID_STAR_USE_AVE; - fit = fitAve; + fit = fitPos; break; case FIT_PM_ONLY: @@ -471,5 +317,5 @@ case FIT_PM_AND_PAR: average[0].flags |= ID_STAR_USE_PAR; - fit = fitPAR; + fit = fitPar; break; } @@ -503,19 +349,19 @@ // what is the offset relative to the mean fit position? - coords.crval1 = average[0].R; - coords.crval2 = average[0].D; - if (isnan(coords.crval1)) { + fitStats->coords.crval1 = average[0].R; + fitStats->coords.crval2 = average[0].D; + if (isnan(fitStats->coords.crval1)) { return (FALSE); } - if (isnan(coords.crval2)) { + if (isnan(fitStats->coords.crval2)) { return (FALSE); } double dXoff, dYoff; - RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords); + RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &fitStats->coords); float dPos = hypot (dXoff, dYoff); if (dPos > MaxMeanOffset) { if (fitStats->Noffset < 100) { - fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff); + fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fit.Ro, fit.Do, dXoff, dYoff); } fitStats->Noffset ++; @@ -531,5 +377,5 @@ average[0].uR, average[0].uD, - fitAve.chisq, fitPM.chisq, fitPAR.chisq); + fitPos.chisq, fitPM.chisq, fitPar.chisq); average[0].R = fit.Ro; // RA in degrees @@ -546,9 +392,9 @@ average[0].dP = fit.dp; // parallax error in arcsec - average[0].ChiSqAve = fitAve.chisq; + average[0].ChiSqAve = fitPos.chisq; average[0].ChiSqPM = fitPM.chisq; - average[0].ChiSqPar = fitPAR.chisq; - - average[0].Tmean = (Tmean * 86400 * 365.25) + T2000; + average[0].ChiSqPar = fitPar.chisq; + + average[0].Tmean = (Tmean * 86400 * 365.25) + fitStats->T2000; average[0].Trange = (Trange * 86400 * 365.25); average[0].Npos = fit.Nfit; @@ -563,6 +409,4 @@ // be called with just MeasureTiny set and Measure == NULL int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) { - - off_t k; // set the default values @@ -573,11 +417,8 @@ /* calculate the average value of R,D for a single star */ - PMFit fitAve; - memset (&fitAve, 0, sizeof(fitAve)); - fitAve.chisq = NAN; + FitAstromResult fitPos; + FitAstromResultInit (&fitPos); if (average[0].Nmeasure == 0) return TRUE; - - int N = 0; int XVERB = FALSE; @@ -585,17 +426,107 @@ XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST); + // select the measurements to be used in this analysis + UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE); + + // too few measurements for average position (require 2 values) + if (fitStats->Npoints < 1) return FALSE; // XXX ?? + + double Tmean, Trange, parRange; + FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange); + + FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints); + FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE); + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); + + // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE; + fitStats->Nave ++; + + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", + average[0].R, + average[0].D, + fitPos.Ro, fitPos.Do, + 3600*(average[0].R - fitPos.Ro), + 3600*(average[0].D - fitPos.Do)); + + // make sure that the fit succeeded + int status = TRUE; + status &= finite(fitPos.Ro); + status &= finite(fitPos.Do); + status &= finite(fitPos.dRo); + status &= finite(fitPos.dDo); + if (!status) { + fitStats->Nskip ++; + return FALSE; + } + + // what is the offset relative to the mean fit position? + fitStats->coords.crval1 = average[0].R; + fitStats->coords.crval2 = average[0].D; + + double dXoff, dYoff; + RD_to_XY (&dXoff, &dYoff, fitPos.Ro, fitPos.Do, &fitStats->coords); + float dPos = hypot (dXoff, dYoff); + if (dPos > MaxMeanOffset) { + if (fitStats->Noffset < 100) { + fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitPos.Ro, fitPos.Do, dXoff, dYoff); + } + fitStats->Noffset ++; + return FALSE; + } + + // set the stack position values + average[0].Rstk = fitPos.Ro; // RA in degrees + average[0].Dstk = fitPos.Do; // DEC in degrees + average[0].dRstk = fitPos.dRo; // RA scatter in arcsec + average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec + + return (TRUE); +} + +int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack) { + + // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries + + int has2MASS = FALSE; + + int Npoints = fit->Npoints = 0; + FitAstromPoint *points = fit->points; + + int TESTPT2 = FALSE; + TESTPT2 |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC); + TESTPT2 |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST); + if (TESTPT2) { + fprintf (stderr, "got test det\n"); + } + // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) + off_t k; for (k = 0; k < average[0].Nmeasure; k++) { - if (XVERB) { + if (0) { char *date = ohana_sec_to_date (measure[k].t); int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; - fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); + fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); free (date); } - // SKIP everything except gpc1 stack data - if (!isGPC1stack(measure[k].photcode)) continue; - + // SKIP gpc1 forced-warp data + if (isGPC1warp(measure[k].photcode)) continue; + + // SKIP gpc1 stack data + if (isStack) { + if (!isGPC1stack(measure[k].photcode)) continue; + } else { + if ( isGPC1stack(measure[k].photcode)) continue; + } + + // reset the bit to note that a detection was used (or not) + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } + + // does the measurement pass the supplied filtering constraints? + // MeasFilterTestTiny does not test psfQF // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); @@ -604,121 +535,304 @@ } - R[N] = getMeanR (&measure[k], average, secfilt); - D[N] = getMeanD (&measure[k], average, secfilt); - - // dX, dY : error in arcsec -- - dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); - dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); + double Ri = measure[k].R; + double Di = measure[k].D; + + // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius + if (ExcludeBogus) { + fit->coords.crval1 = average[0].R; + fit->coords.crval2 = average[0].D; + double Xi, Yi; + RD_to_XY (&Xi, &Yi, Ri, Di, &fit->coords); + double radius = hypot(Xi, Yi); + if (radius > ExcludeBogusRadius) { + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; } + continue; + } + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; } + } + + // outlier rejection + if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { + continue; + } + + FitAstromPointInit (&points[Npoints]); + + points[Npoints].R = Ri; + points[Npoints].D = Di; + + // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000. + // T[] is time in years since J2000 (jd = 2451545) + points[Npoints].T = (measure[k].t - fit->T2000) / (86400*365.25) ; // time relative to J2000 in years + + // add measured systematic error in quadrature? only do this after the fit has + // converged (or you will never improve the poor images) + + // dX,dY are the X and Y direction errors in arcseconds. dR, dD are the errors in + // those directions in degrees. IF we have non-circular errors (different values for + // X and Y), then dR and dD will be incorrect: they would need to be rotated to take + // out the position angle + + // dX, dY : error in arcsec: + points[Npoints].dX = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); + points[Npoints].dY = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); // allow a given photcode or measurement to be // ignored if the error is NAN (for photcode, set astromErrSys to NaN) - if (isnan(dX[N])) continue; - if (isnan(dY[N])) continue; - - // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in - // arcseconds. dR, dD are the errors in those directions in degrees. IF we have - // non-circular errors (different values for X and Y), then dR and dD will be - // incorrect: they would need to be rotated to take out the position angle - dR[N] = dX[N] / 3600.0; - dD[N] = dY[N] / 3600.0; - - // XXX use a different flag for stack measurements? - // measure[k].dbFlags |= ID_MEAS_USED_OBJ; - // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } - - N++; + if (isnan(points[Npoints].dX)) continue; + if (isnan(points[Npoints].dY)) continue; + + points[Npoints].dT = measure[k].dt; + + points[Npoints].measure = k; + Npoints++; + + if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) { + has2MASS = TRUE; + } + + myAssert (Npoints <= fit->NpointsAlloc, "oops"); } // loop over measurements : average[0].Nmeasure - // if we have too few good detections for the desired fit, or too limited a - // baseline, use a fit with fewer parameters. XXX if we have too few measurements - // for even the average position, consider including the lower-quality detections? - - // too few measurements for average position (require 2 values) - if (N < 1) return FALSE; // XXX ?? - - // find the mean position - StatType statsR, statsD; - liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) - liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) - - fitAve.Ro = statsR.mean; - fitAve.dRo = 3600.0*statsR.sigma; - - fitAve.Do = statsD.mean; - fitAve.dDo = 3600.0*statsD.sigma; - - fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq); - fitAve.Nfit = N; - - // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE; - fitStats->Nave ++; - - if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", - average[0].R, - average[0].D, - fitAve.Ro, fitAve.Do, - 3600*(average[0].R - fitAve.Ro), - 3600*(average[0].D - fitAve.Do)); - - // make sure that the fit succeeded - int status = TRUE; - status &= finite(fitAve.Ro); - status &= finite(fitAve.Do); - status &= finite(fitAve.dRo); - status &= finite(fitAve.dDo); - if (!status) { - fitStats->Nskip ++; - return FALSE; - } - - // what is the offset relative to the mean fit position? - coords.crval1 = average[0].R; - coords.crval2 = average[0].D; - - double dXoff, dYoff; - RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords); - float dPos = hypot (dXoff, dYoff); - if (dPos > MaxMeanOffset) { - if (fitStats->Noffset < 100) { - fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff); - } - fitStats->Noffset ++; - return FALSE; - } - - // set the stack position values - average[0].Rstk = fitAve.Ro; // RA in degrees - average[0].Dstk = fitAve.Do; // DEC in degrees - average[0].dRstk = fitAve.dRo; // RA scatter in arcsec - average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec - + int TESTPT = FALSE; + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST); + if (TESTPT) { + fprintf (stderr, "got test det\n"); + } + + // XXX flag measurements from stars with 2MASS + for (k = 0; k < average[0].Nmeasure; k++) { + // reset the bit to note that a detection was used (or not) + if (has2MASS) { + measure[k].dbFlags |= ID_MEAS_OBJECT_HAS_2MASS; + } else { + measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_2MASS; + } + } + + fit->Npoints = Npoints; + return TRUE; +} + +int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) { + + int k; + + int Npoints = fitStats->Npoints; + FitAstromPoint *points = fitStats->points; + + // find Tmin & Tmax from the list of accepted measurements + double Tmin = points[0].T; + double Tmax = points[0].T; + double pRmin = +2.0; + double pRmax = -2.0; + double pDmin = +2.0; + double pDmax = -2.0; + + *Tmean = 0.0; + + double Tsum = 0.0; + double Wsum = 0.0; + for (k = 0; k < Npoints; k++) { + Tmin = MIN(Tmin, points[k].T); + Tmax = MAX(Tmax, points[k].T); + + float wx = 1.0 / SQ(points[k].dX); + + Tsum += points[k].T * wx; + Wsum += wx; + + // at this point, T is in years since J2000 + ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T); + pRmin = MIN (pRmin, points[k].pR); + pRmax = MAX (pRmax, points[k].pR); + pDmin = MIN (pDmin, points[k].pD); + pDmax = MAX (pDmax, points[k].pD); + } + *Trange = Tmax - Tmin; + + // mean epoch + *Tmean = Tsum / Wsum; + + // for HIGH_SPEED, just use the center of the range + if (RELASTRO_OP == OP_HIGH_SPEED) { + *Tmean = 0.5*(Tmax - Tmin); + } + + *parRange = hypot (pRmax - pRmin, pDmax - pDmin); + + /* we need to do the fit in a locally linear space; choose a ref coordinate */ + fitStats->coords.crval1 = points[0].R; + fitStats->coords.crval2 = points[0].D; + + // project all of the R,D coordinates to a plane centered on this coordinate. set + // the times to be relative to Tmean + for (k = 0; k < Npoints; k++) { + RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords); + points[k].T -= *Tmean; + } + return TRUE; +} + +int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) { + + int i, j; + + int Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); + } + } + return Nmax; +} + +int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) { + int i; + + // I need to draw Npoints random entries from 'points' with replacement: + for (i = 0; i < Npoints; i++) { + int N = Npoints * drand48(); + sample[i] = points[N]; + } + return TRUE; +} + +// calculate mean and sigma points for the 5 fit parameter +int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) { + + // generate a histogram for the selected element + double *values = NULL; + ALLOCATE (values, double, Nfit); + + int i; + + for (i = 0; i < Nfit; i++) { + switch (mode) { + case FIT_RESULT_RA: + values[i] = fit[i].Ro; + break; + case FIT_RESULT_DEC: + values[i] = fit[i].Do; + break; + case FIT_RESULT_uR: + values[i] = fit[i].uR; + break; + case FIT_RESULT_uD: + values[i] = fit[i].uD; + break; + case FIT_RESULT_PLX: + values[i] = fit[i].p; + break; + default: + myAbort ("invalid option"); + } + } + + dsort (values, Nfit); + + double median; + if (Nfit % 2) { + int Ncenter = Nfit / 2; + median = values[Ncenter]; + } else { + int Ncenter = Nfit / 2 - 1; + median = 0.5*(values[Ncenter] + values[Ncenter + 1]); + } + + double Slo = VectorFractionInterpolate (values, 0.158655, Nfit); + double Shi = VectorFractionInterpolate (values, 0.841345, Nfit); + double sigma = (Shi - Slo) / 2.0; + + switch (mode) { + case FIT_RESULT_RA: + result->Ro = median; + result->dRo = sigma; + break; + case FIT_RESULT_DEC: + result->Do = median; + result->dDo = sigma; + break; + case FIT_RESULT_uR: + result->uR = median; + result->duR = sigma; + break; + case FIT_RESULT_uD: + result->uD = median; + result->duD = sigma; + break; + case FIT_RESULT_PLX: + result->p = median; + result->dp = sigma; + break; + default: + myAbort ("invalid option"); + } + + return TRUE; +} + +double VectorFractionInterpolate (double *values, float fraction, int Npts) { + + float F = fraction * Npts; + int N = fraction * Npts; + + if (N < 0 ) return NAN; + if (N >= Npts - 2) return NAN; + + // interpolate between N,N+1 + + double S = (F - N) * (values[N+1] - values[N]) + values[N]; + return S; +} + +int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) { + + int i; + + // add up the chi square for the fit + double chisq = 0.0; + for (i = 0; i < Npoints; i++) { + double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR; + double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD; + chisq += SQ(points[i].X - Xf) / SQ(points[i].dX); + chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY); + } + switch (mode) { + case FIT_AVERAGE: + fit->chisq = chisq / (2.0*Npoints - 2.0); + break; + case FIT_PM_ONLY: + fit->chisq = chisq / (2.0*Npoints - 4.0); + break; + case FIT_PM_AND_PAR: + fit->chisq = chisq / (2.0*Npoints - 5.0); + break; + default: + myAbort ("invalid mode"); + } + fit->Nfit = Npoints; return (TRUE); } - - -/* fitting proper-motion and parallax: - - given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) - - let x,y be the coordinate in the linearized frame with y parallel to DEC lines - - L,B are the ecliptic longitude and latitude of the object, - dL and dB are the offsets in the L and B directions - - dL = sin(t - topp) - dB = cos(t - topp)*sin(B) - - these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion - for the source in the x,y frame is: - - x = Ro + uR * (t - to) + p * pR - y = Do + uD * (t - to) + p * pD - - the unknowns in these equations are Ro, uR, Do, uD, and p - - XXX think through the concepts for the pole a bit better. all objects near the pole - move the same way with the same phase. choose a projection center and define dL,dB relative - to that center point coordinate system? - -*/ +int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average) { + + int i; + + if (USE_GALAXY_MODEL) { + for (i = 0; i < Nfit; i++) { + fit->uR = average->uRgal; + fit->uD = average->uDgal; + } + } else { + for (i = 0; i < Nfit; i++) { + fit->uR = 0.0; + fit->uD = 0.0; + } + } + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 38601) +++ trunk/Ohana/src/relastro/src/args.c (revision 38986) @@ -253,4 +253,13 @@ } + // If we are looking at the whole sky (or whole relevant sky), use the full image table + // -- this save substantial memory. this could be automatic if the skyregion covers + // more than 2pi. + USE_ALL_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-use-all-images"))) { + remove_argument (N, &argc, argv); + USE_ALL_IMAGES = TRUE; + } + USE_BASIC_CHECK = FALSE; if ((N = get_argument (argc, argv, "-basic-image-search"))) { @@ -271,4 +280,15 @@ remove_argument (N, &argc, argv); APPLY_OFFSETS = TRUE; + } + + N_BOOTSTRAP_SAMPLES = 1; + if ((N = get_argument (argc, argv, "-bootstrap-samples"))) { + remove_argument (N, &argc, argv); + N_BOOTSTRAP_SAMPLES = atoi (argv[N]); + remove_argument (N, &argc, argv); + if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) { + fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n"); + exit (2); + } } @@ -685,4 +705,13 @@ } + // If we are looking at the whole sky (or whole relevant sky), use the full image table + // -- this save substantial memory. this could be automatic if the skyregion covers + // more than 2pi. + USE_ALL_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-use-all-images"))) { + remove_argument (N, &argc, argv); + USE_ALL_IMAGES = TRUE; + } + USE_BASIC_CHECK = FALSE; if ((N = get_argument (argc, argv, "-basic-image-search"))) { @@ -697,4 +726,15 @@ remove_argument (N, &argc, argv); MaxDensityUse = TRUE; + } + + N_BOOTSTRAP_SAMPLES = 1; + if ((N = get_argument (argc, argv, "-bootstrap-samples"))) { + remove_argument (N, &argc, argv); + N_BOOTSTRAP_SAMPLES = atoi (argv[N]); + remove_argument (N, &argc, argv); + if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) { + fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n"); + exit (2); + } } Index: trunk/Ohana/src/relastro/src/assign_images.c =================================================================== --- trunk/Ohana/src/relastro/src/assign_images.c (revision 38601) +++ trunk/Ohana/src/relastro/src/assign_images.c (revision 38986) @@ -26,5 +26,5 @@ char mapfile[DVO_MAX_PATH]; snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR); - AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE); + AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE); // assign images.coords.offsetMap -> table->map[i] Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 38601) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 38986) @@ -112,7 +112,16 @@ // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly - if (USE_GALAXY_MODEL && !catalog[0].average[i].Nstarpar) { - // we may want to add some starpar quality filters... - continue; + if (USE_GALAXY_MODEL) { + // if we use the galaxy model, we need to require the model pm to exist + if (!catalog[0].average[i].Nstarpar) continue; + if (!isfinite(catalog[0].average[i].uRgal)) continue; + if (!isfinite(catalog[0].average[i].uDgal)) continue; + } + + int TESTPT = FALSE; + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (catalog[0].average[i].catID == CAT_ID_SRC) && (catalog[0].average[i].objID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (catalog[0].average[i].catID == CAT_ID_DST) && (catalog[0].average[i].objID == OBJ_ID_DST); + if (TESTPT) { + fprintf (stderr, "got test det\n"); } Index: trunk/Ohana/src/relastro/src/high_speed_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 38601) +++ trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 38986) @@ -153,4 +153,10 @@ PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + strextend (&command, "relastro_client -high-speed %s %s %f %s", PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, outputDir); + strextend (&command, "-hostID %d", table->hosts[i].hostID); + strextend (&command, "-D CATDIR %s", CATDIR); + strextend (&command, "-hostdir %s", table->hosts[i].pathname); + strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + free (outputDir); @@ -166,4 +172,5 @@ if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); } + if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } if (FlagOutlier) { strextend (&command, "-clip %d", CLIP_THRESH); } @@ -176,6 +183,7 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (MinBadQF > 0.0) { strextend (&command, "-min-bad-psfqf %f", MinBadQF); } - if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset %f", MaxMeanOffset); } + if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (WHERE_A[0]) { strextend (&command, "-D WHERE_A \"%s\"", WHERE_A); } Index: trunk/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 38601) +++ trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 38986) @@ -270,6 +270,10 @@ fprintf (stderr, "found %d matches\n", Nmatch); - dvo_catalog_save (&catalogOut, VERBOSE2); - dvo_catalog_unlock (&catalogOut); + + SetProtect (TRUE); + if (!dvo_catalog_save (&catalogOut, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalogOut.filename); exit (1); } + if (!dvo_catalog_unlock (&catalogOut)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalogOut.filename); exit (1); } + SetProtect (FALSE); + dvo_catalog_free (&catalogOut); free (slowMoving); Index: trunk/Ohana/src/relastro/src/hpm_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/hpm_catalogs.c (revision 38601) +++ trunk/Ohana/src/relastro/src/hpm_catalogs.c (revision 38986) @@ -148,6 +148,9 @@ char *command = NULL; - strextend (&command, "relastro_client -hpm %f %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", - RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + strextend (&command, "relastro_client -hpm %f %s", RADIUS, outputDir); + strextend (&command, " -hostID %d", table->hosts[i].hostID); + strextend (&command, " -D CATDIR %s", CATDIR); + strextend (&command, " -hostdir %s", table->hosts[i].pathname); + strextend (&command, " -region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); free (outputDir); @@ -164,4 +167,5 @@ if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); } + if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } if (FlagOutlier) { strextend (&command, "-clip %d", CLIP_THRESH); } @@ -174,6 +178,7 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (MinBadQF > 0.0) { strextend (&command, "-min-bad-psfqf %f", MinBadQF); } - if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset %f", MaxMeanOffset); } + if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (TimeSelect) { Index: trunk/Ohana/src/relastro/src/hpm_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/hpm_objects.c (revision 38601) +++ trunk/Ohana/src/relastro/src/hpm_objects.c (revision 38986) @@ -110,6 +110,10 @@ if (catalog[0].Naverage == Nslow) { fprintf (stderr, "no possible fast objects, skipping this catalog\n"); - dvo_catalog_save (&catalogOut, VERBOSE2); - dvo_catalog_unlock (&catalogOut); + + SetProtect (TRUE); + if (!dvo_catalog_save (&catalogOut, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalogOut.filename); exit (1); } + if (!dvo_catalog_unlock (&catalogOut)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalogOut.filename); exit (1); } + SetProtect (FALSE); + dvo_catalog_free (&catalogOut); free (slowMoving); Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 38601) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 38986) @@ -6,5 +6,5 @@ int i; - // do not remote the sync and fits files if we do a manual run -- user must clear if needed + // do not remove the sync and fits files if we do a manual run -- user must clear if needed if (!PARALLEL_REGIONS_MANUAL) { // clear the I/O files Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 38601) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 38986) @@ -163,6 +163,12 @@ char *command = NULL; - strextend (&command, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f -D RELASTRO_SIGMA_LIM %f", - table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM); + strextend (&command, "relastro_client -load-objects %s", table->hosts[i].results); + strextend (&command, " -hostID %d", table->hosts[i].hostID); + strextend (&command, " -D CATDIR %s", CATDIR); + strextend (&command, " -hostdir %s", table->hosts[i].pathname); + strextend (&command, " -region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + strextend (&command, " -statmode %s", STATMODE); + strextend (&command, " -minerror %f", MIN_ERROR); + strextend (&command, " -D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm"); @@ -188,4 +194,5 @@ } + if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 38601) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 38986) @@ -7,9 +7,9 @@ */ -int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) { +int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap, int UseAllImages) { Image *image, *subset; off_t Nimage, Nsubset; - off_t *LineNumber; + off_t *LineNumber, i; INITTIME; @@ -29,5 +29,5 @@ char mapfile[DVO_MAX_PATH]; snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR); - table = AstromOffsetMapLoad (mapfile, VERBOSE); + table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE); // assign images.coords.offsetMap -> table->map[i] @@ -39,15 +39,22 @@ // select the images which overlap the selected sky regions - subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap); - MARKTIME(" select images: %f sec\n", dtime); + if (UseAllImages) { + ALLOCATE (LineNumber, off_t, Nimage); + for (i = 0; i < Nimage; i++) LineNumber[i] = i; + subset = image; + Nsubset = Nimage; + } else { + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap); + MARKTIME(" select images: %f sec\n", dtime); - if (Nsubset == Nimage) { - free (subset); - subset = image; + if (Nsubset == Nimage) { + free (subset); + subset = image; + } } - + initImages (subset, LineNumber, Nsubset); MARKTIME(" init images: %f sec\n", dtime); - + initMosaics (subset, Nsubset); MARKTIME(" init mosaics: %f sec\n", dtime); Index: trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro.c (revision 38986) @@ -4,4 +4,5 @@ /* get configuration info, args */ + SetSignals (); initialize (argc, argv); Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 38986) @@ -22,4 +22,5 @@ // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) + SetSignals (); initialize_client (argc, argv); @@ -78,5 +79,5 @@ /* load regions and images based on specified sky patch (default depth) */ - load_images (&db, skylist, FALSE); + load_images (&db, skylist, FALSE, USE_ALL_IMAGES); // I can free the database after I have loaded the images... Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 38986) @@ -21,5 +21,5 @@ /* load regions and images based on specified sky patch (default depth) (require full overlap) */ - load_images (&db, skylist, TRUE); + load_images (&db, skylist, TRUE, USE_ALL_IMAGES); MARKTIME("load images: %f sec\n", dtime); @@ -73,11 +73,6 @@ if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog); for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog, (i > 0)); // calculate ,; if (i > 0), apply Galaxy Motion Model (if desired) - if ((i > 1) || !USE_GALAXY_MODEL) { - // if GALAXY_MODEL is selected, we want to delay the frame correction until we have - // applied the galaxy model a couple of times. - FrameCorrectionSerial (catalog, Ncatalog); - } - UpdateChips (catalog, Ncatalog); // measure.X,Y -> R,D, fit image.coords + UpdateObjects (catalog, Ncatalog, i); // calculate ,; if (i > 0), apply Galaxy Motion Model (if desired) + UpdateChips (catalog, Ncatalog, i); // measure.X,Y -> R,D, fit image.coords MARKTIME("update chips: %f sec\n", dtime); } @@ -103,6 +98,11 @@ freeImageBins (1); - // XXX why do I do this exactly? - reload_images (&db); + // If we did NOT use all images, then we applied the measured corrections to a subset of + // images. we now need to set the images associated with db to have those values so + // they will be written out by the dvo_image_update() operations. If we USE_ALL_IMAGES, + // then we do not need to do this, and we should call dvo_image_save (not dvo_image_update) + if (!USE_ALL_IMAGES) { + reload_images (&db); + } if (PARALLEL) { @@ -110,15 +110,34 @@ // need to also save the image map table... save_astrom_table (); - dvo_image_update (&db, VERBOSE); + if (USE_ALL_IMAGES) { + dvo_image_save (&db, VERBOSE); + } else { + dvo_image_update (&db, VERBOSE); + } dvo_image_unlock (&db); } + // NOTE: if we have parallel partitions, then we need to save the Images.dat and + // AstroMap.fits tables BEFORE the remote relastro_clients are launched (they load + // Images.dat and AstroMap.fits) + + // if we do NOT have parallel partitions, we must NOT write them out yet: the act of + // writing out the files byte-swaps the data and makes the values invalide for the + // following calls to UpdateObjectOffsets (which attempt to apply the image values from + // the structure in memory) + // iterate over catalogs to make detection coordinates consistant - UpdateObjectOffsets (skylist, 0, NULL); + if (APPLY_OFFSETS) { + UpdateObjectOffsets (skylist, 0, NULL); + } if (!PARALLEL) { // save the updated image parameters save_astrom_table (); - dvo_image_update (&db, VERBOSE); + if (USE_ALL_IMAGES) { + dvo_image_save (&db, VERBOSE); + } else { + dvo_image_update (&db, VERBOSE); + } dvo_image_unlock (&db); } Index: trunk/Ohana/src/relastro/src/relastro_merge_source.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_merge_source.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro_merge_source.c (revision 38986) @@ -122,6 +122,8 @@ resort_catalog (&catalog_src); - dvo_catalog_save (&catalog_src, VERBOSE2); - dvo_catalog_unlock (&catalog_src); + SetProtect (TRUE); + if (!dvo_catalog_save (&catalog_src, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog_src.filename); exit (1); } + if (!dvo_catalog_unlock (&catalog_src)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog_src.filename); exit (1); } + SetProtect (FALSE); dvo_catalog_free (&catalog_src); Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 38986) @@ -150,6 +150,10 @@ char *command = NULL; - strextend (&command, "relastro_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s", - table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE); + strextend (&command, "relastro_client -update-objects"); + strextend (&command, "-hostID %d", table->hosts[i].hostID); + strextend (&command, "-D CATDIR %s", CATDIR); + strextend (&command, "-hostdir %s", table->hosts[i].pathname); + strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + strextend (&command, "-statmode %s", STATMODE); if (FIT_MODE == FIT_PM_ONLY) { strextend (&command, "-pm"); } @@ -164,4 +168,6 @@ if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); } if (FlagOutlier) { strextend (&command, "-clip %d", CLIP_THRESH); } + + if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } if (USE_FIXED_PIXCOORDS) { strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); } if (PHOTCODE_KEEP_LIST) { strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); } @@ -172,6 +178,7 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (MinBadQF > 0.0) { strextend (&command, "-min-bad-psfqf %f", MinBadQF); } - if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset %f", MaxMeanOffset); } + if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (TimeSelect) { Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 38601) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 38986) @@ -40,5 +40,5 @@ char mapfile[DVO_MAX_PATH]; snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.%d.fits", CATDIR, REGION_HOST_ID); - AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE); + AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE); // assign images.coords.offsetMap -> table->map[i] @@ -120,12 +120,7 @@ } for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog, (i > 0)); + UpdateObjects (catalog, Ncatalog, i); LOGRTIME("UpdateObjects loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime); - if ((i > 1) || !USE_GALAXY_MODEL) { - // if GALAXY_MODEL is selected, we want to delay the frame correction until we have - // applied the galaxy model a couple of times. - FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i); - LOGRTIME("FrameCorrection loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime); - } + share_mean_pos (catalog, Ncatalog, regionHosts, i); LOGRTIME("share_mean_pos loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime); @@ -133,5 +128,5 @@ LOGRTIME("slurp_mean_pos loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime); - UpdateChips (catalog, Ncatalog); + UpdateChips (catalog, Ncatalog, i); LOGRTIME("UpdateChips loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime); Index: trunk/Ohana/src/relastro/src/save_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/save_catalogs.c (revision 38601) +++ trunk/Ohana/src/relastro/src/save_catalogs.c (revision 38986) @@ -9,6 +9,8 @@ if (VERBOSE2) fprintf (stderr, "saving catalog %s\n", catalog[i].filename); - dvo_catalog_save (&catalog[i], VERBOSE2); - dvo_catalog_unlock (&catalog[i]); + SetProtect (TRUE); + if (!dvo_catalog_save (&catalog[i], VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog[i].filename); exit (1); } + if (!dvo_catalog_unlock (&catalog[i])) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog[i].filename); exit (1); } + SetProtect (FALSE); dvo_catalog_free (&catalog[i]); } Index: trunk/Ohana/src/relastro/src/share_images_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_images_pos.c (revision 38601) +++ trunk/Ohana/src/relastro/src/share_images_pos.c (revision 38986) @@ -115,5 +115,5 @@ snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, regionHosts->hosts[i].hostID); - AstromOffsetTable *table = AstromOffsetMapLoad (mapname, VERBOSE); + AstromOffsetTable *table = AstromOffsetMapLoad (mapname, 100000, VERBOSE); LOGRTIME("image_maps_load host %d loop %d on %s, host %d: %f sec\n", i, nloop, myHostName, REGION_HOST_ID, dtime);