Index: /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/average.d =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/average.d (revision 12205) +++ /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/average.d (revision 12206) @@ -46,2 +46,4 @@ # the DVO object IDs are generated internally and are not equivalent to the PSPS object IDs # probably need to add position chisq + +# XXX include the number of measurements used to determine the positional information? Index: /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/measure.d =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/measure.d (revision 12205) +++ /branches/dvo-mods-2007-02/Ohana/src/libautocode/def/measure.d (revision 12206) @@ -47,4 +47,7 @@ FIELD imageID_lo, IMAGE_ID_LO, unsigned int, reference to image -# we may need 2 extra bytes for padding purposes... -FIELD dummy, dummy, char[6], padding +# note that with airmass = 1.0 / cos(90 - alt), we have full alt/az representation +FIELD az, AZ, float, telescope azimuth + +# we need extra bytes for padding purposes... +FIELD dummy, dummy, char[2], padding Index: /branches/dvo-mods-2007-02/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/libdvo/include/dvo.h (revision 12205) +++ /branches/dvo-mods-2007-02/Ohana/src/libdvo/include/dvo.h (revision 12206) @@ -39,5 +39,5 @@ PROJ_GLS, // pseudocyl PROJ_PAR, // pseudocyl -} OhanaProjections; +} OhanaProjection; typedef enum { @@ -46,5 +46,5 @@ PROJ_MODE_ZENITHAL, PROJ_MODE_PSEUDOCYL, -} OhanaProjectionModes; +} OhanaProjectionMode; /* RegImage.flag values */ @@ -226,5 +226,5 @@ OhanaProjection GetProjection (char *ctype); int SetProjection (char *ctype, OhanaProjection proj); -OhanaProjectionModes GetProjectionMode (OhanaProjections proj); +OhanaProjectionMode GetProjectionMode (OhanaProjection proj); char *libdvo_version (); Index: /branches/dvo-mods-2007-02/Ohana/src/libdvo/src/coordops.c =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/libdvo/src/coordops.c (revision 12205) +++ /branches/dvo-mods-2007-02/Ohana/src/libdvo/src/coordops.c (revision 12206) @@ -7,16 +7,35 @@ } -int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) { - - int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl; - char *type; - double L, M, X, Y, T, Z, Z2; - double R, sphi, cphi, stht, ctht; +int XY_to_LM (double *L, double *M, double x, double y, Coords *coords) { + + double X, Y; + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + + *L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + *M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + /** extra polynomial terms **/ + if (coords[0].Npolyterms > 1) { + *L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; + *M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; + } + if (coords[0].Npolyterms > 2) { + *L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; + *M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; + } + + return (TRUE); +} + +int LM_to_RD (double *ra, double *dec, double L, double M, Coords *coords) { + + OhanaProjection proj; + OhanaProjectionMode mode; + double R, T, Z, Z2, sphi, cphi, stht, ctht; double alpha, delta, salp, calp, sdel, sdp, cdp; - - *ra = 0; - *dec = 0; - stht = ctht = 1; - + proj = GetProjection (coords[0].ctype); mode = GetProjectionMode (proj); @@ -24,20 +43,5 @@ if (proj == PROJ_MODE_NONE) return (FALSE); - /** convert pixel coordinates to cartesian system **/ - X = coords[0].cdelt1*(x - coords[0].crpix1); - Y = coords[0].cdelt2*(y - coords[0].crpix2); - - L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); - M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); - - /** extra polynomial terms **/ - if (coords[0].Npolyterms > 1) { - L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; - M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; - } - if (coords[0].Npolyterms > 2) { - L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; - M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; - } + stht = ctht = 1; /** Locally Cartesian Projections **/ @@ -77,5 +81,5 @@ break; case PROJ_STG: - stht = (4 - RAD_DEG*R) / (4 + RAD_DEC*R); + stht = (4 - RAD_DEG*R) / (4 + RAD_DEG*R); ctht = sqrt (1 - stht*stht); break; @@ -94,4 +98,6 @@ ctht = sqrt (1 - stht*stht); break; + default: + return (FALSE); } sdp = sin(RAD_DEG*coords[0].crval2); @@ -115,5 +121,5 @@ /**** Other Conventional Projections ****/ - if (mode == PROJ_MODE_PSEUDOCLY) { + if (mode == PROJ_MODE_PSEUDOCYL) { switch (proj) { case PROJ_AIT: @@ -135,4 +141,6 @@ delta = 3 * DEG_RAD * asin (M/180.0); break; + default: + return (FALSE); } *ra = alpha + coords[0].crval1; @@ -148,20 +156,15 @@ } -int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) { - - char *type; - int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl; +int RD_to_LM (double *L, double *M, double ra, double dec, Coords *coords) { + double phi, theta; - double determ; - double X, Y, L, M, Lo, Mo, dL, dM; + double Lo, Mo; double sphi, cphi, stht; double salp, calp, sdel, cdel, sdp, cdp; double P, A, Rc; - - status = TRUE; - *x = 0; - *y = 0; - type = &coords[0].ctype[4]; - L = M = 0; + OhanaProjection proj; + OhanaProjectionMode mode; + + *L = *M = 0; proj = GetProjection (coords[0].ctype); @@ -175,10 +178,11 @@ if (mosaic == NULL) return (FALSE); RD_to_XY (&Lo, &Mo, ra, dec, mosaic); - L = (Lo - coords[0].crval1); - M = (Mo - coords[0].crval2); - } else { - L = (ra - coords[0].crval1); - M = (dec - coords[0].crval2); - } + *L = (Lo - coords[0].crval1); + *M = (Mo - coords[0].crval2); + return (TRUE); + } + *L = (ra - coords[0].crval1); + *M = (dec - coords[0].crval2); + return (TRUE); } @@ -195,24 +199,26 @@ sphi = cdel*salp; /* = cos(theta)*sin(phi) */ cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ - if (stht < 0) status = FALSE; + if (stht < 0) return (FALSE); switch (proj) { case PROJ_TAN: case PROJ_DIS: - L = DEG_RAD * sphi / stht; - M = -DEG_RAD * cphi / stht; - break; + *L = DEG_RAD * sphi / stht; + *M = -DEG_RAD * cphi / stht; + return (TRUE); case PROJ_SIN: - L = DEG_RAD * sphi; - M = -DEG_RAD * cphi; - break; + *L = DEG_RAD * sphi; + *M = -DEG_RAD * cphi; + return (TRUE); case PROJ_ZEA: case PROJ_ZPL: Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht); - L = Rc * sphi; - M = -Rc * cphi; - status = TRUE; - break; - } + *L = Rc * sphi; + *M = -Rc * cphi; + return (TRUE); + default: + return (FALSE); + } + return (FALSE); } @@ -224,27 +230,41 @@ theta = RAD_DEG*(dec - coords[0].crval2); P = 1.0 + cos (theta) * cos (0.5*phi); - if (P != 0.0) { - A = DEG_RAD * sqrt (2.0 / P); - L = 2.0 * A * cos (theta) * sin (0.5*phi); - M = A * sin (theta); - } else { - L = 0.0; - M = 0.0; - } - break; + if (P == 0.0) { + *L = 0.0; + *M = 0.0; + return (TRUE); + } + A = DEG_RAD * sqrt (2.0 / P); + *L = 2.0 * A * cos (theta) * sin (0.5*phi); + *M = A * sin (theta); + return (TRUE); case PROJ_GLS: phi = ra - coords[0].crval1; theta = dec - coords[0].crval2; - L = phi * cos(RAD_DEG * theta); - M = theta; - break; + *L = phi * cos(RAD_DEG * theta); + *M = theta; + return (TRUE); case PROJ_PAR: phi = ra - coords[0].crval1; theta = dec - coords[0].crval2; - L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1); - M = 180.0 * sin (RAD_DEG*theta/3.0); - break; - } - } + *L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1); + *M = 180.0 * sin (RAD_DEG*theta/3.0); + return (TRUE); + default: + return (FALSE); + } + return (FALSE); + } + return (FALSE); +} + +int LM_to_XY (double *x, double *y, double L, double M, Coords *coords) { + + int i; + double determ; + double X, Y, Lo, Mo, dL, dM; + + *x = 0; + *y = 0; /* convert L,M to X,Y */ @@ -278,4 +298,28 @@ *y = Y / coords[0].cdelt2 + coords[0].crpix2; + return (TRUE); +} + +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) { + + double L, M; + int status; + + status = XY_to_LM (&L, &M, x, y, coords); + if (!status) return FALSE; + + status = LM_to_RD (ra, dec, L, M, coords); + return (status); +} + +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) { + + double L, M; + int status; + + status = RD_to_LM (&L, &M, ra, dec, coords); + if (!status) return FALSE; + + status = LM_to_XY (x, y, L, M, coords); return (status); } @@ -601,5 +645,5 @@ WRP is equiv to PLY, with implied mosaic */ -int GetProjection (char *ctype) { +OhanaProjection GetProjection (char *ctype) { if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA; if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL; @@ -620,24 +664,25 @@ } -int SetProjection (char *ctype, int proj) { +int SetProjection (char *ctype, OhanaProjection proj) { switch (proj) { - case PROJ_ZEA: strcpy(&ctype[4], "-ZEA") return TRUE; - case PROJ_ZPL: strcpy(&ctype[4], "-ZPL") return TRUE; - case PROJ_ARC: strcpy(&ctype[4], "-ARC") return TRUE; - case PROJ_STG: strcpy(&ctype[4], "-STG") return TRUE; - case PROJ_SIN: strcpy(&ctype[4], "-SIN") return TRUE; - case PROJ_TAN: strcpy(&ctype[4], "-TAN") return TRUE; - case PROJ_DIS: strcpy(&ctype[4], "-DIS") return TRUE; - case PROJ_LIN: strcpy(&ctype[4], "-LIN") return TRUE; - case PROJ_PLY: strcpy(&ctype[4], "-PLY") return TRUE; - case PROJ_WRP: strcpy(&ctype[4], "-WRP") return TRUE; - case PROJ_AIT: strcpy(&ctype[4], "-AIT") return TRUE; - case PROJ_GLS: strcpy(&ctype[4], "-GLS") return TRUE; - case PROJ_PAR: strcpy(&ctype[4], "-PAR") return TRUE; + case PROJ_ZEA: strcpy(&ctype[4], "-ZEA"); return TRUE; + case PROJ_ZPL: strcpy(&ctype[4], "-ZPL"); return TRUE; + case PROJ_ARC: strcpy(&ctype[4], "-ARC"); return TRUE; + case PROJ_STG: strcpy(&ctype[4], "-STG"); return TRUE; + case PROJ_SIN: strcpy(&ctype[4], "-SIN"); return TRUE; + case PROJ_TAN: strcpy(&ctype[4], "-TAN"); return TRUE; + case PROJ_DIS: strcpy(&ctype[4], "-DIS"); return TRUE; + case PROJ_LIN: strcpy(&ctype[4], "-LIN"); return TRUE; + case PROJ_PLY: strcpy(&ctype[4], "-PLY"); return TRUE; + case PROJ_WRP: strcpy(&ctype[4], "-WRP"); return TRUE; + case PROJ_AIT: strcpy(&ctype[4], "-AIT"); return TRUE; + case PROJ_GLS: strcpy(&ctype[4], "-GLS"); return TRUE; + case PROJ_PAR: strcpy(&ctype[4], "-PAR"); return TRUE; + case PROJ_NONE: return FALSE; } return FALSE; } -int GetProjectionMode (int proj) { +OhanaProjectionMode GetProjectionMode (OhanaProjection proj) { switch (proj) { case PROJ_ZEA: @@ -656,5 +701,6 @@ case PROJ_GLS: case PROJ_PAR: - return PROJ_MODE_PSEUDOCLY; + return PROJ_MODE_PSEUDOCYL; + default: PROJ_MODE_NONE; } return PROJ_MODE_NONE; Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h (revision 12205) +++ /branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h (revision 12206) @@ -3,4 +3,17 @@ # include # include + +// structure to hold coordinate fitting terms +typedef struct { + int Npts; + int Nterms; + int Npower; + int Norder; + double **xsum; + double **ysum; + + double **vector; + double **matrix; +} CoordFit; typedef struct { @@ -12,4 +25,7 @@ double p, dp; + + double chisq; + int Nfit; } PMfit;