Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 3388) +++ /trunk/Ohana/src/addstar/Makefile (revision 3389) @@ -1,2 +1,4 @@ +#install: +# @echo "addstar under re-construction" default: addstar help: @@ -110,5 +112,5 @@ # utilities ################################################# -install: +# install: for i in $(INSTALL); do make $$i.install; done Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 3388) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 3389) @@ -62,4 +62,5 @@ PhotCode *thiscode; GSCRegion UserPatch; +Coords *MOSAIC; char DateKeyword[64]; Index: /trunk/Ohana/src/addstar/src/args.c =================================================================== --- /trunk/Ohana/src/addstar/src/args.c (revision 3388) +++ /trunk/Ohana/src/addstar/src/args.c (revision 3389) @@ -14,5 +14,5 @@ int args (int argc, char **argv) { - int i, N; + int N; /* check for help request */ @@ -61,11 +61,33 @@ /* provide a time for dataset */ TIMEREF = 0; - if ((i = get_argument (argc, argv, "-time"))) { - remove_argument (i, &argc, argv); - if (!str_to_time (argv[i], &TIMEREF)) { + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!str_to_time (argv[N], &TIMEREF)) { fprintf (stderr, "syntax error in time\n"); exit (1); } - remove_argument (i, &argc, argv); + remove_argument (N, &argc, argv); + } + /* provide a mosaic for distortion */ + MOSAIC = NULL; + if ((N = get_argument (argc, argv, "-mosaic"))) { + Header header; + ALLOCATE (MOSAIC, Coords, 1); + + remove_argument (N, &argc, argv); + if (!fits_read_header (argv[N], &header)) { + fprintf (stderr, "ERROR: can't read header for mosaic %s\n", argv[N]); + exit (1); + } + if (!GetCoords (MOSAIC, &header)) { + fprintf (stderr, "ERROR: no astrometric solution in header\n"); + exit (1); + } + if (strcmp(&MOSAIC[0].ctype[4], "-DIS")) { + fprintf (stderr, "ERROR: not a mosaic distortion header\n"); + exit (1); + } + remove_argument (N, &argc, argv); + fits_free_header (&header); } @@ -73,7 +95,7 @@ /* only add to existing objects */ ONLY_MATCH = FALSE; - if ((i = get_argument (argc, argv, "-only-match"))) { + if ((N = get_argument (argc, argv, "-only-match"))) { ONLY_MATCH = TRUE; - remove_argument (i, &argc, argv); + remove_argument (N, &argc, argv); } /* don't add missed pts to Missed table (image only) */ @@ -85,7 +107,7 @@ /* replace measurement, don't duplicate (ref/cat only) */ REPLACE = FALSE; - if ((i = get_argument (argc, argv, "-replace"))) { + if ((N = get_argument (argc, argv, "-replace"))) { REPLACE = TRUE; - remove_argument (i, &argc, argv); + remove_argument (N, &argc, argv); } /* only add image potion to image table */ @@ -133,5 +155,5 @@ fprintf (stderr, "USAGE: addstar -ref (filename)\n"); fprintf (stderr, " cur: %s ", argv[0]); - for (i = 0; i < argc; i++) { fprintf (stderr, "%s ", argv[i]); } + for (N = 0; N < argc; N++) { fprintf (stderr, "%s ", argv[N]); } fprintf (stderr, "\n"); exit (2); Index: /trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches.c (revision 3388) +++ /trunk/Ohana/src/addstar/src/find_matches.c (revision 3389) @@ -49,12 +49,22 @@ REALLOCATE (catalog[0].missing, Missing, NMISS); - /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ - /* reference for coords is this image */ - tcoords = image[0].coords; - tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - strcpy (tcoords.ctype, "RA---TAN"); - + /* project onto rectilinear grid with 1 arcsec pixels */ + /* we keep the original crpix1,2 and crref1,2 */ + /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + tcoords = MOSAIC[0]; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 1; + strcpy (tcoords.ctype, "RA---TAN"); + } else { + tcoords = image[0].coords; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + } + for (i = 0; i < Nstars; i++) { fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); @@ -134,5 +144,4 @@ /* calculate accurate per-star airmass */ secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude); - /** add measurements for this star **/ @@ -162,7 +171,7 @@ if (*Mval == NO_MAG) *Mval = Mrel; + /** this is a little tricky: the new measures are not in the measure sequence + this call can run into the next star if we have had two measures found at this location. **/ if (CALIBRATE) { - /** this is a little tricky: the new measures are not in the measure sequence - this call can run into the next star if we have had two measures found at this location. **/ AddToCalibration (&catalog[0].average[n], &catalog[0].measure[m], &catalog[0].measure[Nmeas], N); } @@ -202,4 +211,5 @@ /* add reference for undetected catalog stars */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) RegisterMosaic (MOSAIC); for (j = 0; (j < Nave) && !SKIP_MISSED; j++) { n = N2[j]; @@ -287,4 +297,5 @@ /** now add references from all previous non-detection observations of this spot on the sky */ for (j = 0; (j < Noverlap) && !SKIP_MISSED; j++) { + if (!FindMosaicForImage (overlap, Noverlap, j)) continue; if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue; if (catalog[0].average[Nave].Nn < 1) { Index: /trunk/Ohana/src/addstar/src/gimages.c =================================================================== --- /trunk/Ohana/src/addstar/src/gimages.c (revision 3388) +++ /trunk/Ohana/src/addstar/src/gimages.c (revision 3389) @@ -38,9 +38,21 @@ /* project onto rectilinear grid with 1 arcsec pixels */ /* we keep the original crpix1,2 and crref1,2 */ - tcoords = image[0].coords; - tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - strcpy (tcoords.ctype, "RA---TAN"); + /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + tcoords = MOSAIC[0]; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 1; + strcpy (tcoords.ctype, "RA---TAN"); + RegisterMosaic (MOSAIC); + /* register so image->sky conversions below have correct mosaic */ + } else { + tcoords = image[0].coords; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + } /* define original corners */ @@ -51,4 +63,5 @@ Xo[4] = 0; Yo[4] = 0; /* so we can make a loop easily */ for (j = 0; j < 5; j++) { + /* XY-to_RD is two-level if ctype == WRP */ XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords); RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords); @@ -77,5 +90,6 @@ } - NTIMAGE = 100; + /* kind of dangerous, see note below */ + NTIMAGE = MAX (MIN (Nimage, 1000), 1); ALLOCATE (timage, Image, NTIMAGE); @@ -87,5 +101,15 @@ for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) { ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image"); + if (ntimage == 0) Shutdown ("expected %d images, only found %d", Nimage, Ntimage); + BuildChipMatch (timage, ntimage); + for (i = 0; i < ntimage; i++) { + + /* if any of these images are WRP images, need to find matching DIS */ + if (!FindMosaicForImage (timage, ntimage, i)) { + fprintf (stderr, "!"); + continue; + } + addtolist = FALSE; /* define image corners */ @@ -116,9 +140,7 @@ } } + } - if (ntimage == 0) { - Shutdown ("expected %d images, only found %d", Nimage, Ntimage); - } - } + BuildChipMatch (pimage, npimage); if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage); @@ -129,2 +151,9 @@ } +/* this function tries to find overlapping images with the current image. + it loads the image table in blocks (up to 1000 at a time) and looks + for overlaps within the set. the search is not very robust. in + addition, it needs to look for mosaic DIS images matching any WRP images. + if the WRP image is loaded near the beginning of a block, then the corresponding DIS + may be in the previous block, and thus missed +*/ Index: /trunk/Ohana/src/addstar/src/gregion_image.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_image.c (revision 3388) +++ /trunk/Ohana/src/addstar/src/gregion_image.c (revision 3389) @@ -12,4 +12,9 @@ impatch.DEC[0] = +90; impatch.DEC[1] = -90; + + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + if (MOSAIC == NULL) Shutdown ("no mosaic for WRP image (use -mosaic)"); + RegisterMosaic (MOSAIC); + } for (x = 0; x <= 1; x++) { Index: /trunk/Ohana/src/addstar/src/gstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/gstars.c (revision 3388) +++ /trunk/Ohana/src/addstar/src/gstars.c (revision 3389) @@ -52,4 +52,11 @@ exit (1); } + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + if (MOSAIC == NULL) { + fprintf (stderr, "ERROR: no mosaic for WRP image (use -mosaic)\n"); + exit (1); + } + RegisterMosaic (MOSAIC); + } /* get ST (used for airmass calculation) */ Index: /trunk/Ohana/src/libohana/Makefile =================================================================== --- /trunk/Ohana/src/libohana/Makefile (revision 3388) +++ /trunk/Ohana/src/libohana/Makefile (revision 3389) @@ -28,5 +28,6 @@ $(SRC)/Fread.$(ARCH).o $(SRC)/glockfile.$(ARCH).o \ $(SRC)/LoadPhotcodes.$(ARCH).o $(SRC)/photfits.$(ARCH).o \ -$(SRC)/phot_catalog.$(ARCH).o $(SRC)/imreg_datatypes.$(ARCH).o +$(SRC)/phot_catalog.$(ARCH).o $(SRC)/imreg_datatypes.$(ARCH).o \ +$(SRC)/mosaic_astrom.$(ARCH).o $(DESTLIB)/libohana.a: $(LIB)/libohana.$(ARCH).a Index: /trunk/Ohana/src/libohana/doc/ChangeLog.txt =================================================================== --- /trunk/Ohana/src/libohana/doc/ChangeLog.txt (revision 3388) +++ /trunk/Ohana/src/libohana/doc/ChangeLog.txt (revision 3389) @@ -1,2 +1,10 @@ + + changed meaning of PLY / WRP / DIS (see coordtrans.txt) + added check on ALLOCATE/REALLOCATE to avoid failures on Nvalues < 1 + +libohana-1-2 + added 'STG' to allowed coordinate transforms (inaccurate & stop-gap) + set default Npolyterms to 1 + initialized filemode in glockfile.c libohana-1-1 Index: /trunk/Ohana/src/libohana/doc/coordtrans.txt =================================================================== --- /trunk/Ohana/src/libohana/doc/coordtrans.txt (revision 3388) +++ /trunk/Ohana/src/libohana/doc/coordtrans.txt (revision 3389) @@ -1,2 +1,15 @@ + +I have changed the meaning of coord type PLY from the historical +meaning. The only previous use that is in a DVO-style database is in +the loneos data (all cfht data needed to go to linear terms). these +can be easily converted by changing from PLY to DIS. + +There are now three types of higher-order transformation which may be +used depending on the circumstance. The new naming scheme is: + +PLY : a cartesian transformation (ie, no projection implied) +DIS : equivalent to TAN projection with up-to-3rd order polynomial +WRP : a second-level cartesian transformation, implying the presence + of a DIS coord frame (must be registered with libohana function) { Index: /trunk/Ohana/src/libohana/include/loneos.h =================================================================== --- /trunk/Ohana/src/libohana/include/loneos.h (revision 3388) +++ /trunk/Ohana/src/libohana/include/loneos.h (revision 3389) @@ -410,4 +410,8 @@ int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords); int GetCoords (Coords *coords, Header *header); +int PutCoords (Coords *coords, Header *header); +void RegisterMosaic (Coords *coords); +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); + short int putMi (double value); double getMi (short int value); @@ -458,6 +462,4 @@ void ParseColorTerms (char *terms, float *X, int *N); -void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); - int get_image_type (char *name); char *get_type_name (int type); Index: /trunk/Ohana/src/libohana/src/coordops.c =================================================================== --- /trunk/Ohana/src/libohana/src/coordops.c (revision 3388) +++ /trunk/Ohana/src/libohana/src/coordops.c (revision 3389) @@ -2,7 +2,13 @@ # include +static Coords *mosaic = NULL; + +void RegisterMosaic (Coords *coords) { + mosaic = coords; +} + int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) { - int Polynomial, Zenithal, Cartesian, PseudoCyl; + int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl; char *type; double L, M, X, Y, T, Z; @@ -15,9 +21,14 @@ type = &coords[0].ctype[4]; - /* DIS and PLY are identical & apply projection; WRP is cartesian */ + /* DIS is equiv to TAN with higher order terms + PLY is cartesian with higher order terms + WRP is equiv to PLY, with implied mosaic */ + Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP"); - Zenithal = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM"); - Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); + Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR"); + Zenith1 = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG"); + Zenith2 = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM"); + Zenithal = !strcmp(type, "-ZEA") || Zenith1 || Zenith2; if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE); @@ -45,4 +56,10 @@ *ra = L + coords[0].crval1; *dec = M + coords[0].crval2; + + /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */ + if (!strcmp(type, "-WRP")) { + if (mosaic == NULL) return (FALSE); + XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic); + } return (TRUE); } @@ -60,5 +77,5 @@ /* this is wrong : STG is not TAN - need to put in correct relationships. but is a close approx */ - if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) { + if (Zenith1) { if (R == 0) { stht = 1.0; @@ -70,5 +87,5 @@ } } - if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + if (Zenith2) { ctht = RAD_DEG * R; stht = sqrt (1 - ctht*ctht); @@ -124,5 +141,5 @@ char *type; - int i, status, Polynomial, Zenithal, Cartesian, PseudoCyl; + int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl; double phi, theta; double determ; @@ -139,13 +156,22 @@ Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP"); - Zenithal = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM"); - Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); + Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR"); + Zenith1 = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG"); + Zenith2 = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM"); + Zenithal = !strcmp(type, "-ZEA") || Zenith1 || Zenith2; if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE); /**** Locally Cartesian Projections ****/ if (Cartesian) { - L = (ra - coords[0].crval1); - M = (dec - coords[0].crval2); + if (!strcmp(type, "-WRP")) { + 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); + } } @@ -164,9 +190,9 @@ if (stht < 0) status = FALSE; - if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) { + if (Zenith1) { L = DEG_RAD * sphi / stht; M = -DEG_RAD * cphi / stht; } - if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + if (Zenith2) { L = DEG_RAD * sphi; M = -DEG_RAD * cphi; Index: /trunk/Ohana/src/libohana/src/mosaic_astrom.c =================================================================== --- /trunk/Ohana/src/libohana/src/mosaic_astrom.c (revision 3389) +++ /trunk/Ohana/src/libohana/src/mosaic_astrom.c (revision 3389) @@ -0,0 +1,122 @@ +# include +# include + +/* chip-match table: j = ChipMatch[i], images[j] is DIS for images[i] WRP */ + +static int Ndis = 0; +static int *DISentry = NULL; +static e_time *DIStzero = NULL; +static int *ChipMatch = NULL; + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage (Image *images, int Nimages, int entry) { + + int status; + + if (ChipMatch == NULL) { + status = FindMosaicForImage_TableSearch (images, Nimages, entry); + return (status); + } + status = FindMosaicForImage_MatchSearch (images, Nimages, entry); + return (status); +} + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage_TableSearch (Image *images, int Nimages, int entry) { + + e_time tzero; + int i; + + /* search backwards for image with same time and type DIS */ + + if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { + /* not a wrp image, do nothing */ + return (TRUE); /* error or not */ + } + + tzero = images[entry].tzero; + + for (i = entry - 1; i >= 0; i--) { + if (images[i].tzero != tzero) continue; + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + + /* found a valid match */ + RegisterMosaic (&images[i].coords); + } + + for (i = entry + 1; i < Nimages; i++) { + if (images[i].tzero != tzero) continue; + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + + /* found a valid match */ + RegisterMosaic (&images[i].coords); + return (TRUE); + } + return (FALSE); +} + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage_MatchSearch (Image *images, int Nimages, int entry) { + + int N; + + if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { + /* not a wrp image, do nothing */ + return (TRUE); /* error or not */ + } + + if (ChipMatch == NULL) return (FALSE); + + N = ChipMatch[entry]; + if (N < 0) return (FALSE); + if (N >= Nimages) return (FALSE); + + RegisterMosaic (&images[N].coords); + return (TRUE); +} + +int BuildChipMatch (Image *images, int Nimages) { + + int i, j, NDIS; + + if (DISentry != NULL) free (DISentry); + if (DIStzero != NULL) free (DIStzero); + if (ChipMatch != NULL) free (ChipMatch); + + /* find all DIS images, save tzero (& instrument) */ + Ndis = 0; + NDIS = 100; + ALLOCATE (DISentry, int, NDIS); + ALLOCATE (DIStzero, e_time, NDIS); + + for (i = 0; i < Nimages; i++) { + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + DISentry[Ndis] = i; + DIStzero[Ndis] = images[i].tzero; + Ndis ++; + if (Ndis >= NDIS) { + NDIS += 100; + REALLOCATE (DISentry, int, NDIS); + REALLOCATE (DIStzero, e_time, NDIS); + } + } + + /* find all matched WRP images */ + ALLOCATE (ChipMatch, int, Nimages); + for (i = 0; i < Nimages; i++) { + ChipMatch[i] = -1; + if (strcmp(&images[i].coords.ctype[4], "-WRP")) continue; + for (j = 0; j < Ndis; j++) { + if (DIStzero[j] != images[i].tzero) continue; + ChipMatch[i] = DISentry[j]; + goto got_match; + } + + fprintf (stderr, "WARNING: can't find matching mosaic \n"); + ChipMatch[i] = -2; + + got_match: + ChipMatch[i] = j; + } + return (TRUE); +}