IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16798


Ignore:
Timestamp:
Mar 4, 2008, 12:18:13 PM (18 years ago)
Author:
eugene
Message:

substantial work to fix FitChips

Location:
branches/eam_branch_20080223/Ohana/src/relastro
Files:
2 added
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20080223/Ohana/src/relastro/Makefile

    r16040 r16798  
    3737$(SRC)/UpdateObjects.$(ARCH).o   \
    3838$(SRC)/UpdateSimple.$(ARCH).o    \
     39$(SRC)/UpdateMeasures.$(ARCH).o  \
     40$(SRC)/GetAstromError.$(ARCH).o  \
    3941$(SRC)/args.$(ARCH).o            \
    4042$(SRC)/bcatalog.$(ARCH).o        \

  • branches/eam_branch_20080223/Ohana/src/relastro/include/relastro.h

    r16060 r16798  
    1010} CoordMode;
    1111
     12typedef enum {ERROR_MODE_RA, ERROR_MODE_DEC, ERROR_MODE_POS} ErrorMode;
     13
    1214typedef struct {
    1315  double R, D;  /* Sky Coords    - degrees */
    … …  
    1517  double L, M;  /* Focal Plane   - pixels  */
    1618  double X, Y;  /* Chip Coords   - pixels  */
    17   double Mag, dMag;
     19  double Mag;
     20  double dMag;
     21  double dPos;
    1822  int mask;
    1923} StarData;
    … …  
    263267void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
    264268void fit_eval (CoordFit *fit);
     269void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1);
    265270double **poly2d_dx (double **poly, int Nx, int Ny);
    266271double **poly2d_dy (double **poly, int Nx, int Ny);
    267272double **poly2d_copy (double **poly, int Nx, int Ny);
    268273double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
    269 void fit_apply_coords (CoordFit *fit, Coords *coords);
     274CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords);
    270275int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
    271276CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
    … …  
    278283int UpdateChips (Catalog *catalog, int Ncatalog);
    279284int UpdateMosaic (Catalog *catalog, int Ncatalog);
     285int UpdateMeasures (Catalog *catalog, int Ncatalog);
     286void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
     287
    280288int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
    281289int ParFactor (double *pR, double *pD, double R, double D, time_t T);

  • branches/eam_branch_20080223/Ohana/src/relastro/src/FitMosaic.c

    r15238 r16798  
    55  int i;
    66  CoordFit *fit;
     7  double dP, dQ, dR;
    78
    89  fit = fit_init (coords[0].Npolyterms);
    910  for (i = 0; i < Nmatch; i++) {
    1011    if (raw[i].mask) continue;
     12
     13    // require radius of XXX arcsec
     14    dP = raw[i].P - ref[i].P;
     15    dQ = raw[i].Q - ref[i].Q;
     16    dR = 3600.0 * hypot (dP, dQ);
     17
     18    // XXX the value needs to be set in a more intelligent way
     19    if (dR > 0.15) continue;
     20   
    1121    fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0);
     22  }
     23  if (fit[0].Npts == 0) {
     24    fit_free (fit);
     25    return;
    1226  }
    1327  fit_eval (fit);

  • branches/eam_branch_20080223/Ohana/src/relastro/src/FitSimple.c

    r16785 r16798  
    55  int i;
    66  CoordFit *fit;
     7  double dP, dQ, dR;
    78
    89  fit = fit_init (coords[0].Npolyterms);
    910  for (i = 0; i < Nmatch; i++) {
    1011    if (raw[i].mask) continue;
     12
     13    // require radius of XXX arcsec
     14    dP = raw[i].P - ref[i].P;
     15    dQ = raw[i].Q - ref[i].Q;
     16    dR = 3600.0 * hypot (dP, dQ);
     17
     18    // XXX the value needs to be set in a more intelligent way
     19    if (dR > 0.15) continue;
     20   
    1121    fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);
     22  }
     23  if (fit[0].Npts == 0) {
     24    fit_free (fit);
     25    return;
    1226  }
    1327  fit_eval (fit);
    … …  
    3549
    3650/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
     51
     52/* XXX See notes in FitChips.c */

  • branches/eam_branch_20080223/Ohana/src/relastro/src/ImageOps.c

    r16060 r16798  
    3333
    3434  for (i = 0; i < Nimage; i++) {
    35     start[i] = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
    36     stop[i]  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
     35    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
     36    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    3737  }
    3838}
    … …  
    8181
    8282  int i, j;
     83  char *name;
    8384
    8485  for (i = 0; i < Ncatalog; i++) {
    … …  
    8788    }
    8889  }
     90
     91  for (i = 0; VERBOSE && (i < Nimage); i++) {
     92    name = GetPhotcodeNamebyCode (image[i].photcode);
     93    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i],
     94             ohana_sec_to_date(image[i].tzero), name);
     95  }
    8996}
    9097
    … …  
    182189// return StarData values for detections in the specified image, converting coordinates from the
    183190// chip positions: X,Y -> L,M -> P,Q -> R,D
     191void fixImageRaw (Catalog *catalog, int Ncatalog, int im) {
     192
     193  int i, m, c, n;
     194  double X, Y, L, M, P, Q, R, D, dR, dD;
     195 
     196  Mosaic *mosaic;
     197  Coords *moscoords, *imcoords;
     198 
     199  moscoords = NULL;
     200  mosaic = getMosaicForImage (im);
     201  if (mosaic != NULL) {
     202      moscoords = &mosaic[0].coords;
     203  }
     204  imcoords = &image[im].coords;
     205
     206  for (i = 0; i < Nlist[im]; i++) {
     207    m = mlist[im][i];
     208    c = clist[im][i];
     209
     210    X = catalog[c].measure[m].Xccd;
     211    Y = catalog[c].measure[m].Yccd;
     212    n = catalog[c].measure[m].averef;
     213
     214    if (moscoords == NULL) {
     215      // this is a Simple image (not a mosaic)
     216      // note that for a Simple image, L,M = P,Q
     217      XY_to_LM (&L, &M, X, Y, imcoords);
     218      LM_to_RD (&R, &D, L, M, imcoords);
     219    } else {
     220      XY_to_LM (&L, &M, X, Y, imcoords);
     221      XY_to_LM (&P, &Q, L, M, moscoords);
     222      LM_to_RD (&R, &D, P, Q, moscoords);
     223    }
     224
     225    // new dR, dD : test
     226    dR = 3600.0*(catalog[c].average[n].R - R);
     227    dD = 3600.0*(catalog[c].average[n].D - D);
     228   
     229    if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
     230      // XXXXX running into this still for last megacam exposure: wrong mosaic?
     231      // ???? inconsistently hitting this????
     232      fprintf (stderr, "!");
     233      abort ();
     234    }
     235    if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
     236      fprintf (stderr, "*");
     237      abort ();
     238    }
     239
     240    catalog[c].measure[m].dR = dR;
     241    catalog[c].measure[m].dD = dD;
     242
     243    if (catalog[c].measure[m].dR > +180.0*3600.0) {
     244      // average on high end of boundary, move star up
     245      R += 360.0;
     246      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     247    }
     248    if (catalog[c].measure[m].dR < -180.0*3600.0) {
     249      // average on low end of boundary, move star down
     250      R -= 360.0;
     251      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     252    }
     253  } 
     254  return;
     255}
     256
     257// return StarData values for detections in the specified image, converting coordinates from the
     258// chip positions: X,Y -> L,M -> P,Q -> R,D
    184259StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
    185260
    186   int i, m, c;
     261  int i, m, c, n;
    187262 
    188263  Mosaic *mosaic;
    … …  
    213288    raw[i].Mag  = catalog[c].measure[m].M;
    214289    raw[i].dMag = catalog[c].measure[m].dM;
     290    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
     291
     292    n = catalog[c].measure[m].averef;
    215293
    216294    raw[i].mask = FALSE;
     295    if (catalog[c].average[n].Nmeasure < 2) {
     296      raw[i].mask = TRUE;
     297    }
    217298
    218299    switch (mode) {
    … …  
    273354    ref[i].Mag  = catalog[c].measure[m].M;
    274355    ref[i].dMag = catalog[c].measure[m].dM;
     356    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
    275357
    276358    ref[i].mask = FALSE;
    … …  
    286368      case MODE_MOSAIC:
    287369      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
    288       LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords);
     370      LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
    289371      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
    290372      break;

  • branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateObjects.c

    r16787 r16798  
    3737  ALLOCATE (pY, double, MAX (1, Nmax));
    3838} 
    39 
    40 enum {ERROR_MODE_RA, ERROR_MODE_DEC};
    41 
    42 float GetAstromError (Measure *measure, int mode) {
    43 
    44   PhotCode *code;
    45   float dPobs, dPsys, dPtotal, dM, AS, MS;
    46 
    47   switch (mode) {
    48     case ERROR_MODE_RA:
    49       dPobs = measure[0].dXccd;  // need to redefine this as RAerr
    50       break;
    51     case ERROR_MODE_DEC:
    52       dPobs = measure[0].dYccd;  // need to redefine this as RAerr
    53       break;
    54     default:
    55       abort();
    56   }
    57 
    58   /* the astrometric errors are not being carried yet (but should be!) */
    59   /* we use the photometric mag error as a weighting term */
    60 
    61   code  = GetPhotcodebyCode (measure[0].photcode);
    62   AS    = code[0].astromErrScale;
    63   MS    = code[0].astromErrMagScale;
    64   dPsys = code[0].astromErrSys;
    65   dM    = measure[0].dM;
    66  
    67   dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
    68   dPtotal = MAX (dPtotal, MIN_ERROR);
    69 
    70   return (dPtotal);
    71 }
    7239
    7340int UpdateObjects (Catalog *catalog, int Ncatalog) {

  • branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c

    r16060 r16798  
    128128      matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
    129129    }
     130
     131    // mask the terms not represented by the Coords terms
     132    if (ix + iy > fit[0].Norder) {
     133      for (j = 0; j < fit[0].Nelems; j++) {
     134        matrix[i][j] = 0.0;
     135      }
     136      vector[i][0] = 0.0;
     137      vector[i][1] = 0.0;
     138      matrix[i][i] = 1.0;
     139    }     
    130140  }
    131141
    … …  
    133143    ix = i % fit[0].Nterms;
    134144    iy = i / fit[0].Nterms;
    135     fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n",
    136             ix, iy, vector[i][0], ix, iy, vector[i][1]);
     145    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n",
     146    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
    137147  }     
    138148
    … …  
    142152    ix = i % fit[0].Nterms;
    143153    iy = i / fit[0].Nterms;
    144     fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n",
    145             ix, iy, vector[i][0], ix, iy, vector[i][1]);
     154    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n",
     155    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
    146156  }     
    147157
    … …  
    158168}
    159169
     170void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) {
     171 
     172  int ix, iy;
     173  double xterm, yterm, term;
     174  double x, y;
     175
     176  x = 0.0;
     177  y = 0.0;
     178
     179  xterm = 1;
     180  for (ix = 0; ix < fit[0].Nterms; ix++) {
     181    yterm = 1;
     182    for (iy = 0; iy < fit[0].Nterms; iy++) {
     183      term = xterm*yterm;
     184      x += fit[0].xfit[ix][iy]*term;
     185      y += fit[0].yfit[ix][iy]*term;
     186      yterm *= y1;
     187    }
     188    xterm *= x1;
     189  }
     190  *x2 = x;
     191  *y2 = y;
     192}
     193
    160194// Nx, Ny is the number of terms in x and in y
    161195double **poly2d_dx (double **poly, int Nx, int Ny) {
    … …  
    209243  for (i = 0; i < Nx; i++) {
    210244    for (j = 0; j < Ny; j++) {
    211       out[i][j] = poly[i][j+1];
     245      out[i][j] = poly[i][j];
    212246    }
    213247  }
    … …  
    237271/* this should only apply to the polynomial, not the projection terms */
    238272/* compare with psastro supporting code */
    239 void fit_apply_coords (CoordFit *fit, Coords *coords) {
     273CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) {
    240274
    241275  double c11, c12;
    242276  double c21, c22;
    243   double Xo, Yo, R;
     277  double Xo, Yo, R1, R2;
    244278  CoordFit *modfit;
    245279
    … …  
    256290  modfit = CoordsSetCenter (fit, Xo, Yo);
    257291
     292  /* we do not modify crval1,2: these are kept at the default values */
     293
     294  // set cdelt1, cdelt2
     295  coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]);
     296  coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]);
     297  R1 = 1 / coords[0].cdelt1;
     298  R2 = 1 / coords[0].cdelt2;
     299
    258300  // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
    259   coords[0].pc1_1 = modfit[0].xfit[1][0];
    260   coords[0].pc1_2 = modfit[0].xfit[0][1];
    261   coords[0].pc2_1 = modfit[0].yfit[1][0];
    262   coords[0].pc2_2 = modfit[0].yfit[0][1];
     301  coords[0].pc1_1 = modfit[0].xfit[1][0] * R1;
     302  coords[0].pc1_2 = modfit[0].xfit[0][1] * R2;
     303  coords[0].pc2_1 = modfit[0].yfit[1][0] * R1;
     304  coords[0].pc2_2 = modfit[0].yfit[0][1] * R2;
    263305
    264306  // set the polyterm elements
    265307  if (coords->Npolyterms > 1) {
    266     coords[0].polyterms[0][0] = modfit[0].xfit[2][0];
    267     coords[0].polyterms[1][0] = modfit[0].xfit[1][1];
    268     coords[0].polyterms[2][0] = modfit[0].xfit[0][2];
    269 
    270     coords[0].polyterms[0][1] = modfit[0].yfit[2][0];
    271     coords[0].polyterms[1][1] = modfit[0].yfit[1][1];
    272     coords[0].polyterms[2][1] = modfit[0].yfit[0][2];
     308    coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1;
     309    coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2;
     310    coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2;
     311
     312    coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1;
     313    coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2;
     314    coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2;
    273315  }
    274316
    275317  // I need to validate Norder
    276318  if (coords->Npolyterms > 2) {
    277     coords[0].polyterms[3][0] = modfit[0].xfit[3][0];
    278     coords[0].polyterms[4][0] = modfit[0].xfit[2][1];
    279     coords[0].polyterms[5][0] = modfit[0].xfit[1][2];
    280     coords[0].polyterms[6][0] = modfit[0].xfit[0][3];
    281 
    282     coords[0].polyterms[3][1] = modfit[0].yfit[3][0];
    283     coords[0].polyterms[4][1] = modfit[0].yfit[2][1];
    284     coords[0].polyterms[5][1] = modfit[0].yfit[1][2];
    285     coords[0].polyterms[6][1] = modfit[0].yfit[0][3];
    286   }
    287 
    288   /* we do not modify crval1,2: these are kept at the default values */
    289 
    290   /* normalize pc11,etc */
    291   c11 = coords[0].pc1_1;
    292   c21 = coords[0].pc1_1;
    293   c12 = coords[0].pc1_1;
    294   c22 = coords[0].pc1_1;
    295   coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
    296   R = 1 / coords[0].cdelt1;
    297 
    298   coords[0].pc1_1  = c11*R;
    299   coords[0].pc2_1  = c21*R;
    300   coords[0].pc1_2  = c12*R;
    301   coords[0].pc2_2  = c22*R;
     319    coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1;
     320    coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2;
     321    coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2;
     322    coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2;
     323
     324    coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1;
     325    coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2;
     326    coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2;
     327    coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2;
     328  }
    302329
    303330  fit_free (modfit);

  • branches/eam_branch_20080223/Ohana/src/relastro/src/mkpolyterm.c

    r16060 r16798  
    2020
    2121  int i, Nx, Ny;
    22   double R, Xo, Yo, dPos;
     22  double R, Xo, Yo, dPos, dPosRef;
    2323  double **XdX, **XdY, **YdX, **YdY, **alpha, **beta;
    2424  double **xfit, **yfit;
    … …  
    5151     */
    5252    dPos = tol + 1;
     53    dPosRef = dPos;
    5354    for (i = 0; (dPos > tol) && (i < 20); i++) {
    5455      // NOTE: order for alpha is: [y][x]
    … …  
    6667      Yo -= beta[1][0];
    6768      dPos = hypot(beta[0][0], beta[1][0]);
     69      if (i == 0) {
     70        dPosRef = dPos;
     71      }
     72    }
     73    if (dPos > dPosRef) {
     74      fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n");
    6875    }
    6976    array_free (alpha, 2);

  • branches/eam_branch_20080223/Ohana/src/relastro/src/relastro.c

    r16788 r16798  
    6161  if (!UPDATE) exit (0);
    6262
    63   save_catalogs (catalog, Ncatalog);
     63  if (FIT_TARGET == TARGET_OBJECTS) {
     64    save_catalogs (catalog, Ncatalog);
     65  } else {
     66    // XXX for now, reload all catalogs at once; as memory use gets large, reload one-at-a-time
    6467
    65   if (FIT_TARGET != TARGET_OBJECTS) {
     68    // load catalog data from region files : do not subselect
     69    catalog = load_catalogs (skylist, &Ncatalog, FALSE);
     70   
     71    // match measurements with images
     72    initImageBins (catalog, Ncatalog);
     73    findImages (catalog, Ncatalog);
     74
     75    // update the detection coordinates using the new image parameters
     76    UpdateMeasures (catalog, Ncatalog);
     77
     78    // write the updated detections to disk
     79    save_catalogs (catalog, Ncatalog);
     80
     81    // save the updated image parameters
    6682    dvo_image_update (&db, VERBOSE);
    6783    dvo_image_unlock (&db);
Note: See TracChangeset for help on using the changeset viewer.