IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10937


Ignore:
Timestamp:
Jan 5, 2007, 10:51:17 AM (19 years ago)
Author:
eugene
Message:

extensive work to support MEF images and SPLIT multichip process

Location:
trunk/Ohana/src/addstar
Files:
4 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/Makefile

    r10342 r10937  
    4141$(SRC)/fakeimage.$(ARCH).o \
    4242$(SRC)/find_matches.$(ARCH).o \
    43 $(SRC)/find_matches_closest.$(ARCH).o \
    4443$(SRC)/find_matches_refstars.$(ARCH).o \
    4544$(SRC)/find_subset.$(ARCH).o \
    … …  
    5352$(SRC)/getusno.$(ARCH).o \
    5453$(SRC)/getusnob.$(ARCH).o \
    55 $(SRC)/gimages.$(ARCH).o \
    5654$(SRC)/greference.$(ARCH).o \
    5755$(SRC)/grefstars.$(ARCH).o \
    … …  
    6361$(SRC)/replace_match.$(ARCH).o \
    6462$(SRC)/resort_catalog.$(ARCH).o \
    65 $(SRC)/rfits.$(ARCH).o \
    66 $(SRC)/rtext.$(ARCH).o \
    67 $(SRC)/sort_lists.$(ARCH).o \
     63$(SRC)/ReadStarsFITS.$(ARCH).o \
     64$(SRC)/ReadStarsTEXT.$(ARCH).o \
     65$(SRC)/FilterStars.$(ARCH).o \
     66$(SRC)/sort_lists.$(ARCH).o \
     67$(SRC)/ImageOptions.$(ARCH).o \
     68$(SRC)/GetFileMode.$(ARCH).o \
     69$(SRC)/ReadImageHeader.$(ARCH).o \
    6870$(SRC)/update_coords.$(ARCH).o
     71
     72# $(SRC)/gimages.$(ARCH).o \
     73# $(SRC)/find_matches_closest.$(ARCH).o \
    6974
    7075ADDSTARD = \

  • trunk/Ohana/src/addstar/doc/notes.txt

    r6763 r10937  
     1
     22007.01.04
     3
     4  Use of options.photcode:
     5    - M_REFLIST: used to assign a photcode value to the loaded data
     6    - M_REFCAT: used to select data from the reference of the given photcode
     7    - M_FAKEIMAGE: used to assign a photcode value to the generated data
     8    - M_RESORT: unused
     9    - M_IMAGE: is used to override or supply the image header photcode
     10               * only valid for single chip runs
     11
     12  Re: MOSAIC_PHU
     13      - the mosaic phu needs to have the time to enable the image/mosaic match
     14
     15      - NX,NY are used for plotting the mosaic outline.  the center of
     16        the mosaic is 0,0, and the mosaic outline is drawn from
     17        -NX/2,-NY/2 to NX/2,NY/2.  A problem: NX,NY in Image is using
     18        a short int.  If the units of the FPA/TPA coords are microns, this overflows.
     19        * change NX,NY to int?
     20        * require the use of mm not microns?
     21        * for GPC, NX,NY ~ 384000 um excluding gaps
     22
     23      - only the astrometry, time, and NX,NY are loaded for MOSAIC_PHU headers.
     24
     252007.01.03
     26
     27  I am working on adding the ability to read MEF files containing
     28  multiple tables and multiple image headers.  Here are the details:
     29
     30  SIMPLE CMP:
     31    PHU:
     32      NAXIS = 2 or
     33      SIMPLE = FALSE or
     34      TEXTMODE = TRUE (use option; global variable)
     35
     36  SIMPLE CMF:
     37    PHU:
     38      NAXIS = 0
     39      NAXIS1 exists, > 0
     40      NAXIS2 exists, > 0
     41    EXT:
     42      EXTNAME = SMPDATA, PS1DATA, etc (defines layout)
     43
     44  SIMPLE MEF:
     45    PHU:
     46      NAXIS = 0
     47      NAXIS1 does not exist
     48      NAXIS2 does not exist
     49      &CTYPE[4] is not "-DIS"
     50      EXTEND is true
     51    EXT HEAD:
     52      EXTTYPE = IMAGE
     53      EXTNAME = image.name
     54      EXTDATA = table.name
     55    EXT DATA:
     56      EXTTYPE = SMPDATA, PS1DATA, etc (defines layout)
     57      EXTNAME = table.name
     58      EXTHEAD = image.name
     59
     60  MOSAIC MEF:
     61    PHU:
     62      NAXIS = 0
     63      NAXIS1 does not exist
     64      NAXIS2 does not exist
     65      &CTYPE[4] is "-DIS"
     66      EXTEND is true
     67    EXT HEAD:
     68      EXTTYPE = IMAGE
     69      EXTNAME = image.name
     70      EXTDATA = table.name
     71    EXT DATA:
     72      EXTTYPE = SMPDATA, PS1DATA, etc (defines layout)
     73      EXTNAME = table.name
     74      EXTHEAD = image.name
     75
     76  MOSAIC PHU:
     77    PHU:
     78      NAXIS = 0
     79      NAXIS1 does not exist
     80      NAXIS2 does not exist
     81      &CTYPE[4] is "-DIS"
     82      NEXTEND = 0 or 
     83      EXTEND is false
     84
     85   
    186
    2872006.04.02

  • trunk/Ohana/src/addstar/include/addstar.h

    r10889 r10937  
    1111
    1212/* used in find_matches, find_matches_refstars */
    13 # define IN_CATALOG(R,D) ( \
     13# define IN_REGION(R,D) ( \
    1414((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \
    1515((R) >= region[0].Rmin)  && ((R) < region[0].Rmax))
    … …  
    3535
    3636enum {M_IMAGE, M_REFLIST, M_REFCAT, M_FAKEIMAGE, M_RESORT};
    37 
    38 /* global which define database info / data sources (KEEP) */
     37enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU};
     38/* note: MEF implies CMF */
     39
     40/* globals which define database info / data sources (KEEP) */
    3941char   ImageCat[256];
    4042char   GSCFILE[256];
    … …  
    108110
    109111/*** addstar prototypes ***/
     112
    110113void       AddToCalibration       PROTO((Average *average, Measure *measure, Measure *new, int *next, int Nstar));
    111114AddstarClientOptions ConfigInit   PROTO((int *argc, char **argv));
    … …  
    126129int        edge_check             PROTO((double *x1, double *y1, double *x2, double *y2));
    127130Image     *fakeimage              PROTO((char *rootname, int *Nimage, int photcode));
    128 int        find_matches           PROTO((SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options));
    129 int        find_matches_closest   PROTO((SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options));
     131int        find_matches           PROTO((SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
     132int        find_matches_closest   PROTO((SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
    130133int        find_matches_refstars  PROTO((SkyRegion *region, Stars **stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
    131134Stars    **find_subset            PROTO((SkyRegion *region, Stars *stars, int Nstars, int *NSTARS));
    … …  
    140143Stars     *grefcat                PROTO((char *Refcat, SkyRegion *catstats, int photcode, int *nstars));
    141144Stars     *grefstars              PROTO((char *file, int photcode, int *Nstars));
    142 Stars     *gstars                 PROTO((char *file, int *NSTARS, int photcode, Image *image));
     145Stars     *gstars                 PROTO((char *file, int *Nstars, Image **images, int *Nimages, int photcode));
    143146int        in_image               PROTO((double r, double d, Image *image));
    144147int        load_pt_catalog        PROTO((Catalog *catalog, SkyRegion *region));  /*** choose new name ***/
    … …  
    156159Stars     *rfits                  PROTO((FILE *f, unsigned int *nstars));
    157160Stars     *rtext                  PROTO((FILE *f, unsigned int *nstars));
     161Stars     *ReadStarsFITS          PROTO((FILE *f, Header *header, Header *in_theader, unsigned int *nstars));
     162Stars     *ReadStarsTEXT          PROTO((FILE *f, unsigned int *nstars));
     163int        ReadImageHeader        PROTO((Header *header, Image *image, int photcode));
     164Stars     *FilterStars            PROTO((Stars *instars, Image *image));
     165Stars     *MergeStars             PROTO((Stars *stars, int *Nstars, Stars *instars, int Ninstars));
    158166void       save_pt_catalog        PROTO((Catalog *catalog));  /*** choose new name ***/
    159167double     scat_subpix            PROTO((double x, double y));
    … …  
    169177void       create_image_db        PROTO((FITS_DB *db));
    170178void       set_db                 PROTO((FITS_DB *in));
    171 
    172179void       uppercase              PROTO((char *string));
    173180void       fsort                  PROTO((float *X, int N));
    174181void       fsort2                 PROTO((float *X, float *Y, int N));
    175 
    176 int *init_measure_links (Average *average, int Naverage, Measure *measure, int Nmeasure);
    177 int *init_missing_links (Average *average, int Naverage, Missing *missing, int Nmissing);
    178 int add_meas_link (Average *average, int *next, int Nmeasure, int NMEASURE);
    179 int add_miss_link (Average *average, int *next, int Nmissing);
    180 int *build_measure_links (Average *average, int Naverage, Measure *measure, int Nmeasure);
    181 Measure *sort_measure (Average *average, int Naverage, Measure *measure, int Nmeasure, int *next);
    182 Missing *sort_missing (Average *average, int Naverage, Missing *missing, int Nmissing, int *next_miss);
     182int       *init_measure_links     PROTO((Average *average, int Naverage, Measure *measure, int Nmeasure));
     183int       *init_missing_links     PROTO((Average *average, int Naverage, Missing *missing, int Nmissing));
     184int        add_meas_link          PROTO((Average *average, int *next, int Nmeasure, int NMEASURE));
     185int        add_miss_link          PROTO((Average *average, int *next, int Nmissing));
     186int       *build_measure_links    PROTO((Average *average, int Naverage, Measure *measure, int Nmeasure));
     187Measure   *sort_measure           PROTO((Average *average, int Naverage, Measure *measure, int Nmeasure, int *next));
     188Missing   *sort_missing           PROTO((Average *average, int Naverage, Missing *missing, int Nmissing, int *next_miss));
     189Stars     *ConvertSMPDATA         PROTO((FTable *table, int *nstars));
     190// Stars *ConvertPS1DATA          PROTO((FTable *table, int *nstars));
     191int        ImageOptions           PROTO((AddstarClientOptions *options, Image *images, int Nimages));
     192int        GetFileMode            PROTO((Header *header));
     193AddstarClientOptions args_client  PROTO((int argc, char **argv, AddstarClientOptions options));
     194AddstarClientOptions args_load2mass PROTO((int argc, char **argv, AddstarClientOptions options));
     195AddstarClientOptions args_sedstar PROTO((int argc, char **argv, AddstarClientOptions options));
     196void       args_server            PROTO((int argc, char **argv));
     197int        CheckPassword          PROTO((int BindSocket));
     198int        NewImage               PROTO((int BindSocket));
     199int        NewReflist             PROTO((int BindSocket));
     200int        NewRefcat              PROTO((int BindSocket));
     201int        InitServerSocket       PROTO((SockAddress *Address));
     202int        WaitServerSocket       PROTO((int InitSocket, SockAddress *Address, int *validIP, int Nvalid));
     203int        GetClientSocket        PROTO((char *hostname));
     204int        UpdateDatabase_Image   PROTO((AddstarClientOptions *options, Image *image, Coords *mosaic, Stars *stars, int Nstars));
     205int        UpdateDatabase_Reflist PROTO((AddstarClientOptions *options, Stars *stars, int Nstars));
     206int        UpdateDatabase_Refcat  PROTO((AddstarClientOptions *options, SkyRegion *UserPatch, char *refcat));
     207SkyList   *SkyListForStars        PROTO((SkyTable *table, int depth, Stars *stars, int Nstars));
     208SkyList   *SkyListExistingSubset  PROTO((SkyList *input, char *path));
     209int        SkyListSetPath         PROTO((SkyList *list, char *path));
     210SkyTable  *SkyTableLoadOptimal    PROTO(());
     211int        InitDataset            PROTO(());
     212int        PushDataset            PROTO((DVO_DATA *data));
     213DVO_DATA  *PopDataset             PROTO((void));
     214void      *ListenClients_Thread   PROTO((void *data));
     215int        NewImage_Thread        PROTO((int BindSocket));
     216int        NewRefcat_Thread       PROTO((int BindSocket));
     217int        NewReflist_Thread      PROTO((int BindSocket));
     218
     219// this is a gnu extension?? caution!
     220void *memrchr(const void *s, int c, size_t n);
    183221
    184222/**
    … …  
    191229**/
    192230
     231
    193232/** function for client / server **/
    194233
    195 AddstarClientOptions args_client (int argc, char **argv, AddstarClientOptions options);
    196 AddstarClientOptions args_load2mass (int argc, char **argv, AddstarClientOptions options);
    197 AddstarClientOptions args_sedstar (int argc, char **argv, AddstarClientOptions options);
    198 
    199 void args_server (int argc, char **argv);
    200 
    201 int CheckPassword (int BindSocket);
    202 int NewImage (int BindSocket);
    203 int NewReflist (int BindSocket);
    204 int NewRefcat (int BindSocket);
    205 
    206 int InitServerSocket (SockAddress *Address);
    207 int WaitServerSocket (int InitSocket, SockAddress *Address, int *validIP, int Nvalid);
    208 int GetClientSocket (char *hostname);
    209 
    210 int UpdateDatabase_Image (AddstarClientOptions *options, Image *image, Coords *mosaic, Stars *stars, int Nstars);
    211 int UpdateDatabase_Reflist (AddstarClientOptions *options, Stars *stars, int Nstars);
    212 int UpdateDatabase_Refcat (AddstarClientOptions *options, SkyRegion *UserPatch, char *refcat);
    213 
    214 SkyList *SkyListForStars (SkyTable *table, int depth, Stars *stars, int Nstars);
    215 SkyList *SkyListExistingSubset (SkyList *input, char *path);
    216 int SkyListSetPath (SkyList *list, char *path);
    217 SkyTable *SkyTableLoadOptimal ();
    218 
    219 int InitDataset ();
    220 int PushDataset (DVO_DATA *data);
    221 DVO_DATA *PopDataset (void);
    222 void *ListenClients_Thread (void *data);
    223 int NewImage_Thread (int BindSocket);
    224 int NewRefcat_Thread (int BindSocket);
    225 int NewReflist_Thread (int BindSocket);
    226 
    227 // this is a gnu extension?? caution!
    228 void *memrchr(const void *s, int c, size_t n);

  • trunk/Ohana/src/addstar/src/FilterStars.c

    r10897 r10937  
    44
    55  int j, N;
     6  float MTIME, dMs;
     7  Stars *stars;
     8
     9  /* correct instrumental mags for exposure time */
     10  MTIME = (image[0].exptime > 0) ? 2.500*log10(image[0].exptime) : 0.0;
    611
    712  /* modify resulting star list */
    8   ALLOCATE (stars, Stars, image[0].nstars);
    9   for (N = j = 0; j < image[0].nstars; j++) {
     13  ALLOCATE (stars, Stars, image[0].nstar);
     14  for (N = j = 0; j < image[0].nstar; j++) {
    1015    /* allow for some dynamic filtering of star list */
    1116    if (SNLIMIT && instars[j].dM > SNLIMIT) continue;
    … …  
    2025    while (stars[N].R >= 360.0) stars[N].R -= 360.0;
    2126    stars[N].found = -1;
    22     stars[N].code = photcode;
     27    stars[N].code = image[0].source;
    2328
     29    /** additional quantities to supply to Stars based on the image data **/
     30
     31    /* calculate accurate per-star airmass */
     32    stars[N].airmass = airmass (image[0].secz_PS, stars[N].R, stars[N].D, image[0].sidtime, image[0].latitude);
     33    stars[N].Mcal    = image[0].Mcal_PS;
     34    stars[N].t       = image[0].tzero + 1e-4*stars[N].Y*image[0].trate;  /* trate is in 0.1 msec / row */
     35
     36    stars[N].M       = MIN (stars[N].M + MTIME, NO_MAG);
     37    stars[N].dM      = MIN (stars[N].dM, NO_ERR);
     38
     39    stars[N].Mgal    = MIN (stars[N].Mgal + MTIME, NO_MAG);
     40   
    2441    if (SUBPIX) {
    2542      dMs = get_subpix (stars[N].X, stars[N].Y);
    … …  
    4057}
    4158
    42 Stars *MergeStars (Stars *instars, int Ninstars, Stars *stars, int *Nstars) {
     59Stars *MergeStars (Stars *stars, int *Nstars, Stars *instars, int Ninstars) {
     60
     61  int i, j;
    4362
    4463  if (stars == NULL) {
    … …  
    4867  }
    4968
    50   for (j = 0, i = *Nstars; i < *Nstars + Ninstars; i++) {
     69  for (j = 0, i = *Nstars; i < *Nstars + Ninstars; i++, j++) {
    5170    stars[i] = instars[j];
    5271  }

  • trunk/Ohana/src/addstar/src/ReadStarsFITS.c

    r10897 r10937  
    55Stars *ReadStarsFITS (FILE *f, Header *header, Header *in_theader, unsigned int *nstars) {
    66
    7   int i, Nstars;
     7  int i, Nskip, Nstars;
     8  char type[80];
    89  Header theader;
    910  FTable table;
    … …  
    3334  }
    3435
    35   if (!gfits_scan (header, "EXTTYPE", "%s", 1, type)) {
     36  if (!gfits_scan (&theader, "EXTTYPE", "%s", 1, type)) {
    3637    strcpy (type, "SMPDATA");
    3738  }
    … …  
    4142    stars = ConvertSMPDATA (&table, &Nstars);
    4243  }
     44# if (0)
    4345  if (!strcmp (type, "PS1DATA")) {
    4446    stars = ConvertPS1DATA (&table, &Nstars);
    4547  }
     48# endif
    4649  if (stars == NULL) {
    4750    fprintf (stderr, "ERROR: invalid table type %s\n", type);
    … …  
    4952  }
    5053  if (*nstars != Nstars) {
    51     fprintf (stderr, "ERROR: inconsistent number of stars? %d vs %d\n", *nstars, Nstars);
    52     exit (1);
     54    fprintf (stderr, "WARNING: inconsistent number of stars? %d vs %d\n", *nstars, Nstars);
    5355  }
     56  *nstars = Nstars;
    5457  return stars;
    5558}
    … …  
    5760Stars *ConvertSMPDATA (FTable *table, int *nstars) {
    5861
    59   int i, Nstars;
    60   Stars *stars;
    61   SMPData *smpdata;
     62  int i, Nstars, swapped;
     63  Stars *stars = NULL;
     64  SMPData *smpdata = NULL;
    6265
    63   smpdata = gfits_table_get_SMPData (table, &Nstars, NULL);
     66  swapped = FALSE;
     67  smpdata = gfits_table_get_SMPData (table, &Nstars, &swapped);
     68  /* XXX we need to check at least the size of the loaded table */
    6469
    6570  ALLOCATE (stars, Stars, Nstars);
    … …  
    6873    stars[i].Y      = smpdata[i].Y;
    6974    stars[i].M      = smpdata[i].M;
    70     stars[i].dM     = smpdata[i].dM;
     75    stars[i].dM     = smpdata[i].dM*0.001;
    7176    stars[i].dophot = smpdata[i].dophot;
    7277
    … …  
    8186}
    8287
     88# if (0)
    8389Stars *ConvertPS1DATA (FTable *table, int *nstars) {
    8490
    … …  
    106112  return (stars);
    107113}
     114# endif

  • trunk/Ohana/src/addstar/src/Shutdown.c

    r8458 r10937  
    11# include "addstar.h"
    2 
    3 static FITS_DB *db;
    4 
    5 void set_db (FITS_DB *in) {
    6   db = in;
    7 }
    82
    93/* clean up open / locked ImageCat before shutting down */
    … …  
    2115  va_end (argp);
    2216
    23   SetProtect (TRUE);
    24   if (db != NULL) gfits_db_close (db);
    2517  fprintf (stderr, "ERROR: addstar halted\n");
    2618  exit (1);

  • trunk/Ohana/src/addstar/src/addstar.c

    r10889 r10937  
    44
    55  int Nmatch, status;
    6   int i, Nstars, Nimages, Noverlap, Nsubset;
     6  int i, Nstars, Nimages, Nsubset;
    77  unsigned long long Naverage, Nmeasure;
    88  Stars *stars, **subset;
    9   Image *images, *overlap;
     9  Image *images;
    1010  Catalog catalog;
    1111  FITS_DB db;
    … …  
    2929 
    3030  stars = NULL;
    31   overlap = NULL;
    32   set_db (&db);
    3331
    34   /* we use the image table to lock db access -- perhaps this is not necessary? */
    35   db.mode   = dvo_catalog_catmode (CATMODE);
    36   db.format = dvo_catalog_catformat (CATFORMAT);
    37   status    = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD);
    38   if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
    39 
     32  /*** load in the new data (images, stars) ***/
    4033  switch (options.mode) {
    4134    case M_IMAGE:
    42       stars = gstars (argv[1], &Nstars, options.photcode, &images, &Nimages);
     35      stars = gstars (argv[1], &Nstars, &images, &Nimages, options.photcode);
    4336      if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars (stars, Nstars);
    4437      RegisterMosaic (MOSAIC);
    4538      for (i = 0; i < Nimages; i++) {
    46         newlist = SkyListByImage (sky, -1, &image);
    47         SkyListMergeLists (&skylist, &newlist);
     39        newlist = SkyListByImage (sky, -1, &images[i]);
     40        SkyListMerge (&skylist, newlist);
     41        SkyListFree (newlist, FALSE);
    4842      }
    49       overlap = gimages (&db, &image, MOSAIC, &Noverlap);
     43      ImageOptions (&options, images, Nimages);
    5044      break;
    5145    case M_REFLIST:
    … …  
    6054      images = fakeimage (argv[1], &Nimages, options.photcode);
    6155      RegisterMosaic (MOSAIC);
     56      ALLOCATE (skylist, SkyList, 1);
     57      skylist[0].Nregions = 0;
     58      break;
    6259
    63       if (db.dbstate == LCK_EMPTY) {
    64         dvo_image_create (&db, ZeroPt);
    65       } else {
    66         if (!dvo_image_load (&db, VERBOSE, FORCE_READ)) {
    67           Shutdown ("can't read image catalog %s", db.filename);
    68         }
    69       }   
     60    default:
     61      fprintf (stderr, "ERROR: invalid mode\n");
     62      exit (2);
     63  }
    7064
    71       dvo_image_addrows (&db, imageSet, Nimages);
    72       SetProtect (TRUE);
    73       dvo_image_update (&db, VERBOSE);
    74       SetProtect (FALSE);
    75 
    76       ohana_memcheck (FALSE);
    77       dvo_image_unlock (&db); /* unlock? */
    78 
    79       gettimeofday (&stop, NULL);
    80       dtime = DTIME (stop, start);
    81       fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for fake image\n", dtime);
    82       exit (0);
    83   }
     65  // in these cases, limit the sky catalogs to an existing subset
    8466  if (options.only_match || options.existing_regions) {
    8567    SkyList *tmp;
    … …  
    9274  /* XXX clean this up a bit : for only_image, I don't need to do this loop
    9375     unless we are getting the calibration. */
     76  /*** match stars to existing catalog data (or otherwise manipulate catalog data) ***/
    9477  Nmatch = Naverage = Nmeasure = 0;
    9578  for (i = 0; i < skylist[0].Nregions; i++) {
    … …  
    11699    }
    117100
    118     Naverage += catalog.Naverage;
    119     Nmeasure += catalog.Nmeasure;
    120 
    121101    switch (options.mode) {
    122102      case M_IMAGE:
    123103        Nsubset = Nstars;
    124104        if (options.closest) {
    125           Nmatch += find_matches_closest (skylist[0].regions[i], stars, Nstars, &catalog, &image, overlap, Noverlap, MOSAIC, options);
     105          // Nmatch += find_matches_closest (skylist[0].regions[i], stars, Nstars, &catalog, options);
    126106        } else {
    127           Nmatch += find_matches (skylist[0].regions[i], stars, Nstars, &catalog, &image, overlap, Noverlap, MOSAIC, options);
     107          Nmatch += find_matches (skylist[0].regions[i], stars, Nstars, &catalog, options);
    128108        }
    129109        break;
    … …  
    141121        break;
    142122    }
     123    /* report total updated values */
     124    Naverage += catalog.Naverage;
     125    Nmeasure += catalog.Nmeasure;
    143126
    144127    // write out catalog, if appropriate
    … …  
    154137  }
    155138
    156   // XXX is it necessary to lock the image catalog during this entire process?
    157   if (options.calibrate) { FindCalibration (&image); }
     139  // We only measure a single value for the entire mosaic (add all images to this function)
     140  if (options.calibrate) { FindCalibration (&images[0]); }
    158141
     142  /*** update the image table ***/
     143  /* setup image table format and lock */
     144  db.mode   = dvo_catalog_catmode (CATMODE);
     145  db.format = dvo_catalog_catformat (CATFORMAT);
     146  status    = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD);  // shorter timeout?
     147  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
     148
     149  /* load or create the image table */
    159150  if (db.dbstate == LCK_EMPTY) {
    160151    if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat);
    161152    dvo_image_create (&db, ZeroPt);
     153  } else {
     154    if (!dvo_image_load (&db, VERBOSE, FORCE_READ)) {
     155      Shutdown ("can't read image catalog %s", db.filename);
     156    }
    162157  }
    163  
    164   ohana_memcheck (FALSE);
    165158
    166   /* write out new image */
     159  /* add the new images and save */
    167160  if (options.mode == M_IMAGE) {
    168     dvo_image_addrows (&db, &image, 1);
     161    dvo_image_addrows (&db, images, Nimages);
    169162    SetProtect (TRUE);
    170163    dvo_image_update (&db, VERBOSE);
    171164    SetProtect (FALSE);
    172165  }
    173   ohana_memcheck (FALSE);
    174166  dvo_image_unlock (&db); /* unlock? */
    175167
    … …  
    177169  dtime = DTIME (stop, start);
    178170  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for %5d stars (%5d matches), %6lld average, %7lld measure\n", dtime, Nstars, Nmatch, Naverage, Nmeasure);
     171
    179172  exit (0);
    180173}
    … …  
    185178   patch   - RA,DEC bounded portion of sky
    186179*/
     180
     181// add in case of failures:
     182// ohana_memcheck (FALSE);

  • trunk/Ohana/src/addstar/src/args.c

    r7691 r10937  
    137137  }
    138138  /* don't add missed pts to Missed table (image only) */
    139   options.skip_missed = FALSE;
     139  options.skip_missed = TRUE;
    140140  if ((N = get_argument (argc, argv, "-missed"))) {
    141141    options.skip_missed = TRUE;
    142142    remove_argument (N, &argc, argv);
     143    fprintf (stderr, "ERROR: addstar no longer supports -missed\n");
     144    exit (2);
    143145  }
    144146  /* replace measurement, don't duplicate (ref/cat only) */

  • trunk/Ohana/src/addstar/src/find_matches.c

    r10889 r10937  
    11# include "addstar.h"
    22
    3 int find_matches (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options) {
     3int find_matches (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, AddstarClientOptions options) {
    44
    55  int i, j, n, N, J, status, Nstars;
    6   double X, Y, RADIUS, RADIUS2, secz;
     6  double X, Y, RADIUS, RADIUS2;
    77  float *X1, *Y1, *X2, *Y2;
    88  float dX, dY, dR;
    9   int *N1, *N2,  *next_meas, *next_miss;
    10   int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch;
     9  int *N1, *N2,  *next_meas;
     10  int Nave, NAVE, Nmeas, NMEAS, Nmatch;
    1111  int Nsecfilt, Nsec;
    1212  float Mcat, *Mval, MTIME;
    … …  
    1414  Coords tcoords;
    1515
    16   /* changes needed to handle mosaic MEF images:
    17      - options should not carry photcode: this is set per star in gstars
    18      - validate that all images have the same equiv code value (same Nsec)
    19      - add dt to Stars.d and set in gstars, not here
    20      - add airmass to Stars.d and set in gstars, not here
    21      - set stars.t in gstars
    22      - disallow mosaic and missed?  otherwise, need to loop over all images
    23        for each missed source to see if it should have been found
    24      - Mcal_PS is supplied by image, but is not correctly set: it is
    25        calculated *after* find_matches in addstar
    26   */
    27 
    28   /* XXX EAM : options.photcode overridden by image.source.... */
    29   /* XXX this function could be modified to handle multiple photcodes as long
    30      as they match the same equiv code value, resulting in the same value of Nsec */
    31   code = GetPhotcodebyCode (image[0].source);
    32 
    33   /* photcode data - must by of type DEP, (PRI, SEC) - probably should restrict to DEP */
    34   // XXX : rely on catalog[0].Nsecfilt or GetPhotcodeNsecfilt??
     16  /* photcode data - must by of type DEP; options.photcode is equiv PRI/SEC photcode */
     17  /* this function requires incoming stars to have the same photcode.equiv value */
    3518  Nsecfilt = GetPhotcodeNsecfilt ();
    36   Nsec = (code[0].type == PHOT_DEP) ? GetPhotcodeNsec (code[0].equiv) : GetPhotcodeNsec (code[0].code);
    37   /* this function requires incoming stars to have the same photcode */
     19  Nsec     = GetPhotcodeNsec (options.photcode);
    3820
    3921  /** allocate local arrays (stars) **/
    … …  
    5335  Nmatch = 0;
    5436  NMEAS = Nmeas = catalog[0].Nmeasure;
    55   NMISS = Nmiss = catalog[0].Nmissing;
    5637 
    5738  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ZEA projection has the
    58      advantage that every point in R,D has a mapping to a unique X,Y.  However, note that not all
    59      possible X,Y points map back to R,D and the local plate scale changes substantially far from
    60      the projection pole.  a better mapping might be ARC, not yet implemented (see
    61      coordops.update.c).  We keep the original crpix1,2 and crref1,2.  For mosaic astrometry, the
    62      grid should be w.r.t. the tangent-plane, not chip coords */
    63 
    64   if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
    65     tcoords = mosaic[0];
    66     tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    67     tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    68     tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    69     tcoords.Npolyterms = 1;
    70     strcpy (tcoords.ctype, "RA---ZEA");
     39   * advantage that every point in R,D has a mapping to a unique X,Y.  However, note that not all
     40   * possible X,Y points map back to R,D and the local plate scale changes substantially far from
     41   * the projection pole.  a better mapping might be ARC, not yet implemented (see
     42   * coordops.update.c).  We use the center of the region (catalog) for crval1,2.
     43   */
     44  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
     45  if (region[0].Rmax < 90) {
     46    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
    7147  } else {
    72     tcoords = image[0].coords;
    73     tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    74     tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    75     tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    76     tcoords.Npolyterms = 1;
    77     strcpy (tcoords.ctype, "RA---ZEA");
    78   }
     48    tcoords.crval2 = 90.0;
     49  }
     50  tcoords.crpix1 = 0;
     51  tcoords.crpix2 = 0;
     52  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
     53  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
     54  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
     55  tcoords.Npolyterms = 1;
     56  strcpy (tcoords.ctype, "RA---ZEA");
    7957
    8058  /* build spatial index (RA sort) */
    … …  
    9573    free (Y2);
    9674    free (N2);
    97     return;
     75    return (0);
    9876  }
    9977  if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars);
    … …  
    10785  if (Nave > 1) sort_lists (X2, Y2, N2, Nave);
    10886
    109   /* set up pointers for linked list of measure, missing */
     87  /* set up pointers for linked list of measure */
    11088  if (catalog[0].sorted) {
    11189    next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
    … …  
    11391    next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
    11492  }   
    115   next_miss = init_missing_links (catalog[0].average, Nave, catalog[0].missing, Nmiss);
    116   /* missing MUST be written 'sorted', or not at all */
    117 
    118   /* choose a radius for matches */
    119   if (options.radius == 0) {
    120     RADIUS = options.Nsigma * 0.02 * image[0].cerror;  /* 0.02 corrects cerror to arcsec from storage units */
    121   } else {
    122     RADIUS = options.radius; /* provided by config */
    123   }
     93
     94  /* choose a radius for matches (defined in args.c or ImageOptions.c) */
     95  RADIUS = options.radius;
    12496  RADIUS2 = RADIUS*RADIUS;
    125 
    126   /* correct instrumental mags for exposure time */
    127   MTIME = (image[0].exptime > 0) ? 2.500*log10(image[0].exptime) : 0.0;
    12897
    12998  /** find matched stars **/
    … …  
    163132      add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
    164133
    165       /* calculate accurate per-star airmass */
    166       secz = airmass (image[0].secz_PS, stars[N].R, stars[N].D, image[0].sidtime, image[0].latitude);
    167      
    168134      /** add measurements for this star **/
    169135      /** dR,dD now represent arcsec **/
    … …  
    187153      }
    188154      catalog[0].measure[Nmeas].dD_PS       = 3600.0*(catalog[0].average[n].D - stars[N].D);
    189       catalog[0].measure[Nmeas].M_PS        = MIN (stars[N].M + MTIME, NO_MAG);
    190       catalog[0].measure[Nmeas].dM_PS       = MIN (stars[N].dM, NO_ERR);  /* error in input files stored in thousandths of mag */
    191       catalog[0].measure[Nmeas].Mcal_PS     = image[0].Mcal_PS;
    192       catalog[0].measure[Nmeas].t           = image[0].tzero + 1e-4*stars[N].Y*image[0].trate;  /* trate is in 0.1 msec / row */
     155      catalog[0].measure[Nmeas].M_PS        = stars[N].M;
     156      catalog[0].measure[Nmeas].dM_PS       = stars[N].dM;  /* error in input files stored in thousandths of mag */
     157      catalog[0].measure[Nmeas].Mcal_PS     = stars[N].Mcal;
     158      catalog[0].measure[Nmeas].t           = stars[N].t;
    193159      catalog[0].measure[Nmeas].averef      = n;              /* this must be an absolute sequence number, if partial average is loaded */
    194160      catalog[0].measure[Nmeas].source      = stars[N].code;  /* photcode */
    195161      catalog[0].measure[Nmeas].dophot      = stars[N].dophot; 
    196162      catalog[0].measure[Nmeas].flags       = 0;
    197       catalog[0].measure[Nmeas].dt_PS       = MTIME;
    198       catalog[0].measure[Nmeas].airmass_PS  = secz;
    199 
    200       catalog[0].measure[Nmeas].Mgal_PS     = MIN (stars[N].Mgal + MTIME, NO_MAG);
     163      catalog[0].measure[Nmeas].dt_PS       = stars[N].dt;
     164      catalog[0].measure[Nmeas].airmass_PS  = stars[N].airmass;
     165
     166      catalog[0].measure[Nmeas].Mgal_PS     = stars[N].Mgal;
    201167      catalog[0].measure[Nmeas].FWx         = MIN (100*stars[N].fx, NO_MAG);
    202168      catalog[0].measure[Nmeas].FWy         = MIN (100*stars[N].fy, NO_MAG);
    … …  
    247213  }
    248214
    249   /* add reference for undetected catalog stars */
    250   if (!strcmp (&image[0].coords.ctype[4], "-WRP")) RegisterMosaic (mosaic);
    251   for (j = 0; (j < Nave) && !options.skip_missed; j++) {
    252     n = N2[j];
    253     if (catalog[0].found[n] < 0) {
    254       /* make sure there is space for next entry */
    255       if (Nmiss >= NMISS) {
    256         NMISS = Nmiss + 1000;
    257         REALLOCATE (next_miss, int, NMISS);
    258         REALLOCATE (catalog[0].missing, Missing, NMISS);
    259       }
    260 
    261       /* should the catalog star be on this image? project into image coords */
    262       if (!in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) continue;
    263       add_miss_link (&catalog[0].average[n], next_miss, Nmiss);
    264 
    265       /* calculate time of exposure for this coordinate in the image */
    266       RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords);     
    267       catalog[0].missing[Nmiss].t  = image[0].tzero + 1e-4*Y*image[0].trate;  /* trate is in 0.1 msec / row */
    268       catalog[0].average[n].Nn ++;
    269       Nmiss ++;
    270     }
    271   }
    272 
    273215  /* incorporate unmatched image stars, if this star is in field of this catalog */
    274216  /* these new entries are all written out in UPDATE mode */
    … …  
    288230    N = N1[i];
    289231    if (stars[N].found >= 0) continue;
    290     if (!IN_CATALOG (stars[N].R, stars[N].D)) continue;
    291 
    292     secz = airmass (image[0].secz_PS, stars[N].R, stars[N].D, image[0].sidtime, image[0].latitude);
     232    if (!IN_REGION (stars[N].R, stars[N].D)) continue;
    293233
    294234    catalog[0].average[Nave].R         = stars[N].R;
    … …  
    322262    catalog[0].measure[Nmeas].dR_PS       = 0.0;
    323263    catalog[0].measure[Nmeas].dD_PS       = 0.0;
    324     catalog[0].measure[Nmeas].M_PS        = MIN (stars[N].M + MTIME, NO_MAG);
    325     catalog[0].measure[Nmeas].dM_PS       = MIN (stars[N].dM, NO_ERR);
    326     catalog[0].measure[Nmeas].Mcal_PS     = image[0].Mcal_PS;
    327     catalog[0].measure[Nmeas].t           = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */
     264    catalog[0].measure[Nmeas].M_PS        = stars[N].M;
     265    catalog[0].measure[Nmeas].dM_PS       = stars[N].dM;
     266    catalog[0].measure[Nmeas].Mcal_PS     = stars[N].Mcal;
     267    catalog[0].measure[Nmeas].t           = stars[N].t;
    328268    catalog[0].measure[Nmeas].averef      = Nave;           /* XXX EAM : must be absolute Nave if partial read */
    329269    catalog[0].measure[Nmeas].source      = stars[N].code;  /* photcode */
    330270    catalog[0].measure[Nmeas].dophot      = stars[N].dophot; 
    331271    catalog[0].measure[Nmeas].flags       = 0;
    332     catalog[0].measure[Nmeas].dt_PS       = MTIME;
    333     catalog[0].measure[Nmeas].airmass_PS  = secz;
     272    catalog[0].measure[Nmeas].dt_PS       = stars[N].dt;
     273    catalog[0].measure[Nmeas].airmass_PS  = stars[N].airmass;
    334274
    335275    catalog[0].measure[Nmeas].Mgal_PS     = MIN (stars[N].Mgal + MTIME, NO_MAG);
    … …  
    341281    Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M_PS;
    342282    if (*Mval == NO_MAG) *Mval = Mcat;
    343 
    344     /** now add references from all previous non-detection observations of this spot on the sky */
    345     for (j = 0; (j < Noverlap) && !options.skip_missed; j++) {
    346       /* make sure there is space for next entry */
    347       if (Nmiss >= NMISS) {
    348         NMISS = Nmiss + 1000;
    349         REALLOCATE (next_miss, int, NMISS);
    350         REALLOCATE (catalog[0].missing, Missing, NMISS);
    351       }
    352       if (!FindMosaicForImage (overlap, Noverlap, j)) continue;
    353       if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue;
    354       add_miss_link (&catalog[0].average[Nave], next_miss, Nmiss);
    355 
    356       /* get time of exposure of this portion of the image */
    357       RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j].coords);     
    358       catalog[0].missing[Nmiss].t  = overlap[j].tzero + 1e-4*Y*overlap[j].trate;  /* rough guess at time */
    359       catalog[0].average[Nave].Nn ++;
    360       Nmiss ++;
    361     }
    362283
    363284    /* next[Nmeas] should always be -1 in this context (it is always the only
    … …  
    371292  REALLOCATE (catalog[0].average, Average, Nave);
    372293  REALLOCATE (catalog[0].measure, Measure, Nmeas);
    373   REALLOCATE (catalog[0].missing, Missing, Nmiss);
    374294 
    375295  if (options.nosort) {
    … …  
    379299    catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
    380300  }
    381   catalog[0].missing = sort_missing (catalog[0].average, Nave, catalog[0].missing, Nmiss, next_miss);
    382   /* missing is REQUIRED to be sorted */
    383301
    384302  /* note stars which have been found in this catalog */
    … …  
    394312  catalog[0].Naverage = Nave;
    395313  catalog[0].Nmeasure = Nmeas;
    396   catalog[0].Nmissing = Nmiss;
    397   if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
     314  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmatch);
    398315
    399316  free (catalog[0].found);

  • trunk/Ohana/src/addstar/src/gstars.c

    r10897 r10937  
    11# include "addstar.h"
    22
    3 enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU};
    4 /* note: MEF implies CMF */
    5 
    6 Stars *gstars (char *filename, int *NSTARS, int photcode, Image **imageList, int Nimages) {
    7 
    8   int i, Nfile, Nimage;
     3// XXX where and when should I call RegisterMosaic??
     4
     5Stars *gstars (char *filename, int *Nstars, Image **images, int *Nimages, int photcode) {
     6
     7  int i, j, Nfile, Nheader, NHEADER, Nimage, NIMAGE;
     8  int Nskip, Nhead, Ndata, done, status, mode, NinStars;
    99  char **file;
    1010  FILE *f;
    1111  glob_t globList;
    12   char **exthead, **extdata;
     12  char **exthead, **extdata, **exttype, tmpword[80];
    1313  int *extnum_head, *extnum_data, *extsize;
    1414  Header *header, **headers;
    … …  
    1818  // parse the filename as a glob
    1919  globList.gl_offs = 0;
    20   glob (file, 0, NULL, &globList);
     20  glob (filename, 0, NULL, &globList);
    2121
    2222  // if the glob does not match, save the literal word:
    … …  
    4646
    4747  stars = NULL;
    48   Nstars = 0;
    49 
    50   // valid for all but MEF files
     48  *Nstars = 0;
     49
     50  /*** load data from a single PHU or a collection of PHU files ***/
    5151  if ((mode != SIMPLE_MEF) && (mode != MOSAIC_MEF)) {
    5252    Nimage = Nfile;
    … …  
    6262      }
    6363
    64       ReadImageHeader (header, &image[i]);
     64      ReadImageHeader (header, &image[i], photcode);
    6565   
    6666      switch (mode) {
    6767        case SIMPLE_CMP:
    6868        case MOSAIC_CMP:
    69           inStars = ReadStarsTEXT (f, header, NULL, &image[i].nstars);
     69          inStars = ReadStarsTEXT (f, &image[i].nstar);
    7070          inStars = FilterStars (inStars, &image[i]);
    71           stars = MergeStars (inStars, NinStars, stars, &Nstars);
     71          stars = MergeStars (stars, Nstars, inStars, image[i].nstar);
    7272          break;
    7373
    7474        case SIMPLE_CMF:
    7575        case MOSAIC_CMF:
    76           inStars = ReadStarsFITS (f, header, NULL, &image[i].nstars);
     76          inStars = ReadStarsFITS (f, header, NULL, &image[i].nstar);
    7777          inStars = FilterStars (inStars, &image[i]);
    78           stars = MergeStars (inStars, NinStars, stars, &Nstars);
     78          stars = MergeStars (stars, Nstars, inStars, image[i].nstar);
    7979          break;
    8080
    … …  
    8787      gfits_free_header (header);
    8888    }
    89   } else {
    90    
    91     /* we need to examine the extensions to determine the headers and the data */
    92     NHEADER = 10;
    93     ALLOCATE (headers, Header *, NHEADER);
    94 
    95     // the first header is already loaded
    96     headers[0] = header;
    97     Nskip = gfits_matrix_size (header);
    98     fseek (f, Nskip, SEEK_CUR);
    99 
    100     // load all headers into memory
    101     done = FALSE;
    102     for (i = 1; !done; i++) {
     89    *Nimages = Nimage;
     90    *images = image;
     91   
     92    return stars;
     93  }
     94   
     95  /* we have a multi-chip image */
     96
     97  /* supplied photcode is incompatible with multi-chip images */
     98  if (photcode) {
     99    fprintf (stderr, "ERROR: photcode cannot be supplied to multi-chip images -- manually adjust the headers\n");
     100    exit (1);
     101  }
     102
     103  /* we need to examine the extensions to determine the headers and the data */
     104  NHEADER = 10;
     105  ALLOCATE (headers, Header *, NHEADER);
     106
     107  // the first header is already loaded
     108  headers[0] = header;
     109  Nskip = gfits_matrix_size (header);
     110  fseek (f, Nskip, SEEK_CUR);
     111
     112  // load all headers into memory
     113  done = FALSE;
     114  for (i = 1; !done; i++) {
    103115      ALLOCATE (header, Header, 1);
    104116     
    105117      status = gfits_fread_header (f, header);
    106118      if (!status) {
    107         done = TRUE;
     119          done = TRUE;
    108120      } else {
    109         headers[i] = header;
    110         Nskip = gfits_matrix_size (header);
    111         fseek (f, Nskip, SEEK_CUR);
     121          headers[i] = header;
     122          Nskip = gfits_matrix_size (header);
     123          fseek (f, Nskip, SEEK_CUR);
    112124      }
    113125      if (i == NHEADER - 1) {
    114         NHEADER += 10;
    115         REALLOCATE (headers, Header *, NHEADER);
    116       }
    117     }
    118     Nheader = i;
    119    
    120     // space to store the images, indexes to the matching headers
    121     Nimage = 0;
    122     NIMAGE = Nheader;
    123     ALLOCATE (image, Image, NIMAGE);
    124     ALLOCATE (exthead, char **, NIMAGE);
    125     ALLOCATE (extdata, char **, NIMAGE);
    126     ALLOCATE (extnum_head, int, NIMAGE);
    127     ALLOCATE (extnum_data, int, NIMAGE);
    128     ALLOCATE (extsize, int, Nheader);
    129 
    130     if (mode == MOSAIC_MEF) {
     126          NHEADER += 10;
     127          REALLOCATE (headers, Header *, NHEADER);
     128      }
     129  }
     130  Nheader = i - 1; /* we failed on the last loop */
     131   
     132  // space to store the images, indexes to the matching headers
     133  Nimage = 0;
     134  NIMAGE = Nheader;
     135  ALLOCATE (image, Image, NIMAGE);
     136  ALLOCATE (exthead, char *, NIMAGE);
     137  ALLOCATE (extdata, char *, NIMAGE);
     138  ALLOCATE (extnum_head, int, NIMAGE);
     139  ALLOCATE (extnum_data, int, NIMAGE);
     140  ALLOCATE (extsize, int, Nheader);
     141
     142  if (mode == MOSAIC_MEF) {
    131143      exthead[Nimage] = strcreate ("PHU");
    132144      extdata[Nimage] = strcreate ("NONE");
    … …  
    134146      extnum_head[Nimage] = 0;
    135147      Nimage ++;
    136     }
    137 
    138     // now examine the headers, count the table entries, find corresponding headers
    139     for (i = 0; i < Nheader; i++) {
     148  }
     149
     150  // now examine the headers, count the table entries, find corresponding headers
     151  for (i = 0; i < Nheader; i++) {
    140152      extsize[i] = headers[i][0].size + gfits_matrix_size (headers[i]);
    141153      gfits_scan (headers[i], "EXTTYPE", "%s", 1, tmpword);
    142154      if (!strcmp (tmpword, "SMPDATA") || 
    143155          !strcmp (tmpword, "PS1DATA")) {
    144         exttype[Nimage] = strcreate (tmpword);
    145         gfits_scan (headers[i], "EXTNAME", "%s", 1, tmpword);
    146         extdata[Nimage] = strcreate (tmpword);
    147         gfits_scan (header[i], "EXTHEAD", "%s", 1, tmpword);
    148         exthead[Nimage] = strcreate (tmpword);
    149         extnum_data[Nimage] = i;
    150         extnum_head[Nimage] = -1;
    151         // find the matching exthead entry
    152         for (j = 0; j < Nheader; j++) {
    153           gfits_scan (header[j], "EXTNAME", "%s", 1, tmpword);
    154           if (!strcmp (tmpword, exthead[Nimage])) {
    155             extnum_head[Nimage] = j;
     156          exttype[Nimage] = strcreate (tmpword);
     157          gfits_scan (headers[i], "EXTNAME", "%s", 1, tmpword);
     158          extdata[Nimage] = strcreate (tmpword);
     159          gfits_scan (headers[i], "EXTHEAD", "%s", 1, tmpword);
     160          exthead[Nimage] = strcreate (tmpword);
     161          extnum_data[Nimage] = i;
     162          extnum_head[Nimage] = -1;
     163          // find the matching exthead entry
     164          for (j = 0; j < Nheader; j++) {
     165            if (gfits_scan (headers[j], "EXTNAME", "%s", 1, tmpword)) {
     166              if (!strcmp (tmpword, exthead[Nimage])) {
     167                extnum_head[Nimage] = j;
     168              }
     169            }
    156170          }
    157         }
    158         // skip or crash on table with missing matching header?
    159         if (extnum_head[Nimage] == -1) {
    160           fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
    161           exit (1);
    162         }
    163         Nimage ++;
    164       }
    165     }
    166 
    167     // now run through the images, interpret the headers and read the stars
    168     for (i = 0; i < Nimage; i++) {
     171          // skip or crash on table with missing matching header?
     172          if (extnum_head[Nimage] == -1) {
     173              fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
     174              exit (1);
     175          }
     176          Nimage ++;
     177      }
     178  }
     179  if (VERBOSE) fprintf (stderr, "file %s has %d headers, including %d images\n", file[0], Nheader, Nimage);
     180
     181  // now run through the images, interpret the headers and read the stars
     182  for (i = 0; i < Nimage; i++) {
    169183      Nhead = extnum_head[i];
    170       ReadImageHeader (header[Nhead], &image[i]);
     184
     185      if (VERBOSE) fprintf (stderr, "reading header for %s (%s)\n", exthead[i], extdata[i]);
     186      ReadImageHeader (headers[Nhead], &image[i], 0);
    171187
    172188      if (!strcmp(exthead[i], "PHU")) continue;
    … …  
    176192      Nskip = 0;
    177193      for (j = 0; j < Ndata; j++) {
    178         Nskip += extsize[j];
     194          Nskip += extsize[j];
    179195      }
    180196      fseek (f, Nskip, SEEK_SET);
    181197         
    182       inStars = ReadStarsFITS (f, header[Nhead], header[Ndata], &NinStars);
    183       inStars = FilterStars (inStars, &NinStars, &image[i]);
    184       stars = MergeStars (inStars, NinStars, stars, &Nstars);
    185     }
     198      inStars = ReadStarsFITS (f, headers[Nhead], headers[Ndata], &image[i].nstar);
     199      inStars = FilterStars (inStars, &image[i]);
     200      stars = MergeStars (stars, Nstars, inStars, image[i].nstar);
     201  }
     202  *Nimages = Nimage;
     203  *images = image;
     204
     205  return stars;
     206}

  • trunk/Ohana/src/addstar/src/parse_time.c

    r8428 r10937  
    1414  if (strcasecmp (JDKeyword, "NONE")) {
    1515    uppercase (JDKeyword);
    16     gfits_scan (header, JDKeyword, "%lf", 1, &jd);
     16    if (!gfits_scan (header, JDKeyword, "%lf", 1, &jd)) {
     17      fprintf (stderr, "ERROR: missing JD Keyword %s\n", JDKeyword);
     18      exit (1);
     19    }
    1720    Nsec = (jd - 2440587.5)*86400;
    1821    return (Nsec);
    … …  
    2225  if (strcasecmp (MJDKeyword, "NONE")) {
    2326    uppercase (MJDKeyword);
    24     gfits_scan (header, MJDKeyword, "%lf", 1, &jd);
     27    if (!gfits_scan (header, MJDKeyword, "%lf", 1, &jd)) {
     28      fprintf (stderr, "ERROR: missing MJD Keyword %s\n", MJDKeyword);
     29      exit (1);
     30    }
    2531    Nsec = (jd - 40587.0)*86400;
    2632    return (Nsec);
    … …  
    4248  /* get UT and DATE */
    4349  uppercase (UTKeyword);
    44   gfits_scan (header, UTKeyword, "%s", 1, line);
     50  if (!gfits_scan (header, UTKeyword, "%s", 1, line)) {
     51      fprintf (stderr, "ERROR: missing UT Keyword %s\n", UTKeyword);
     52      exit (1);
     53    }
    4554  /* remove ':' characters */
    4655  for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; }
    … …  
    8190  /* parse date entry */
    8291  uppercase (DateKeyword);
    83   gfits_scan (header, DateKeyword, "%s",  1, line);
     92  if (!gfits_scan (header, DateKeyword, "%s",  1, line)) {
     93    fprintf (stderr, "ERROR: missing DATE Keyword %s\n", DateKeyword);
     94    exit (1);
     95  }
    8496  /* remove possible separators: ':', '/' '.', '-' */
    8597  for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; }
Note: See TracChangeset for help on using the changeset viewer.