IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 8, 2007, 12:21:18 PM (19 years ago)
Author:
eugene
Message:

complete merge from branch dvo-mods-2007-02

Location:
trunk/Ohana/src/relastro
Files:
14 added
1 deleted
26 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r11880 r12332  
    1515FULL_CFLAGS   = $(BASE_CFLAGS)
    1616FULL_CPPFLAGS = $(BASE_CPPFLAGS)
    17 FULL_LDFLAGS  = -ldvo -lFITS -lohana $(BASE_LDFLAGS)
     17FULL_LDFLAGS  = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS)
    1818
    1919relastro: $(BIN)/relastro.$(ARCH)
    … …  
    2121
    2222RELASTRO = \
    23 $(SRC)/relastro.$(ARCH).o        \
    24 $(SRC)/initialize.$(ARCH).o      \
    2523$(SRC)/ConfigInit.$(ARCH).o      \
    26 $(SRC)/args.$(ARCH).o            \
    27 $(SRC)/liststats.$(ARCH).o       \
    28 $(SRC)/load_images.$(ARCH).o     \
    29 $(SRC)/name_region.$(ARCH).o     \
    30 $(SRC)/find_images.$(ARCH).o     \
    31 $(SRC)/find_regions.$(ARCH).o    \
    32 $(SRC)/get_regions.$(ARCH).o     \
    33 $(SRC)/load_catalogs.$(ARCH).o   \
    34 $(SRC)/gcatalog.$(ARCH).o        \
    35 $(SRC)/bcatalog.$(ARCH).o        \
    36 $(SRC)/GridOps.v2.$(ARCH).o      \
    37 $(SRC)/StarOps.$(ARCH).o         \
     24$(SRC)/FitChip.$(ARCH).o         \
     25$(SRC)/FitMosaic.$(ARCH).o       \
     26$(SRC)/FitPM.$(ARCH).o           \
     27$(SRC)/FitPMandPar.$(ARCH).o     \
     28$(SRC)/FitSimple.$(ARCH).o       \
    3829$(SRC)/ImageOps.$(ARCH).o        \
    3930$(SRC)/MosaicOps.$(ARCH).o       \
    40 $(SRC)/global_stats.$(ARCH).o    \
     31$(SRC)/ParFactor.$(ARCH).o       \
     32$(SRC)/SetSignals.$(ARCH).o      \
     33$(SRC)/Shutdown.$(ARCH).o        \
     34$(SRC)/UpdateChips.$(ARCH).o     \
     35$(SRC)/UpdateMosaic.$(ARCH).o    \
     36$(SRC)/UpdateObjects.$(ARCH).o   \
     37$(SRC)/UpdateSimple.$(ARCH).o    \
     38$(SRC)/args.$(ARCH).o            \
     39$(SRC)/bcatalog.$(ARCH).o        \
     40$(SRC)/dvo_astrom_ops.$(ARCH).o  \
     41$(SRC)/fitpoly.$(ARCH).o         \
     42$(SRC)/initialize.$(ARCH).o      \
     43$(SRC)/liststats.$(ARCH).o       \
     44$(SRC)/load_catalogs.$(ARCH).o   \
     45$(SRC)/load_images.$(ARCH).o     \
     46$(SRC)/mkpolyterm.$(ARCH).o      \
    4147$(SRC)/plot_scatter.$(ARCH).o    \
    4248$(SRC)/plotstuff.$(ARCH).o       \
    43 $(SRC)/misc.$(ARCH).o            \
     49$(SRC)/select_images.$(ARCH).o   \
     50$(SRC)/sort.$(ARCH).o            \
     51$(SRC)/relastro.$(ARCH).o        \
    4452$(SRC)/reload_catalogs.$(ARCH).o \
    45 $(SRC)/image-db.$(ARCH).o \
    46 $(SRC)/free_catalogs.$(ARCH).o   \
    47 $(SRC)/check_permissions.$(ARCH).o   \
     53$(SRC)/save_catalogs.$(ARCH).o   \
     54$(SRC)/setExclusions.$(ARCH).o   \
     55$(SRC)/write_coords.$(ARCH).o
     56
     57DROP = \
    4858$(SRC)/setMrelFinal.$(ARCH).o    \
    49 $(SRC)/write_coords.$(ARCH).o    \
    50 $(SRC)/setExclusions.$(ARCH).o   \
    51 $(SRC)/Shutdown.$(ARCH).o        \
    52 $(SRC)/SetSignals.$(ARCH).o      \
    53 $(SRC)/wcatalog.$(ARCH).o
     59$(SRC)/StarOps.$(ARCH).o         
     60
    5461
    5562$(RELASTRO): $(INC)/relastro.h

  • trunk/Ohana/src/relastro/doc/parallax.pro

    r6809 r12332  
    332332end
    333333
     334
     335macro parfactor
     336
     337  $RA = $1
     338  $Dec = $2
     339  $e = 23 + 27/60
     340 
     341  create s 0 360
     342  set pR = dcos($e)*d
     343
     344  set pR =  +(dcos($e)*dsin(s)*dcos($RA) - dcos(s)*dsin($RA))
     345  set pD =  -(dcos($e)*dsin(s)*dsin($RA) + dcos(s)*dcos($RA))*dsin($Dec) + dsin($e)*dsin(s)*dcos($Dec)
     346
     347  lim -1.1 1.1 -1.1 1.1; clear; box; plot pR pD
     348end

  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void ConfigInit (int *argc, char **argv) {
    … …  
    55  double ZERO_POINT;
    66  char  *config, *file;
    7   char   PhotCodeFile[256];
     7  char CatdirPhotcodeFile[256];
     8  char MasterPhotcodeFile[256];
    89
    910  /*** load configuration info ***/
    … …  
    3233  ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
    3334  ScanConfig (config, "CATFORMAT",              "%s",  0, CATFORMAT);
     35  ScanConfig (config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
     36
    3437  sprintf (ImageCat, "%s/Images.dat", CATDIR);
    3538
    … …  
    4144  }
    4245
    43   GetConfig (config, "PHOTCODE_FILE",          "%s",  0, PhotCodeFile);
    4446  GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
    4547
    … …  
    5254  if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
    5355
    54   if (!LoadPhotcodes (PhotCodeFile)) {
    55     fprintf (stderr, "error loading photcodes\n");
     56  /* XXX this does not yet write out the master photcode table */
     57  sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR);
     58  if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) {
     59    fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile);
    5660    exit (1);
    5761  }

  • trunk/Ohana/src/relastro/src/ImageOps.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33static unsigned int *start;
    … …  
    1313static int          Nimage;
    1414
     15Image *getimages (int *N) {
     16  *N = Nimage;
     17  return (image);
     18}
     19
     20Image *getimage (int N) {
     21  return (&image[N]);
     22}
     23
    1524void initImages (Image *input, int N) {
    1625
    … …  
    6877}
    6978
    70 /* select all image equivalent to the current photcode */
     79/* match measurements to images */
    7180void findImages (Catalog *catalog, int Ncatalog) {
    7281
    73   int i, j, ecode;
     82  int i, j;
    7483
    7584  for (i = 0; i < Ncatalog; i++) {
    7685    for (j = 0; j < catalog[i].Nmeasure; j++) {
    77       ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
    78       if (photcode[0].code != ecode) continue;
    7986      matchImage (catalog, j, i);
    8087    }
    … …  
    8289}
    8390
     91/* modify this function to use the measure->imageID field */
    8492void matchImage (Catalog *catalog, int meas, int cat) {
    8593
    86   int i, ave, ccdnum;
    87   double ra, dec, X, Y;
    88   char *pname, *filter, *p, base[256];
     94  int i;
    8995  Measure *measure;
    9096 
    9197  measure = &catalog[cat].measure[meas];
    9298  for (i = 0; i < Nimage; i++) {
    93     if (image[0].source == -1) continue;
    94     if (measure[0].source != image[i].source) continue;
     99    if (image[0].photcode == -1) continue;
     100    if (measure[0].photcode != image[i].photcode) continue;
    95101    if (measure[0].t < start[i]) continue;
    96102    if (measure[0].t > stop[i]) continue;
    97103   
    98 # ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
    99     if (USE_GRID) {
    100 
    101       /* identify the ccd on the basis of the photcode name */
    102       pname = GetPhotcodeNamebyCode (image[i].source);
    103       filter = photcode[0].name;
    104       sprintf (base, "%s.%s.", MOSAICNAME, filter);
    105       if (strncmp (pname, base, strlen (base))) continue;
    106       p = pname + strlen(base);
    107       if (*p == 0) continue;
    108       ccdnum = atoi (p);
    109       /* some ambiguity here between seq number and id number */
    110 
    111       /* add this measurement to the grid cell for this chip */
    112       ave = measure[0].averef;
    113       ra  = catalog[cat].average[ave].R - measure[0].dR_PS / 3600.0;
    114       dec = catalog[cat].average[ave].D - measure[0].dD_PS / 3600.0;
    115        
    116       /* X,Y always positive-definite in range 0,0 - dX, dY */
    117       RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    118       setGridMeasure (meas, cat, X, Y, ccdnum);
    119     }
    120 # endif
    121 
    122104    bin[cat][meas] = i;
    123105
    … …  
    133115    return;
    134116  }
    135   /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].source); */
    136 }
    137 
    138 float getMcal (int meas, int cat) {
    139 
    140   int i;
    141   float value;
    142 
    143   i = bin[cat][meas];
    144   if (i == -1) return (NO_MAG);
    145 
    146   if (image[i].code & IMAGE_BAD)  return (NO_MAG); 
    147   value = image[i].Mcal_PS;
    148   return (value);
    149117}
    150118
    … …  
    158126}
    159127
    160 /* determine Mcal values for all images */
    161 void setMcal (Catalog *catalog, int PoorImages) {
    162 
    163   int i, j, m, c, n, N, Nmax, mark, bad;
    164   float Msys, Mrel, Mmos, Mgrid;
    165   double *list, *dlist;
    166   StatType stats;
    167 
    168   if (FREEZE_IMAGES) return;
    169 
    170   if (PoorImages) {
    171     IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    172   }
    173 
    174   Nmax = 0;
    175   for (i = 0; i < Nimage; i++) {
    176     Nmax = MAX (Nmax, Nlist[i]);
    177   }
    178   ALLOCATE (list, double, Nmax);
    179   ALLOCATE (dlist, double, Nmax);
    180 
    181   for (i = 0; i < Nimage; i++) {
    182    
    183     /* on PoorImages run, skip good images */
    184     if (PoorImages) {
    185       bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
    186       if (!bad) continue;
    187     }     
    188 
    189     N = 0;
    190     for (j = 0; j < Nlist[i]; j++) {
    191      
    192       m = mlist[i][j];
    193       c = clist[i][j];
    194      
    195       if (catalog[c].measure[m].flags & MEAS_BAD) continue;
    196       if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
    197       if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
    198       if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
    199      
    200       n = catalog[c].measure[m].averef;
    201       Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
    202       list[N] = Msys - Mrel - Mmos - Mgrid;
    203       dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
    204       N++;
    205     }
    206     /* Nlist[i] is all measurements, N is good measurements */
    207 
    208     /* too few good measurements or too many bad measurements */
    209     if (!PoorImages) {
    210       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
    211       if (mark) {
    212         image[i].code |= ID_IMAGE_FEW;
    213       } else {
    214         image[i].code &= ~ID_IMAGE_FEW;
    215       }     
    216     }
    217    
    218     liststats (list, dlist, N, &stats);
    219     image[i].Mcal_PS  = stats.mean;
    220     image[i].dMcal_PS = stats.sigma;
    221     image[i].Xm    = 100.0*log10(stats.chisq);
    222   }
    223   free (list);
    224   free (dlist);
    225   if (PoorImages) {
    226     IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
    227     STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    228     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
    229   }
    230   return;
    231 }
    232 
    233 /* mark image if: abs(Mcal) too large, dMcal too large */
    234 void clean_images () {
    235 
    236   int i, N, mark, Nmark;
    237   double *mlist, *slist, *dlist;
    238   double MaxOffset, MaxScatter, MedOffset;
    239   StatType stats;
    240 
    241   if (FREEZE_IMAGES) return;
    242 
    243   if (VERBOSE) fprintf (stderr, "marking poor images\n");
    244 
    245   ALLOCATE (mlist, double, Nimage);
    246   ALLOCATE (slist, double, Nimage);
    247   ALLOCATE (dlist, double, Nimage);
    248 
    249   for (i = N = 0; i < Nimage; i++) {
    250     if (image[i].code & IMAGE_BAD) continue;
    251     mlist[N] = fabs (image[i].Mcal_PS);
    252     slist[N] = image[i].dMcal_PS;
    253     dlist[N] = 1;
    254     N++;
    255   }
    256   initstats ("MEAN");
    257   liststats (mlist, dlist, N, &stats);
    258   MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
    259   MedOffset = stats.median;
    260   liststats (slist, dlist, N, &stats);
    261   MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    262   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    263  
    264   Nmark = 0;
    265   for (i = 0; i < Nimage; i++) {
    266     mark = FALSE;
    267     image[i].code &= ~ID_IMAGE_POOR;
    268     mark = (image[i].dMcal_PS > MaxScatter) || (fabs(image[i].Mcal_PS - MedOffset) > MaxOffset);
    269     if (mark) {
    270       Nmark ++;
    271       image[i].code |= ID_IMAGE_POOR;
    272     } else {
    273       image[i].code &= ~ID_IMAGE_POOR;
    274     }
    275   }
    276 
    277   fprintf (stderr, "%d images marked poor\n", Nmark);
    278   initstats (STATMODE);
    279   free (mlist);
    280   free (slist);
    281   free (dlist);
    282 }
    283 
    284128void plot_images () {
    285129
    … …  
    296140  /**** dMcal vs airmass ****/
    297141  for (i = 0; i < Nimage; i++) {
    298     Mlist[i] = image[i].Mcal_PS;
    299     dlist[i] = image[i].dMcal_PS;
    300     xlist[i] = image[i].secz_PS;
     142    Mlist[i] = image[i].Mcal;
     143    dlist[i] = image[i].dMcal;
     144    xlist[i] = image[i].secz;
    301145  }
    302146
    … …  
    316160  bzero (Mlist, NBIN*sizeof(double));
    317161  for (i = 0; i < Nimage; i++) {
    318     bin = image[i].dMcal_PS / 0.00025;
     162    bin = image[i].dMcal / 0.00025;
    319163    bin = MAX (0, MIN (NBIN - 1, bin));
    320164    Mlist[bin] += 1.0;
    … …  
    330174}
    331175
    332 StatType statsImageN (Catalog *catalog) {
    333 
    334   int i, j, m, c, n, N;
    335   double *list, *dlist;
    336   StatType stats;
    337 
    338   bzero (&stats, sizeof (StatType));
    339   if (FREEZE_IMAGES) return (stats);
    340 
    341   ALLOCATE (list, double, Nimage);
    342   ALLOCATE (dlist, double, Nimage);
    343 
    344   n = 0;
    345   for (i = 0; i < Nimage; i++) {
    346     if (image[i].code & IMAGE_BAD)  continue;
    347 
    348     N = 0;
    349     for (j = 0; j < Nlist[i]; j++) {
    350 
    351       m = mlist[i][j];
    352       c = clist[i][j];
    353 
    354       if (getMcal  (m, c) == NO_MAG) continue;
    355       if (getMmos  (m, c) == NO_MAG) continue;
    356       if (getMgrid (m, c) == NO_MAG) continue;
    357       N++;
     176/* XXX re-write to use image->nstars if that is possible */
     177StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
     178
     179  int i, m, c;
     180
     181  Image *mosaic;
     182  Coords *moscoords;
     183  StarData *raw;
     184 
     185  ALLOCATE (raw, StarData, Nlist[im]);
     186
     187  mosaic = NULL;
     188  moscoords = NULL;
     189  if (isMosaic) {
     190    mosaic = getMosaicForImage (im);
     191    moscoords = &mosaic[0].coords;
     192  }
     193
     194  for (i = 0; i < Nlist[im]; i++) {
     195    m = mlist[im][i];
     196    c = clist[im][i];
     197
     198    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     199    raw[i].X = catalog[c].measure[m].Xccd;
     200    raw[i].Y = catalog[c].measure[m].Yccd;
     201   
     202    raw[i].Mag  = catalog[c].measure[m].M;
     203    raw[i].dMag = catalog[c].measure[m].dM;
     204
     205    /* note that for a Simple image, L,M = P,Q */
     206    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     207    if (isMosaic) {
     208      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
     209      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
     210    } else {
     211      raw[i].P = raw[i].L;
     212      raw[i].Q = raw[i].M;
     213      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
    358214    }
    359     list[n] = N;
    360     dlist[n] = 1;
    361     n++;
    362   }
    363 
    364   liststats (list, dlist, n, &stats);
    365   free (list);
    366   free (dlist);
    367   return (stats);
    368 }
    369 
    370 StatType statsImageX (Catalog *catalog) {
    371 
    372   int i, n;
    373   double *list, *dlist;
    374   StatType stats;
    375 
    376   bzero (&stats, sizeof (StatType));
    377   if (FREEZE_IMAGES) return (stats);
    378 
    379   ALLOCATE (list, double, Nimage);
    380   ALLOCATE (dlist, double, Nimage);
    381 
    382   n = 0;
    383   for (i = 0; i < Nimage; i++) {
    384 
    385     if (image[i].code & IMAGE_BAD)  continue;
    386 
    387     list[n] = pow (10.0, 0.01*image[i].Xm);
    388     dlist[n] = 1;
    389     n++;
    390   }
    391 
    392   liststats (list, dlist, n, &stats);
    393   free (list);
    394   free (dlist);
    395   return (stats);
    396 }
    397 
    398 StatType statsImageM (Catalog *catalog) {
    399 
    400   int i, n;
    401   double *list, *dlist;
    402   StatType stats;
    403 
    404   bzero (&stats, sizeof (StatType));
    405   if (FREEZE_IMAGES) return (stats);
    406 
    407   ALLOCATE (list, double, Nimage);
    408   ALLOCATE (dlist, double, Nimage);
    409 
    410   n = 0;
    411   for (i = 0; i < Nimage; i++) {
    412 
    413     if (image[i].code & IMAGE_BAD)  continue;
    414 
    415     list[n] = image[i].Mcal_PS;
    416     dlist[n] = 1;
    417     n++;
    418   }
    419 
    420   liststats (list, dlist, n, &stats);
    421   free (list);
    422   free (dlist);
    423   return (stats);
    424 }
    425 
    426 StatType statsImagedM (Catalog *catalog) {
    427 
    428   int i, n;
    429   double *list, *dlist;
    430   StatType stats;
    431 
    432   bzero (&stats, sizeof (StatType));
    433   if (FREEZE_IMAGES) return (stats);
    434 
    435   ALLOCATE (list, double, Nimage);
    436   ALLOCATE (dlist, double, Nimage);
    437 
    438   n = 0;
    439   for (i = 0; i < Nimage; i++) {
    440 
    441     if (image[i].code & IMAGE_BAD)  continue;
    442 
    443     list[n] = image[i].dMcal_PS;
    444     dlist[n] = 1;
    445     n++;
    446   }
    447 
    448   liststats (list, dlist, n, &stats);
    449   free (list);
    450   free (dlist);
    451   return (stats);
    452 }
    453 
    454 Image *getimages (int *N) {
    455 
    456   *N = Nimage;
    457   return (image);
    458 }
    459 
    460 Image *getimage (int N) {
    461   return (&image[N]);
    462 }
    463 
     215  }
     216 
     217  *Nstars = Nlist[im];
     218  return (raw);
     219}
     220
     221StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
     222
     223  int i, m, c, n;
     224
     225  Image *mosaic;
     226  Coords *moscoords;
     227  StarData *ref;
     228 
     229  ALLOCATE (ref, StarData, Nlist[im]);
     230
     231  mosaic = NULL;
     232  moscoords = NULL;
     233  if (isMosaic) {
     234    mosaic = getMosaicForImage (im);
     235    moscoords = &mosaic[0].coords;
     236  }
     237
     238  for (i = 0; i < Nlist[im]; i++) {
     239    m = mlist[im][i];
     240    c = clist[im][i];
     241    n = catalog[c].measure[m].averef;
     242
     243    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     244    ref[i].R = catalog[c].average[n].R;
     245    ref[i].D = catalog[c].average[n].D;
     246
     247    /* note that for a Simple image, L,M = P,Q */
     248    RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
     249    if (isMosaic) {
     250      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords);
     251      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, moscoords);
     252    } else {
     253      ref[i].L = ref[i].P;
     254      ref[i].M = ref[i].Q;
     255      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
     256    }
     257  }
     258 
     259  *Nstars = Nlist[im];
     260  return (ref);
     261}

  • trunk/Ohana/src/relastro/src/MosaicOps.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33Image *getimages (int *N);
    … …  
    1616
    1717/* find mosaic frames (unique time periods & photcode name matches mosaic) */
     18/* XXX what is a mosaic?  do we need to use the 'DIS/WRP' info to track? */
    1819void initMosaics (Image *image, int Nimage) {
    1920
    … …  
    3637
    3738    /* select valid mosaic images by photcode */
    38     pname = GetPhotcodeNamebyCode (image[i].source);
     39    pname = GetPhotcodeNamebyCode (image[i].photcode);
    3940    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    4041    if (status) continue;
    … …  
    6970    mosaic[Nmosaic].Xm    = 0.0;
    7071    mosaic[Nmosaic].code  = image[i].code;
    71     mosaic[Nmosaic].secz  = image[i].secz_PS;
     72    mosaic[Nmosaic].secz  = image[i].secz;
    7273
    7374    /* add image to mosaic image list */
    … …  
    8687    }
    8788  }
    88 
    89   initMosaicGrid (image, Nimage);
    9089  return;
    91 }
    92 
    93 void initMosaicGrid (Image *image, int Nimage) {
    94 
    95   /* find max dR, dD range for all mosaics */
    96   /* define mosaic.coords to cover dR, dD */
    97   /* send results to initGridBins */
    98 
    99   int i, j, m, NX, NY;
    100   int dXmax, dYmax;
    101   double dS, dX, dY;
    102   double R, D, Rmin, Rmax, Dmin, Dmax;
    103   double Mcal, dMcal, Xm;
    104 
    105   dXmax = dYmax = 0.0;
    106   for (i = 0; i < Nmosaic; i++) {
    107     Dmin = Rmin =  360.0;
    108     Dmax = Rmax = -360.0;
    109     dS = 0.0;
    110     Mcal = dMcal = Xm = 0;
    111     for (j = 0; j < Nimlist[i]; j++) {
    112       m = imlist[i][j];
    113       NX = image[m].NX;
    114       NY = image[m].NY;
    115       dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
    116       XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
    117       Rmin = MIN (Rmin, R);
    118       Rmax = MAX (Rmax, R);
    119       Dmin = MIN (Dmin, D);
    120       Dmax = MAX (Dmax, D);
    121       XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
    122       Rmin = MIN (Rmin, R);
    123       Rmax = MAX (Rmax, R);
    124       Dmin = MIN (Dmin, D);
    125       Dmax = MAX (Dmax, D);
    126       XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
    127       Rmin = MIN (Rmin, R);
    128       Rmax = MAX (Rmax, R);
    129       Dmin = MIN (Dmin, D);
    130       Dmax = MAX (Dmax, D);
    131       XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
    132       Rmin = MIN (Rmin, R);
    133       Rmax = MAX (Rmax, R);
    134       Dmin = MIN (Dmin, D);
    135       Dmax = MAX (Dmax, D);
    136       Mcal += image[m].Mcal_PS;
    137       dMcal += image[m].dMcal_PS;
    138       Xm += image[m].Xm;
    139       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    140       image[m].Mcal_PS = 0.0;
    141       image[m].dMcal_PS = image[m].Xm = NO_MAG;
    142     }
    143     dS /= Nimlist[i];
    144     strcpy (mosaic[i].coords.ctype, "RA---TAN");
    145     mosaic[i].coords.crval1 = Rmin;
    146     mosaic[i].coords.crval2 = Dmin;
    147     mosaic[i].coords.crpix1 = 0.0;
    148     mosaic[i].coords.crpix2 = 0.0;
    149     mosaic[i].coords.cdelt1 = dS;
    150     mosaic[i].coords.cdelt2 = dS;
    151     mosaic[i].coords.pc1_1  = 1.0;
    152     mosaic[i].coords.pc2_2  = 1.0;
    153     mosaic[i].coords.pc1_2  = 0.0;
    154     mosaic[i].coords.pc2_1  = 0.0;
    155     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
    156 
    157     mosaic[i].Mcal  = Mcal / Nimlist[i];
    158     mosaic[i].dMcal = dMcal / Nimlist[i];
    159     mosaic[i].Xm    = Xm / Nimlist[i];
    160   }
    161   if (!USE_GRID) return;
    162 
    163   dXmax = MAX (dXmax, dX);
    164   dYmax = MAX (dYmax, dY);
    165   initGrid (dXmax, dYmax);
    166   return;
    167 }
    168 
    169 void setMcalFinal () {
    170 
    171   int i, j, im, Nimage;
    172   Image *image;
    173 
    174   if (!MOSAICNAME[0]) return;
    175 
    176   image = getimages (&Nimage);
    177 
    178   for (i = 0; i < Nmosaic; i++) {
    179     for (j = 0; j < Nimlist[i]; j++) {
    180       im = imlist[i][j];
    181       image[im].Mcal_PS = mosaic[i].Mcal;
    182       image[im].dMcal_PS = mosaic[i].dMcal;
    183       image[im].Xm = mosaic[i].Xm;
    184       image[im].code |= (mosaic[i].code & ID_IMAGE_FEW);
    185       image[im].code |= (mosaic[i].code & ID_IMAGE_POOR);
    186     }
    187   }     
    18890}
    18991
    … …  
    238140int findMosaics (Catalog *catalog, int Ncatalog) {
    239141 
    240   int i, j, ecode;
     142  int i, j;
    241143
    242144  if (!MOSAICNAME[0]) return (FALSE);
    … …  
    247149        if (catalog[i].measure[j].t > TSTOP) continue;
    248150      }
    249       ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
    250       if (photcode[0].code != ecode) continue;
    251151      matchMosaics (catalog, j, i);
    252152    }
    … …  
    255155}
    256156
     157/* modify this function to use the measure->imageID field ? */
    257158void matchMosaics (Catalog *catalog, int meas, int cat) {
    258159
    … …  
    263164    if (catalog[cat].measure[meas].t > mosaic[i].stop) continue;
    264165   
    265 # ifdef GRID_V1
    266     if (USE_GRID) {
    267       ave = catalog[cat].measure[meas].averef;
    268       ra  = catalog[cat].average[ave].R_PS - catalog[cat].measure[meas].dR_PS / 3600.0;
    269       dec = catalog[cat].average[ave].D_PS - catalog[cat].measure[meas].dD_PS / 3600.0;
    270 
    271       /* X,Y always positive-definite in range 0,0 - dX, dY */
    272       RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords);
    273       setGridMeasure (meas, cat, X, Y);
    274     }
    275 # endif
    276 
    277166    bin[cat][meas] = i;
    278167
    … …  
    288177    return;
    289178  }
    290   fprintf (stderr, "missed measurement\n");
    291179  return;
    292 }
    293 
    294 float getMmos (int meas, int cat) {
    295 
    296   int i;
    297   float value;
    298 
    299   if (!MOSAICNAME[0]) return (0);
    300   i = bin[cat][meas];
    301   if (i == -1) return (NO_MAG);
    302 
    303   if (mosaic[i].code & IMAGE_BAD) return (NO_MAG); 
    304   value = mosaic[i].Mcal;
    305   return (value);
    306 }
    307 
    308 int setMmos (Catalog *catalog, int PoorImages) {
    309 
    310   int i, j, m, c, n, N, Nmax, mark, bad;
    311   float Msys, Mrel, Mcal, Mgrid;
    312   double *list, *dlist, *Mlist, *dMlist;
    313   StatType stats;
    314   Image *image;
    315 
    316   if (!MOSAICNAME[0]) return (FALSE);
    317 
    318   image = getimages (&N);
    319 
    320   if (PoorImages) {
    321     IMAGE_BAD = 0;
    322   }
    323 
    324   Nmax = 0;
    325   for (i = 0; i < Nmosaic; i++) {
    326     Nmax = MAX (Nmax, Nlist[i]);
    327   }
    328   ALLOCATE (list, double, Nmax);
    329   ALLOCATE (dlist, double, Nmax);
    330   ALLOCATE (Mlist, double, Nmax);
    331   ALLOCATE (dMlist, double, Nmax);
    332 
    333   for (i = 0; i < Nmosaic; i++) {
    334    
    335     /* on PoorImages run, skip good images */
    336     if (PoorImages) {
    337       bad = mosaic[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
    338       if (!bad) continue;
    339     }     
    340 
    341     N = 0;
    342     for (j = 0; j < Nlist[i]; j++) {
    343      
    344       m = mlist[i][j];
    345       c = clist[i][j];
    346      
    347       if (catalog[c].measure[m].flags & MEAS_BAD) continue;
    348       if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
    349       if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
    350       if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
    351      
    352       n = catalog[c].measure[m].averef;
    353       Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
    354       list[N]  = Msys - Mrel - Mcal - Mgrid;
    355       dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
    356       Mlist[N] = Msys;
    357       dMlist[N] = list[N];
    358       N++;
    359     }
    360     /* Nlist[i] is all measurements, N is good measurements */
    361 
    362     /* too few good measurements or too many bad measurements (skip in PoorImages run) */
    363 
    364     if (!PoorImages) {
    365       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
    366       if (mark) {
    367         fprintf (stderr, "marked image %s (%d), (%d < %d) || (%d < %f*%d)\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]);
    368         mosaic[i].code |= ID_IMAGE_FEW;
    369       } else {
    370         mosaic[i].code &= ~ID_IMAGE_FEW;
    371       }
    372     }
    373     liststats (list, dlist, N, &stats);
    374     if (PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
    375     mosaic[i].Mcal  = stats.mean;
    376     mosaic[i].dMcal = stats.sigma;
    377     mosaic[i].Xm    = 100.0*log10(stats.chisq);
    378   }
    379   free (list);
    380   free (dlist);
    381   free (Mlist);
    382   free (dMlist);
    383 
    384   if (PoorImages) {
    385     IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
    386     STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    387     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
    388   }
    389   return (TRUE);
    390 }
    391  
    392 StatType statsMosaicM (Catalog *catalog) {
    393 
    394   int i, n;
    395   double *list, *dlist;
    396   StatType stats;
    397 
    398   bzero (&stats, sizeof (StatType));
    399   if (!MOSAICNAME[0]) return (stats);
    400 
    401   ALLOCATE (list, double, Nmosaic);
    402   ALLOCATE (dlist, double, Nmosaic);
    403 
    404   n = 0;
    405   for (i = 0; i < Nmosaic; i++) {
    406     if (mosaic[i].code & IMAGE_BAD) continue;
    407     list[n] = mosaic[i].Mcal;
    408     dlist[n] = 1;
    409     n++;
    410   }
    411 
    412   liststats (list, dlist, n, &stats);
    413   free (list);
    414   free (dlist);
    415   return (stats);
    416 }
    417 
    418 StatType statsMosaicdM (Catalog *catalog) {
    419 
    420   int i, n;
    421   double *list, *dlist;
    422   StatType stats;
    423 
    424   bzero (&stats, sizeof (StatType));
    425   if (!MOSAICNAME[0]) return (stats);
    426 
    427   ALLOCATE (list, double, Nmosaic);
    428   ALLOCATE (dlist, double, Nmosaic);
    429 
    430   n = 0;
    431   for (i = 0; i < Nmosaic; i++) {
    432 
    433     if (mosaic[i].code & IMAGE_BAD) continue;
    434     list[n] = mosaic[i].dMcal;
    435     dlist[n] = 1;
    436     n++;
    437   }
    438 
    439   liststats (list, dlist, n, &stats);
    440   free (list);
    441   free (dlist);
    442   return (stats);
    443 }
    444 
    445 StatType statsMosaicN (Catalog *catalog) {
    446 
    447   int i, j, m, c, n, N;
    448   double *list, *dlist;
    449   StatType stats;
    450 
    451   bzero (&stats, sizeof (StatType));
    452   if (!MOSAICNAME[0]) return (stats);
    453 
    454   ALLOCATE (list, double, Nmosaic);
    455   ALLOCATE (dlist, double, Nmosaic);
    456 
    457   n = 0;
    458   for (i = 0; i < Nmosaic; i++) {
    459     if (mosaic[i].code & IMAGE_BAD)  continue;
    460 
    461     N = 0;
    462     for (j = 0; j < Nlist[i]; j++) {
    463 
    464       m = mlist[i][j];
    465       c = clist[i][j];
    466 
    467       if (getMcal  (m, c) == NO_MAG) continue;
    468       if (getMgrid (m, c) == NO_MAG) continue;
    469       if (getMrel  (catalog, m, c) == NO_MAG) continue;
    470       N++;
    471     }
    472     list[n] = N;
    473     dlist[n] = 1;
    474     n++;
    475   }
    476   fprintf (stderr, "Nmosaic: %d, n: %d\n", Nmosaic, n);
    477 
    478   liststats (list, dlist, n, &stats);
    479   free (list);
    480   free (dlist);
    481   return (stats);
    482 }
    483 
    484 StatType statsMosaicX (Catalog *catalog) {
    485 
    486   int i, n;
    487   double *list, *dlist;
    488   StatType stats;
    489 
    490   bzero (&stats, sizeof (StatType));
    491   if (!MOSAICNAME[0]) return (stats);
    492 
    493   ALLOCATE (list, double, Nmosaic);
    494   ALLOCATE (dlist, double, Nmosaic);
    495 
    496   n = 0;
    497   for (i = 0; i < Nmosaic; i++) {
    498 
    499     if (mosaic[i].code & IMAGE_BAD) continue;
    500     list[n] = pow(10.0, 0.01*mosaic[i].Xm);
    501     dlist[n] = 1;
    502     n++;
    503   }
    504 
    505   liststats (list, dlist, n, &stats);
    506   free (list);
    507   free (dlist);
    508   return (stats);
    509 }
    510 
    511 /* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
    512 void clean_mosaics () {
    513 
    514   int i, N, mark, Nmark;
    515   double *mlist, *slist, *dlist;
    516   double MaxOffset, MedOffset, MaxScatter;
    517   StatType stats;
    518 
    519   if (!MOSAICNAME[0]) return;
    520 
    521   if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
    522 
    523   ALLOCATE (mlist, double, Nmosaic);
    524   ALLOCATE (slist, double, Nmosaic);
    525   ALLOCATE (dlist, double, Nmosaic);
    526 
    527   for (i = N = 0; i < Nmosaic; i++) {
    528     if (mosaic[i].code & IMAGE_BAD) continue;
    529     mlist[N] = mosaic[i].Mcal;
    530     slist[N] = mosaic[i].dMcal;
    531     dlist[N] = 1;
    532     N++;
    533   }
    534   initstats ("MEAN");
    535   liststats (mlist, dlist, N, &stats);
    536   MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
    537   MedOffset = stats.median;
    538   liststats (slist, dlist, N, &stats);
    539   MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    540   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    541  
    542   Nmark = 0;
    543   for (i = 0; i < Nmosaic; i++) {
    544     mark = FALSE;
    545     mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
    546     if (mark) {
    547       Nmark ++;
    548       mosaic[i].code |= ID_IMAGE_POOR;
    549     } else {
    550       mosaic[i].code &= ~ID_IMAGE_POOR;
    551     }
    552   }
    553 
    554   fprintf (stderr, "%d mosaics marked poor\n", Nmark);
    555   initstats (STATMODE);
    556   free (mlist);
    557   free (slist);
    558   free (dlist);
    559180}
    560181
    … …  
    591212
    592213      ave = catalog[c].measure[m].averef;
    593       xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR_PS / 3600.0;
    594       ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD_PS / 3600.0;
     214      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0;
     215      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.0;
    595216      N++;
    596217    }
    … …  
    655276}
    656277
     278Image *getMosaicForImage (int Nim) {
     279
     280    Image *mosaic;
     281
     282    mosaic = NULL;
     283    return mosaic;
     284}
     285
    657286int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage) {
    658287

  • trunk/Ohana/src/relastro/src/SetSignals.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33static int Protect = FALSE;

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

    r7080 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33static FITS_DB *db;
    … …  
    2727}
    2828
     29
     30/* XXX this is probably not needed anymore : just protect the write statements */

  • trunk/Ohana/src/relastro/src/StarOps.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33static int Nmax;
    … …  
    2828  if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG); 
    2929 
    30   value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M_PS;
     30  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
    3131  return (value);
    3232}
    … …  
    5454        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    5555        list[N] = Msys - Mcal - Mmos - Mgrid;
    56         dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     56        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    5757        N++;
    5858      }
    … …  
    6565      liststats (list, dlist, N, &stats);
    6666
    67       if (PhotPrimary) {
    68         catalog[i].average[j].M = stats.mean;
    69         catalog[i].average[j].dM = stats.sigma;
    70         catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    71       } else {
    72         catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
    73         catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = stats.sigma;
    74         catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    75       }     
     67      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
     68      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
     69      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    7670    }
    7771  }
    … …  
    112106        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    113107        list[N] = Msys - Mcal - Mmos - Mgrid;
    114         dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     108        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    115109        N++;
    116110      }
    … …  
    121115
    122116      /* use sigma or error in dM for output? */
    123       if (PhotPrimary) {
    124         catalog[i].average[j].M = stats.mean;
    125         catalog[i].average[j].dM = MAX (stats.sigma, stats.error);
    126         catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    127       } else {
    128         catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
    129         catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = MAX (stats.error, stats.sigma);
    130         catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    131       }     
     117      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
     118      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
     119      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
    132120    }
    133121  }
    … …  
    156144        if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
    157145        if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
    158         catalog[i].measure[m].Mcal_PS = Mcal + Mmos + Mgrid;
     146        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
    159147      }
    160148    }
    … …  
    183171    for (j = 0; j < catalog[i].Naverage; j++) {
    184172      if (catalog[i].average[j].code & STAR_BAD) continue;
    185       Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
     173      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
    186174      if (Xm == -1) continue;
    187175      Chisq = pow (10.0, 0.01*Xm);
    188176      xlist[Ntot] = Chisq;
    189       slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
     177      slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
    190178      dlist[Ntot] = 1;
    191179      Ntot ++;
    … …  
    203191  for (i = 0; i < Ncatalog; i++) {
    204192    for (j = 0; j < catalog[i].Naverage; j++) {
    205       dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
    206       Xm = PhotPrimary ? catalog[i].average[j].Xm    : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
     193      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
     194      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
    207195      Chisq = pow (10.0, 0.01*Xm);
    208196      mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq);
    … …  
    269257        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    270258        list[N] = Msys - Mcal - Mmos - Mgrid;
    271         dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     259        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    272260        N++;
    273261      }
    … …  
    301289        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    302290        list[N] = Msys - Mcal - Mmos - Mgrid;
    303         dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     291        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    304292        ilist[N] = m;
    305293        N++;
    … …  
    387375      if (catalog[i].average[j].code & STAR_BAD) continue; 
    388376
    389       Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
     377      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
    390378      if (Xm == NO_MAG) continue;
    391379      list[n] = pow (10.0, 0.01*Xm);
    … …  
    423411      if (catalog[i].average[j].code & STAR_BAD) continue; 
    424412
    425       dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
     413      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
    426414      list[n] = dM;
    427415      dlist[n] = 1;
    … …  
    452440    for (j = 0; j < catalog[i].Naverage; j++) {
    453441      if (catalog[i].average[j].code & STAR_BAD) continue; 
    454       dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
     442      dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
    455443      bin = dMrel / 0.00025;
    456444      bin = MAX (0, MIN (NBIN-1, bin));
    … …  
    483471    for (j = 0; j < catalog[i].Naverage; j++) {
    484472      if (catalog[i].average[j].code & STAR_BAD) continue;
    485       xlist[N] = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
    486       value    = PhotPrimary ? catalog[i].average[j].Xm   : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
     473      xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
     474      value    = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
    487475      if (value == NO_MAG) continue;
    488476      ylist[N] = 0.01*value;

  • trunk/Ohana/src/relastro/src/UpdateImages.c

    r11742 r12332  
    4242        /* the astrometric errors are not being carried yet (but should be!) */
    4343        /* we use the photometric mag error as a weighting term */
    44         dlistR[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
    45         dlistD[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     44        dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     45        dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    4646        N++;
    4747      }

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r11742 r12332  
    22
    33static int Nmax;
    4 static double *listR, *listD;
    5 static double *dlistR, *dlistD;
     4static double *X, *dX;
     5static double *Y, *dY;
     6static double *R, *dR;
     7static double *D, *dD;
     8static double *pX;
     9static double *pY;
     10static time_t *T;
     11static double *dT;
    612
    713void initObjectData (Catalog *catalog, int Ncatalog) {
    … …  
    1622  }
    1723
    18   ALLOCATE (list, double, MAX (1, Nmax));
    19   ALLOCATE (dlist, double, MAX (1, Nmax));
     24  ALLOCATE (R, double, MAX (1, Nmax));
     25  ALLOCATE (D, double, MAX (1, Nmax));
     26  ALLOCATE (T, time_t, MAX (1, Nmax));
     27  ALLOCATE (X, double, MAX (1, Nmax));
     28  ALLOCATE (Y, double, MAX (1, Nmax));
     29
     30  ALLOCATE (dR, double, MAX (1, Nmax));
     31  ALLOCATE (dD, double, MAX (1, Nmax));
     32  ALLOCATE (dT, double, MAX (1, Nmax));
     33  ALLOCATE (dX, double, MAX (1, Nmax));
     34  ALLOCATE (dY, double, MAX (1, Nmax));
     35
     36  ALLOCATE (pX, double, MAX (1, Nmax));
     37  ALLOCATE (pY, double, MAX (1, Nmax));
    2038} 
    2139
    … …  
    2341
    2442  int i, j, k, m, N;
    25   float chisq;
    2643  StatType statsR, statsD;
     44  Coords coords;
     45  PMFit fit;
     46
     47  initObjectData (catalog, Ncatalog);
     48
     49  coords.crval1 = 0;
     50  coords.crval2 = 0;
     51  coords.crpix1 = 0;
     52  coords.crpix2 = 0;
     53  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
     54  coords.pc1_1  = coords.pc2_2 = 1.0;
     55  coords.pc1_2  = coords.pc2_1 = 0.0;
     56  coords.Npolyterms = 1;
     57  strcpy (coords.ctype, "RA---SIN");
    2758
    2859  for (i = 0; i < Ncatalog; i++) {
    … …  
    3162      /* calculate the average value of R,D for a single star */
    3263      if (catalog[i].average[j].code & STAR_BAD) continue; 
    33       m = catalog[i].average[j].offset;
    3464
    3565      N = 0;
     66      m = catalog[i].average[j].offset;
    3667      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
    3768        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    3869       
    39         listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    40         listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     70        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     71        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     72        T[N] = catalog[i].measure[m].t;
    4173
    4274        /* the astrometric errors are not being carried yet (but should be!) */
    4375        /* we use the photometric mag error as a weighting term */
    44         dlistR[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
    45         dlistD[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
     76        dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     77        dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     78        dT[N] = catalog[i].measure[m].dt;
     79
    4680        N++;
    4781      }
     82
     83      // XXX This criterion needs to be better considered: adjust to match Ndof
    4884      if (N < STAR_TOOFEW) { /* too few measurements */
    4985        catalog[i].average[j].code |= ID_STAR_FEW;
     86        continue;
    5087      } else {
    5188        catalog[i].average[j].code &= ~ID_STAR_FEW;
    5289      }
    5390
    54       /* in here, we should be fitting the parallax and proper-motion components */
     91      /* we need to do the fit in a locally linear space; choose a ref coordinate */
     92      coords.crval1 = R[0];
     93      coords.crval2 = D[0];
     94     
     95      /* project all of the R,D coordinates to a plane centered on this coordinate */
     96      for (k = 0; k < N; k++) {
     97          RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
     98      }   
    5599
    56       liststats (listR, dlistR, N, &statsR);
    57       liststats (listR, dlistR, N, &statsD);
     100      /* fit the model components as needed */
     101      switch (FIT_MODE) {
     102        case FIT_AVERAGE:
     103          liststats (R, dR, N, &statsR);
     104          liststats (D, dD, N, &statsD);
    58105
    59       catalog[i].average[j].R = statsR.mean;
    60       catalog[i].average[j].dR = statsR.sigma;
     106          fit.Ro = statsR.mean;
     107          fit.dRo = statsR.sigma;
    61108
    62       catalog[i].average[j].D = statsD.mean;
    63       catalog[i].average[j].dD = statsD.sigma;
     109          fit.Do = statsD.mean;
     110          fit.dDo = statsD.sigma;
    64111
    65       chisq = 0.5*(statsR.chisq + statsD.chisq);
    66       catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
     112          fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
     113          fit.Nfit = N;
     114          break;
     115
     116        case FIT_PM_ONLY:
     117          FitPM (&fit, X, dX, Y, dY, T, N);
     118          break;
     119
     120        case FIT_PM_AND_PAR:
     121          for (k = 0; k < N; k++) {
     122            ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
     123          }
     124          FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
     125          break;
     126        default:
     127          fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__);
     128          exit (2);
     129      }   
     130
     131      if (0 && (j < 100)) {
     132          fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
     133                   catalog[i].average[j].R,
     134                   catalog[i].average[j].D,
     135                   fit.Ro, fit.Do,
     136                   3600*(catalog[i].average[j].R - fit.Ro),
     137                   3600*(catalog[i].average[j].D - fit.Do));
     138      }
     139
     140      // the measure fields must be updated before the average fields
     141      m = catalog[i].average[j].offset;
     142      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     143        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
     144        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     145        setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     146      }     
     147
     148      catalog[i].average[j].R   = fit.Ro;
     149      catalog[i].average[j].D   = fit.Do;
     150      catalog[i].average[j].dR  = fit.dRo;
     151      catalog[i].average[j].dD  = fit.dDo;
     152
     153      catalog[i].average[j].uR  = fit.uR;
     154      catalog[i].average[j].uD  = fit.uD;
     155      catalog[i].average[j].duR = fit.duR;
     156      catalog[i].average[j].duD = fit.duD;
     157
     158      catalog[i].average[j].P   = fit.p;
     159      catalog[i].average[j].dP  = fit.dp;
     160
     161      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
    67162
    68163    }
    … …  
    71166}
    72167
     168/* fitting proper-motion and parallax:
     169
     170 given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
     171
     172 let x,y be the coordinate in the linearized frame with y parallel to DEC lines
     173
     174 L,B are the ecliptic longitude and latitude of the object,
     175 dL and dB are the offsets in the L and B directions
     176
     177 dL = sin(t - topp)
     178 dB = cos(t - topp)*sin(B)
     179
     180 these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
     181 for the source in the x,y frame is:
     182
     183 x = Ro + uR * (t - to) + p * pR
     184 y = Do + uD * (t - to) + p * pD
     185
     186 the unknowns in these equations are Ro, uR, Do, uD, and p
     187
     188 XXX think through the concepts for the pole a bit better.  all objects near the pole
     189 move the same way with the same phase.  choose a projection center and define dL,dB relative
     190 to that center point coordinate system?
     191
     192*/

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

    r11742 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22void usage (void);
    33
    … …  
    66  int N;
    77  double trange;
     8
     9  /* possible operations */
     10  FIT_TARGET = TARGET_NONE;
     11  FIT_MODE = FIT_AVERAGE;
     12  if ((N = get_argument (argc, argv, "-update-objects"))) {
     13    remove_argument (N, &argc, argv);
     14    FIT_TARGET = TARGET_OBJECTS;
     15
     16    // check for object fitting modes (not valid for images)
     17    if ((N = get_argument (argc, argv, "-pm"))) {
     18        remove_argument (N, &argc, argv);
     19        FIT_MODE = FIT_PM_ONLY;
     20    }
     21    if ((N = get_argument (argc, argv, "-pmpar"))) {
     22        remove_argument (N, &argc, argv);
     23        FIT_MODE = FIT_PM_AND_PAR;
     24    }
     25  }
     26  if ((N = get_argument (argc, argv, "-update-simple"))) {
     27    remove_argument (N, &argc, argv);
     28    FIT_TARGET = TARGET_SIMPLE;
     29  }
     30  if ((N = get_argument (argc, argv, "-update-chips"))) {
     31    remove_argument (N, &argc, argv);
     32    FIT_TARGET = TARGET_CHIPS;
     33  }
     34  if ((N = get_argument (argc, argv, "-update-mosaics"))) {
     35    remove_argument (N, &argc, argv);
     36    FIT_TARGET = TARGET_MOSAICS;
     37  }
     38  if (FIT_TARGET == TARGET_NONE) usage();
     39
     40  /* specify portion of the sky : allow default of all sky? */
     41  UserPatch.Rmin = 0;
     42  UserPatch.Rmax = 360;
     43  UserPatch.Dmin = -90;
     44  UserPatch.Dmax = +90;
     45  if ((N = get_argument (argc, argv, "-region"))) {
     46    remove_argument (N, &argc, argv);
     47    UserPatch.Rmin = atof (argv[N]);
     48    remove_argument (N, &argc, argv);
     49    UserPatch.Rmax = atof (argv[N]);
     50    remove_argument (N, &argc, argv);
     51    UserPatch.Dmin = atof (argv[N]);
     52    remove_argument (N, &argc, argv);
     53    UserPatch.Dmax = atof (argv[N]);
     54    remove_argument (N, &argc, argv);
     55  } else {
     56    usage ();
     57  }
    858
    959  /* define time */
    … …  
    3484  }
    3585
    36   /* specify portion of the sky */
    37   UserPatch.Rmin = 0;
    38   UserPatch.Rmax = 360;
    39   UserPatch.Dmin = -90;
    40   UserPatch.Dmax = +90;
    41   UserPatchSelect = FALSE;
    42   if ((N = get_argument (argc, argv, "-region"))) {
    43     remove_argument (N, &argc, argv);
    44     UserPatch.Rmin = atof (argv[N]);
    45     remove_argument (N, &argc, argv);
    46     UserPatch.Rmax = atof (argv[N]);
    47     remove_argument (N, &argc, argv);
    48     UserPatch.Dmin = atof (argv[N]);
    49     remove_argument (N, &argc, argv);
    50     UserPatch.Dmax = atof (argv[N]);
    51     remove_argument (N, &argc, argv);
    52     UserPatchSelect = TRUE;
     86  PHOTCODE_LIST = NULL;
     87  if ((N = get_argument (argc, argv, "-photcode"))) {
     88    remove_argument (N, &argc, argv);
     89    PHOTCODE_LIST = strcreate(argv[N]);
     90    remove_argument (N, &argc, argv);
    5391  }
    5492
    … …  
    79117  }
    80118
    81   NLOOP = 8;
    82   if ((N = get_argument (argc, argv, "-n"))) {
    83     remove_argument (N, &argc, argv);
    84     NLOOP = atof (argv[N]);
    85     remove_argument (N, &argc, argv);
    86   }
    87 
    88119  RESET = FALSE;
    89120  if ((N = get_argument (argc, argv, "-reset"))) {
    … …  
    104135  }
    105136
     137  /* XXX update these for relevant plots */
    106138  PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0;
    107139  if ((N = get_argument (argc, argv, "-plrange"))) {
    … …  
    117149  }
    118150
    119   /* group images by mosaic, find Mmos */
    120   MOSAICNAME[0] = 0;
    121   if ((N = get_argument (argc, argv, "-mosaic"))) {
    122     remove_argument (N, &argc, argv);
    123     strcpy (MOSAICNAME, argv[N]);
    124     remove_argument (N, &argc, argv);
    125   }
    126 
    127   FREEZE_IMAGES = FALSE;
    128   if ((N = get_argument (argc, argv, "-imfreeze"))) {
    129     remove_argument (N, &argc, argv);
    130     FREEZE_IMAGES = TRUE;
    131   }
    132 
    133   USE_GRID = FALSE;
    134   if ((N = get_argument (argc, argv, "-grid"))) {
    135     remove_argument (N, &argc, argv);
    136     USE_GRID = TRUE;
    137     if (!MOSAICNAME[0]) {
    138       fprintf (stderr, "-grid is only valid with -mosaic\n");
    139       exit (2);
    140     }
    141   }
    142 
     151  /* XXX update this */
    143152  MIN_ERROR = 0.001;
    144153  if ((N = get_argument (argc, argv, "-minerror"))) {
    … …  
    172181    ImagSelect = TRUE;
    173182  }
    174 
    175   /* possible operations */
    176   UpdateObjects = FALSE;
    177   if ((N = get_argument (argc, argv, "-update-objects"))) {
    178     remove_argument (N, &argc, argv);
    179     UpdateObjects = TRUE;
    180   }
    181   UpdateImages = FALSE;
    182   if ((N = get_argument (argc, argv, "-update-images"))) {
    183     remove_argument (N, &argc, argv);
    184     UpdateImages = TRUE;
    185   }
    186 
    187   if ( UserPatchSelect && (argc != 2)) usage ();
    188   if (!UserPatchSelect && (argc != 3)) usage ();
    189 
     183 
     184  /* XXX drop this? */
     185  DophotSelect = FALSE;
     186  if ((N = get_argument (argc, argv, "-dophot"))) {
     187    remove_argument (N, &argc, argv);
     188    DophotValue = atof (argv[N]);
     189    remove_argument (N, &argc, argv);
     190    DophotSelect = TRUE;
     191  }
     192
     193  if (argc != 1) usage ();
    190194  return TRUE;
    191195}
    192196
    193197void usage () {
    194   fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
    195   fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
    196   fprintf (stderr, "  options: \n");
     198  fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
     199  fprintf (stderr, "  working options: \n");
    197200  fprintf (stderr, "  -time (start) (stop)\n");
    198201  fprintf (stderr, "  -v\n");
    199202  fprintf (stderr, "  -plot\n");
    200203  fprintf (stderr, "  -plotdelay (seconds)\n");
    201   fprintf (stderr, "  -statmode (mode)\n");
    202   fprintf (stderr, "  -n (nloop)\n");
    203   fprintf (stderr, "  -reset\n");
    204204  fprintf (stderr, "  -update\n");
    205205  fprintf (stderr, "  -params\n");
    206   fprintf (stderr, "  -mosaic (mosaic)\n");
    207   fprintf (stderr, "  -imfreeze\n");
    208   fprintf (stderr, "  -grid\n");
     206  fprintf (stderr, "  -reset\n");
    209207  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
    210   fprintf (stderr, "  -instmag min max\n");
     208  fprintf (stderr, "  -instmag min max\n\n");
     209  fprintf (stderr, "  planned options: \n");
     210  fprintf (stderr, "  -photcode code[,code,code,..]\n");
    211211  fprintf (stderr, "  \n");
    212212  exit (2);

  • trunk/Ohana/src/relastro/src/bcatalog.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33int bcatalog (Catalog *subcatalog, Catalog *catalog) {
    44 
    5   int i, j, offset, ecode;
     5  int i, j, k, offset, found;
    66  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    77  float mag;
    … …  
    1919    if (catalog[0].average[i].Nm < 2) continue;
    2020
    21     /* XXX this limitation is absurd
    22     if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
    23     if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
    24     if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
    25     if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
    26     */
    27 
    2821    /* start with all stars good */
    2922    subcatalog[0].average[Naverage] = catalog[0].average[i];
    … …  
    3427
    3528    if (RESET) {
    36       float *p;
    37       p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M_PS;
    38       *p = NO_MAG;
    39       p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM_PS;
    40       *p = NO_MAG;
     29      // XXX reset the ra,dec coords?
    4130      subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW;
    4231      subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR;
    … …  
    4837      offset = catalog[0].average[i].offset + j;
    4938
    50       /* select measurements by photcode */
    51       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
    52       if (ecode != photcode[0].code) continue;
     39      /* select measurements by photcode, if specified */
     40      if (Nphotcodes > 0) {
     41        found = FALSE;
     42        for (k = 0; (k < Nphotcodes) && !found; k++) {
     43          if (photcodes[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
     44        }
     45        if (!found) continue;
     46      }
    5347
    5448      /* select measurements by time */
    … …  
    5953
    6054      /* select measurements by quality */
    61       if (catalog[0].measure[offset].dophot != 1) continue;
     55      if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
    6256
    6357      /* select measurements by mag limit */
    … …  
    6660
    6761      /* select measurements by measurement error */
    68       if (catalog[0].measure[offset].dM_PS > SIGMA_LIM) continue;
     62      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue;
    6963
    7064      /* select measurements by mag limit */
    … …  
    7973      subcatalog[0].measure[Nmeasure].averef = Naverage;
    8074      if (RESET) {
    81         subcatalog[0].measure[Nmeasure].Mcal_PS = 0;
     75        subcatalog[0].measure[Nmeasure].Mcal = 0;
    8276        subcatalog[0].measure[Nmeasure].flags &= 0xff00;
    8377        subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_POOR;
    … …  
    111105  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
    112106
    113   if (catalog[0].Naverage > 0) free (catalog[0].average);
    114   if (catalog[0].Nmeasure > 0) free (catalog[0].measure);
    115   if (catalog[0].Naverage > 0) free (catalog[0].secfilt);
    116 
    117107  if (VERBOSE) {
    118108    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,

  • trunk/Ohana/src/relastro/src/global_stats.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void global_stats (Catalog *catalog, int Ncatalog) {

  • trunk/Ohana/src/relastro/src/initialize.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void initialize (int argc, char **argv) {
    44
    5   int N;
     5  int NPHOTCODES;
     6  char *codename, *ptr;
    67
    78  ConfigInit (&argc, argv);
    89  args (argc, argv);
    910
    10   N = UserPatchSelect ? 1 : 2;
    11   if ((photcode = GetPhotcodebyName (argv[N])) == NULL) {
    12     fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]);
    13     exit (1);
     11  /* XXX we need to build a list of accepted photcodes. these will be used by bcatalog to accept or
     12     reject loaded data */
     13
     14  Nphotcodes = 0;
     15  photcodes = NULL;
     16  if (PHOTCODE_LIST != NULL) {
     17    NPHOTCODES = 10;
     18    ALLOCATE (photcodes, PhotCode *, NPHOTCODES);
     19
     20    /* parse the comma-separated list of photcodes */
     21    strtok_r (PHOTCODE_LIST, ",", &ptr);
     22    while ((codename = strtok_r (NULL, ",", &ptr)) != NULL) {
     23      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_LIST);
     24      fprintf (stderr, "codename: %s\n", codename);
     25      if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) {
     26        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
     27        exit (1);
     28      }
     29      Nphotcodes ++;
     30      CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10);
     31    }
    1432  }
    15   if ((photcode[0].type != PHOT_PRI) && (photcode[0].type != PHOT_SEC)) {
    16     fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[N]);
    17     exit (1);
    18   }
    19 
    20   PhotPrimary = (photcode[0].type == PHOT_PRI);
    2133  PhotNsec = GetPhotcodeNsecfilt ();
    22   PhotSec = GetPhotcodeNsec (photcode[0].code);
    2334
    2435  initstats (STATMODE);
    … …  
    2839  MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
    2940
     41  /* XXX drop irrelevant entries */
    3042  if (SHOW_PARAMS) {
    3143    fprintf (stderr, "current parameter settings:\n");

  • trunk/Ohana/src/relastro/src/liststats.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,

  • trunk/Ohana/src/relastro/src/load_catalogs.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    3 Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) {
     3Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
    44
    55  int i, Nstar;
    6   Catalog *catalog, tcatalog;
     6  Catalog *catalog, *pcatalog, tcatalog;
    77
    88  if (VERBOSE) fprintf (stderr, "loading catalog data\n");
    … …  
    1313  for (i = 0; i < skylist[0].Nregions; i++) {
    1414
     15    pcatalog = subselect ? &tcatalog : &catalog[i];
     16
    1517    // set up the basic catalog info
    16     tcatalog.filename = skylist[0].filename[i];
    17     tcatalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    18     tcatalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    19     tcatalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
     18    pcatalog[0].filename = skylist[0].filename[i];
     19    pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     20    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
     21    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
    2022
    21     if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
    22       fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename);
     23    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
     24      fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
    2325      exit (1);
    2426    }
    25     if (VERBOSE && !tcatalog.Nave_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
     27    if (VERBOSE && !pcatalog[0].Nave_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
    2628
    2729    // select only the brighter stars
    28     bcatalog (&catalog[i], &tcatalog);
    29     dvo_catalog_unlock (&tcatalog);
    30     dvo_catalog_free (&tcatalog);
     30    if (subselect) {
     31      bcatalog (&catalog[i], &tcatalog);
     32      dvo_catalog_unlock (&tcatalog);
     33      dvo_catalog_free (&tcatalog);
     34    }
    3135  }
    3236
    … …  
    3741  }
    3842  if (Nstar < 2) {
    39     fprintf (stderr, "insufficient stars %d\n", Nstar);
    40     exit (0);
     43    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
    4144  }
    4245

  • trunk/Ohana/src/relastro/src/load_images.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    3 SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
     3SkyList *load_images (FITS_DB *db, SkyRegion *region) {
    44
    55  Image     *image, *subset;
    6   int        Nimage, Nsubset, Nchar;
     6  int        Nimage, Nsubset;
    77  int       *LineNumber;
    88
    … …  
    1515 
    1616  // determine the populated SkyRegions overlapping the requested area
    17   if (RegionSelect) {
    18     skylist = SkyListByPatch (sky, -1, region);
    19   } else {
    20     Nchar = strlen(regionName);
    21     if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
    22     skylist = SkyListByName (sky, regionName);
    23   }
     17  skylist = SkyListByPatch (sky, -1, region);
    2418
    2519  // convert database table to internal structure
    … …  
    3428  initMosaics (subset, Nsubset);
    3529 
     30  /* unlock, if we can (else, unlocked below) */
     31  if (!UPDATE) dvo_image_unlock (db);
     32
    3633  return (skylist);
    3734}

  • trunk/Ohana/src/relastro/src/plot_scatter.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22   
    33void plot_scatter (Catalog *catalog, int Ncatalog) {
    44
     5# if (0)   
    56  int i, j, k, m, N, Ntot;
    67  float Mrel, Mcal, Mmos, Mgrid;
    … …  
    2829      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
    2930        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    30         if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
    31         if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
    32         if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
    3331
    34         Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
     32        Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
    3533        xlist[N] = Mrel;
    3634        ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;
    … …  
    5553  free (ylist);
    5654  free (ilist);
     55
     56# endif
    5757}
     58
     59/* XXX this should become astrometrically relevant */

  • trunk/Ohana/src/relastro/src/plotstuff.c

    r10843 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22# include <signal.h>
    33
    … …  
    4949  if (Xgraph[N] == 0) return;
    5050
    51   KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename);
     51  KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename, "default");
    5252  return;
    5353}

  • trunk/Ohana/src/relastro/src/relastro.c

    r11742 r12332  
    33int main (int argc, char **argv) {
    44
    5   int i, status, Ncatalog;
     5  int status, Ncatalog;
    66  Catalog *catalog;
    77  FITS_DB db;
    … …  
    2222
    2323  /* load regions and images based on specified sky patch */
    24   // XXX need to mimic old-style load by passing patch name
    25   // XXX need to reduce number of global variables in use.
    26   // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
    27   skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
     24  skylist = load_images (&db, &UserPatch);
    2825
    29   /* unlock, if we can (else, unlocked below) */
    30   if (!UPDATE) dvo_image_unlock (&db);
    31 
    32   /* load catalog data from region files */
    33   catalog = load_catalogs (skylist, &Ncatalog);
     26  /* load catalog data from region files : subselect only if we are not doing the objects */
     27  catalog = load_catalogs (skylist, &Ncatalog, (FIT_TARGET != TARGET_OBJECTS));
    3428
    3529  /* match measurements with images, mosaics */
    3630  initImageBins  (catalog, Ncatalog);
    37   initMosaicBins (catalog, Ncatalog);
    38   initGridBins   (catalog, Ncatalog);
    39   initMrel (catalog, Ncatalog);
     31  // initMosaicBins (catalog, Ncatalog);
    4032
    4133  findImages (catalog, Ncatalog);
    42   findMosaics (catalog, Ncatalog);  /* also sets Grid values */
    43   SAVEPLOT = FALSE;
    44 
    45   setExclusions (catalog, Ncatalog);
     34  // findMosaics (catalog, Ncatalog);  /* also sets Grid values */
    4635
    4736  if (PLOTSTUFF) {
    48     plot_star_coords (catalog, Ncatalog);
     37    // plot_star_coords (catalog, Ncatalog);
    4938    plot_mosaic_fields (catalog);
    5039  }
    5140
    5241  /* major modes */
    53   if (UpdateObjects) {
    54     UpdateObjects (catalog, Ncatalog);
    55   }
    56   if (UpdateImages) {
    57     UpdateImages (catalog, Ncatalog);
     42  switch (FIT_TARGET) {
     43    case TARGET_OBJECTS:
     44      UpdateObjects (catalog, Ncatalog);
     45      break;
     46    case TARGET_SIMPLE:
     47      UpdateSimple (catalog, Ncatalog);
     48      break;
     49
     50    case TARGET_CHIPS:
     51      UpdateChips (catalog, Ncatalog);
     52      break;
     53
     54    case TARGET_MOSAICS:
     55      UpdateMosaic (catalog, Ncatalog);
     56      break;
     57    default:
     58      fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__);
     59      exit (2);
    5860  }
    5961
    … …  
    6971  if (!UPDATE) exit (0);
    7072
    71   /* need to figure out how to update images, etc */
    72   dvo_image_update (&db, VERBOSE);
    73   dvo_image_unlock (&db);
     73  if (FIT_TARGET == TARGET_OBJECTS) {
     74    save_catalogs (catalog, Ncatalog);
     75  } else {
     76    dvo_image_update (&db, VERBOSE);
     77    dvo_image_unlock (&db);
     78  }
    7479
    7580  exit (0);

  • trunk/Ohana/src/relastro/src/reload_catalogs.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void reload_catalogs (SkyList *skylist) {
    … …  
    3838    initImageBins  (&catalog, 1);
    3939    initMosaicBins (&catalog, 1);
    40     initGridBins   (&catalog, 1);
    4140
    4241    findImages (&catalog, 1);
    4342    findMosaics (&catalog, 1);
    4443
    45     setMrelFinal (&catalog);
    4644    dvo_catalog_save (&catalog, VERBOSE);
    4745    dvo_catalog_unlock (&catalog);
    … …  
    5048    freeImageBins (1);
    5149    freeMosaicBins (1);
    52     freeGridBins (1);
    5350  }
    5451}

  • trunk/Ohana/src/relastro/src/select_images.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33/* this function returns a list of all images which overlap the given SkyList (set of
    … …  
    1717  Image *image;
    1818  int i, j, k, m, found, nimage, NIMAGE;
    19   int InRange, ecode;
     19  int InRange;
    2020  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
    2121  int *line_number;
    … …  
    7575  for (i = 0; i < Ntimage; i++) {
    7676     
    77     /* exclude images by photcode */
    78     ecode = GetPhotcodeEquivCodebyCode (timage[i].source);
    79     if (ecode != photcode[0].code) continue;
     77    /* exclude images by photcode, if selected */
     78    if (Nphotcodes > 0) {
     79      found = FALSE;
     80      for (j = 0; (j < Nphotcodes) && !found; j++) {
     81        if (photcodes[j][0].code == timage[i].photcode) found = TRUE;
     82      }
     83      if (!found) continue;
     84    }
    8085
    8186    /* exclude images by time */
    … …  
    135140      if (RESET) {
    136141        assignMcal (&image[nimage], (double *) NULL, -1);
    137         image[nimage].dMcal_PS = NO_MAG;
     142        image[nimage].dMcal = NO_MAG;
    138143        image[nimage].code &= ~ID_IMAGE_POOR;
    139144      }

  • trunk/Ohana/src/relastro/src/setExclusions.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
     3/* XXX I think all of these, except for X,Y selection, are done / can be done in bcatalog */
    34int setExclusions (Catalog *catalog, int Ncatalog) {
    45
    5   int i, j, k, m, Narea, Nnocal, ecode;
     6  int i, j, k, m, Narea, Nnocal;
    67  Coords *coords;
    78  double r, d, x, y;
    … …  
    1415
    1516        /* select measurements by photcode */
    16         ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].source);
    17         if (ecode != photcode[0].code) goto mark_nocal;
     17        // XXXX is this done when loaded, or is it needed for reload?
     18        // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
     19        // if (ecode != photcode[0].code) goto mark_nocal;
    1820       
    1921        /* select measurements by time */
    … …  
    2527        /* select measurements by mag limit */
    2628        if (AreaSelect) {
    27           r = catalog[i].average[j].R + catalog[i].measure[m].dR_PS / 3600.0;
    28           d = catalog[i].average[j].D + catalog[i].measure[m].dD_PS / 3600.0;
     29          r = catalog[i].average[j].R + catalog[i].measure[m].dR / 3600.0;
     30          d = catalog[i].average[j].D + catalog[i].measure[m].dD / 3600.0;
    2931          if ((coords = getCoords (m, i)) == NULL) goto markbad;
    3032          RD_to_XY (&x, &y, r, d, coords);

  • trunk/Ohana/src/relastro/src/setMrelFinal.c

    r6808 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void setMrelFinal (Catalog *catalog) {
    … …  
    1111
    1212    for (i = 0; i < catalog[0].Naverage; i++) {
    13       p = (PhotPrimary) ? &catalog[0].average[i].M  : &catalog[0].secfilt[PhotNsec*i+PhotSec].M_PS;
     13      p = (PhotPrimary) ? &catalog[0].average[i].M  : &catalog[0].secfilt[PhotNsec*i+PhotSec].M;
    1414      *p = NO_MAG;
    15       p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM_PS;
     15      p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM;
    1616      *p = NO_MAG;
    1717      q = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm;
    … …  
    2222       
    2323        /* select measurements by photcode */
    24         ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
     24        ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
    2525        if (ecode != photcode[0].code) continue;
    2626       
    … …  
    3131        }
    3232       
    33         catalog[0].measure[m].Mcal_PS = 0;
     33        catalog[0].measure[m].Mcal = 0;
    3434        catalog[0].measure[m].flags &= 0xff00;
    3535        catalog[0].measure[m].flags &= ~ID_MEAS_POOR;
    … …  
    9797      /** never use these measurements (wrong photcode, bad time range) */
    9898      /* skipped via NOCAL, don't mark as skipped */
    99       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
     99      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
    100100      if (ecode != photcode[0].code) continue;
    101101

  • trunk/Ohana/src/relastro/src/sort.c

    r8390 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33void sortA (double *X, int N) {

  • trunk/Ohana/src/relastro/src/write_coords.c

    r7080 r12332  
    1 # include "relphot.h"
     1# include "relastro.h"
    22
    33# define CD_COORDS 1
Note: See TracChangeset for help on using the changeset viewer.