Index: /trunk/psastro/src/psLibUtils.c =================================================================== --- /trunk/psastro/src/psLibUtils.c (revision 5562) +++ /trunk/psastro/src/psLibUtils.c (revision 5562) @@ -0,0 +1,172 @@ +# include // for strncasecmp +# include +# include "psLibUtils.h" + +// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) +// 2 = default, -v = 3, -vv = 4, -vvv = 5 +psS32 psLogArguments (int *argc, char **argv) { + + int N, level; + + // default log level is 2 + level = 2; + + // set in order, so that -vvv overrides -vv overrides -v + if ((N = psArgumentGet (*argc, argv, "-v"))) { + psArgumentRemove (N, argc, argv); + level = 3; + } + if ((N = psArgumentGet (*argc, argv, "-vv"))) { + psArgumentRemove (N, argc, argv); + level = 4; + } + if ((N = psArgumentGet (*argc, argv, "-vvv"))) { + psArgumentRemove (N, argc, argv); + level = 5; + } + + if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { + if (*argc < N + 2) { + psAbort ("psLogArguments", "USAGE: -logfmt (format)"); + } + psArgumentRemove (N, argc, argv); + psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid + psArgumentRemove (N, argc, argv); + } + + // set the level, return the level + psLogSetLevel (level); + return (level); +} + +// set trace levels by facility +psS32 psTraceArguments (int *argc, char **argv) { + + int N; + + // argument format is: -trace (facil) (level) + while ((N = psArgumentGet (*argc, argv, "-trace"))) { + if (*argc < N + 3) { + psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); + } + psArgumentRemove (N, argc, argv); + psTraceSetLevel (argv[N], atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { + psTracePrintLevels (); + exit (2); + } + return (TRUE); +} + +psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) { + + psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type); + + for (int i = 0; i < input->nX; i++) { + for (int j = 0; j < input->nY; j++) { + output->mask[i][j] = input->mask[i][j]; + output->coeff[i][j] = input->coeff[i][j]; + output->coeffErr[i][j] = input->coeffErr[i][j]; + } + } + return (output); +} + +psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) { + + psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type); + + for (int i = 0; i < input->nX; i++) { + for (int j = 0; j < input->nY; j++) { + for (int k = 0; k < input->nZ; k++) { + for (int m = 0; m < input->nT; m++) { + output->mask[i][j][k][m] = input->mask[i][j][k][m]; + output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; + output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; + } + } + } + } + return (output); +} + +psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input) { + + psPlaneDistort *output = psPlaneDistortAlloc (input->x->nX, input->x->nY, input->x->nZ, input->x->nT); + + for (int i = 0; i < input->x->nX; i++) { + for (int j = 0; j < input->x->nY; j++) { + for (int k = 0; k < input->x->nZ; k++) { + for (int m = 0; m < input->x->nT; m++) { + // x-terms + output->x->mask[i][j][k][m] = input->x->mask[i][j][k][m]; + output->x->coeff[i][j][k][m] = input->x->coeff[i][j][k][m]; + output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m]; + // y-terms + output->y->mask[i][j][k][m] = input->y->mask[i][j][k][m]; + output->y->coeff[i][j][k][m] = input->y->coeff[i][j][k][m]; + output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m]; + } + } + } + } + return (output); +} + +psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) { + + psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY); + + for (int i = 0; i < input->x->nX; i++) { + for (int j = 0; j < input->x->nY; j++) { + // x-terms + output->x->mask[i][j] = input->x->mask[i][j]; + output->x->coeff[i][j] = input->x->coeff[i][j]; + output->x->coeffErr[i][j] = input->x->coeffErr[i][j]; + // y-terms + output->y->mask[i][j] = input->y->mask[i][j]; + output->y->coeff[i][j] = input->y->coeff[i][j]; + output->y->coeffErr[i][j] = input->y->coeffErr[i][j]; + } + } + return (output); +} + +psProjection *psProjectionCopy (psProjection *input) { + + psProjection *output = psProjectionAlloc (input->R, input->D, input->Xs, input->Ys, input->type); + return (output); +} + +// very crude distortion inversion: assumes 0 order in z and t, linear in x and y: +psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) { + PS_ASSERT_PTR_NON_NULL(distort, 0); + PS_ASSERT_PTR_NON_NULL(distort->x, 0); + PS_ASSERT_PTR_NON_NULL(distort->y, 0); + + psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0); + + /* simple matrix inversion code */ + + psF64 r11 = distort->x->coeff[1][0][0][0]; + psF64 r12 = distort->x->coeff[0][1][0][0]; + psF64 r21 = distort->y->coeff[1][0][0][0]; + psF64 r22 = distort->y->coeff[0][1][0][0]; + psF64 xo = distort->x->coeff[0][0][0][0]; + psF64 yo = distort->y->coeff[0][0][0][0]; + + psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant + + out->x->coeff[1][0][0][0] = +invDet * r22; + out->x->coeff[0][1][0][0] = -invDet * r12; + out->y->coeff[1][0][0][0] = -invDet * r21; + out->y->coeff[0][1][0][0] = +invDet * r11; + + out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo); + out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo); + + return(out); +} Index: /trunk/psastro/src/psLibUtils.h =================================================================== --- /trunk/psastro/src/psLibUtils.h (revision 5562) +++ /trunk/psastro/src/psLibUtils.h (revision 5562) @@ -0,0 +1,24 @@ +# ifndef PS_LIB_UTILS +# define PS_LIB_UTILS + +// structure to carry a dynamic string +typedef struct { + int NLINE; + int Nline; + char *line; +} psLine; + +# define psMemCopy(A)(psMemIncrRefCounter((A))) +# define DEG_RAD 57.295779513082322 +# define RAD_DEG 0.017453292519943 + +// psLib extra utilities +psS32 psLogArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) +psS32 psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) + +// psLine functions -- keep out for now? +psLine *psLineAlloc (int Nline); +bool psLineInit (psLine *line); +bool psLineAdd (psLine *line, char *format, ...); + +# endif Index: /trunk/psastro/src/psastro-mktest.c =================================================================== --- /trunk/psastro/src/psastro-mktest.c (revision 5562) +++ /trunk/psastro/src/psastro-mktest.c (revision 5562) @@ -0,0 +1,264 @@ +# include "psastro.h" + +bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); +bool testWriteRef (char *filename, psArray *sources); +bool testWriteRaw (char *filename, psArray *sources); + +int main (int argc, char **argv) { + + bool status; + + // load configuration information + psMetadata *config = testArguments (&argc, argv); + + // fill in the header with astrometry info as we go... + psMetadata *header = psMetadataAlloc(); + + // make fake star set + psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); + + int NX = psMetadataLookupS32 (&status, config, "NX"); + int NY = psMetadataLookupS32 (&status, config, "NY"); + int NSTARS = psMetadataLookupS32 (&status, config, "NSTARS"); + + psMetadataAdd (header, PS_LIST_TAIL, "SIMPLE", PS_DATA_BOOL | PS_META_REPLACE, "foo", true); + psMetadataAdd (header, PS_LIST_TAIL, "BITPIX", PS_DATA_S32 | PS_META_REPLACE, "foo", 16); + psMetadataAdd (header, PS_LIST_TAIL, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "foo", 2); + psMetadataAdd (header, PS_LIST_TAIL, "NAXIS1", PS_DATA_S32 | PS_META_REPLACE, "foo", NX); + psMetadataAdd (header, PS_LIST_TAIL, "NAXIS2", PS_DATA_S32 | PS_META_REPLACE, "foo", NY); + + psArray *rawstars = psArrayAlloc (NSTARS); + for (int i = 0; i < NSTARS; i++) { + pmAstromObj *star = pmAstromObjAlloc (); + star->pix.x = NX * psRandomUniform (rnd); + star->pix.y = NY * psRandomUniform (rnd); + star->Mag = 10.0; + star->dMag = 0.05; + rawstars->data[i] = star; + } + + pmFPA *fpa = pmFPAAlloc (NULL, NULL); + pmChip *chip = pmChipAlloc (fpa); + pmCell *cell = pmCellAlloc (chip); + pmReadout *readout = pmReadoutAlloc (cell); + + { // define readout -> cell transformation + readout->col0 = psMetadataLookupF32 (&status, config, "READOUT.COL.0"); + readout->row0 = psMetadataLookupF32 (&status, config, "READOUT.ROW.0"); + readout->colBins = psMetadataLookupF32 (&status, config, "READOUT.COL.BIN"); + readout->rowBins = psMetadataLookupF32 (&status, config, "READOUT.ROW.BIN"); + } + + { // define cell -> chip transformation + cell->toChip = psPlaneTransformAlloc (1, 1); + cell->toChip->x->coeff[0][0] = psMetadataLookupF32 (&status, config, "CELL.COL.0"); + cell->toChip->y->coeff[0][0] = psMetadataLookupF32 (&status, config, "CELL.ROW.0"); + + cell->toChip->x->coeff[1][0] = 1; + cell->toChip->x->coeff[0][1] = 0; + cell->toChip->x->coeff[1][1] = 0; + + cell->toChip->y->coeff[1][0] = 0; + cell->toChip->y->coeff[0][1] = 1; + cell->toChip->y->coeff[1][1] = 0; + } + + { // define chip -> fpa transformation + chip->toFPA = psPlaneTransformAlloc (1, 1); + chip->toFPA->x->coeff[0][0] = psMetadataLookupF32 (&status, config, "CHIP.COL.0"); + chip->toFPA->y->coeff[0][0] = psMetadataLookupF32 (&status, config, "CHIP.ROW.0"); + + double theta = psMetadataLookupF32 (&status, config, "CHIP.THETA"); + double px = psMetadataLookupF32 (&status, config, "CHIP.PARITY.X"); + double py = psMetadataLookupF32 (&status, config, "CHIP.PARITY.Y"); + + double cs = cos (theta*RAD_DEG); + double sn = sin (theta*RAD_DEG); + + chip->toFPA->x->coeff[1][0] = +px*cs; + chip->toFPA->x->coeff[0][1] = +px*sn; + chip->toFPA->x->coeff[1][1] = 0; + + chip->toFPA->y->coeff[1][0] = -py*sn; + chip->toFPA->y->coeff[0][1] = +py*cs; + chip->toFPA->y->coeff[1][1] = 0; + + psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]); + + chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); + } + + { // define fpa -> tpa transformation + fpa->toTangentPlane = psPlaneDistortAlloc (1, 1, 0, 0); + fpa->toTangentPlane->x->coeff[0][0][0][0] = psMetadataLookupF32 (&status, config, "FPA.COL.0"); + fpa->toTangentPlane->y->coeff[0][0][0][0] = psMetadataLookupF32 (&status, config, "FPA.ROW.0"); + + double theta = psMetadataLookupF32 (&status, config, "FPA.THETA"); + // double px = psMetadataLookupF32 (&status, config, "FPA.PARITY.X"); + // double py = psMetadataLookupF32 (&status, config, "FPA.PARITY.Y"); + + double cs = cos (theta*RAD_DEG); + double sn = sin (theta*RAD_DEG); + + fpa->toTangentPlane->x->coeff[1][0][0][0] = +cs; + fpa->toTangentPlane->x->coeff[0][1][0][0] = +sn; + fpa->toTangentPlane->x->coeff[1][1][0][0] = 0; + + fpa->toTangentPlane->y->coeff[1][0][0][0] = -sn; + fpa->toTangentPlane->y->coeff[0][1][0][0] = +cs; + fpa->toTangentPlane->y->coeff[1][1][0][0] = 0; + + fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane); + } + + { // define tpa -> sky projection, update output header + + double RA = psMetadataLookupF32 (&status, config, "RA"); + double DEC = psMetadataLookupF32 (&status, config, "DEC"); + double PS = psMetadataLookupF32 (&status, config, "PLATE.SCALE"); + // XXX EAM : PS units are inverted in psDeproject (fixed in eam branch) + PS /= 3600.0 * DEG_RAD; // convert from arcsec to radians + RA *= RAD_DEG; // convert from degrees to radians + DEC *= RAD_DEG; // convert from degrees to radians + + fpa->projection = psProjectionAlloc (RA, DEC, PS, PS, PS_PROJ_SIN); + + // find the readout coord for ref pixel: + psSphere tmpSky; + psPlane tmp1, tmp2; + tmpSky.r = RA; + tmpSky.d = DEC; + + psCoordSkyToTP (&tmp1, &tmpSky, fpa->projection); + psCoordTPToFPA (&tmp2, &tmp1, 0.0, 0.0, fpa); + psCoordFPAToChip (&tmp1, &tmp2, chip); + psCoordChipToCell_EAM (&tmp2, &tmp1, cell); + psCoordCellToReadout (&tmp1, &tmp2, readout); + + // write output header WCS info + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "foo", RA*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "foo", DEC*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN"); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.x); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.y); + } + + // project stars from readout up to sky + psastroProjectRawstars (rawstars, readout); + + // generate refstars from rawstars (for test comparison) + psArray *refstars = psArrayAlloc (rawstars->n); + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *oStar = pmAstromObjAlloc (); + pmAstromObj *iStar = rawstars->data[i]; + oStar->sky = iStar->sky; + oStar->Mag = iStar->Mag; + oStar->dMag = iStar->dMag; + refstars->data[i] = oStar; + } + + // project stars from sky down to readout + psastroProjectRefstars (refstars, readout); + + testWriteCMP (header, argv[2], rawstars); + testWriteRef (argv[3], rawstars); + testWriteRaw (argv[4], rawstars); + testWriteRaw (argv[5], refstars); + exit (0); +} + +// elixir-style pseudo FITS table (header + ascii list) +bool testWriteCMP (psMetadata *header, char *filename, psArray *sources) { + + int i; + psMetadataItem *mdi; + + psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "foo", sources->n); + + // create file, write-out header + unlink (filename); + psFits *fits = psFitsAlloc (filename); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + psFitsWriteHeader (header, fits); + psFree (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f\n", + star->pix.x, star->pix.y, star->Mag, (int)(star->dMag*1000), 1, 2.5, 11.0, 12.0, 2.0, 0.5, 15.0); + } + fclose (f); + return true; +} + +// elixir-style pseudo FITS table (header + ascii list) +bool testWriteRef (char *filename, psArray *sources) { + + int i; + + // re-open, add data to end of file + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%10.7f %10.7f %6.3f %6.3f\n", star->sky.r, star->sky.d+10*RAD_DEG/3600, 10.0, 0.05); + } + fclose (f); + return true; +} + +// elixir-style pseudo FITS table (header + ascii list) +bool testWriteRaw (char *filename, psArray *sources) { + + int i; + + // re-open, add data to end of file + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", + star->pix.x, star->pix.y, + star->cell.x, star->cell.y, + star->chip.x, star->chip.y, + star->FP.x, star->FP.y, + star->TP.x, star->TP.y, + star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, + star->Mag, star->dMag); + } + fclose (f); + return true; +} Index: /trunk/psastro/src/psastroIO.c =================================================================== --- /trunk/psastro/src/psastroIO.c (revision 5562) +++ /trunk/psastro/src/psastroIO.c (revision 5562) @@ -0,0 +1,125 @@ +# include "psastro.h" + +// crude function to load CMP data (header + ascii) +// XXX EAM : assumes fixed line, expects NSTARS entries +psArray *psastroReadCMP (psMetadata **header, char *filename) { + + bool status; + char line[80]; + double X, Y, Mag, dMag; + psMetadata *myHeader; + + psFits *fits = psFitsAlloc (filename); + myHeader = psFitsReadHeader (NULL, fits); + psFree (fits); + + // how many lines in the header? + // XXX EAM : is this calculation robust? + int nLines = myHeader->list->n; + int nBytes = nLines * 80; + if (nBytes % 2880) { + int nBlock = 1 + (int)(nBytes / 2880); + nBytes = nBlock * 2880; + } + + // re-open, seek to end of header + FILE *f = fopen (filename, "r"); + if (f == NULL) { + psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename); + return NULL; + } + fseek (f, nBytes, SEEK_SET); + + // prepare array to store data + int nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS"); + psArray *stars = psArrayAlloc (nStars); + stars->n = 0; + + // we have fixed bytes / line : use that info + for (int i = 0; i < nStars; i++) { + fread (line, 1, 66, f); + line[66] = 0; + sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag); + // fprintf (stderr, "%f %f %f %f\n", X, Y, Mag, dMag); + + pmAstromObj *obj = pmAstromObjAlloc (); + obj->pix.x = X; + obj->pix.y = Y; + obj->Mag = Mag; + obj->dMag = dMag; + psArrayAdd (stars, 100, obj); + } + fclose (f); + + *header = myHeader; + return (stars); +} + +// elixir-style pseudo FITS table (header + ascii list) +bool psastroWriteCMP (pmFPA *fpa, char *filename) { + + bool status; + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + for (int j = 0; j < chip->cells->n; j++) { + pmCell *cell = chip->cells->data[j]; + + // create file, write-out header + unlink (filename); + psFits *fits = psFitsAlloc (filename); + psFitsWriteHeader (cell->header, fits); + psFree (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("psastroWriteCMP", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + for (int k = 0; k < cell->readouts->n; k++) { + pmReadout *readout = cell->readouts->data[k]; + psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET"); + for (int m = 0; m < stars->n; m++) { + pmAstromObj *star = stars->data[m]; + fprintf (f, "%6.1f %6.1f\n", star->pix.x, star->pix.y); + } + } + fclose (f); + } + } + return true; +} + +// fake this by loading from a text file ?? +psArray *psastroLoadReference (psMetadata *header, psMetadata *config) { + + // which reference source? + // ra & dec, area from header + + double RA, DEC, Mag, dMag; + + FILE *f = fopen ("test.cat", "r"); + if (f == NULL) { + psLogMsg ("psastroLoadReference", 3, "can't open dummy reference file test.ref\n"); + return NULL; + } + + psArray *stars = psArrayAlloc (100); + stars->n = 0; + + // we have fixed bytes / line : use that info + while (fscanf (f, "%lf %lf %lf %lf", &RA, &DEC, &Mag, &dMag) != EOF) { + pmAstromObj *obj = pmAstromObjAlloc (); + obj->sky.r = RA; + obj->sky.d = DEC; + obj->Mag = Mag; + obj->dMag = dMag; + psArrayAdd (stars, 100, obj); + } + fclose (f); + return (stars); +} +