IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7696


Ignore:
Timestamp:
Jun 26, 2006, 4:35:44 PM (20 years ago)
Author:
eugene
Message:

finished sedstar

Location:
trunk/Ohana/src/addstar
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/Makefile

    r7691 r7696  
    1717# 
    1818INCS    =       -I$(INC) -I$(LINC) -I$(XINC)
    19 LIBS    =       -L$(LLIB) -ldvo -lFITS -lohana -lz -lm
     19LIBS    =       -L$(LLIB) -ldvo -lkapa -lFITS -lohana -lX11 -lz -lm
    2020CFLAGS  =       $(INCS)
    2121LFLAGS  =       $(LIBS)
    … …  
    101101$(SRC)/replace_match.$(ARCH).o \
    102102$(SRC)/update_coords.$(ARCH).o \
     103$(SRC)/gcatalog.$(ARCH).o \
     104$(SRC)/mkcatalog.$(ARCH).o \
     105$(SRC)/wcatalog.$(ARCH).o \
     106$(SRC)/ConfigInit.$(ARCH).o \
     107$(SRC)/Shutdown.$(ARCH).o \
     108$(SRC)/SetSignals.$(ARCH).o
     109
     110SEDSTAR = \
     111$(SRC)/sedstar.$(ARCH).o \
     112$(SRC)/SEDtableLoad.$(ARCH).o \
     113$(SRC)/SEDfit.$(ARCH).o \
     114$(SRC)/SEDops.$(ARCH).o \
     115$(SRC)/args_sedstar.$(ARCH).o \
    103116$(SRC)/gcatalog.$(ARCH).o \
    104117$(SRC)/mkcatalog.$(ARCH).o \
    … …  
    256269load2mass.install        : $(DESTBIN)/load2mass
    257270
     271sedstar                  : $(BIN)/sedstar.$(ARCH)
     272$(BIN)/sedstar.$(ARCH)   : $(SEDSTAR)
     273sedstar.install          : $(DESTBIN)/sedstar
     274
    258275all: addstar addstarc addstard addstart
    259276

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

    r7691 r7696  
    9191double CAL_INSTMAG_MIN;
    9292int    VERBOSE;
     93int    PLOT;
    9394
    9495/* modify server behavior (make this an addstar cleanup mode?) */
    … …  
    191192AddstarClientOptions args_client (int argc, char **argv, AddstarClientOptions options);
    192193AddstarClientOptions args_load2mass (int argc, char **argv, AddstarClientOptions options);
     194AddstarClientOptions args_sedstar (int argc, char **argv, AddstarClientOptions options);
    193195
    194196void args_server (int argc, char **argv);

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

    r7693 r7696  
     1# include "addstar.h"
     2# include "kapa.h"
    13
    2 typedef structu {
    3   int Nfilter;
    4   float *wavecode;
    5   float *vegaToAB;
    6   int *hashcode;
    7   int codeP;
    8   int codeM;
    9   SEDtableRow **row;
    10   int Nrow;
    11 } SEDtable;
     4typedef enum {
     5  SED_FIT,
     6  SED_REQ,
     7  SED_MODEL,
     8  SED_SAMPLE,
     9} SEDtableModes;
    1210
    1311typedef struct {
    … …  
    2422} SEDfit;
    2523
     24typedef struct {
     25  int Nfilter;
     26  float *wavecode;
     27  float *vegaToAB;
     28  int *mode;
     29  int *hashcode;
     30  int *code;
     31  int codeP;
     32  int codeM;
     33  SEDtableRow **row;
     34  int Nrow;
     35} SEDtable;
     36
    2637SEDtableRow **sort_SEDtable (SEDtableRow *raw, int N);
    27 int SEDcolorBracket (SEDtableRow **table, int Ntable, float color);
    2838SEDfit SEDchisq (SEDtableRow *ref, SEDtableRow *data, SEDtableRow *error, int Nfilter);
    29 
     39SEDtable *SEDtableLoad (char *filename);
     40int SEDcolorBracket (SEDtable *table, float color);
     41int SEDfitInit (SEDtable *table);
     42int SEDfitPlot (SEDtable *table, double R, double D, SEDfit *minFit, SEDtableRow *sourceValue, SEDtableRow *sourceError);
     43int SEDfitCatalog (Catalog *outcat, Catalog *incat, SEDtable *table);
     44void SetLimitsRaw (float *xvec, float *yvec, int Nelements, Graphdata *graphmode);

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

    r7693 r7696  
    1 # include "addstar.h"
    21# include "sedstar.h"
    32
    4 int SEDfit (Catalog *outcat, Catalog *incat, SEDtable *table) {
     3int SEDfitCatalog (Catalog *outcat, Catalog *incat, SEDtable *table) {
    54 
    6   Graphdata graphdata;
    7   KapaSection magSection, resSection;
     5  int i, j, m, idx, start, done, row, Nsec, Nfit, Nphot, fd;
     6  int Nave, Nmeas, NAVE, NMEAS;
     7  unsigned short USNOred, USNOblu;
     8  float color;
     9  int *found, valid;
    810
    911  SEDtableRow sourceValue, sourceError;
    1012  SEDfit minFit, testFit;
    1113
    12   /* defaults */
    13   outcat.average = NULL;
    14   outcat.measure = NULL;
    15   outcat.secfilt = NULL;
    16 
    17   SEDtable = NULL;
    18   SEDtableRaw = NULL;
    1914  sourceValue.mags = NULL;
    2015  sourceError.mags = NULL;
    21   wavecode = NULL;
    22   hashcode = NULL;
    23   RegionName = NULL;
    24   RegionList = NULL;
    25   magSection.name = NULL;
    26   resSection.name = NULL;
    2716
    28   oldsignal = signal (SIGINT, handle_interrupt);
    29   interrupt = FALSE;
     17  Nsec = GetPhotcodeNsecfilt ();
     18  Nave = outcat[0].Naverage;
     19  Nmeas = outcat[0].Nmeasure;
    3020
    31   /* load photcode information */
    32   if (!InitPhotcodes ()) goto escape;
    33   Nsec = GetPhotcodeNsecfilt ();
    34 
    35   /* interpret command-line options */
    36   if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
    37   if (!SetPhotSelections (&argc, argv, 4)) goto usage;
    38 
    39   PLOT = FALSE;
    40   if ((N = get_argument (argc, argv, "-plot"))) {
    41     remove_argument (N, &argc, argv);
    42     PLOT = TRUE;
    43   }
    44 
    45   SAVEDIR = NULL;
    46   if ((N = get_argument (argc, argv, "-save"))) {
    47     remove_argument (N, &argc, argv);
    48     SAVEDIR = strcreate (argv[N]);
    49     remove_argument (N, &argc, argv);
    50   }
    51 
    52   /* interpret command-line options */
    53   if (argc != 6) goto usage;
    54 
    55   Nfit = 0;
    56   colorP = GetPhotcodeCodebyName (argv[3]);
    57   colorM = GetPhotcodeCodebyName (argv[5]);
    58   if (!colorP || !colorM) goto color_undefined;
     21  NAVE = 100;
     22  NMEAS = 100;
     23  REALLOCATE (outcat[0].average, Average, NAVE);
     24  REALLOCATE (outcat[0].secfilt, SecFilt, NAVE*Nsec);
     25  REALLOCATE (outcat[0].measure, Measure, NMEAS);
    5926
    6027  // artificially set USNOred and blu errors to 0.3
    … …  
    6229  USNOblu = GetPhotcodeCodebyName ("USNO_BLUE");
    6330
    64   // load SED table
    65   f = fopen (argv[1], "r");
    66   if (f == NULL) goto table_missing;
     31  // create holder for the source data
     32  ALLOCATE (sourceValue.mags, float, table[0].Nfilter);
     33  ALLOCATE (sourceError.mags, float, table[0].Nfilter);
     34  ALLOCATE (found, int, table[0].Nfilter);
    6735
    68   // XXX add error checks for header data
    69   scan_line (f, line);
    70   sscanf (line, "%*s %*s %d", &Nfilter);
     36  if (PLOT) SEDfitInit (table);
    7137
    72   // load SED table photcodes, generate the photcode hashtable
    73   ALLOCATE (hashcode, int, 0x10000);
    74   ALLOCATE (wavecode, float, Nfilter);
    75   ALLOCATE (vegaToAB, float, Nfilter);
     38  // perform the fit to all sources
     39  for (i = 0; i < incat[0].Naverage; i++) {
    7640
    77   for (i = 0; i < 0x10000; i++) hashcode[i] = -1;
    78   for (i = 0; i < Nfilter; i++) {
    79     scan_line (f, line);
    80     sscanf (line, "%*s %s %f %f", name, &wavecode[i], &vegaToAB[i]);
    81     code = GetPhotcodeCodebyName (name);
    82     if (code == 0) goto code_missing;
    83     hashcode[code] = i;
    84   }
    85   codeP = hashcode[colorP];
    86   codeM = hashcode[colorM];
    87   if ((codeP == -1) || (codeM == -1)) goto color_missing;
    88    
    89   // skip remaining header lines
    90   scan_line (f, line);
    91   scan_line (f, line);
    92   scan_line (f, line);
    93   scan_line (f, line);
    94  
    95   // load the SED table data
    96   Nrow = 0;
    97   NROW = 100;
    98   ALLOCATE (SEDtableRaw, SEDtableRow, NROW);
    99   while (scan_line(f, line) != EOF) {
    100     fparse (&SEDtableRaw[Nrow].Temp, 1, line);
    101     fparse (&SEDtableRaw[Nrow].Av, 2, line);
    102     ALLOCATE (SEDtableRaw[Nrow].mags, float, Nfilter);
    103     for (i = 0; i < Nfilter; i++) {
    104       fparse (&SEDtableRaw[Nrow].mags[i], i + 3, line);
     41    // blank out the source array
     42    for (j = 0; j < table[0].Nfilter; j++) {
     43      sourceValue.mags[j] = 100;
     44      found[j] = FALSE;
     45    }   
     46
     47    // load the measurements for this source
     48    m = incat[0].average[i].offset;
     49    Nphot = 0;
     50    for (j = 0; j < incat[0].average[i].Nm; j++) {
     51      idx = table[0].hashcode[incat[0].measure[m+j].source];
     52      if (idx == -1) continue;
     53      // only fit the selected photcodes (mode == "fit")
     54      if (table[0].mode[idx] == SED_MODEL) continue;
     55      if (table[0].mode[idx] == SED_SAMPLE) continue;
     56      // XXX do something more clever if more than one value exists per photcode
     57      sourceValue.mags[idx] = incat[0].measure[m+j].M_PS + table[0].vegaToAB[idx];
     58      sourceError.mags[idx] = incat[0].measure[m+j].dM_PS;
     59      if (incat[0].measure[m+j].source == USNOred) sourceError.mags[idx] = 0.3;
     60      if (incat[0].measure[m+j].source == USNOblu) sourceError.mags[idx] = 0.3;
     61      found[idx] = TRUE;
     62      Nphot ++;
    10563    }
    106     SEDtableRaw[Nrow].color = SEDtableRaw[Nrow].mags[codeP] - SEDtableRaw[Nrow].mags[codeM];
    107     Nrow ++;
    108     CHECK_REALLOCATE (SEDtableRaw, SEDtableRow, NROW, Nrow, 100);
    109   }     
     64    if (Nphot < 3) continue;
    11065
    111   // sort the SEDtable by the reference colors
    112   SEDtable = sort_SEDtable (SEDtableRaw, Nrow);
     66    valid = TRUE;
     67    for (j = 0; valid && (j < table[0].Nfilter); j++) {
     68      if ((table[0].mode[j] == SED_REQ) && !found[j]) valid = FALSE;
     69    }
     70    if (!valid) continue;
    11371
    114   // create holder for the source data
    115   ALLOCATE (sourceValue.mags, float, Nfilter);
    116   ALLOCATE (sourceError.mags, float, Nfilter);
     72    // skip sources without ref color
     73    if (sourceValue.mags[table[0].codeP] > 50) continue;
     74    if (sourceValue.mags[table[0].codeM] > 50) continue;
     75    color = sourceValue.mags[table[0].codeP] - sourceValue.mags[table[0].codeM];
    11776
    118   if (PLOT) {
    119     if (!GetGraph (&graphdata, &fd, NULL)) return (FALSE);
    120     SetLimitsRaw (wavecode, NULL, Nfilter, &graphdata);
    121     graphdata.style = 2;
    122     graphdata.ptype = 2;
    123     KapaClear (fd, TRUE);
    124     magSection.name = strcreate ("mag");
    125     magSection.x  = 0;
    126     magSection.dx = 1;
    127     magSection.y  = 0.5;
    128     magSection.dy = 0.5;
    129     resSection.name = strcreate ("res");
    130     resSection.x  = 0.0;
    131     resSection.dx = 1.0;
    132     resSection.y  = 0.0;
    133     resSection.dy = 0.5;
    134    
    135     KapaSetFont (fd, "helvetica", 14);
    136     ALLOCATE (fitmags, float, Nfilter);
    137     ALLOCATE (fiterrs, float, Nfilter);
    138   }
     77    // find tableRow within 0.1 mag of color
     78    // XXX : check on the delta value
     79    start = SEDcolorBracket (table, color);
     80    minFit = SEDchisq (table[0].row[start], &sourceValue, &sourceError, table[0].Nfilter);
     81    minFit.row = start;
    13982
    140   /* load region corresponding to selection above */
    141   if ((skylist = SelectRegions (RegionName, RegionList)) == NULL) goto escape;
    142 
    143   /* loop over regions, extract data for each region */
    144   // XXX add interrupt checks
    145   fprintf (stderr, "using %d possible regions\n", skylist[0].Nregions);
    146   for (k = 0; k < skylist[0].Nregions; k++) {
    147     /* lock, load, unlock catalog */
    148     catalog.filename = skylist[0].filename[k];
    149     switch (lock_catalog (&catalog, LCK_SOFT)) {
    150       case 2:
    151         unlock_catalog (&catalog);
    152       case 0:
    153         catalog.Naverage = 0;
    154         continue;
    155     }
    156     catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
    157     if (!load_catalog (&catalog, TRUE)) {
    158       catalog.Naverage = 0;
    159       unlock_catalog (&catalog);
    160       continue;
     83    // search for min chisq backwards
     84    // XXX : check on the delta value
     85    done = FALSE;
     86    row = start - 1;
     87    while (!done && (row > 0)) {
     88      testFit = SEDchisq (table[0].row[row], &sourceValue, &sourceError, table[0].Nfilter);
     89      if (testFit.chisq < minFit.chisq) {
     90        minFit = testFit;
     91        minFit.row = row;
     92      }
     93      if (fabs(table[0].row[row][0].color - color) > 0.5) done = TRUE;
     94      row --;
    16195    }
    16296
    163     // perform the fit to all sources
    164     for (i = 0; i < catalog.Naverage; i++) {
     97    // search for min chisq forwards
     98    // XXX : check on the delta value
     99    done = FALSE;
     100    row = start + 1;
     101    while (!done && (row < table[0].Nrow)) {
     102      testFit = SEDchisq (table[0].row[row], &sourceValue, &sourceError, table[0].Nfilter);
     103      if (testFit.chisq < minFit.chisq) {
     104        minFit = testFit;
     105        minFit.row = row;
     106      }
     107      if (fabs(table[0].row[row][0].color - color) > 0.5) done = TRUE;
     108      row ++;
     109    }
    165110
    166       // blank out the source array
    167       for (j = 0; j < Nfilter; j++) {
    168         sourceValue.mags[j] = 100;
    169       }
     111    Nfit ++;
     112    // create the vectors for the example plots
     113    if (PLOT) SEDfitPlot (table, incat[0].average[i].R, incat[0].average[i].D, &minFit, &sourceValue, &sourceError);
    170114
    171       // load the measurements for this source
    172       m = catalog.average[i].offset;
    173       for (j = 0; j < catalog.average[i].Nm; j++) {
    174         idx = hashcode[catalog.measure[m+j].source];
    175         if (idx == -1) continue;
    176         // XXX do something more clever if more than one value exists per photcode
    177         sourceValue.mags[idx] = catalog.measure[m+j].M_PS + vegaToAB[idx];
    178         sourceError.mags[idx] = catalog.measure[m+j].dM_PS;
    179         if ((catalog.measure[m+j].source == USNOred) || (catalog.measure[m+j].source == USNOblu)) {
    180           sourceError.mags[idx] = 0.3;
    181         }
    182       }
     115    // construct an average object for this object
     116    // XXX for now, the output objects will have limited astrometric interpretation...
     117    outcat[0].average[Nave].R         = incat[0].average[i].R;
     118    outcat[0].average[Nave].D         = incat[0].average[i].D;
     119    outcat[0].average[Nave].dR    = 0;
     120    outcat[0].average[Nave].dD    = 0;
     121    outcat[0].average[Nave].uR    = 0;
     122    outcat[0].average[Nave].uD    = 0;
     123    outcat[0].average[Nave].duR   = 0;
     124    outcat[0].average[Nave].duD   = 0;
     125    outcat[0].average[Nave].P     = 0;
     126    outcat[0].average[Nave].dP    = 0;
    183127
    184       // XXX for the moment, skip sources without ref color
    185       if (sourceValue.mags[codeP] > 50) continue;
    186       if (sourceValue.mags[codeM] > 50) continue;
    187       color = sourceValue.mags[codeP] - sourceValue.mags[codeM];
     128    // XXX for now, set the average mag data to NULL
     129    outcat[0].average[Nave].M         = NO_MAG;
     130    outcat[0].average[Nave].dM        = NO_MAG;
     131    outcat[0].average[Nave].Nm        = 0;
     132    outcat[0].average[Nave].Nn        = 0;
     133    outcat[0].average[Nave].Xp        = NO_MAG;
     134    outcat[0].average[Nave].Xm        = NO_MAG;
     135    outcat[0].average[Nave].Xg        = NO_MAG;
     136    outcat[0].average[Nave].offset    = Nmeas;
     137    outcat[0].average[Nave].missing   = -1;
     138    outcat[0].average[Nave].code      = 0;
    188139
    189       // XXX find tableRow within 0.1 mag of color
    190       start = SEDcolorBracket (SEDtable, Nrow, color);
    191       minFit = SEDchisq (SEDtable[start], &sourceValue, &sourceError, Nfilter);
    192       minFit.row = start;
     140    for (j = 0; j < Nsec; j++) {
     141      outcat[0].secfilt[Nave*Nsec+j].M_PS  = NO_MAG;
     142      outcat[0].secfilt[Nave*Nsec+j].dM_PS = NO_MAG;
     143      outcat[0].secfilt[Nave*Nsec+j].Xm    = NO_MAG;
     144    }
    193145
    194       // search for min chisq backwards
    195       done = FALSE;
    196       row = start - 1;
    197       while (!done && (row > 0)) {
    198         testFit = SEDchisq (SEDtable[row], &sourceValue, &sourceError, Nfilter);
    199         if (testFit.chisq < minFit.chisq) {
    200           minFit = testFit;
    201           minFit.row = row;
    202         }
    203         if (fabs(SEDtable[row][0].color - color) > 0.5) done = TRUE;
    204         row --;
    205       }
     146    // we now have the min chisq row. use this to supply the other filter values....
     147    for (j = 0; valid && (j < table[0].Nfilter); j++) {
     148      if (table[0].mode[j] != SED_MODEL) continue;
     149      outcat[0].measure[Nmeas].dR_PS       = 0.0;
     150      outcat[0].measure[Nmeas].dD_PS       = 0.0;
     151      outcat[0].measure[Nmeas].M_PS        = MIN (table[0].row[minFit.row][0].mags[j] + minFit.Md,  NO_MAG);
     152      outcat[0].measure[Nmeas].dM_PS       = 0.0;
     153      outcat[0].measure[Nmeas].Mcal_PS     = 0;
     154      outcat[0].measure[Nmeas].t           = TIMEREF;
     155      outcat[0].measure[Nmeas].averef      = Nave;
     156      outcat[0].measure[Nmeas].source      = table[0].code[j];
     157      outcat[0].measure[Nmeas].dophot      = 0;
     158      outcat[0].measure[Nmeas].flags       = 0;
     159      outcat[0].measure[Nmeas].dt_PS       = 0xffff;
    206160
    207       // search for min chisq forwards
    208       done = FALSE;
    209       row = start + 1;
    210       while (!done && (row < Nrow)) {
    211         testFit = SEDchisq (SEDtable[row], &sourceValue, &sourceError, Nfilter);
    212         if (testFit.chisq < minFit.chisq) {
    213           minFit = testFit;
    214           minFit.row = row;
    215         }
    216         if (fabs(SEDtable[row][0].color - color) > 0.5) done = TRUE;
    217         row ++;
    218       }
     161      outcat[0].measure[Nmeas].Mgal_PS     = NO_MAG;
     162      outcat[0].measure[Nmeas].airmass_PS  = 0;
     163      outcat[0].measure[Nmeas].FWx         = NO_MAG;
     164      outcat[0].measure[Nmeas].FWy         = NO_ERR;
     165      outcat[0].measure[Nmeas].theta       = NO_ERR;
     166      Nmeas ++;
     167      CHECK_REALLOCATE (outcat[0].measure, Measure, NMEAS, Nmeas, 100);
     168    }
    219169
    220       Nfit ++;
    221       // create the vectors for the example plots
    222       if (PLOT) {
    223         double tmp;
    224         // find plot range
    225         SetLimitsRaw (NULL, SEDtable[minFit.row][0].mags, Nfilter, &graphdata);
    226         SWAP (graphdata.ymin, graphdata.ymax);
    227 
    228         KapaClear (fd, TRUE);
    229         KapaSetSection (fd, &magSection);
    230         KapaSetLimits (fd, &graphdata);
    231         KapaBox (fd, &graphdata);
    232         graphdata.color = KapaColorByName ("blue");
    233         graphdata.etype = 0;
    234         KapaPrepPlot (fd, Nfilter, &graphdata);
    235         KapaPlotVector (fd, Nfilter, wavecode);
    236         KapaPlotVector (fd, Nfilter, SEDtable[minFit.row][0].mags);
    237         graphdata.color = KapaColorByName ("red");
    238         graphdata.etype = 1;
    239         for (j = 0; j < Nfilter; j++) {
    240           fitmags[j] = 100;
    241           fiterrs[j] = 0;
    242           if (sourceValue.mags[j] > 50) continue;
    243           fitmags[j] = sourceValue.mags[j] - minFit.Md;
    244           fiterrs[j] = sourceError.mags[j];
    245         }
    246         KapaPrepPlot (fd, Nfilter, &graphdata);
    247         KapaPlotVector (fd, Nfilter, wavecode);
    248         KapaPlotVector (fd, Nfilter, fitmags);
    249         KapaPlotVector (fd, Nfilter, fiterrs);
    250         KapaPlotVector (fd, Nfilter, fiterrs);
    251         KapaSendLabel (fd, "model,fit (mags)", 1);
    252 
    253         sprintf (line, "star: %10.6f %10.6f  T: %5.0fK  A_V|: %4.2f  M_D|: %5.2f  &sc&h^2|: %5.2f",
    254                  catalog.average[i].R, catalog.average[i].D,
    255                  SEDtable[minFit.row][0].Temp, SEDtable[minFit.row][0].Av, minFit.Md, minFit.chisq);
    256         KapaSendLabel (fd, line, 2);
    257         KapaSendLabel (fd, "model,fit (mags)", 1);
    258 
    259         KapaSetSection (fd, &resSection);
    260         graphdata.ymin = -1.0;
    261         graphdata.ymax = +1.0;
    262         KapaSetLimits (fd, &graphdata);
    263         KapaBox (fd, &graphdata);
    264         graphdata.color = KapaColorByName ("red");
    265         graphdata.etype = 1;
    266 
    267         for (j = 0; j < Nfilter; j++) {
    268           fitmags[j] = 100;
    269           fiterrs[j] = 0;
    270           if (sourceValue.mags[j] > 50) continue;
    271           fitmags[j] = sourceValue.mags[j] - minFit.Md - SEDtable[minFit.row][0].mags[j];
    272           fiterrs[j] = sourceError.mags[j];
    273         }
    274         KapaPrepPlot (fd, Nfilter, &graphdata);
    275         KapaPlotVector (fd, Nfilter, wavecode);
    276         KapaPlotVector (fd, Nfilter, fitmags);
    277         KapaPlotVector (fd, Nfilter, fiterrs);
    278         KapaPlotVector (fd, Nfilter, fiterrs);
    279         KapaSendLabel (fd, "wavelength (nm)", 0);
    280         KapaSendLabel (fd, "resid (mags)", 1);
    281 
    282         KiiCursorOn (fd);
    283         while (KiiCursorRead (fd, &X, &Y, key)) {
    284           fprintf (stderr, "window: %f %f (%s)\n", X, Y, key);
    285           if (!strcasecmp (key, "Q")) {
    286             KiiCursorOff (fd);
    287             break;
    288           }
    289           if (!strcasecmp (key, "ESCAPE")) {
    290             KiiCursorOff (fd);
    291             unlock_catalog (&catalog);
    292             goto escape;
    293           }
    294         }
    295       }
    296 
    297       // we now have the min chisq row. use this to supply the other filter values....
     170    Nave ++;
     171    if (Nave >= NAVE) {
     172      NAVE += 100;
     173      REALLOCATE (outcat[0].average, Average, NAVE);
     174      REALLOCATE (outcat[0].secfilt, SecFilt, NAVE*Nsec);
    298175    }
    299     unlock_catalog (&catalog);
    300176  }
    301   fprintf (stderr, "fitted %d stars\n", Nfit);
    302   status = TRUE;
    303   goto finish;
     177  outcat[0].Naverage = Nave;
     178  outcat[0].Nmeasure = Nmeas;
    304179 
    305 usage:
    306   fprintf (stderr, "USAGE: fitset (sedtable) : (F) - (F)\n");
    307   goto escape;
    308 
    309 table_missing:
    310   fprintf (stderr, "ERROR: can't open the SED table\n");
    311   goto escape;
    312 
    313 color_missing:
    314   fprintf (stderr, "ERROR: reference color not in SED table\n");
    315   goto escape;
    316 
    317 color_undefined:
    318   fprintf (stderr, "ERROR: undefined photcode in reference color\n");
    319   goto escape;
    320 
    321 code_missing:
    322   fprintf (stderr, "ERROR: undefined photcode in SED table\n");
    323   goto escape;
    324 
    325 escape:
    326   status = FALSE;
    327   goto finish;
    328 
    329 finish:
    330   if (skylist != NULL) SkyListFree (skylist, ((RegionName != NULL) || (RegionList != NULL)));
    331   if (catalog.average != NULL) free (catalog.average);
    332   if (catalog.secfilt != NULL) free (catalog.secfilt);
    333   if (catalog.measure != NULL) free (catalog.measure);
    334 
    335   if (RegionName != NULL) free (RegionName);
    336   if (RegionList != NULL) free (RegionList);
    337   if (wavecode != NULL) free (wavecode);
    338   if (hashcode != NULL) free (hashcode);
    339   if (SEDtableRaw != NULL) {
    340     for (i = 0; i < Nrow; i++) {
    341       free (SEDtableRaw[i].mags);
    342     }
    343     free (SEDtableRaw);
    344   }
    345   if (SEDtable != NULL) free (SEDtable);
    346   if (sourceValue.mags != NULL) free (sourceValue.mags);
    347   if (sourceError.mags != NULL) free (sourceError.mags);
    348 
    349   signal (SIGINT, oldsignal);
    350   return (status);
     180  SEDfitClear ();
    351181}
    352 

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

    r7693 r7696  
    1 # include "addstar.h"
    21# include "sedstar.h"
    32
    … …  
    3534
    3635// find the first table row within 0.1 mag of the requested color (or within 10)
    37 int SEDcolorBracket (SEDtableRow **table, int Ntable, float color) {
     36int SEDcolorBracket (SEDtable *table, float color) {
    3837
    3938  int Nlo, Nhi, N;
    4039  float tcolor;
    4140
    42   N = Nlo = 0; Nhi = Ntable;
    43   tcolor = table[Nlo][0].color;
     41  N = Nlo = 0; Nhi = table[0].Nrow;
     42  tcolor = table[0].row[Nlo][0].color;
    4443  while ((Nhi - Nlo > 10) && (fabs(tcolor-color) > 0.1)) {
    4544    N = 0.5*(Nlo + Nhi);
    4645    N = MAX (N, 0);
    47     N = MIN (N, Ntable - 1);
    48     tcolor = table[N][0].color;
     46    N = MIN (N, table[0].Nrow - 1);
     47    tcolor = table[0].row[N][0].color;
    4948    if (tcolor < color) {
    5049      Nlo = N;
    … …  
    10099}
    101100
     101static Graphdata graphdata;
     102static KapaSection magSection, resSection;
     103static int Xgraph;
     104static float *fitmags, *fiterrs;
     105
     106int SEDfitInit (SEDtable *table) {
     107
     108  Xgraph = KiiOpen ("kapa", "sedstar");
     109  KapaInitGraph (&graphdata);
     110  SetLimitsRaw (table[0].wavecode, NULL, table[0].Nfilter, &graphdata);
     111  graphdata.style = 2;
     112  graphdata.ptype = 2;
     113  KapaClear (Xgraph, TRUE);
     114  magSection.name = strcreate ("mag");
     115  magSection.x  = 0;
     116  magSection.dx = 1;
     117  magSection.y  = 0.5;
     118  magSection.dy = 0.5;
     119  resSection.name = strcreate ("res");
     120  resSection.x  = 0.0;
     121  resSection.dx = 1.0;
     122  resSection.y  = 0.0;
     123  resSection.dy = 0.5;
     124   
     125  KiiResize (Xgraph, 900, 500);
     126  KapaSetFont (Xgraph, "helvetica", 14);
     127  ALLOCATE (fitmags, float, table[0].Nfilter);
     128  ALLOCATE (fiterrs, float, table[0].Nfilter);
     129  return (TRUE);
     130}
     131
     132int SEDfitClear () {
     133
     134  free (fitmags);
     135  free (fiterrs);
     136  KiiClose (Xgraph);
     137  return (TRUE);
     138}
     139
     140int SEDfitPlot (SEDtable *table, double R, double D, SEDfit *minFit, SEDtableRow *sourceValue, SEDtableRow *sourceError) {
     141
     142  int j, minRow, Nfilter;
     143  double tmp, X, Y;
     144  char line[1024], key[20];
     145
     146  minRow = minFit[0].row;
     147  Nfilter = table[0].Nfilter;
     148
     149  // we want to plot the OBSERVED magnitudes
     150  for (j = 0; j < Nfilter; j++) {
     151    fitmags[j] = table[0].row[minRow][0].mags[j] + minFit[0].Md;
     152  }
     153
     154  // find plot range
     155  SetLimitsRaw (NULL, fitmags, Nfilter, &graphdata);
     156  SWAP (graphdata.ymin, graphdata.ymax);
     157
     158  KapaClear (Xgraph, TRUE);
     159  KapaSetSection (Xgraph, &magSection);
     160  KapaSetLimits (Xgraph, &graphdata);
     161  KapaBox (Xgraph, &graphdata);
     162  graphdata.color = KapaColorByName ("blue");
     163  graphdata.etype = 0;
     164  graphdata.ptype = 7;
     165  KapaPrepPlot (Xgraph, Nfilter, &graphdata);
     166  KapaPlotVector (Xgraph, Nfilter, table[0].wavecode);
     167  KapaPlotVector (Xgraph, Nfilter, fitmags);
     168
     169  graphdata.color = KapaColorByName ("red");
     170  graphdata.etype = 1;
     171  graphdata.ptype = 2;
     172  for (j = 0; j < Nfilter; j++) {
     173    fitmags[j] = 100;
     174    fiterrs[j] = 0;
     175    if (sourceValue[0].mags[j] > 50) continue;
     176    fitmags[j] = sourceValue[0].mags[j];
     177    fiterrs[j] = sourceError[0].mags[j];
     178  }
     179  KapaPrepPlot (Xgraph, Nfilter, &graphdata);
     180  KapaPlotVector (Xgraph, Nfilter, table[0].wavecode);
     181  KapaPlotVector (Xgraph, Nfilter, fitmags);
     182  KapaPlotVector (Xgraph, Nfilter, fiterrs);
     183  KapaPlotVector (Xgraph, Nfilter, fiterrs);
     184  KapaSendLabel (Xgraph, "model,fit (mags)", 1);
     185
     186  sprintf (line, "star: %10.6f %10.6f  T: %5.0fK  A_V|: %4.2f  M_D|: %5.2f  &sc&h^2|: %5.2f",
     187           R, D,
     188           table[0].row[minRow][0].Temp,
     189           table[0].row[minRow][0].Av,
     190           minFit[0].Md, minFit[0].chisq);
     191  KapaSendLabel (Xgraph, line, 2);
     192  KapaSendLabel (Xgraph, "model,fit (mags)", 1);
     193
     194  KapaSetSection (Xgraph, &resSection);
     195  graphdata.ymin = -1.0;
     196  graphdata.ymax = +1.0;
     197  KapaSetLimits (Xgraph, &graphdata);
     198  KapaBox (Xgraph, &graphdata);
     199  graphdata.color = KapaColorByName ("red");
     200  graphdata.etype = 1;
     201
     202  for (j = 0; j < Nfilter; j++) {
     203    fitmags[j] = 100;
     204    fiterrs[j] = 0;
     205    if (sourceValue[0].mags[j] > 50) continue;
     206    fitmags[j] = sourceValue[0].mags[j] - minFit[0].Md - table[0].row[minRow][0].mags[j];
     207    fiterrs[j] = sourceError[0].mags[j];
     208  }
     209  KapaPrepPlot (Xgraph, Nfilter, &graphdata);
     210  KapaPlotVector (Xgraph, Nfilter, table[0].wavecode);
     211  KapaPlotVector (Xgraph, Nfilter, fitmags);
     212  KapaPlotVector (Xgraph, Nfilter, fiterrs);
     213  KapaPlotVector (Xgraph, Nfilter, fiterrs);
     214  KapaSendLabel (Xgraph, "wavelength (nm)", 0);
     215  KapaSendLabel (Xgraph, "resid (mags)", 1);
     216
     217  KiiCursorOn (Xgraph);
     218  while (KiiCursorRead (Xgraph, &X, &Y, key)) {
     219    // fprintf (stderr, "window: %f %f (%s)\n", X, Y, key);
     220    if (!strcasecmp (key, "Q")) {
     221      KiiCursorOff (Xgraph);
     222      break;
     223    }
     224    if (!strcasecmp (key, "ESCAPE")) {
     225      KiiCursorOff (Xgraph);
     226      PLOT = FALSE;
     227      return (TRUE);
     228    }
     229    if (!strcasecmp (key, "X")) {
     230      KiiCursorOff (Xgraph);
     231      Shutdown ("quitting sedstar");
     232    }
     233  }
     234  return (TRUE);
     235}
     236
     237void SetLimitsRaw (float *xvec, float *yvec, int Nelements, Graphdata *graphmode) {
     238
     239  double maxX, minX, maxY, minY, range;
     240  int i;
     241
     242  if (xvec != NULL) {
     243    maxX = minX = xvec[0];
     244    for (i = 1; i < Nelements; i++) {
     245      if (!finite(xvec[i])) continue;
     246      maxX = MAX (maxX, xvec[i]);
     247      minX = MIN (minX, xvec[i]);
     248    }
     249    range = maxX - minX;
     250    if (range == 0) range = 0.001 * maxX;
     251    if (range == 0) range = 0.001;
     252    graphmode[0].xmin = minX - 0.05*range;
     253    graphmode[0].xmax = maxX + 0.05*range;
     254  }
     255
     256  if (yvec != NULL) {
     257    maxY = minY = yvec[0];
     258    for (i = 1; i < Nelements; i++) {
     259      if (!finite(yvec[i])) continue;
     260      maxY = MAX (maxY, yvec[i]);
     261      minY = MIN (minY, yvec[i]);
     262    }
     263    range = maxY - minY;
     264    if (range == 0) range = 0.0011 * maxY;
     265    if (range == 0) range = 0.0011;
     266    graphmode[0].ymin = minY - 0.05*range;
     267    graphmode[0].ymax = maxY + 0.05*range;
     268  }
     269}

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

    r7693 r7696  
    1 # include "addstar.h"
    21# include "sedstar.h"
    32
    … …  
    65
    76  FILE *f;
    8   char line[1024], name[64];
     7  char line[1024], name[64], mode[64];
    98  char colornameP[64], colornameM[64];
    10   int colorP, colorM;
     9  int colorP, colorM, code;
    1110  int i, Nrow, NROW;
     11  SEDtable *table;
     12  SEDtableRow *raw;
    1213
    1314  ALLOCATE (table, SEDtable, 1);
    … …  
    1516  // load SED table
    1617  f = fopen (filename, "r");
    17   if (f == NULL)
     18  if (f == NULL) Shutdown ("failure to open SED table");
     19
    1820
    1921  // XXX add error checks for header data
    … …  
    2527  ALLOCATE (table[0].wavecode, float, table[0].Nfilter);
    2628  ALLOCATE (table[0].vegaToAB, float, table[0].Nfilter);
     29  ALLOCATE (table[0].mode,     int,   table[0].Nfilter);
     30  ALLOCATE (table[0].code,     int,   table[0].Nfilter);
    2731
    2832  for (i = 0; i < 0x10000; i++) table[0].hashcode[i] = -1;
    2933  for (i = 0; i < table[0].Nfilter; i++) {
    3034    scan_line (f, line);
    31     sscanf (line, "%*s %s %f %f", name, &table[0].wavecode[i], &table[0].vegaToAB[i]);
     35    sscanf (line, "%*s %s %f %f %s", name, &table[0].wavecode[i], &table[0].vegaToAB[i], mode);
    3236    code = GetPhotcodeCodebyName (name);
    33     if (code == 0) goto code_missing;
     37    table[0].code[i] = code;
     38    if (code == 0) Shutdown ("undefined photcode in SED table");
    3439    table[0].hashcode[code] = i;
     40
     41    table[0].mode[i] = -1;     
     42    if (!strcasecmp(mode, "fit")) table[0].mode[i] = SED_FIT;
     43    if (!strcasecmp(mode, "req")) table[0].mode[i] = SED_REQ;
     44    if (!strcasecmp(mode, "model")) table[0].mode[i] = SED_MODEL;
     45    if (!strcasecmp(mode, "sample")) table[0].mode[i] = SED_SAMPLE;
     46    if (table[0].mode[i] == -1) Shutdown ("invalid photcode mode in SED table");
    3547  }
     48
     49  // load color key
     50  scan_line (f, line);
    3651
    3752  // define the selection color codes
    … …  
    4459  if (table[0].codeM == -1) Shutdown ("missing positive color filter");
    4560   
    46   // skip remaining header lines
    47   scan_line (f, line);
    48   scan_line (f, line);
    49   scan_line (f, line);
    50   scan_line (f, line);
    51  
    5261  // load the SED raw table rows
    5362  Nrow = 0;
    5463  NROW = 100;
    55   ALLOCATE (SEDtableRaw, SEDtableRow, NROW);
     64  ALLOCATE (raw, SEDtableRow, NROW);
    5665  while (scan_line(f, line) != EOF) {
    57     fparse (&SEDtableRaw[Nrow].Temp, 1, line);
    58     fparse (&SEDtableRaw[Nrow].Av, 2, line);
    59     ALLOCATE (SEDtableRaw[Nrow].mags, float, table[0].Nfilter);
     66    stripwhite (line);
     67    if (line[0] == '#') continue;
     68    fparse (&raw[Nrow].Temp, 1, line);
     69    fparse (&raw[Nrow].Av, 2, line);
     70    ALLOCATE (raw[Nrow].mags, float, table[0].Nfilter);
    6071    for (i = 0; i < table[0].Nfilter; i++) {
    61       fparse (&SEDtableRaw[Nrow].mags[i], i + 3, line);
     72      fparse (&raw[Nrow].mags[i], i + 3, line);
    6273    }
    63     SEDtableRaw[Nrow].color = SEDtableRaw[Nrow].mags[codeP] - SEDtableRaw[Nrow].mags[codeM];
     74    raw[Nrow].color = raw[Nrow].mags[table[0].codeP] - raw[Nrow].mags[table[0].codeM];
    6475    Nrow ++;
    65     CHECK_REALLOCATE (SEDtableRaw, SEDtableRow, NROW, Nrow, 100);
     76    CHECK_REALLOCATE (raw, SEDtableRow, NROW, Nrow, 100);
    6677  }     
    6778
    6879  // sort the SEDtable by the reference colors
    69   table[0].row = sort_SEDtable (SEDtableRaw, Nrow);
     80  table[0].row = sort_SEDtable (raw, Nrow);
    7081  table[0].Nrow = Nrow;
    7182  return (table);

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

    r7693 r7696  
    4242  }
    4343
     44  /* extra error messages */
     45  PLOT = FALSE;
     46  if ((N = get_argument (argc, argv, "-plot"))) {
     47    PLOT = TRUE;
     48    remove_argument (N, &argc, argv);
     49  }
     50
    4451  /* other defaults */
    4552  options.timeref = 0;

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

    r7693 r7696  
    1 # include "addstar.h"
     1# include "sedstar.h"
    22
    33int main (int argc, char **argv) {
    44
    5   char *path;
    6   int i, N, Nrefcat;
     5  char *root, *ext;
     6  int i, N, Nrefcat, status;
     7  SkyList *skylist;
    78  SkyTable *sky, *sky2mass;
    89  AddstarClientOptions options;
    910  Catalog incatalog, outcatalog;
     11  SEDtable *sedtable;
    1012
    1113  // need to construct these options with args_load2mass...
    … …  
    2426  for (i = 0; i < skylist[0].Nregions; i++) {
    2527    incatalog.filename = skylist[0].filename[i];
    26     load_pt_catalog (&incatalog, skylist[0].regions[i]);
     28    incatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
     29
     30    // lock and open input catalog
     31    status = lock_catalog (&incatalog, LCK_XCLD);
     32    if (status != 1) {
     33        if (VERBOSE) fprintf (stderr, "skipping empty region\n");
     34        continue;
     35    }
     36    gcatalog (&incatalog);
    2737
    2838    // create output catalog filename
     39    // XXX make the output catalog optional
    2940    root = strstr (incatalog.filename, CATDIR);
    3041    if (root == NULL) Shutdown ("error with input catalog name");
    3142
    32     extension = incatalog.filename + strlen(CATDIR);
    33     ALLOCATE (outcatalog.filename, char, strlen(argv[2]) + strlen(extension) + 2);
    34     sprintf (outcatalog.filename, "%s/%s", argv[2], extension);
     43    ext = incatalog.filename + strlen(CATDIR);
     44    ALLOCATE (outcatalog.filename, char, strlen(argv[2]) + strlen(ext) + 2);
     45    sprintf (outcatalog.filename, "%s/%s", argv[2], ext);
    3546
    36     status = lock_catalog (outcatalog, LCK_XCLD);
     47    if (!check_file_access (outcatalog.filename, TRUE, TRUE)) Shutdown ("can't create output directory");
     48
     49    status = lock_catalog (&outcatalog, LCK_XCLD);
    3750    if (status != 2) Shutdown ("error: output catalog exists");
    38     mkcatalog (region, catalog); /* fills in new header info */
     51    mkcatalog (skylist[0].regions[i], &outcatalog); /* fills in new header info */
    3952
    40     SEDfit (&outcatalog, &incatalog, sedtable);
     53    SEDfitCatalog (&outcatalog, &incatalog, sedtable);
    4154   
    42     wcatalog (outcatalog);
    43     unlock_catalog (outcatalog);
     55    wcatalog (&outcatalog);
     56    unlock_catalog (&outcatalog);
    4457    free (outcatalog.filename);
     58
     59    unlock_catalog (&incatalog);
    4560  }
    4661
    4762  exit (0);
    4863} 
    49 
    5064
    5165/**  sedstar:
Note: See TracChangeset for help on using the changeset viewer.