IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37622


Ignore:
Timestamp:
Nov 18, 2014, 1:58:02 AM (12 years ago)
Author:
eugene
Message:

working on ICRF tie

Location:
branches/eam_branches/ipp-20140904/Ohana/src/relastro
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r37606 r37622  
    583583int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords);
    584584int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE);
     585
     586int dump_stardata_pts (StarData *raw, int Npts, char *filename);
     587void printNcatTotal ();
     588

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

    r37606 r37622  
    323323  return TRUE;
    324324}
     325
     326int dump_stardata_pts (StarData *raw, int Npts, char *filename) {
     327
     328  FILE *f = fopen (filename, "w");
     329  if (!f) {
     330    fprintf (stderr, "failed to open file %s\n", filename);
     331    return FALSE;
     332  }
     333
     334  int i;
     335  for (i = 0; i < Npts; i++) {
     336    fprintf (f, "%4d %10.6f %10.6f  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
     337             i, raw[i].R, raw[i].D, raw[i].P, raw[i].Q, raw[i].L, raw[i].M, raw[i].X, raw[i].Y);
     338  }
     339  fclose (f);
     340
     341  return TRUE;
     342}

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

    r37604 r37622  
    11# include "relastro.h"
     2int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection);
    23
    34/* The sequence of calibration is something like this:
    … …  
    2728
    2829  // binD = (DEC + 89.0) / scale
    29   ALLOCATE (frame->Nra, int, frame->Ndec);
     30  ALLOCATE (frame->Nra,  int, frame->Ndec);
     31  ALLOCATE (frame->dR,   double, frame->Ndec);
    3032  ALLOCATE (frame->Roff, double *, frame->Ndec);
    3133  ALLOCATE (frame->Doff, double *, frame->Ndec);
    … …  
    8688  FrameCorrectionType *frame = NULL;
    8789
    88 # define LOCAL_FRAME 1
    89 # define SH_FRAME 0
     90# define LOCAL_FRAME 0
     91# define SH_FRAME 1
    9092
    9193  int POLE = FALSE;
    … …  
    144146
    145147    // XXX write out an image to represent the correction
     148    FrameCorrectionWriteImageSH (frame, "map.dR.fits", TRUE);
     149    FrameCorrectionWriteImageSH (frame, "map.dD.fits", FALSE);
    146150
    147151    SHtermsFree (dRc);
    … …  
    166170
    167171  int Nicrf = ICRFmax();
     172  fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, Nicrf);
    168173
    169174  int i, j;
    … …  
    186191    double dR = average->R - measure->R;
    187192    double dD = average->D - measure->D;
     193
     194    if (isnan(dR)) continue;
     195    if (isnan(dD)) continue;
    188196
    189197    // get the value of Ylm at this coordinate
    … …  
    238246}
    239247
     248// write out the maps as images
     249int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection) {
     250
     251  Header header;
     252  Matrix matrix;
     253  Coords coords;
     254   
     255  gfits_init_header (&header);
     256  header.bitpix = -32;
     257  header.Naxes = 2;
     258  header.Naxis[0] = 370.0 / frame->scale;
     259  header.Naxis[1] = 190.0 / frame->scale;
     260
     261  int Nx = header.Naxis[0];
     262  int Ny = header.Naxis[1];
     263
     264  gfits_create_header (&header);
     265  gfits_create_matrix (&header, &matrix);
     266
     267  InitCoords (&coords, "DEC--AIT");
     268  coords.cdelt1 = coords.cdelt2 = frame->scale;
     269  coords.crval1 = 0.0; // SAS center
     270  coords.crval2 = 0.0;
     271  coords.crpix1 = 0.5*Nx; // middle of projection is middle of map
     272  coords.crpix2 = 0.5*Ny; // middle of projection is middle of map
     273   
     274  PutCoords (&coords, &header);
     275
     276  float *buffer = (float *)matrix.buffer;
     277
     278  int ix, iy;
     279  for (ix = 0; ix < Nx; ix++) {
     280    for (iy = 0; iy < Ny; iy++) {
     281      double R, D;
     282      int status = XY_to_RD (&R, &D, ix, iy, &coords);
     283      if (!status) continue;
     284
     285      int iD = (D + 89.0) / frame->scale;
     286      int iR = R / frame->dR[iD];
     287
     288      float value = NAN;
     289      if (raDirection) {
     290        value = frame->Roff[iD][iR];
     291      } else {
     292        value = frame->Doff[iD][iR];
     293      }
     294      buffer[ix + Nx*iy] = value;
     295    }
     296  }
     297  gfits_write_header (filename, &header);
     298  gfits_write_matrix (filename, &matrix);
     299  return TRUE;
     300}
     301
    240302int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords) {
    241303

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

    r37604 r37622  
    1818int ICRFsave (int cat, int ave, int meas) {
    1919 
    20   fprintf (stderr, "ICRF: %d %d %d\n", cat, ave, meas);
     20  // fprintf (stderr, "ICRF: %d %d %d\n", cat, ave, meas);
    2121
    2222  ICRFtoCatalog[Nicrf] = cat;

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

    r37606 r37622  
    816816    // XXX subtract off dR,dD GAL here
    817817    if (USE_GALAXY_MODEL) {
     818      myAssert (!isnan(measure[0].RoffGAL), "oops");
     819      myAssert (!isnan(measure[0].DoffGAL), "oops");
    818820      ref[i].R += measure[0].RoffGAL / 3600.0;
    819821      ref[i].D += measure[0].DoffGAL / 3600.0;

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

    r35763 r37622  
    11# include "relastro.h"
    2 
    3 void printNcatTotal ();
    42
    53// this function operates on Measure, not MeasureTiny
    … …  
    3230  }
    3331
    34   printNcatTotal ();
     32  // printNcatTotal ();
    3533
    3634  return (TRUE);

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

    r37606 r37622  
    7474    save_catalogs (&catalog, 1);
    7575  }
    76  
     76  printNcatTotal();
    7777  return (TRUE);
    7878}

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

    r37606 r37622  
    132132      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop);
    133133    }
    134     if (VERBOSE) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
    135     if (VERBOSE) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
     134    if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
     135    if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
    136136    sumFitStats (&fitStatsChips, &sumStatsChips);
    137137    sumFitStats (&fitStatsStack, &sumStatsStack);
    … …  
    139139  freeObjectData ();
    140140
    141   if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
    142   if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
     141  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
     142  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
    143143  return (TRUE);
    144144}
    … …  
    501501  coords.crval1 = average[0].R;
    502502  coords.crval2 = average[0].D;
     503  if (isnan(coords.crval1)) {
     504    return (FALSE);
     505  }
     506  if (isnan(coords.crval2)) {
     507    return (FALSE);
     508  }
    503509
    504510  double dXoff, dYoff;

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

    r37604 r37622  
    247247  }
    248248
    249   fprintf (stderr, "added %d ICRF QSO\n", Nicrf);
     249  // fprintf (stderr, "added %d ICRF QSO\n", Nicrf);
    250250 
    251251  return (TRUE);

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

    r37438 r37622  
    353353  /* keep the order and type from initial values */
    354354 
     355  if (isnan(coords[0].crval1)) {
     356    return FALSE;
     357  }
     358  if (isnan(coords[0].crval2)) {
     359    return FALSE;
     360  }
     361  if (isnan(coords[0].crpix1)) {
     362    return FALSE;
     363  }
     364  if (isnan(coords[0].crpix2)) {
     365    return FALSE;
     366  }
     367
    355368  return (TRUE);
    356369}

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

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