Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/Makefile =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/Makefile (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/Makefile (revision 12220) @@ -21,35 +21,40 @@ RELASTRO = \ -$(SRC)/relastro.$(ARCH).o \ -$(SRC)/initialize.$(ARCH).o \ $(SRC)/ConfigInit.$(ARCH).o \ -$(SRC)/args.$(ARCH).o \ -$(SRC)/liststats.$(ARCH).o \ -$(SRC)/load_images.$(ARCH).o \ -$(SRC)/name_region.$(ARCH).o \ -$(SRC)/find_images.$(ARCH).o \ -$(SRC)/find_regions.$(ARCH).o \ -$(SRC)/get_regions.$(ARCH).o \ -$(SRC)/load_catalogs.$(ARCH).o \ -$(SRC)/gcatalog.$(ARCH).o \ -$(SRC)/bcatalog.$(ARCH).o \ -$(SRC)/GridOps.v2.$(ARCH).o \ -$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/FitChip.$(ARCH).o \ +$(SRC)/FitMosaic.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/UpdateChips.$(ARCH).o \ +$(SRC)/UpdateMosaic.$(ARCH).o \ +$(SRC)/UpdateObjects.$(ARCH).o \ +$(SRC)/UpdateSimple.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/dvo_astrom_ops.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ $(SRC)/global_stats.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ $(SRC)/plot_scatter.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ -$(SRC)/misc.$(ARCH).o \ +$(SRC)/relastro.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ -$(SRC)/image-db.$(ARCH).o \ -$(SRC)/free_catalogs.$(ARCH).o \ -$(SRC)/check_permissions.$(ARCH).o \ +$(SRC)/setExclusions.$(ARCH).o \ +$(SRC)/write_coords.$(ARCH).o + +DROP = \ $(SRC)/setMrelFinal.$(ARCH).o \ -$(SRC)/write_coords.$(ARCH).o \ -$(SRC)/setExclusions.$(ARCH).o \ -$(SRC)/Shutdown.$(ARCH).o \ -$(SRC)/SetSignals.$(ARCH).o \ -$(SRC)/wcatalog.$(ARCH).o + $(RELASTRO): $(INC)/relastro.h Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h (revision 12220) @@ -4,15 +4,25 @@ # include +typedef struct { + double R, D; /* Sky Coords - degrees */ + double P, Q; /* Tangent Plane - pixels */ + double L, M; /* Focal Plane - pixels */ + double X, Y; /* Chip Coords - pixels */ + double Mag, dMag; + int mask; +} StarData; + // structure to hold coordinate fitting terms typedef struct { int Npts; int Nterms; - int Npower; int Norder; + int Nsums; + int Nelems; + double **sum; double **xsum; double **ysum; - - double **vector; - double **matrix; + double **xfit; + double **yfit; } CoordFit; @@ -28,5 +38,5 @@ double chisq; int Nfit; -} PMfit; +} PMFit; typedef struct { @@ -62,4 +72,6 @@ char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ char CameraConfig[256]; +char SKY_TABLE[256]; +int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ double MAG_LIM; @@ -108,14 +120,19 @@ enum {black, white, red, orange, yellow, green, blue, indigo, violet}; +int DophotSelect, DophotValue; + int TimeSelect; time_t TSTART, TSTOP; -# ifdef GRID_V1 -int setGridMeasure (int meas, int cat, double X, double Y); -# endif - -# ifdef GRID_V2 -int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum); -# endif +int FITTING_MODE; +enum {FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR}; + +SkyRegion UserPatch; +int UserPatchSelect; + +int DoUpdateObjects; +int DoUpdateSimple; +int DoUpdateChips; +int DoUpdateMosaics; /*** relphot prototypes ***/ @@ -126,5 +143,5 @@ int *SelectRefMosaic PROTO((Mosaic **refmosaic, int *Nimage)); int args PROTO((int argc, char **argv)); -int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion)); +int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); void clean_images PROTO(()); void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final)); @@ -171,6 +188,7 @@ void initstats PROTO((char *mode)); int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); -Catalog *load_catalogs PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion)); -GSCRegion *load_images PROTO((FITS_DB *db, char *seed, int *nregion, GSCRegion *fullregion)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog)); +SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect)); +Image *select_images PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage)); void check_permissions (char *basefile); @@ -195,5 +213,5 @@ void plot_star_coords PROTO((Catalog *catalog, int Ncatalog)); void plot_stars PROTO((Catalog *catalog, int Ncatalog)); -void reload_catalogs PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion)); +void reload_catalogs PROTO((SkyList *skylist)); int setExclusions PROTO((Catalog *catalog, int Ncatalog)); void setMcal PROTO((Catalog *catalog, int Poor)); @@ -226,2 +244,37 @@ void wimages PROTO(()); void write_coords PROTO((Header *header, Coords *coords)); + + +double **array_init (int Nx, int Ny); +void array_free (double **array, int Nx); +CoordFit *fit_init (int order); +void fit_free (CoordFit *fit); +void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt); +void fit_eval (CoordFit *fit); +double **poly2d_dx (double **poly, int Nx, int Ny); +double **poly2d_dy (double **poly, int Nx, int Ny); +double **poly2d_copy (double **poly, int Nx, int Ny); +double poly2d_eval (double **poly, int Nx, int Ny, double x, double y); +void fit_apply_coords (CoordFit *fit, Coords *coords); +int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo); +CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo); +void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void initObjectData (Catalog *catalog, int Ncatalog); +int UpdateObjects (Catalog *catalog, int Ncatalog); +int UpdateSimple (Catalog *catalog, int Ncatalog); +int UpdateChips (Catalog *catalog, int Ncatalog); +int UpdateMosaic (Catalog *catalog, int Ncatalog); +int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); +int ParFactor (double *pR, double *pD, double R, double D, time_t T); +int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts); +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts); + +Coords *getMosaicForImage (int N); + +StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic); +StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic); + +double getMeanR (Measure *measure, Average *average, SecFilt *secfilt); +double getMeanD (Measure *measure, Average *average, SecFilt *secfilt); Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c (revision 12220) @@ -4,12 +4,14 @@ int i; + CoordFit *fit; - fit_init (coords[0].Npolyterms); + fit = fit_init (coords[0].Npolyterms); for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - fit_add (raw[i].L, raw[i].M, ref[i].P, ref[i].Q); + fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0); } - fit_eval (); - fit_apply_coords (coords); + fit_eval (fit); + fit_apply_coords (fit, coords); + fit_free (fit); /* FitChip and FitSimple update the coords in different ways? maybe not... */ } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c (revision 12220) @@ -2,18 +2,16 @@ /* 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 FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts) { + + int i; double **A, **B; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; /* do I need to do this as 2 2x2 matrix equations? */ - ALLOCATE (A, double *, 4); - ALLOCATE (B, double *, 4); - for (i = 0; i < 4; i++) { - ALLOCATE (A[i], double, 4); - ALLOCATE (B[i], double, 1); - bzero (A[i], 4*sizeof(double)); - bzero (B[i], 1*sizeof(double)); - } + A = array_init (4, 4); + B = array_init (4, 1); + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; for (i = 0; i < Npts; i++) { /* handle case where dX or dY = 0.0 */ @@ -56,8 +54,8 @@ dgaussj (A, 4, B, 1); - fit[0].Ro = B[0]; - fit[0].uR = B[1]; - fit[0].Do = B[2]; - fit[0].uD = B[3]; + 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; @@ -68,10 +66,6 @@ fit[0].dp = 0.0; - for (i = 0; i < 4; i++) { - free (A[i]); - free (B[i]); - } - free (A); - free (B); + array_free (A, 4); + array_free (B, 4); return (TRUE); Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c (revision 12220) @@ -2,17 +2,17 @@ /* 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, double *pR, double *pD, int Npts) { +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts) { + + 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; - ALLOCATE (A, double *, 5); - ALLOCATE (B, double *, 5); - for (i = 0; i < 5; i++) { - ALLOCATE (A[i], double, 5); - ALLOCATE (B[i], double, 1); - bzero (A[i], 5*sizeof(double)); - bzero (B[i], 1*sizeof(double)); - } + A = array_init (5, 5); + B = array_init (5, 1); + 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++) { /* handle case where dX or dY = 0.0 */ @@ -78,9 +78,9 @@ dgaussj (A, 5, B, 1); - fit[0].Ro = B[0]; - fit[0].uR = B[1]; - fit[0].Do = B[2]; - fit[0].uD = B[3]; - fit[0].p = B[4]; + 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[0].dRo = sqrt(A[0][0]); @@ -90,10 +90,6 @@ fit[0].dp = sqrt(A[4][4]); - for (i = 0; i < 5; i++) { - free (A[i]); - free (B[i]); - } - free (A); - free (B); + array_free (A, 5); + array_free (B, 5); /* get the chisq from the matrix values */ Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c (revision 12220) @@ -4,12 +4,14 @@ int i; + CoordFit *fit; - fit_init (coords[0].Npolyterms); + fit = fit_init (coords[0].Npolyterms); for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M); + fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0); } - fit_eval (); - fit_apply_coords (coords); + fit_eval (fit); + fit_apply_coords (fit, coords); + fit_free (fit); } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c (revision 12220) @@ -94,7 +94,5 @@ void matchImage (Catalog *catalog, int meas, int cat) { - int i, ave, ccdnum; - double ra, dec, X, Y; - char *pname, *filter, *p, base[256]; + int i; Measure *measure; @@ -178,10 +176,15 @@ } -StarData *getImageRaw (int im, int *Nstars, int isMosaic) { - +/* XXX re-write to use image->nstars if that is possible */ +StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) { + + int i, m, c; + + Coords *mosaic; StarData *raw; ALLOCATE (raw, StarData, Nlist[im]); + mosaic = NULL; if (isMosaic) { mosaic = getMosaicForImage (im); @@ -202,6 +205,6 @@ XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); if (isMosaic) { - XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, &mosaic); - LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &mosaic); + XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, mosaic); + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, mosaic); } else { raw[i].P = raw[i].L; @@ -215,10 +218,14 @@ } -StarData *getImageRef (int im, int *Nstars, int isMosaic) { - +StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) { + + int i, m, c, n; + + Coords *mosaic; StarData *ref; ALLOCATE (ref, StarData, Nlist[im]); + mosaic = NULL; if (isMosaic) { mosaic = getMosaicForImage (im); @@ -237,6 +244,6 @@ RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords); if (isMosaic) { - LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, &mosaic); - LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &mosaic); + LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, mosaic); + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, mosaic); } else { ref[i].L = ref[i].P; Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c (revision 12220) @@ -28,5 +28,7 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) { - double n, L, g, lambda,epsilon,alpha,delta,x,y,z; + double n, L, g; + +# define J2000 2451545. /* Julian date at standard epoch */ n = jd - J2000; @@ -43,5 +45,5 @@ double jd; - double L, B; + double L, B, E, e, s, r, d; /* given a time T in UNIX seconds, determine the solar longitude S */ Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/StarOps.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/StarOps.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/StarOps.c (revision 12220) @@ -28,5 +28,5 @@ if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG); - value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M; + value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M; return (value); } @@ -65,13 +65,7 @@ liststats (list, dlist, N, &stats); - if (PhotPrimary) { - catalog[i].average[j].M = stats.mean; - catalog[i].average[j].dM = stats.sigma; - catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; - } else { - catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; - catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma; - catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; - } + catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma; + catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; } } @@ -121,13 +115,7 @@ /* use sigma or error in dM for output? */ - if (PhotPrimary) { - catalog[i].average[j].M = stats.mean; - catalog[i].average[j].dM = MAX (stats.sigma, stats.error); - catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; - } else { - catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; - catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma); - catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; - } + catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma); + catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; } } @@ -183,9 +171,9 @@ for (j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].average[j].code & STAR_BAD) continue; - Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; if (Xm == -1) continue; Chisq = pow (10.0, 0.01*Xm); xlist[Ntot] = Chisq; - slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; dlist[Ntot] = 1; Ntot ++; @@ -203,6 +191,6 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM; - Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; Chisq = pow (10.0, 0.01*Xm); mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq); @@ -387,5 +375,5 @@ if (catalog[i].average[j].code & STAR_BAD) continue; - Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; if (Xm == NO_MAG) continue; list[n] = pow (10.0, 0.01*Xm); @@ -423,5 +411,5 @@ if (catalog[i].average[j].code & STAR_BAD) continue; - dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; list[n] = dM; dlist[n] = 1; @@ -452,5 +440,5 @@ for (j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].average[j].code & STAR_BAD) continue; - dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; bin = dMrel / 0.00025; bin = MAX (0, MIN (NBIN-1, bin)); @@ -483,6 +471,6 @@ for (j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].average[j].code & STAR_BAD) continue; - xlist[N] = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M; - value = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M; + value = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; if (value == NO_MAG) continue; ylist[N] = 0.01*value; Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c (revision 12220) @@ -1,10 +1,10 @@ # include "relastro.h" -int UpdateChips () { +int UpdateChips (Catalog *catalog, int Ncatalog) { /* we can measure new image parameters for each non-mosaic chip independently */ - + int i, Nimage, Nstars; Image *image; - int Nimage; + StarData *raw, *ref; image = getimages (&Nimage); @@ -18,8 +18,8 @@ /* convert measure coordinates to raw entries */ - raw = getImageRaw (&image[i], &Nstars, TRUE); + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE); /* convert average coordinates to ref entries */ - ref = getImageRef (&image[i], &Nstars, TRUE); + ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE); /* XXX grab this code from mosastro */ Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c (revision 12220) @@ -1,10 +1,10 @@ # include "relastro.h" -int UpdateMosaic () { +int UpdateMosaic (Catalog *catalog, int Ncatalog) { /* we can measure new image parameters for each mosaic independently */ - + int i, Nimage, Nstars; Image *image; - int Nimage; + StarData *raw, *ref; image = getimages (&Nimage); @@ -18,8 +18,8 @@ /* convert measure coordinates to raw entries */ - raw = getImageRaw (&image[i], &Nstars, TRUE); + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE); /* convert average coordinates to ref entries */ - ref = getImageRef (&image[i], &Nstars, TRUE); + ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE); /* XXX grab this code from mosastro */ Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c (revision 12220) @@ -6,4 +6,6 @@ static double *R, *dR; static double *D, *dD; +static double *pX; +static double *pY; static time_t *T; static double *dT; @@ -25,6 +27,4 @@ ALLOCATE (X, double, MAX (1, Nmax)); ALLOCATE (Y, double, MAX (1, Nmax)); - ALLOCATE (pR, double, MAX (1, Nmax)); - ALLOCATE (pD, double, MAX (1, Nmax)); ALLOCATE (dR, double, MAX (1, Nmax)); @@ -33,4 +33,7 @@ ALLOCATE (dX, double, MAX (1, Nmax)); ALLOCATE (dY, double, MAX (1, Nmax)); + + ALLOCATE (pX, double, MAX (1, Nmax)); + ALLOCATE (pY, double, MAX (1, Nmax)); } @@ -38,5 +41,4 @@ int i, j, k, m, N; - float chisq; StatType statsR, statsD; Coords coords; @@ -88,5 +90,4 @@ /* project all of the R,D coordinates to a plane centered on this coordinate */ - /* calculate pR[i], pD[i] for each point */ for (k = 0; k < N; k++) { RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); @@ -100,8 +101,8 @@ fit.Ro = statsR.mean; - fit.dR = statsR.sigma; + fit.dRo = statsR.sigma; fit.Do = statsD.mean; - fit.dD = statsD.sigma; + fit.dDo = statsD.sigma; fit.chisq = 0.5*(statsR.chisq + statsD.chisq); @@ -117,5 +118,5 @@ ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]); } - FitPMandPar (&fit, X, dX, Y, dY, T, pR, pD, N); + FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N); break; default: @@ -137,5 +138,5 @@ catalog[i].average[j].dP = fit.dp; - catalog[i].average[j].Xp = (fit.Npts > 1) ? 100.0*log10(fit.chisq) : NO_MAG; + catalog[i].average[j].Xp = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG; } } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c (revision 12220) @@ -1,10 +1,10 @@ # include "relastro.h" -int UpdateSimple () { +int UpdateSimple (Catalog *catalog, int Ncatalog) { /* we can measure new image parameters for each non-mosaic chip independently */ - + int i, Nimage, Nstars; Image *image; - int Nimage; + StarData *raw, *ref; image = getimages (&Nimage); @@ -19,8 +19,8 @@ /* convert measure coordinates to raw entries */ - raw = getImageRaw (&image[i], &Nstars, FALSE); + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, FALSE); /* convert average coordinates to ref entries */ - ref = getImageRef (&image[i], &Nstars, FALSE); + ref = getImageRef (catalog, Ncatalog, i, &Nstars, FALSE); /* XXX grab this code from mosastro */ Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/args.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/args.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/args.c (revision 12220) @@ -174,13 +174,31 @@ /* possible operations */ - UpdateObjects = FALSE; + DoUpdateObjects = FALSE; if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - UpdateObjects = TRUE; - } - UpdateImages = FALSE; - if ((N = get_argument (argc, argv, "-update-images"))) { - remove_argument (N, &argc, argv); - UpdateImages = TRUE; + DoUpdateObjects = TRUE; + } + DoUpdateSimple = FALSE; + if ((N = get_argument (argc, argv, "-update-simple"))) { + remove_argument (N, &argc, argv); + DoUpdateSimple = TRUE; + } + DoUpdateChips = FALSE; + if ((N = get_argument (argc, argv, "-update-chips"))) { + remove_argument (N, &argc, argv); + DoUpdateChips = TRUE; + } + DoUpdateMosaics = FALSE; + if ((N = get_argument (argc, argv, "-update-mosaics"))) { + remove_argument (N, &argc, argv); + DoUpdateMosaics = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/bcatalog.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/bcatalog.c (revision 12220) @@ -19,11 +19,4 @@ if (catalog[0].average[i].Nm < 2) continue; - /* XXX this limitation is absurd - if (catalog[0].average[i].R < fullregion[0].RA[0]) continue; - if (catalog[0].average[i].R > fullregion[0].RA[1]) continue; - if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue; - if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue; - */ - /* start with all stars good */ subcatalog[0].average[Naverage] = catalog[0].average[i]; @@ -34,9 +27,5 @@ if (RESET) { - float *p; - p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M; - *p = NO_MAG; - p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM; - *p = NO_MAG; + // XXX reset the ra,dec coords? subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW; subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR; @@ -59,5 +48,5 @@ /* select measurements by quality */ - if (catalog[0].measure[offset].dophot != 1) continue; + if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue; /* select measurements by mag limit */ @@ -66,5 +55,5 @@ /* select measurements by measurement error */ - if (catalog[0].measure[offset].dM > SIGMA_LIM) continue; + if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue; /* select measurements by mag limit */ @@ -111,4 +100,5 @@ subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + // XXX really do this or not??? if (catalog[0].Naverage > 0) free (catalog[0].average); if (catalog[0].Nmeasure > 0) free (catalog[0].measure); Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c (revision 12220) @@ -12,8 +12,31 @@ */ +double **array_init (int Nx, int Ny) { + + int i; + double **array; + + ALLOCATE (array, double *, Nx); + for (i = 0; i < Nx; i++) { + ALLOCATE (array[i], double, Ny); + memset (array[i], 0, Ny*sizeof(double)); + } + return (array); +} + +void array_free (double **array, int Nx) { + + int i; + + for (i = 0; i < Nx; i++) { + free (array[i]); + } + free (array); +} + // XXX define a fit structure and drop the file static variables? CoordFit *fit_init (int order) { - int i; + CoordFit *fit; ALLOCATE (fit, CoordFit, 1); @@ -29,32 +52,30 @@ // xsum[i][j] holds \sum (x2 wt x1^i y1^j) // ysum[i][j] holds \sum (y2 wt x1^i y1^j) - ALLOCATE (fit[0].xsum, double *, fit[0].Nterms); - ALLOCATE (fit[0].ysum, double *, fit[0].Nterms); - for (i = 0; i < fit[0].Nterms; i++) { - ALLOCATE (fit[0].xsum[i], double, fit[0].Nterms); - ALLOCATE (fit[0].ysum[i], double, fit[0].Nterms); - bzero (fit[0].xsum[i], fit[0].Nterms*sizeof(double)); - bzero (fit[0].ysum[i], fit[0].Nterms*sizeof(double)); - } + fit[0].xsum = array_init (fit[0].Nterms, fit[0].Nterms); + fit[0].ysum = array_init (fit[0].Nterms, fit[0].Nterms); + + // xfit[i][j] holds x2 coeff for x1^i y1^j + // yfit[i][j] holds y2 coeff for x1^i y1^j + fit[0].xfit = array_init (fit[0].Nterms, fit[0].Nterms); + fit[0].yfit = array_init (fit[0].Nterms, fit[0].Nterms); // sum[i][j] holds \sum (wt x1^i y1^j) - ALLOCATE (fit[0].sum, double *, fit[0].Nsums); - for (i = 0; i < fit[0].Nsums; i++) { - ALLOCATE (fit[0].sum[i], double, fit[0].Nsums); - bzero (fit[0].sum[i], fit[0].Nsums*sizeof(double)); - } - - // matrix, vector hold the final linear system - ALLOCATE (fit[0].matrix, double *, fit[0].Nelems); - ALLOCATE (fit[0].vector, double *, fit[0].Nelems); - for (i = 0; i < fit[0].Nelems; i++) { - ALLOCATE (fit[0].matrix[i], double, fit[0].Nelems); - ALLOCATE (fit[0].vector[i], double, fit[0].Nelems); - bzero (fit[0].matrix[i], fit[0].Nelems*sizeof(double)); - bzero (fit[0].vector[i], 2*sizeof(double)); - } + fit[0].sum = array_init (fit[0].Nsums, fit[0].Nsums); + return (fit); } +void fit_free (CoordFit *fit) { + + array_free (fit[0].xfit, fit[0].Nterms); + array_free (fit[0].yfit, fit[0].Nterms); + + array_free (fit[0].sum, fit[0].Nsums); + array_free (fit[0].xsum, fit[0].Nterms); + array_free (fit[0].ysum, fit[0].Nterms); + + free (fit); +} + // XXX use implicit masks as below or explicit masks (with function to set?) // XXX eg, add a global mask to this file and @@ -84,9 +105,14 @@ void fit_eval (CoordFit *fit) { - int i, j, n, m, M, N; + int i, j, ix, iy, jx, jy; + double **matrix, **vector; if (fit[0].Npts == 0) { fprintf (stderr, "warning: no valid pts\n"); } + + // matrix, vector hold the final linear system + matrix = array_init (fit[0].Nelems, fit[0].Nelems); + vector = array_init (fit[0].Nelems, 2); /* remap the xsum,ysum terms into the vector */ @@ -94,15 +120,15 @@ ix = i % fit[0].Nterms; iy = i / fit[0].Nterms; - vector[i][0] = xsum[ix][iy]; - vector[i][1] = ysum[ix][iy]; + vector[i][0] = fit[0].xsum[ix][iy]; + vector[i][1] = fit[0].ysum[ix][iy]; for (j = 0; j < fit[0].Nelems; j++) { jx = j % fit[0].Nterms; jy = j / fit[0].Nterms; - matrix[i][j] = sum[ix+jx][iy+jy]; - } - } - - dgaussj (matrix, NPARS, vector, 2); + matrix[i][j] = fit[0].sum[ix+jx][iy+jy]; + } + } + + dgaussj (matrix, fit[0].Nelems, vector, 2); for (i = 0; i < fit[0].Nelems; i++) { @@ -112,4 +138,91 @@ ix, iy, vector[i][0], ix, iy, vector[i][1]); } + + /* remap the vector terms into xfit,yfit */ + for (i = 0; i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + fit[0].xfit[ix][iy] = vector[i][0]; + fit[0].yfit[ix][iy] = vector[i][1]; + } + + array_free (matrix, fit[0].Nelems); + array_free (vector, fit[0].Nelems); +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_dx (double **poly, int Nx, int Ny) { + + int i, j, Nxout, Nyout; + double **dpoly; + + Nxout = Nx - 1; + Nyout = Ny; + + // poly[i][j] holds coeff for x1^i y1^j + dpoly = array_init (Nxout, Nyout); + + for (i = 0; i < Nxout; i++) { + for (j = 0; j < Nyout; j++) { + dpoly[i][j] = poly[i+1][j] * (i+1); + } + } + return dpoly; +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_dy (double **poly, int Nx, int Ny) { + + int i, j, Nxout, Nyout; + double **dpoly; + + Nxout = Nx; + Nyout = Ny - 1; + + // poly[i][j] holds coeff for x1^i y1^j + dpoly = array_init (Nxout, Nyout); + + for (i = 0; i < Nxout; i++) { + for (j = 0; j < Nyout; j++) { + dpoly[i][j] = poly[i][j+1] * (j+1); + } + } + return dpoly; +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_copy (double **poly, int Nx, int Ny) { + + int i, j; + double **out; + + // poly[i][j] holds coeff for x1^i y1^j + out = array_init (Nx, Ny); + + for (i = 0; i < Nx; i++) { + for (j = 0; j < Ny; j++) { + out[i][j] = poly[i][j+1]; + } + } + return out; +} + +double poly2d_eval (double **poly, int Nx, int Ny, double x, double y) { + + int i, j; + double xsum, ysum, sum; + + sum = 0; + xsum = ysum = 1.0; + + for (i = 0; i < Nx; i++) { + ysum = xsum; + for (j = 0; j < Ny; j++) { + sum += poly[i][j] * ysum; + ysum *= y; + } + xsum *= x; + } + return (sum); } @@ -119,8 +232,8 @@ void fit_apply_coords (CoordFit *fit, Coords *coords) { - int i, j, Np, Nv, N; double c11, c12; double c21, c22; - double R; + double Xo, Yo, R; + CoordFit *modfit; /* I have L,M = fit(X,Y). set corresponding terms for coords */ @@ -129,23 +242,36 @@ // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... - coords[0].crpix1,crpix2 = CoordsGetCenter (fit); - modfit = CoordsSetCenter (fit, xo, yo); - - // set crpix1,2 - fit_to_coordterms (coords, modfit, 0, 0); + CoordsGetCenter (fit, 0.001, &Xo, &Yo); + coords[0].crpix1 = Xo; + coords[0].crpix2 = Yo; + + // resulting fit should have zero constant terms + modfit = CoordsSetCenter (fit, Xo, Yo); // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0) - fit_to_coordterms (coords, modfit, 0, 1); - fit_to_coordterms (coords, modfit, 1, 0); + coords[0].pc1_1 = modfit[0].xfit[1][0]; + coords[0].pc1_2 = modfit[0].xfit[0][1]; + coords[0].pc2_1 = modfit[0].yfit[1][0]; + coords[0].pc2_2 = modfit[0].yfit[0][1]; // set the polyterm elements - fit_to_coordterms (coords, modfit, 0, 2); - fit_to_coordterms (coords, modfit, 1, 1); - fit_to_coordterms (coords, modfit, 2, 0); - - fit_to_coordterms (coords, modfit, 0, 3); - fit_to_coordterms (coords, modfit, 1, 2); - fit_to_coordterms (coords, modfit, 2, 1); - fit_to_coordterms (coords, modfit, 3, 0); + coords[0].polyterms[0][0] = modfit[0].xfit[2][0]; + coords[0].polyterms[1][0] = modfit[0].xfit[1][1]; + coords[0].polyterms[2][0] = modfit[0].xfit[0][2]; + + coords[0].polyterms[0][1] = modfit[0].yfit[2][0]; + coords[0].polyterms[1][1] = modfit[0].yfit[1][1]; + coords[0].polyterms[2][1] = modfit[0].yfit[0][2]; + + // I need to validate Norder + coords[0].polyterms[3][0] = modfit[0].xfit[3][0]; + coords[0].polyterms[4][0] = modfit[0].xfit[2][1]; + coords[0].polyterms[5][0] = modfit[0].xfit[1][2]; + coords[0].polyterms[6][0] = modfit[0].xfit[0][3]; + + coords[0].polyterms[3][1] = modfit[0].yfit[3][0]; + coords[0].polyterms[4][1] = modfit[0].yfit[2][1]; + coords[0].polyterms[5][1] = modfit[0].yfit[1][2]; + coords[0].polyterms[6][1] = modfit[0].yfit[0][3]; /* we do not modify crval1,2: these are kept at the default values */ @@ -164,31 +290,6 @@ coords[0].pc2_2 = c22*R; + fit_free (modfit); /* keep the order and type from initial values */ } -void fit_free (CoordFit *fit) { - - int i; - - for (i = 0; i < fit[0].Nterms; i++) { - free (fit[0].xsum[i]); - free (fit[0].ysum[i]); - } - free (fit[0].xsum); - free (fit[0].ysum); - - for (i = 0; i < fit[0].Nsums; i++) { - free (fit[0].sum[i]); - } - free (fit[0].sum); - - for (i = 0; i < fit[0].Nelems; i++) { - free (fit[0].matrix[i]); - free (fit[0].vector[i]); - } - free (fit[0].matrix); - free (fit[0].vector); - - free (fit); -} - Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c (revision 12220) @@ -1,166 +1,146 @@ -# include "mosastro.h" +# include "relastro.h" -fit_to_coordterms (Coords *coords, CoordFit *fit, int nx, int ny) { +double factorial (int N) { - int nsum; + int i; + double F; - nsum = nx + ny; - - switch (nsum) { - case 0: - xterm = &coords[0].crpix1; - yterm = &coords[0].crpix2; - return; - case 1: - if (ny == 0) { - xterm = &coords[0].pc1_1; - yterm = &coords[0].pc2_1; - } else { - xterm = &coords[0].pc1_2; - yterm = &coords[0].pc2_2; - } - return; - case 2: - xterm = &coords[0].polyterm[ny][0]; - yterm = &coords[0].polyterm[ny][1]; - return; - case 3: - xterm = &coords[0].polyterm[ny+3][0]; - yterm = &coords[0].polyterm[ny+3][1]; - return; - default: - fprintf (stderr, "programming error %s:%d\n", __func__, __line__); - exit (2); + F = 1; + for (i = N; i > 1; i--) { + F *= i; } - fprintf (stderr, "programming error %s:%d\n", __func__, __line__); - exit (2); + return F; } - -// XXXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ? - -CoordsGetCenter () { +// XXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ? // given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y // for which L,M = 0,0. tol is the allowed error on x,y. -psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, double tol) -{ +int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo) { - // crpix1,2 = X,Y(crval1,2) - // start with linear solution for Xo,Yo: - double R = (trans->x->coeff[1][0]*trans->y->coeff[0][1] - trans->x->coeff[0][1]*trans->y->coeff[1][0]); - double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R; - double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R; + int i, Nx, Ny; + double R, Xo, Yo, dPos; + double **XdX, **XdY, **YdX, **YdY, **alpha, **beta; + double **xfit, **yfit; - // iterate to actual solution: requires small non-linear terms - if (trans->x->nX > 1) { - psPolynomial2D *XdX = psPolynomial2D_dX(NULL, trans->x); - psPolynomial2D *XdY = psPolynomial2D_dY(NULL, trans->x); + xfit = fit[0].xfit; + yfit = fit[0].yfit; - psPolynomial2D *YdX = psPolynomial2D_dX(NULL, trans->y); - psPolynomial2D *YdY = psPolynomial2D_dY(NULL, trans->y); + // solve for Xo,Yo where L,M(Xo,Yo) = 0,0 + // start with linear solution for Xo,Yo: + R = (xfit[1][0]*yfit[0][1] - xfit[0][1]*yfit[1][0]); + Xo = (yfit[0][0]*xfit[0][1] - xfit[0][0]*yfit[0][1])/R; + Yo = (xfit[0][0]*yfit[1][0] - yfit[0][0]*xfit[1][0])/R; - psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32); - psVector *Beta = psVectorAlloc (2, PS_DATA_F32); + Nx = fit[0].Nterms; + Ny = fit[0].Nterms; - /* this loop uses the Newton-Raphson method to solve for Xo,Yo - * it needs the high order terms to be small - * Xo,Yo are in pixels; - */ - double dPos = tol + 1; - for (int i = 0; (dPos > tol) && (i < 20); i++) { - // NOTE: order for Alpha is: [y][x] - Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo); - Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo); - Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo); - Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo); + // iterate to actual solution: requires small non-linear terms + if (fit[0].Norder > 1) { + XdX = poly2d_dx (xfit, Nx, Ny); + XdY = poly2d_dy (xfit, Nx, Ny); + YdX = poly2d_dx (yfit, Nx, Ny); + YdY = poly2d_dy (yfit, Nx, Ny); - Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo); - Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo); + alpha = array_init (2, 2); + beta = array_init (2, 1); - psMatrixGJSolveF32 (Alpha, Beta); + /* this loop uses the Newton-Raphson method to solve for Xo,Yo + * it needs the high order terms to be small + * Xo,Yo are in pixels; + */ + dPos = tol + 1; + for (i = 0; (dPos > tol) && (i < 20); i++) { + // NOTE: order for alpha is: [y][x] + alpha[0][0] = poly2d_eval (XdX, Nx-1, Ny, Xo, Yo); + alpha[1][0] = poly2d_eval (XdY, Nx, Ny-1, Xo, Yo); + alpha[0][1] = poly2d_eval (YdX, Nx-1, Ny, Xo, Yo); + alpha[1][1] = poly2d_eval (YdY, Nx, Ny-1, Xo, Yo); - Xo -= Beta->data.F32[0]; - Yo -= Beta->data.F32[1]; - dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]); - } - psFree (Alpha); - psFree (Beta); - psFree (XdX); - psFree (XdY); - psFree (YdX); - psFree (YdY); + beta[0][0] = poly2d_eval (xfit, Nx, Ny, Xo, Yo); + beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo); + + dgaussj (alpha, 2, beta, 1); + + Xo -= beta[0][0]; + Yo -= beta[1][0]; + dPos = hypot(beta[0][0], beta[1][0]); } - psPlane *center = psPlaneAlloc (); - center->x = Xo; - center->y = Yo; - - return center; + array_free (alpha, 2); + array_free (beta, 2); + array_free (XdX, Nx - 1); + array_free (XdY, Nx); + array_free (YdX, Nx - 1); + array_free (YdY, Nx); + } + *xo = Xo; + *yo = Yo; + return TRUE; } // convert a transformation L(x,y) to L'(x-xo,y-yo) -psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo) -{ +CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo) { - // validate fit order + int i, j, Nx, Ny; + double **xPx, **yPx, **xPy, **yPy, **tmp; - if (output == NULL) { - output = psPlaneTransformAlloc(input->x->nX, input->x->nY); + CoordFit *output; + output = fit_init (input->Norder); + + /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j + * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and + * noting that the constant term in each is the coefficient in the case of L(x,y) and is the + * value of L'(-xo,-yo) in the second case. + */ + + Nx = input->Nterms; + Ny = input->Nterms; + + xPx = poly2d_copy (input->xfit, Nx, Ny); + yPx = poly2d_copy (input->yfit, Nx, Ny); + + for (i = 0; i <= input->Nterms; i++) { + xPy = poly2d_copy (xPx, Nx, Ny); + yPy = poly2d_copy (yPx, Nx, Ny); + for (j = 0; j <= input->Nterms; j++) { + output->xfit[i][j] = poly2d_eval (xPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j); + output->yfit[i][j] = poly2d_eval (yPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j); + + // take the next derivative wrt y, catch output (is NULL on last pass) + if (Ny > 0) { + tmp = poly2d_dy(xPy, Nx, Ny); + array_free (xPy, Nx); + xPy = tmp; + + tmp = poly2d_dy(yPy, Nx, Ny); + array_free (yPy, Nx); + yPy = tmp; + } else { + array_free (xPy, Nx); + array_free (yPy, Nx); + } + + Ny --; } - - /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j - * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and - * noting that the constant term in each is the coefficient in the case of L(x,y) and is the - * value of L'(-xo,-yo) in the second case. - */ - - psPolynomial2D *tmp; - - psPolynomial2D *xPx = psPolynomial2DCopy (NULL, input->x); - psPolynomial2D *yPx = psPolynomial2DCopy (NULL, input->y); - - for (int i = 0; i <= input->x->nX; i++) { - psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx); - psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx); - for (int j = 0; j <= input->x->nY; j++) { - output->x->mask[i][j] = input->x->mask[i][j]; - output->y->mask[i][j] = input->y->mask[i][j]; - output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1); - output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1); - - // take the next derivative wrt y, catch output (is NULL on last pass) - tmp = psPolynomial2D_dY(NULL, xPy); - psFree (xPy); - xPy = tmp; - tmp = psPolynomial2D_dY(NULL, yPy); - psFree (yPy); - yPy = tmp; - } - // take the next derivative wrt x, catch output (is NULL on last pass) - tmp = psPolynomial2D_dX(NULL, xPx); - psFree (xPx); - xPx = tmp; - tmp = psPolynomial2D_dX(NULL, yPx); - psFree (yPx); - yPx = tmp; + Ny = input->Nterms; + // take the next derivative wrt x, catch output (is NULL on last pass) + if (Nx > 0) { + tmp = poly2d_dx(xPx, Nx, Ny); + array_free (xPx, Nx); + xPx = tmp; + tmp = poly2d_dx(yPx, Nx, Ny); + array_free (yPx, Nx); + yPx = tmp; + Nx --; + } else { + array_free (xPx, Nx); + array_free (yPx, Nx); } - return output; -} - -int mkpolyterm (int n, int m) { - - int i, nt, N; - - N = 0; - nt = n + m; - for (i = 2; i < nt; i++) { - N += i + 1; } - N += m; - return (N); + return output; } /* - coords uses an odd, rigid sequence for the coefficients: + Coords uses a rigid sequence for the coefficients: x^0 y^0 : crpix1,2 (but note this is applied before higher order terms) Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plot_scatter.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plot_scatter.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plot_scatter.c (revision 12220) @@ -32,5 +32,5 @@ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; - Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M; + Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M; xlist[N] = Mrel; ylist[N] = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel; Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plotstuff.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plotstuff.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/plotstuff.c (revision 12220) @@ -49,5 +49,5 @@ if (Xgraph[N] == 0) return; - KiiPS (Xgraph[N], TRUE, KAPA_NEWPLOT, filename); + KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename, "default"); return; } Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c (revision 12219) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c (revision 12220) @@ -3,5 +3,5 @@ int main (int argc, char **argv) { - int i, status, Ncatalog; + int status, Ncatalog; Catalog *catalog; FITS_DB db; @@ -43,15 +43,15 @@ /* major modes */ - if (UpdateObjects) { + if (DoUpdateObjects) { UpdateObjects (catalog, Ncatalog); } - if (UpdateSimple) { - UpdateSimple (); + if (DoUpdateSimple) { + UpdateSimple (catalog, Ncatalog); } - if (UpdateChips) { - UpdateChips (); + if (DoUpdateChips) { + UpdateChips (catalog, Ncatalog); } - if (UpdateMosaics) { - UpdateMosaics (); + if (DoUpdateMosaics) { + UpdateMosaic (catalog, Ncatalog); }