IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

merge changes from eam_branch_20080223 : updates to DVO output format, addstar (SDSS and 2MASS flags), relastro (FitChip)

Location:
trunk/Ohana/src/relastro
Files:
2 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r16040 r16810  
    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        \

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

    r16060 r16810  
    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);

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

    r15600 r16810  
    11# include "relastro.h"
     2# define SCATTER_MAX_ERROR 0.05
    23
    34void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
    45
    5   int i;
     6  int i, Nscatter, Niter;
    67  CoordFit *fit;
     8  double dL, dM, dR, dRmax, *values;
    79
    8   fit = fit_init (coords[0].Npolyterms);
    9   for (i = 0; i < Nmatch; i++) {
    10     if (raw[i].mask) continue;
    11     fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
    12   }
    13   fit_eval (fit);
    14   fit_apply_coords (fit, coords);
    15   fit_free (fit);
     10  ALLOCATE (values, double, Nmatch);
    1611
    17   // apply new coords to raw (X,Y -> L,M)
    18   for (i = 0; i < Nmatch; i++) {
    19     XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
     12  for (Niter = 0; Niter < 3; Niter ++) {
     13    // measure the scatter distribution (use only the bright end detections)
     14    for (i = Nscatter = 0; i < Nmatch; i++) {
     15      if (raw[i].mask) continue;
     16      if (raw[i].dMag > SCATTER_MAX_ERROR) continue;
     17
     18      dL = raw[i].L - ref[i].L;
     19      dM = raw[i].M - ref[i].M;
     20      dR = hypot (dL, dM);
     21   
     22      values[Nscatter] = dR;
     23      Nscatter++;
     24    }
     25
     26    // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
     27    dsort (values, Nscatter);
     28    dRmax = 3.0*values[(int)(0.40*Nscatter)];
     29
     30    fit = fit_init (coords[0].Npolyterms);
     31
     32    for (i = 0; i < Nmatch; i++) {
     33      if (raw[i].mask) continue;
     34
     35      // require radius of XXX arcsec
     36      dL = raw[i].L - ref[i].L;
     37      dM = raw[i].M - ref[i].M;
     38      dR = hypot (dL, dM);
     39
     40      if (dR > dRmax) {
     41        continue;
     42        raw[i].mask = TRUE;
     43      }
     44   
     45      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
     46    }
     47    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n",
     48             dRmax, Nscatter, fit[0].Npts, Nmatch);
     49
     50    if (fit[0].Npts < 25) {
     51      fit_free (fit);
     52      free (values);
     53      return;
     54    }
     55    fit_eval (fit);
     56    fit_apply_coords (fit, coords);
     57    fit_free (fit);
     58
     59    for (i = 0; i < Nmatch; i++) {
     60      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
     61    }
    2062  }
    2163
    22   fit_free (fit);
     64  free (values);
     65  return;
    2366}
    2467
    … …  
    3477*/
    3578
    36 /* XXX I'm not using the errors at all : this could at least be done with the dMag values */
     79/* example using fit_apply() :
     80
     81  f = fopen ("test3.dat", "w");
     82
     83  // apply new coords to raw (X,Y -> L,M)
     84  for (i = 0; i < Nmatch; i++) {
     85    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
     86    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
     87    fprintf (f, "%f %f\n", L1, M1);
     88  }
     89  fclose (f);
     90*/
     91
     92/*** XXX this function can be improved in 4 important ways:
     93
     941) the initial value of the clipping radius could be set based on the distribution of the
     95   value (eg, inner 50%) (better to use dL, dM rather than dP, dQ)
     96
     972) multiple clipping passes could be performed
     98
     993) the per-star astrometric errors could be included in the fit (non-systematic terms at
     100   least, or scaled values of the mag error, or simply the mag error itself)
     101
     1024) we could ignore input detections on the basis of the photFlags
     103
     104*/

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

    r15238 r16810  
    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);

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

    r15600 r16810  
    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);
    … …  
    2135    raw[i].Q = raw[i].M;
    2236  }
    23 
    24   fit_free (fit);
    2537}
    2638
    … …  
    3749
    3850/* 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 */

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

    r16060 r16810  
    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;

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

    r16060 r16810  
    3838
    3939    /* set image time range */
    40     start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
    41     stop  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
     40    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
     41    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    4242
    4343    /* a new mosaic, define ranges */
    … …  
    7777
    7878    /* set image time range */
    79     start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
    80     stop  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
     79    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
     80    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    8181
    8282    /* find existing mosaic with this time range */

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

    r15600 r16810  
    2121    ref = getImageRef (catalog, Ncatalog, i, &Nstars, MODE_MOSAIC);
    2222
     23    fprintf (stderr, "image %d : Nstars: %d\n", i, Nstars);
    2324    FitChip (raw, ref, Nstars, &image[i].coords);
    2425

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

    r16060 r16810  
    1818  for (i = 0; i < Ncatalog; i++) {
    1919    for (j = 0; j < catalog[i].Naverage; j++) {
    20       Nmax = MAX (Nmax, catalog[i].average[j].Nm);
     20      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
    2121    }
    2222  }
    … …  
    4848  double Tmin, Tmax;
    4949  float errorScale;
    50   PhotCode *code;
    5150
    5251  initObjectData (catalog, Ncatalog);
    … …  
    8180
    8281      N = 0;
    83       m = catalog[i].average[j].offset;
     82      m = catalog[i].average[j].measureOffset;
    8483
    8584      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
    8685      mode = FIT_MODE;
    8786
    88       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     87      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    8988        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    9089       
    … …  
    9695        Tmax = MAX(Tmax, T[N]);
    9796
    98         /* the astrometric errors are not being carried yet (but should be!) */
    99         /* we use the photometric mag error as a weighting term */
    100 
    101         code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
    102         errorScale = code[0].astromScale;
    103         dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
    104         dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
     97        dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
     98        dD[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
    10599        dT[N] = catalog[i].measure[m].dt;
    106100
    … …  
    193187
    194188      // the measure fields must be updated before the average fields
    195       m = catalog[i].average[j].offset;
    196       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     189      m = catalog[i].average[j].measureOffset;
     190      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    197191        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    198192        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);

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

    r15600 r16810  
    202202  fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
    203203  fprintf (stderr, "  working options: \n");
    204   fprintf (stderr, "  -time (start) (stop)\n");
    205   fprintf (stderr, "  -v\n");
     204  fprintf (stderr, "  -update-objects\n");
     205  fprintf (stderr, "    -pm\n");
     206  fprintf (stderr, "    -par\n");
     207  fprintf (stderr, "    -pmpar\n");
     208  fprintf (stderr, "  -update-simple\n");
     209  fprintf (stderr, "  -update-chips\n");
     210  fprintf (stderr, "  -update-mosaics\n");
     211  fprintf (stderr, "  -time (start)(stop)\n");
     212  fprintf (stderr, "  -photcode (code)[,code,code...]\n");
    206213  fprintf (stderr, "  -plot\n");
    207214  fprintf (stderr, "  -plotdelay (seconds)\n");
    208   fprintf (stderr, "  -update\n");
     215  fprintf (stderr, "  -statmode (mode)\n");
     216  fprintf (stderr, "  -reset");
     217  fprintf (stderr, "  -update : apply new fit to database\n");
    209218  fprintf (stderr, "  -params\n");
    210   fprintf (stderr, "  -reset\n");
     219  fprintf (stderr, "  -plrange\n");
     220  fprintf (stderr, "  -minerror\n");
     221  fprintf (stderr, "  -area\n");
    211222  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
    212223  fprintf (stderr, "  -instmag min max\n\n");
    213   fprintf (stderr, "  planned options: \n");
    214   fprintf (stderr, "  -photcode code[,code,code,..]\n");
     224  fprintf (stderr, "  -v\n");
    215225  fprintf (stderr, "  \n");
    216226  exit (2);
    217227}
     228

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

    r15743 r16810  
    2222  /* exclude stars not in range or with too few measurements */
    2323  for (i = 0; i < catalog[0].Naverage; i++) {
    24     if (catalog[0].average[i].Nm < 2) continue;
     24    if (catalog[0].average[i].Nmeasure < 2) continue;
    2525
    2626    /* start with all stars good */
    2727    subcatalog[0].average[Naverage] = catalog[0].average[i];
    28     subcatalog[0].average[Naverage].offset = Nmeasure;
     28    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
    2929    for (j = 0; j < PhotNsec; j++) {
    3030      subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j];
    … …  
    3838
    3939    Nm = 0;
    40     for (j = 0; j < catalog[0].average[i].Nm; j++) {
     40    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
    4141
    42       offset = catalog[0].average[i].offset + j;
     42      offset = catalog[0].average[i].measureOffset + j;
    4343
    4444      /* select measurements by photcode, if specified */
    … …  
    9595      continue;
    9696    }
    97     subcatalog[0].average[Naverage].Nm = Nm;
     97    subcatalog[0].average[Naverage].Nmeasure = Nm;
    9898    Naverage ++;
    9999    if (Naverage == NAVERAGE) {

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

    r16060 r16810  
    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);

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

    r15743 r16810  
    3939        for (j = 0; j < catalog[i].Naverage; j++) {
    4040          catalog[i].average[j].code = 0;
    41           m = catalog[i].average[j].offset;
    42           for (k = 0; k < catalog[i].average[j].Nm; k++) {
     41          m = catalog[i].average[j].measureOffset;
     42          for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
    4343            catalog[i].measure[m+k].dbFlags = 0;
    4444          }

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

    r16060 r16810  
    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);

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

    r12332 r16810  
    1212  for (i = 0; i < Ncatalog; i++) {
    1313    for (j = 0; j < catalog[i].Naverage; j++) {
    14       Ntot += catalog[i].average[j].Nm;
     14      Ntot += catalog[i].average[j].Nmeasure;
    1515    }
    1616  }
    … …  
    2525      /* calculate the average value for a single star */
    2626      if (catalog[i].average[j].code & STAR_BAD) continue; 
    27       m = catalog[i].average[j].offset;
     27      m = catalog[i].average[j].measureOffset;
    2828
    29       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     29      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    3030        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    3131

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

    r16060 r16810  
    6464    save_catalogs (catalog, Ncatalog);
    6565  } else {
     66    // XXX for now, reload all catalogs at once; as memory use gets large, reload one-at-a-time
     67
     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);

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

    r15130 r16810  
    1111  for (i = 0; i < Ncatalog; i++) {
    1212    for (j = 0; j < catalog[i].Naverage; j++) {
    13       m = catalog[i].average[j].offset;
    14       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     13      m = catalog[i].average[j].measureOffset;
     14      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    1515
    1616        /* select measurements by photcode */
Note: See TracChangeset for help on using the changeset viewer.