IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38518


Ignore:
Timestamp:
Jun 23, 2015, 6:50:35 AM (11 years ago)
Author:
eugene
Message:

loadgalphot mostly done

Location:
branches/eam_branches/ipp-20150616/Ohana/src/addstar
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/Makefile

    r38496 r38518  
    256256$(SRC)/loadgalphot_readstars.$(ARCH).o \
    257257$(SRC)/loadgalphot_fit2d.$(ARCH).o \
     258$(SRC)/loadgalphot_join.$(ARCH).o \
    258259$(SRC)/loadgalphot_table.$(ARCH).o \
     260$(SRC)/strhash.$(ARCH).o \
     261$(SRC)/sortIDs.$(ARCH).o \
    259262$(SRC)/psps_ids.$(ARCH).o
    260263

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/include/loadgalphot.h

    r38501 r38518  
    3030} GalPhot_Stars;
    3131
     32typedef struct {
     33  int *ID;
     34  int *type;
     35  int  N;
     36} GalPhotIDset;
     37
    3238AddstarClientOptions args_loadgalphot (int argc, char **argv, AddstarClientOptions options);
    3339
    … …  
    5359int fit2d (Fit2D *fit, float *xval, float *yval, float *zval, float *zfit, char *mask, int Npts);
    5460int fit2d_reset (Fit2D *fit);
     61
     62int *join_IDs (GalPhotIDset *gal, GalPhotIDset *fit);
     63int strhash (char *string, int module);
     64
     65void sort_IDs_seqonly (int *X, int *S, int N);

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_join.c

    r38511 r38518  
    11# include "addstar.h"
    22# include "loadgalphot.h"
    3 
    4 typedef struct {
    5   int *ID;
    6   int *type;
    7   int  N;
    8 } GalPhotIDset;
    93
    104// find the entries in 'fit' which match entries in 'gal'
    … …  
    126int *join_IDs (GalPhotIDset *gal, GalPhotIDset *fit) {
    137 
    14   off_t ifit, igal, first_ifit, Ifit, Igal, *seq_fit, *seq_gal, Nmatch, NMATCH;
    15   opihi_int dID;
     8  int Nmatch = 0;
     9  int NMATCH = gal->N;
     10  int *idx_gal = NULL;
     11  ALLOCATE (idx_gal, int, gal->N);
    1612
    17   NMATCH = gal->N;
    18   ALLOCATE (idx_gal, off_t, gal->N);
    19 
    20   ALLOCATE (seq_gal, off_t, gal->N);
    21   ALLOCATE (seq_fit, off_t, fit->N);
    22 
    23   for (i = 0; i < gal->N; i++) { seq_gal[i] = i; }
    24   for (i = 0; i < fit->N; i++) { seq_fit[i] = i; }
     13  int *seq_fit, *seq_gal;
     14  ALLOCATE (seq_gal, int, gal->N);
     15  ALLOCATE (seq_fit, int, fit->N);
     16  {
     17    int i;
     18    for (i = 0; i < gal->N; i++) { seq_gal[i] = i; }
     19    for (i = 0; i < fit->N; i++) { seq_fit[i] = i; }
     20  }
    2521
    2622  // sort the sequences
    … …  
    2824  sort_IDs_seqonly (fit->ID, seq_fit, fit->N);
    2925
    30   Nmatch = 0;
     26  int ifit, igal;
    3127  for (ifit = igal = 0; (igal < gal->N) && (ifit < fit->N);) {
    32     Igal = seq_gal[i];
    33     Ifit = seq_fit[j];
     28    int Igal = seq_gal[igal];
     29    int Ifit = seq_fit[ifit];
    3430
    35     dID = gal->ID[Igal] - fit->ID[Ifit];
     31    int dID = gal->ID[Igal] - fit->ID[Ifit];
    3632
    3733    if (dID < 0) { igal++; continue; }
    … …  
    4036    // look for all matches of list2() to list1(i)
    4137    // this allows for multiple values of ID1 or ID2
    42     first_ifit = ifit;
    43     for (ifit = first_ifit; (dID == 0) && (ifit < fit->N); ifit++) {
    44       Ifit = seq_fit[j];
     38    int ifit_first = ifit;
     39    int found = FALSE;
     40    for (ifit = ifit_first; (dID == 0) && (ifit < fit->N); ifit++) {
     41      Ifit = seq_fit[ifit];
    4542      dID = gal->ID[Igal] - fit->ID[Ifit];
    46       if (dID == 0) {
    47         // XXX find the matching model types
    48         index1->elements.Int[Nmatch] = I;
    49         index2->elements.Int[Nmatch] = J;
    50 
     43      if ((dID == 0) && (gal->type[Igal] == fit->type[Ifit])){
     44        idx_gal[Igal] = Ifit;
     45        found = TRUE;
    5146        Nmatch ++;
    5247        if (Nmatch >= NMATCH) {
    53           NMATCH += DMATCH;
    54           REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
    55           REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
     48          NMATCH += 1000;
     49          REALLOCATE (idx_gal, int, NMATCH);
    5650        }
     51        break;
    5752      }
    5853    }
    59     j = first_j;
    60     i++;
     54    myAssert (found, "gal entry not found?");
     55    ifit = ifit_first;
     56    igal ++;
    6157  }
    62   index1->Nelements = Nmatch;
    63   index2->Nelements = Nmatch;
     58  free (seq_gal);
     59  free (seq_fit);
    6460
    65   free (N1);
    66   free (N2);
     61  myAssert (Nmatch == gal->N, "did not find matches for all galaxies");
     62 
     63  for (igal = 0; igal < gal->N; igal++) {
     64    ifit = idx_gal[igal];
     65    myAssert (fit->ID[ifit]   == gal->ID[igal],   "ID mis-match");
     66    myAssert (fit->type[ifit] == gal->type[igal], "type mis-match");
     67  }
    6768
    68   return (TRUE);
     69  return (idx_gal);
    6970}

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_readstars.c

    r38511 r38518  
    2222  if (f == NULL) Shutdown ("can't read stellar parameter file: %s", filename);
    2323
    24   int i, Ncol, Nsample;
     24  int i, j, Ncol, Nsample;
    2525  off_t Nrow = 0;
    2626  off_t NrowNew;
    … …  
    3535    fclose (f);
    3636    return NULL;
    37   }
    38 
    39   // in xgal, we have numerical values to specify the model types
    40   // in xfit, we have strings.  these are listed with their numerical
    41   // match in the PHU.  I need to read the list of names and set up a string hash
    42 
    43   if (!gfits_scan (PHU, "MTNUM", "%d", 1, &ModelTypeN)) { myAbort ("fail"); }
    44   ALLOCATE (ModelTypeName, char *, ModelTypeN);
    45   ALLOCATE (ModelTypeNum,  int, ModelTypeN);
    46   ALLOCATE (ModelTypeHash, int, ModelTypeN);
    47 
    48   int tmpnumber, namehash[256][3];
    49   for (i = 0; i < 256; i++) {
    50     namehash[i][0] = -1;
    51     namehash[i][1] = -1;
    52     namehash[i][2] = -1;
    53   }
    54 
    55   char name[256], tmpname[256];
    56   for (i = 0; i < ModelTypeN; i++) {
    57     snprintf (name, 64, "MTNAM%02d", i);
    58     if (!gfits_scan (PHU, name, "%s", 1, tmpname)) { myAbort ("fail"); }
    59     ModelTypeName[i] = strcreate (tmpname);
    60     int found = FALSE;
    61     int myHash = strhash (tmpname);
    62     ModelTypeHash[i] = myHash;
    63     for (j = 0; (j < 3) && !found; j++) {
    64       if (namehash[myHash][j] != -1) continue;
    65       found = TRUE;
    66       namehash[myHash][j] = i;
    67     }
    68     myAssert (found, "hash conflict");
    69     snprintf (name, 64, "MTVAL%02d", i);
    70     if (!gfits_scan (PHU, name, "%d", 1, &tmpnumber)) { myAbort ("fail"); }
    71     ModelTypeNum[i] = tmpnumber;
    7237  }
    7338
    … …  
    8146    return NULL;
    8247  }
    83 
    8448  ftable_xfit.header = &header_xfit;
     49
     50  // in xgal, we have numerical values to specify the model types
     51  // in xfit, we have strings.  these are listed with their numerical
     52  // match in the PHU.  I need to read the list of names and set up a string hash
     53
     54  int ModelTypeN;
     55  char **ModelTypeName;
     56  int *ModelTypeNum;
     57  int *ModelTypeHash;
     58  if (!gfits_scan (&header_xfit, "MTNUM", "%d", 1, &ModelTypeN)) { myAbort ("fail"); }
     59  ALLOCATE (ModelTypeName, char *, ModelTypeN);
     60  ALLOCATE (ModelTypeNum,  int, ModelTypeN);
     61  ALLOCATE (ModelTypeHash, int, ModelTypeN);
     62
     63  int tmpnumber, namehash[256][3];
     64  for (i = 0; i < 256; i++) {
     65    namehash[i][0] = -1;
     66    namehash[i][1] = -1;
     67    namehash[i][2] = -1;
     68  }
     69
     70  char name[256], tmpname[256];
     71  for (i = 0; i < ModelTypeN; i++) {
     72    snprintf (name, 64, "MTNAM%02d", i);
     73    if (!gfits_scan (&header_xfit, name, "%s", 1, tmpname)) { myAbort ("fail"); }
     74    ModelTypeName[i] = strcreate (tmpname);
     75    int found = FALSE;
     76    int myHash = strhash (tmpname, 0xff);
     77    fprintf (stderr, "model %s, hash: %d\n", tmpname, myHash);
     78    ModelTypeHash[i] = myHash;
     79    for (j = 0; (j < 3) && !found; j++) {
     80      if (namehash[myHash][j] != -1) continue;
     81      found = TRUE;
     82      namehash[myHash][j] = i;
     83    }
     84    myAssert (found, "hash conflict");
     85    snprintf (name, 64, "MTVAL%02d", i);
     86    if (!gfits_scan (&header_xfit, name, "%d", 1, &tmpnumber)) { myAbort ("fail"); }
     87    ModelTypeNum[i] = tmpnumber;
     88  }
    8589
    8690  if (!gfits_fread_ftable_data (f, &ftable_xfit, FALSE)) {
    … …  
    9498  firstCol = TRUE;
    9599  GET_COLUMN (xfit, ID_fit,         "IPP_IDET",      int);
    96   GET_COLUMN (xfit, MODEL_TYPE_str, "MODEL_TYPE",    char); NcharModel = Ncol;
     100  GET_COLUMN (xfit, MODEL_TYPE_str, "MODEL_TYPE",    char);  int NcharModel = Ncol;
    97101  GET_COLUMN (xfit, EXT_WIDTH_MAJ,  "EXT_WIDTH_MAJ", float);
    98102  GET_COLUMN (xfit, EXT_WIDTH_MIN,  "EXT_WIDTH_MIN", float);
    … …  
    102106  int Nfit = Nrow;
    103107
     108  myAssert (NcharModel < 256, "max model type name is very long: %d", NcharModel);
     109
    104110  // free the memory associated with the FITS files
    105111  gfits_free_header (&header_xfit);
    106112  gfits_free_table (&ftable_xfit);
    107113
     114  char string[256];
     115
    108116  // convert MODEL_TYPE_str to MODEL_TYPE_fit (transforming strings to values)
     117  int *MODEL_TYPE_fit;
    109118  ALLOCATE (MODEL_TYPE_fit, int, Nfit);
    110119  for (i = 0; i < Nfit; i++) {
    111     int myHash = strhash (MODEL_TYPE_str[i]);
     120    memset (string, 0, 256);
     121    memcpy (string, &MODEL_TYPE_str[i*NcharModel], NcharModel); string[NcharModel] = 0;
     122    stripwhite (string);
     123    int myHash = strhash (string, 0xff);
    112124    int found = FALSE;
    113125    for (j = 0; !found && (j < 3); j++) {
    114126      int n = namehash[myHash][j];
    115       if (strcmp(ModelTypeName[n], MODEL_TYPE_str[i])) continue;
     127      if (strcmp(ModelTypeName[n], string)) continue;
    116128      MODEL_TYPE_fit[i] = ModelTypeNum[n];
    117129      found = TRUE;
    118130    }
    119     myAssert (found, "unknown model name? %s", MODEL_TYPE_str[i]);
     131    myAssert (found, "unknown model name? %s", string);
    120132  }
    121133
    … …  
    155167  GET_COLUMN (xgal, FR_MINOR_MAX,  "FR_MINOR_MAX", float);
    156168  GET_COLUMN (xgal, FR_MINOR_DEL,  "FR_MINOR_DEL", float);
    157   Ngal = Nrow;
     169  int Ngal = Nrow;
    158170
    159171  // i need to match the two lists based on (ID_gal == ID_fit), (MODEL_TYPE_gal == MODEL_TYPE_fit)
    160   GalPhotIDset fit, gal;
    161   fit.ID   = ID_fit;
    162   fit.type = MODEL_TYPE_fit;
    163   fit.N    = Nfit;
    164   gal.ID   = ID_gal;
    165   gal.type = MODEL_TYPE_gal;
    166   gal.N    = Ngal;
    167   int *index = join_IDs (&gal, &fit);
     172  GalPhotIDset fitSet, galSet;
     173  fitSet.ID   = ID_fit;
     174  fitSet.type = MODEL_TYPE_fit;
     175  fitSet.N    = Nfit;
     176  galSet.ID   = ID_gal;
     177  galSet.type = MODEL_TYPE_gal;
     178  galSet.N    = Ngal;
     179  int *idx_gal = join_IDs (&galSet, &fitSet);
    168180
    169181  // free the memory associated with the FITS files
    … …  
    192204  for (i = 0; i < Nrow; i++) {
    193205
     206    int ifit = idx_gal[i];
     207
    194208    double R, D;
    195209    XY_to_RD (&R, &D, X_FIT[i], Y_FIT[i], &coords);
    … …  
    208222    // theta, theta_err, index,
    209223    // note that FR_MAJOR, etc are the fractions of the stack fit value
    210     stars[i].galphot.theta = NAN;
    211     stars[i].galphot.thetaErr = NAN;
     224    stars[i].galphot.theta     = EXT_THETA[ifit];
     225    stars[i].galphot.thetaErr  = EXT_THETA_ERR[ifit];
     226    stars[i].galphot.index     = INDEX[ifit];
     227    stars[i].galphot.Npix      = NPIX[i];
     228    stars[i].galphot.modelType = MODEL_TYPE_gal[i];
     229    stars[i].galphot.detID     = ID_gal[i];
    212230
    213231    // I have a grid of measurements with (Flux, dFlux, Chisq) at each point
    … …  
    217235      Nbad ++;
    218236    }
     237
     238    // I could either multiply FR_MAJOR_MIN, etc above or the fitted values below
     239    stars[i].galphot.majorAxis    *= EXT_WIDTH_MAJ[ifit];
     240    stars[i].galphot.majorAxisErr *= EXT_WIDTH_MAJ[ifit];
     241    stars[i].galphot.minorAxis    *= EXT_WIDTH_MIN[ifit];
     242    stars[i].galphot.minorAxisErr *= EXT_WIDTH_MIN[ifit];
    219243  }
    220244
    … …  
    347371  // use bilinear interpolation to choose Flux @ (Xmin,Ymin)?
    348372  galphot->mag    = -2.5*log10(flux[iMin]); // correct for exptime?
    349   galphot->magErr = -2.5*log10(fluxErr[iMin]); // correct for exptime?
     373  galphot->magErr = fluxErr[iMin] / flux[iMin];
    350374  galphot->chisq  = chisqMin;
    351375 
Note: See TracChangeset for help on using the changeset viewer.