Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37477) @@ -42,4 +42,7 @@ $(SRC)/load_template_images.$(ARCH).o \ $(SRC)/make_fake_images.$(ARCH).o \ +$(SRC)/fit_fake_stars.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/AstromOffsetMapOps.$(ARCH).o \ $(SRC)/airmass.$(ARCH).o \ $(SRC)/get_image_patch.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37477) @@ -9,4 +9,14 @@ typedef struct { + int NfakeImage; + int NFAKEIMAGE; + Image *fakeImage; + + int NtrueImage; + int NTRUEIMAGE; + Image *trueImage; +} ImageInfo; + +typedef struct { double R, D; StarPar starpar; @@ -16,4 +26,6 @@ typedef struct { + double Rref; // "reference" coordinate of the fake stars + double Dref; Average average; Measure measure; @@ -21,4 +33,18 @@ int found; } Stars; + +// structure to hold coordinate fitting terms +typedef struct { + int Npts; + int Nterms; + int Norder; + int Nsums; + int Nelems; + double **sum; + double **xsum; + double **ysum; + double **xfit; + double **yfit; +} CoordFit; /* used in find_matches, find_matches_refstars */ @@ -33,5 +59,4 @@ # define DVO_MAX_PATH 1024 -char ImageCat[DVO_MAX_PATH]; char GSCFILE[DVO_MAX_PATH]; char CATDIR[DVO_MAX_PATH]; @@ -116,5 +141,5 @@ Image *make_fake_images (Image *image, int *nfakeImage); -Image *fakeastro_images_region (Image *image, int *nimage, int *NIMAGE, Image *refImage, int NrefImage, SkyTable *skyTableInput, SkyTable *skyTableOutput, SkyRegion *region); +int fakeastro_images_region (ImageInfo *imageInfo, Image *refImage, int NrefImage, SkyTable *skyTableInput, SkyTable *skyTableOutput, SkyRegion *innerRegion); SkyRegion *get_image_patch (Image *image); @@ -135,2 +160,20 @@ float airmass (float secz_image, double ra, double dec, double st, double latitude); float azimuth (double ha, double dec, double latitude); + +int fit_fake_stars (Stars *stars, int Nstars, Image *image); + +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); +int fit_eval (CoordFit *fit); +void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1); +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); +int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef); + +int save_astrom_table (); +AstromOffsetTable *get_astrom_table (); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/AstromOffsetMapOps.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/AstromOffsetMapOps.c (revision 37477) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/AstromOffsetMapOps.c (revision 37477) @@ -0,0 +1,16 @@ +# include "fakeastro.h" + +static AstromOffsetTable *table = NULL; + +int save_astrom_table () { + + char mapfile[DVO_MAX_PATH]; + snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR); + AstromOffsetMapSave (table, mapfile); + + return TRUE; +} + +AstromOffsetTable *get_astrom_table () { + return table; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37477) @@ -25,6 +25,4 @@ strcpy (CATDIR, tmpcatdir); free (tmpcatdir); - - sprintf (ImageCat, "%s/Images.dat", CATDIR); } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c (revision 37477) @@ -1,27 +1,12 @@ # include "fakeastro.h" +int open_image_db (FITS_DB *db, char *ImageName); +int save_image_db (FITS_DB *db, Image *image, int Nimage); int fakeastro_images () { + FITS_DB dbfake, dbtrue; - FITS_DB db; - - // we are creating a new image table (what about db ID?) - snprintf (ImageCat, DVO_MAX_PATH, "%s/Images.dat", CATDIR_OUTPUT); - - /* setup image table format and lock */ - db.mode = dvo_catalog_catmode (CATMODE); - db.format = dvo_catalog_catformat (CATFORMAT); - int status = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? - if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); - - /* load or create the image table */ - if (db.dbstate == LCK_EMPTY) { - if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); - dvo_image_create (&db, GetZeroPoint()); - } else { - if (!dvo_image_load (&db, VERBOSE, FALSE)) { - Shutdown ("can't read image catalog %s", db.filename); - } - } + open_image_db (&dbfake, "Images.dat"); + open_image_db (&dbtrue, "Images.true.dat"); SkyTable *skyTableInput = SkyTableLoadOptimal (CATDIR_INPUT, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); @@ -31,25 +16,25 @@ SkyTableSetFilenames (skyTableOutput, CATDIR_OUTPUT, "cpt"); - // subset out the DIS exposures here: - IMAGE_ID = 1; int NrefImage; - Image *refImage = load_template_images (&NrefImage); + Image *refImage = load_template_images (&NrefImage); // subset of the DIS exposures - int Nimage = 0; - int NIMAGE = 1000; - Image *image = NULL; - - ALLOCATE (image, Image, NIMAGE); + ImageInfo imageInfo; - image = fakeastro_images_region (image, &Nimage, &NIMAGE, refImage, NrefImage, skyTableInput, skyTableOutput, &UserPatch); + imageInfo.NfakeImage = 0; + imageInfo.NFAKEIMAGE = 1000; + imageInfo.fakeImage = NULL; + ALLOCATE (imageInfo.fakeImage, Image, imageInfo.NFAKEIMAGE); - /* add the new image and save */ - dvo_image_addrows (&db, image, Nimage); - SetProtect (TRUE); - dvo_image_update (&db, VERBOSE); - SetProtect (FALSE); - dvo_image_unlock (&db); /* unlock? */ + imageInfo.NtrueImage = 0; + imageInfo.NTRUEIMAGE = 1000; + imageInfo.trueImage = NULL; + ALLOCATE (imageInfo.trueImage, Image, imageInfo.NTRUEIMAGE); + + fakeastro_images_region (&imageInfo, refImage, NrefImage, skyTableInput, skyTableOutput, &UserPatch); + + save_image_db (&dbfake, imageInfo.fakeImage, imageInfo.NfakeImage); + save_image_db (&dbtrue, imageInfo.trueImage, imageInfo.NtrueImage); exit (0); @@ -69,2 +54,37 @@ // -input catdir -output catdir -images images.fits + +int open_image_db (FITS_DB *db, char *ImageName) { + + char ImageCat[DVO_MAX_PATH]; + + // we are creating a new image table (what about db ID?) + snprintf (ImageCat, DVO_MAX_PATH, "%s/%s", CATDIR_OUTPUT, ImageName); + + /* setup image table format and lock */ + db->mode = dvo_catalog_catmode (CATMODE); + db->format = dvo_catalog_catformat (CATFORMAT); + int status = dvo_image_lock (db, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db->filename); + + /* load or create the image table */ + if (db->dbstate != LCK_EMPTY) { + Shutdown ("image table %s exists, exiting", db->filename); + } + dvo_image_create (db, GetZeroPoint()); + + return TRUE; +} + +int save_image_db (FITS_DB *db, Image *image, int Nimage) { + + /* add the new image and save */ + dvo_image_addrows (db, image, Nimage); + SetProtect (TRUE); + dvo_image_update (db, VERBOSE); + SetProtect (FALSE); + dvo_image_unlock (db); /* unlock? */ + + return TRUE; +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images_region.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images_region.c (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images_region.c (revision 37477) @@ -1,7 +1,5 @@ # include "fakeastro.h" -Image *fakeastro_images_region (Image *image, int *nimage, int *NIMAGE, Image *refImage, int NrefImage, SkyTable *skyTableInput, SkyTable *skyTableOutput, SkyRegion *innerRegion) { - - int Nimage = *nimage; +int fakeastro_images_region (ImageInfo *imageInfo, Image *refImage, int NrefImage, SkyTable *skyTableInput, SkyTable *skyTableOutput, SkyRegion *innerRegion) { // extend the outer region by 2.0 degrees @@ -14,5 +12,5 @@ Catalog *catalog = load_fake_stars (skyListInput, &Ncatalog); - int i; + int i, j; for (i = 0; i < NrefImage; i++) { @@ -32,8 +30,17 @@ Image *fakeImage = make_fake_images (&refImage[i], &NfakeImage); - int j; + // save the new fake images with their true image parameters (not yet fitted to the data) + if (imageInfo->NtrueImage + NfakeImage >= imageInfo->NTRUEIMAGE) { + imageInfo->NTRUEIMAGE += 1000 + NfakeImage; + REALLOCATE (imageInfo->trueImage, Image, imageInfo->NTRUEIMAGE); + } + for (j = 0; j < NfakeImage; j++) { + memcpy (&imageInfo->trueImage[imageInfo->NtrueImage], &fakeImage[j], sizeof(Image)); + imageInfo->NtrueImage ++; + } + for (j = 0; j < NfakeImage; j++) { - // we only want to make fake stars for the fake chips + // we only want to make fake stars for the fake chips (not the PHU entries) if (strcmp(&fakeImage[j].coords.ctype[4], "-WRP")) continue; @@ -41,4 +48,6 @@ Stars *fakeStars = make_fake_stars (catalog, Ncatalog, skyListInput, &fakeImage[j], &NfakeStars); + fit_fake_stars (fakeStars, NfakeStars, &fakeImage[j]); + save_fake_stars (skyTableOutput, &fakeImage[j], fakeStars, NfakeStars); @@ -46,12 +55,12 @@ } - if (Nimage + NfakeImage >= *NIMAGE) { - *NIMAGE += 1000 + NfakeImage; - REALLOCATE (image, Image, *NIMAGE); + // append the new fake images to the end of the full set: + if (imageInfo->NfakeImage + NfakeImage >= imageInfo->NFAKEIMAGE) { + imageInfo->NFAKEIMAGE += 1000 + NfakeImage; + REALLOCATE (imageInfo->fakeImage, Image, imageInfo->NFAKEIMAGE); } - for (j = 0; j < NfakeImage; j++) { - memcpy (&image[Nimage], &fakeImage[j], sizeof(Image)); - Nimage ++; + memcpy (&imageInfo->fakeImage[imageInfo->NfakeImage], &fakeImage[j], sizeof(Image)); + imageInfo->NfakeImage ++; } @@ -61,11 +70,9 @@ SkyListFree (skyListInput); - for (i = 0; i < Ncatalog; i++) { dvo_catalog_free (&catalog[i]); } - *nimage = Nimage; - return (image); + return TRUE; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c (revision 37477) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c (revision 37477) @@ -0,0 +1,100 @@ +# include "fakeastro.h" + +int fit_fake_stars (Stars *stars, int Nstars, Image *image) { + + // XXX I should set this based on Nstars + int order_use = 6; + + // I have Nstars with "reference" positions stars[].Rref,Dref and + // "observed" chip coordinates stars[].measure.Xccd,Yccd. + + // fit Xccd,Yccd to Rref,Dref + + int CHIPMAP = (image[0].coords.Npolyterms == -1); + if (CHIPMAP) image[0].coords.Npolyterms = -1; + + CoordFit *fit = fit_init (image[0].coords.Npolyterms); + + double *L, *M; + ALLOCATE (L, double, Nstars); + ALLOCATE (M, double, Nstars); + + // generate the fit matches + int i; + for (i = 0; i < Nstars; i++) { + RD_to_XY (&L[i], &M[i], stars[i].Rref, stars[i].Dref, image[0].coords.mosaic); + fit_add (fit, stars[i].measure.Xccd, stars[i].measure.Yccd, L[i], M[i], 0.05); + } + + // measure the fit, update the coords & object coordinates + // XXX do something more drastic on failure? + if (!fit_eval (fit)) { + fprintf (stderr, "failed to fit new model\n"); + image[0].flags |= ID_IMAGE_ASTROM_FAIL; + if (fit) fit_free (fit); + free (L); + free (M); + return FALSE; + } + + if (!fit_apply_coords (fit, &image[0].coords, FALSE)) { + fprintf (stderr, "failed to fit new model\n"); + image[0].flags |= ID_IMAGE_ASTROM_FAIL; + if (fit) fit_free (fit); + free (L); + free (M); + return FALSE; + } + + if (CHIPMAP) { + + // apply fit to get the fitted X,Y coordinates. we need these to fit the residual map below + float *dX, *dY, *Xref, *Yref; + ALLOCATE (dX, float, Nstars); + ALLOCATE (dY, float, Nstars); + ALLOCATE (Xref, float, Nstars); + ALLOCATE (Yref, float, Nstars); + + for (i = 0; i < Nstars; i++) { + // we have not yet fitted the map, so Npolyterms needs to be 1 here: + double tmpX, tmpY; + LM_to_XY (&tmpX, &tmpY, L[i], M[i], &image[0].coords); + + // fit the linear terms as above + // calculate dX (raw.X - ref.X) and dY + Xref[i] = tmpX; + Yref[i] = tmpY; + dX[i] = Xref[i] - stars[i].measure.Xccd; + dY[i] = Yref[i] - stars[i].measure.Yccd; + } + + // in coordsops.c:XY_to_LM, the map is defined to carry dX,dY so that: + // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY) + + if (image[0].coords.offsetMap == NULL) { + // allocate a new table and assign to this image + // need to allocate table? + AstromOffsetTable *table = get_astrom_table (); + AstromOffsetTableNewMap(table, order_use, image); + } + AstromOffsetMapFit (image[0].coords.offsetMap, Xref, Yref, dX, Nstars, TRUE); + AstromOffsetMapFit (image[0].coords.offsetMap, Xref, Yref, dY, Nstars, FALSE); + + image[0].coords.Npolyterms = -1; + + free (dX); + free (dY); + free (Xref); + free (Yref); + } + + for (i = 0; i < Nstars; i++) { + // if we have fitted the map above, Npolyterms needs to be -1 here: + XY_to_RD (&stars[i].average.R, &stars[i].average.D, stars[i].measure.Xccd, stars[i].measure.Yccd, &image[0].coords); + } + + free (L); + free (M); + + return TRUE; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fitpoly.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fitpoly.c (revision 37477) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fitpoly.c (revision 37477) @@ -0,0 +1,360 @@ +# include "fakeastro.h" + +/* these functions support simultaneous 2D fits to + x2 = \sum a_i,j x1^i y1^j + y2 = \sum b_i,j x1^i y1^j + + the order of the fit (largest coefficient) is fixed to a single + value for both x1,y1 terms and for x2,y2 fits + + the code is currently confusing because we limit to i+j <= order. + this could be cleaner if we used masks and allowed i <= order, j <= order +*/ + +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) { + + CoordFit *fit; + + ALLOCATE (fit, CoordFit, 1); + + fit[0].Npts = 0; + fit[0].Norder = order; + fit[0].Nterms = order + 1; + fit[0].Nsums = 2*order + 1; + fit[0].Nelems = SQ(order + 1); + + /* summing arrays for fit solution */ + + // xsum[i][j] holds \sum (x2 wt x1^i y1^j) + // ysum[i][j] holds \sum (y2 wt x1^i y1^j) + 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) + 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 +void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) { + + int ix, iy; + double xterm, yterm, term; + + xterm = 1; + for (ix = 0; ix < fit[0].Nsums; ix++) { + yterm = 1; + for (iy = 0; iy < fit[0].Nsums; iy++) { + term = xterm*yterm*wt; + fit[0].sum[ix][iy] += term; + if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) { + fit[0].xsum[ix][iy] += x2*term; + fit[0].ysum[ix][iy] += y2*term; + } + yterm *= y1; + } + xterm *= x1; + } + fit[0].Npts ++; +} + +/* convert the xsum,ysum,sum terms into vector,matrix and solve */ +int fit_eval (CoordFit *fit) { + + int i, j, ix, iy, jx, jy; + double **matrix, **vector; + + if (fit[0].Npts == 0) { + fprintf (stderr, "warning: no valid pts\n"); + return (FALSE); + } + + // 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 */ + for (i = 0; i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + 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] = fit[0].sum[ix+jx][iy+jy]; + } + + // mask the terms not represented by the Coords terms + if (ix + iy > fit[0].Norder) { + for (j = 0; j < fit[0].Nelems; j++) { + matrix[i][j] = 0.0; + } + vector[i][0] = 0.0; + vector[i][1] = 0.0; + matrix[i][i] = 1.0; + } + } + + for (i = 0; FALSE && (i < fit[0].Nelems); i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + ix, iy, vector[i][0], ix, iy, vector[i][1]); + } + + if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) { + array_free (matrix, fit[0].Nelems); + array_free (vector, fit[0].Nelems); + return (FALSE); + } + + for (i = 0; FALSE && i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + 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); + return (TRUE); +} + +void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) { + + int ix, iy; + double xterm, yterm, term; + double x, y; + + x = 0.0; + y = 0.0; + + xterm = 1; + for (ix = 0; ix < fit[0].Nterms; ix++) { + yterm = 1; + for (iy = 0; iy < fit[0].Nterms; iy++) { + term = xterm*yterm; + x += fit[0].xfit[ix][iy]*term; + y += fit[0].yfit[ix][iy]*term; + yterm *= y1; + } + xterm *= x1; + } + *x2 = x; + *y2 = y; +} + +// 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]; + } + } + 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); +} + +/* linear portion of fit : NORDER is 1 */ +/* this should only apply to the polynomial, not the projection terms */ +/* compare with psastro supporting code */ +// this code will work for linear (Npolyterm == 0 or 1) and linear + map fits +int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef) { + + double R1, R2; + CoordFit *modfit; + +# if (0) + if (keepRef) { + // adjust crpix1,2 as needed: + // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... + // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... + + double Xo, Yo; + if (!CoordsGetCenter (fit, 0.001, &Xo, &Yo)) { + fprintf (stderr, "failed to modify model\n"); + return (FALSE); + } + coords[0].crpix1 = Xo; + coords[0].crpix2 = Yo; + /* we do not modify crval1,2: these are kept at the default values */ + + // resulting fit should have zero constant terms + modfit = CoordsSetCenter (fit, Xo, Yo); + } else { +# endif + // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... + // P = cd1*x, Q = cd2*y + // L = pc1_1*P + pc1_2*Q + ... + + /* modify crval1,2: these are kept at the default values */ + + coords[0].crpix1 = 0.0; + coords[0].crpix2 = 0.0; + modfit = fit; +// } + + // set cdelt1, cdelt2 + coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]); + coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]); + R1 = 1 / coords[0].cdelt1; + R2 = 1 / coords[0].cdelt2; + + // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0) + coords[0].pc1_1 = modfit[0].xfit[1][0] * R1; + coords[0].pc1_2 = modfit[0].xfit[0][1] * R2; + coords[0].pc2_1 = modfit[0].yfit[1][0] * R1; + coords[0].pc2_2 = modfit[0].yfit[0][1] * R2; + + // set the polyterm elements + if (coords->Npolyterms > 1) { + coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1; + coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2; + coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2; + + coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1; + coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2; + coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2; + } + + // I need to validate Norder + if (coords->Npolyterms > 2) { + coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1; + coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2; + coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2; + coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2; + + coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1; + coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2; + coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2; + coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2; + } + + if (keepRef) { + fit_free (modfit); + } else { + fit_apply (fit, &coords[0].crval1, &coords[0].crval2, coords[0].crpix1, coords[0].crpix2); + } + /* keep the order and type from initial values */ + + return (TRUE); +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 37476) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 37477) @@ -37,5 +37,5 @@ float Mtime = 2.5*log10(image->exptime); - // XXX fix this!! + // XXX hard-wired plateScale for now? double plateScale = 0.257; @@ -70,4 +70,8 @@ double Rtru = average[i].R; double Dtru = average[i].D; + + // XXX what is the epoch for the src catalog? J2000? + stars[Nstars].Rref = Rtru; // add in any deviations we want here + stars[Nstars].Dref = Dtru; // observed position is scattered from true position by: