IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37545


Ignore:
Timestamp:
Nov 2, 2014, 3:35:12 PM (12 years ago)
Author:
eugene
Message:

added ICRF QSOs and Spherical Harmonics-based FrameCorrection code

Location:
branches/eam_branches/ipp-20140904/Ohana/src/relastro
Files:
2 added
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/Makefile

    r37261 r37545  
    7777$(SRC)/share_mean_pos.$(ARCH).o \
    7878$(SRC)/share_images_pos.$(ARCH).o \
     79$(SRC)/ICRF.$(ARCH).o \
     80$(SRC)/FrameCorrection.$(ARCH).o \
    7981$(SRC)/ImagePosIO.$(ARCH).o \
    8082$(SRC)/ImageTable.$(ARCH).o \
    … …  
    128130$(SRC)/plot_scatter.$(ARCH).o        \
    129131$(SRC)/plotstuff.$(ARCH).o           \
     132$(SRC)/ICRF.$(ARCH).o \
     133$(SRC)/FrameCorrection.$(ARCH).o \
    130134$(SRC)/relastroVisual.$(ARCH).o \
    131135$(SRC)/syncfile.$(ARCH).o \

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/doc/pv2.txt

    r37038 r37545  
     1
     22014.10.28
     3
     4Going over the sequence of steps in relastro.  I'm concerned that the
     5iteration is not modifying the measure.R,D values, preventing any
     6actual iterative improvements.
     7
     8* UpdateObjects : takes measure.R,D to calculate average.R,D and
     9  motion terms
     10
     11* ImageOps :
     12  * getImageRaw : takes measure.X,Y, applies image.coords to get (L,M),
     13                  (P,Q), etc, calculates
     14  * getImageRef : takes average.R,D, applies image.coords to get
     15                  (P,Q), (L,M), etc
     16  * FitChip : fit raw.X,Y to ref.L,M, plus residual map (dX,dY)
     17 
     18  *** modified raw.L,M values are then dropped and not used to update measure.R.D
    119
    2202014.07.03

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h

    r37498 r37545  
    142142} FitStats;
    143143
     144typedef struct {
     145  double scale;
     146  double **Roff;
     147  double **Doff;
     148  double  *dR;
     149  int *Nra;
     150  int Ndec;
     151} FrameCorrectionType;
     152
    144153/* global variables set in parameter file */
    145154# define DVO_MAX_PATH 1024
    … …  
    188197
    189198int    USE_FIXED_PIXCOORDS;
    190 int    USE_GAL_MODEL;
     199int    USE_GALAXY_MODEL;
     200int    USE_ICRF_CORRECT;
    191201
    192202int    RESET;
    … …  
    419429StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
    420430StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
     431int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode);
    421432off_t getImageByID (off_t ID);
    422433
    … …  
    555566int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts);
    556567
     568// ICRF QSOs : these must be marked in the database (flag on average, flag on measure)
     569
     570void ICRFinit ();
     571int ICRFsave (int cat, int ave, int meas);
     572int ICRFdata (int n, int *cat, int *ave, int *meas);
     573int ICRFmax ();
     574
     575FrameCorrectionType *FrameCorrectionInit (double scale);
     576void FrameCorrectionFree (FrameCorrectionType *frame);
     577int FrameCorrection (Catalog *catalog, int Ncatalog);
     578int FrameCorrectionFit (Catalog *catalog, int Ncatalog, SHterms *dR, SHterms *dD);
     579int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD);

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ConfigInit.c

    r37498 r37545  
    3636
    3737  if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS))  USE_FIXED_PIXCOORDS = FALSE;
    38   if (!ScanConfig (config, "USE_GAL_MODEL",       "%d", 0, &USE_GAL_MODEL))        USE_GAL_MODEL = FALSE;
     38  if (!ScanConfig (config, "USE_GALAXY_MODEL",    "%d", 0, &USE_GALAXY_MODEL))     USE_GALAXY_MODEL = FALSE;
     39  if (!ScanConfig (config, "USE_ICRF_CORRECT",    "%d", 0, &USE_ICRF_CORRECT))     USE_ICRF_CORRECT = FALSE;
    3940
    4041  // force CATDIR to be absolute (so parallel mode will work)

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImageOps.c

    r37455 r37545  
    809809
    810810    // XXX subtract off dR,dD GAL here
    811     if (USE_GAL_MODEL) {
     811    if (USE_GALAXY_MODEL) {
    812812      ref[i].R += measure[0].RoffGAL;
    813813      ref[i].D += measure[0].DoffGAL;
    … …  
    848848  *Nstars = N_onImage[im];
    849849  return (ref);
     850}
     851
     852// return StarData values for detections in the specified image, converting coordinates from the
     853// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
     854// which the detections have already been filtered when they were loaded (bcatalog)
     855int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode) {
     856
     857  off_t i, m, c;
     858
     859  Coords *moscoords;
     860
     861  moscoords = NULL;
     862  if (mode == MODE_MOSAIC) {
     863    moscoords = image[im].coords.mosaic;
     864    myAssert (moscoords, "coords.mosaic not defined for image %s (%d)", image[im].name, (int) im);
     865  }
     866
     867  myAssert (Nraw == N_onImage[im], "impossible!");
     868
     869  for (i = 0; i < N_onImage[im]; i++) {
     870    m = ImageToMeasure[im][i];
     871    c = ImageToCatalog[im][i];
     872
     873    // XXX should I use the raw coords or just measure.X,Y -> R,D?
     874
     875    switch (mode) {
     876      case MODE_SIMPLE:
     877        /* note that for a Simple image, L,M = P,Q */
     878        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     879        raw[i].P = raw[i].L;
     880        raw[i].Q = raw[i].M;
     881        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
     882        break;
     883      case MODE_MOSAIC:
     884        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     885        XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
     886        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
     887        break;
     888      default:
     889        fprintf (stderr, "error: invalid mode in getImageRaw");
     890        abort ();
     891    }
     892
     893    MeasureTiny *measure = &catalog[c].measureT[m];
     894    measure->R = raw[i].R;
     895    measure->D = raw[i].D;
     896    if (catalog[c].measure) {
     897      catalog[c].measure[m].R = raw[i].R;
     898      catalog[c].measure[m].D = raw[i].D;
     899    }     
     900  }
     901  return TRUE;
    850902}
    851903

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateChips.c

    r37378 r37545  
    132132      continue;
    133133    }
     134
     135    // apply the modified R,D back to the measures
     136    setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC);
    134137
    135138    saveCenter (image, &Ro[i], &Do[i], i);

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjects.c

    r37531 r37545  
    209209    double Di = getMeanD (&measure[k], average, secfilt);
    210210
    211     if (USE_GAL_MODEL) {
     211    if (USE_GALAXY_MODEL) {
    212212      Ri -= measure[0].RoffGAL;
    213213      Di -= measure[0].DoffGAL;

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/bcatalog.c

    r37386 r37545  
    6161  }
    6262
     63  int Nicrf = 0;
     64
    6365  /* exclude stars not in range or with too few measurements */
    6466  for (i = 0; i < catalog[0].Naverage; i++) {
    6567    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
    6668      continue;
    67       // XXX can't (easily) count this for a given chip
    6869    }
    6970   
    … …  
    7778    if (RESET) {
    7879      subcatalog[0].average[Naverage].flags &= ~averageBits;
     80    }
     81
     82    // if we want to do static correction to ICRF, keep the ICRF QSOs
     83    if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) {
     84      Nm = 0;
     85      for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
     86         
     87        offset = catalog[0].average[i].measureOffset + j;
     88         
     89        CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     90        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
     91        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
     92        subcatalog[0].measureT[Nmeasure].averef   = Naverage;
     93        if (RESET) {
     94          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     95          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     96        }
     97         
     98        // we need to save the location of the ICRF QSOs in the database
     99        if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_ICRF_QSO) {
     100          ICRFsave (i, Naverage, Nmeasure);
     101          Nicrf ++;
     102        }
     103
     104        Nmeasure ++;
     105        Nm ++;
     106        if (Nmeasure == NMEASURE) {
     107          NMEASURE += 1000;
     108          REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
     109        }
     110      }
     111      subcatalog[0].average[Naverage].Nmeasure = Nm;
     112      Naverage ++;
     113      if (Naverage == NAVERAGE) {
     114        NAVERAGE += 50;
     115        REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
     116        REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     117      }
     118      continue;
     119    }
     120
     121    // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly
     122    if (USE_GALAXY_MODEL && !catalog[0].average[i].Nstarpar) {
     123      // we may want to add some starpar quality filters...
     124      continue;
    79125    }
    80126
    … …  
    109155      // exclude bogus
    110156      if (ExcludeBogus) {
    111           double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
    112           double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
    113           coords.crval1 = catalog[0].average[i].R;
    114           coords.crval2 = catalog[0].average[i].D;
    115           double Xi, Yi;
    116           RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
    117           double radius = hypot(Xi, Yi);
    118           if (radius > ExcludeBogusRadius) {
    119               NskipBogus ++;
    120               if (VERBOSE2) {
    121                 FILE *foutput = fbogus ? fbogus : stderr;
    122                 char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
    123                 fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
    124                 free (date);
    125               }
    126               continue;
     157        double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     158        double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     159        coords.crval1 = catalog[0].average[i].R;
     160        coords.crval2 = catalog[0].average[i].D;
     161        double Xi, Yi;
     162        RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
     163        double radius = hypot(Xi, Yi);
     164        if (radius > ExcludeBogusRadius) {
     165          NskipBogus ++;
     166          if (VERBOSE2) {
     167            FILE *foutput = fbogus ? fbogus : stderr;
     168            char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
     169            fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
     170            free (date);
    127171          }
     172          continue;
     173        }
    128174      }
    129175
    … …  
    150196      subcatalog[0].measureT[Nmeasure].averef   = Naverage;
    151197      if (RESET) {
    152         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    153         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     198        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     199        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    154200      }
    155201
    … …  
    157203      Nm ++;
    158204      if (Nmeasure == NMEASURE) {
    159         NMEASURE += 1000;
    160         REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
     205        NMEASURE += 1000;
     206        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    161207      }
    162208    }
    … …  
    179225  assert (Nsecfilt == catalog[0].Nsecfilt);
    180226
    181   // limit the total number of stars in the catalog
     227// limit the total number of stars in the catalog
    182228  if (MaxDensityUse) {
    183229    LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
    … …  
    194240  }
    195241
     242  fprintf (stderr, "added %d ICRF QSO\n", Nicrf);
     243 
    196244  return (TRUE);
    197245}
    … …  
    210258
    211259# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
    212   tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
    213   tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
    214 }
     260    tmpI = X[A]; X[A] = X[B]; X[B] = tmpI;      \
     261    tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT;      \
     262  }
    215263# define COMPARE(A,B)(X[A] > X[B])
    216264
    … …  
    224272/* generate a grid in a locally-projected coordinate system, try to select average entries
    225273   from each grid cell in decending Nmeasure order.
    226  */
     274*/
    227275int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) {
    228276

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_catalogs.c

    r37038 r37545  
    8888
    8989  bcatalog_show_skips();
     90  fprintf (stderr, "included %d ICRF QSOs\n", ICRFmax());
    9091
    9192  Nstar = 0;

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c

    r37528 r37545  
    7474        // XXX is is really possible that I do not update the measure.R,D on each loop??
    7575        UpdateObjects (catalog, Ncatalog); // calculate <R>,<D>
     76        FrameCorrection (catalog, Ncatalog);
    7677        UpdateChips (catalog, Ncatalog);   // measure.X,Y -> R,D, fit image.coords
    77         FrameCorrection (catalog, Ncatalog);
    7878        MARKTIME("update chips: %f sec\n", dtime);
    7979      }
Note: See TracChangeset for help on using the changeset viewer.