IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 6, 2016, 9:44:08 AM (10 years ago)
Author:
eugene
Message:

split galaxy motion model option into two options: galaxy motion model + apply proper motion; remove deprecated chunk of code; add options to skip HSC, CFH, PS1_CHIP, or PS1_STACK images (leaving them static); calculate the error based on the ordinary least squares and regular poisson weights (when doing bootstrap error analysis, choose the MAX); calculate dXpixSys,dYpixSys in arcsec based on plate scale; add option to avoid changing measurements for specific cameras; StarMap tests are now done in real arcsec using plate scales; option to fit the chips without modifying the mean objects; remove excessively verbose memory checks; add tests for HSC and CFH photcodes

Location:
branches/czw_branch/20160809/Ohana/src/relastro
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20160809/Ohana/src/relastro/include/relastro.h

    r39719 r39731  
    3333typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES, OP_REPAIR_STACKS, OP_REPAIR_WARPS, OP_REPAIR_OBJECT_ID} RelastroOp;
    3434
    35 typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
     35typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, SET_CHIPS, TARGET_MOSAICS} FitTarget;
    3636
    3737typedef enum {
    … …  
    144144  int Nfit;
    145145  int converged;
     146  int useWeight;
     147
    146148} FitAstromResult;
    147149
    … …  
    364366int    USE_ALL_IMAGES;
    365367int    CHECK_MEASURE_TO_IMAGE;
     368
     369int    SKIP_PS1_CHIP;
     370int    SKIP_PS1_STACK;
     371int    SKIP_HSC;
     372int    SKIP_CFH;
     373
     374int    UPDATE_PS1_CHIP_MEASURE;
     375int    UPDATE_HSC_MEASURE;
     376int    UPDATE_CFH_MEASURE;
     377
     378int    APPLY_PROPER_MOTION;
    366379
    367380int    RESET;
    … …  
    744757int isGPC1stack (int photcode);
    745758int isGPC1warp (int photcode);
     759int isHSCchip (int photcode);
     760int isCFHchip (int photcode);
    746761
    747762int save_astrom_table ();

  • branches/czw_branch/20160809/Ohana/src/relastro/src/BrightCatalog.c

    r39457 r39731  
    353353  /*** MeasureTiny ***/
    354354  {
    355     ohana_memcheck (1);
     355    // ohana_memcheck (1);
    356356    gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY");
    357357
    358     ohana_memcheck (1);
     358    // ohana_memcheck (1);
    359359
    360360    gfits_define_bintable_column (&theader, "D", "RA",       "ra",                         "degrees", 1.0, 0.0);
    … …  
    429429    gfits_set_bintable_column (&theader, &ftable, "RA",         R,         catalog->Nmeasure);
    430430
    431     fprintf (stderr, "--------------- after set_bintable RA --------------");
    432     ohana_memdump_file (stderr, TRUE);
     431    // fprintf (stderr, "--------------- after set_bintable RA --------------");
     432    // ohana_memdump_file (stderr, TRUE);
    433433   
    434434    gfits_set_bintable_column (&theader, &ftable, "DEC",        D,         catalog->Nmeasure);

  • branches/czw_branch/20160809/Ohana/src/relastro/src/FitAstromOps.c

    r39610 r39731  
    257257  fit->converged = FALSE;
    258258 
     259  // this is an input value
     260  // if true, use the IRLS modified weight
     261  fit->useWeight = FALSE;
     262
    259263  return;
    260264}

  • branches/czw_branch/20160809/Ohana/src/relastro/src/FitChip.c

    r39457 r39731  
    206206  if (VERBOSE2) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);
    207207
    208   image[0].dXpixSys = dLsig;
    209   image[0].dYpixSys = dMsig;
     208  // need to convert dLsig, dMsig back to pixel scale (or up to arcsec)
     209
     210  float plateScale;
     211  if (image[0].coords.mosaic) {
     212    // NOTE: for the full pixel to sky plate scale, use this:
     213    // float plateScaleX = 3600.0*image[0].coords.mosaic->cdelt1*image[0].coords.cdelt1;
     214    // float plateScaleY = 3600.0*image[0].coords.mosaic->cdelt2*image[0].coords.cdelt2;
     215
     216    // since we are compare L,M values, just need to compensate for focal plate to sky:
     217    float plateScaleX = 3600.0*fabs(image[0].coords.mosaic->cdelt1);
     218    float plateScaleY = 3600.0*fabs(image[0].coords.mosaic->cdelt2);
     219    plateScale = 0.5*(plateScaleX + plateScaleY);
     220  } else {
     221    // since we are compare L,M values, just need to compensate for arcsec vs degrees:
     222    plateScale = 3600.0;
     223  }
     224
     225  image[0].dXpixSys = plateScale*dLsig;
     226  image[0].dYpixSys = plateScale*dMsig;
    210227  image[0].nFitAstrom = fit[0].Npts;
    211228
    212   // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom);
     229  fprintf (stderr, "%s | %6.3f %6.3f | %4d %4d | %6.3f %6.3f\n", image[0].name, image[0].refColorRed, image[0].refColorBlue, Ncolor, image[0].nFitAstrom, image[0].dXpixSys, image[0].dYpixSys);
    213230
    214231  if (fit) fit_free (fit);

  • branches/czw_branch/20160809/Ohana/src/relastro/src/FitPM.c

    r39612 r39731  
    2525  }
    2626
     27  fit->useWeight = FALSE; // Ordinary Least Squares
    2728  if (!FitPM_MinChisq (fit, data, points, Npoints)) return FALSE;
    2829  if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE;
    … …  
    4748 
    4849  // Solve OLS equation 
     50  fit->useWeight = FALSE; // Ordinary Least Squares
    4951  if (!FitPM_MinChisq(fit, data, points, Npoints)) {
    5052    return(FALSE);
    … …  
    9092
    9193    // Solve with the new weights
     94    fit->useWeight = TRUE; // Reweighted Least Squares
    9295    if (!FitPM_MinChisq(fit, data, points, Npoints)) {
    9396
    … …  
    159162  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
    160163  if (data->getError) {
    161     double ax = 0.0, ay = 0.0;
    162     double bx = 0.0, by = 0.0;
    163 
    164     for (i = 0; i < Npoints; i++) {
    165       ax += dpsi_cauchy(points[i].rx / points[i].dX);
    166       ay += dpsi_cauchy(points[i].ry / points[i].dY);
    167 
    168       bx += SQ(points[i].Wx);
    169       by += SQ(points[i].Wy);
    170     }
    171     ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    172     ay /= 1.0 * Npoints;
    173     bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
    174     by /= 1.0 * (Npoints - data->Nterms);
    175  
    176     double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
    177     double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
    178  
    179     double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
    180     double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
    181 
    182     // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
    183     double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    184     double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    185 
    186     fit[0].dRo = sqrt(data->Cov[0][0]);
    187     fit[0].duR = sqrt(data->Cov[1][1]);
    188     fit[0].dDo = sqrt(data->Cov[2][2]);
    189     fit[0].duD = sqrt(data->Cov[3][3]);
    190 
    191     fit[0].dRo *= sigma_final_x;
    192     fit[0].duR *= sigma_final_x;
    193     fit[0].dDo *= sigma_final_y;
    194     fit[0].duD *= sigma_final_y;
     164    FitAstromResult fitErrors;
     165    FitAstromResultInit (&fitErrors);
     166    fitErrors.useWeight = FALSE;
     167
     168    FitPM_MinChisq(&fitErrors, data, points, Npoints);
     169
     170    // we use the errors from a simple OLS, ignoring masked points
     171    fit[0].dRo = fitErrors.dRo;
     172    fit[0].duR = fitErrors.duR;
     173    fit[0].dDo = fitErrors.dDo;
     174    fit[0].duD = fitErrors.duD;
    195175  }
    196176
    … …  
    212192    Nfit ++;
    213193
    214     wx = points[i].qx;
    215     wy = points[i].qy;
     194    if (fit->useWeight) {
     195      wx = points[i].qx;
     196      wy = points[i].qy;
     197    } else {
     198      wx = points[i].Qx;
     199      wy = points[i].Qy;
     200    }
    216201
    217202    Wx += wx;

  • branches/czw_branch/20160809/Ohana/src/relastro/src/FitPMandPar.c

    r39612 r39731  
    2525  }
    2626
     27  fit->useWeight = FALSE; // Ordinary Least Squares
    2728  if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE;
    2829  if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE;
    … …  
    4748 
    4849  // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
     50  fit->useWeight = FALSE; // Ordinary Least Squares
    4951  if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
    5052    return(FALSE);
    … …  
    9092
    9193    // Solve with the new weights
     94    fit->useWeight = TRUE; // Reweighted Least Squares
    9295    if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
    9396
    … …  
    157160  }
    158161
    159   // this section calculates the formal error on the weighted fit using the covariance values
    160   // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     162  // this section calculates the formal error on the regular (unweighted) fit using the covariance values
     163  // NOTE 20160929 : use only the unmasked points to calculate the error
    161164  if (data->getError) {
    162     double ax = 0.0, ay = 0.0;
    163     double bx = 0.0, by = 0.0;
    164 
    165     for (i = 0; i < Npoints; i++) {
    166       ax += dpsi_cauchy(points[i].rx / points[i].dX);
    167       ay += dpsi_cauchy(points[i].ry / points[i].dY);
    168 
    169       bx += SQ(points[i].Wx);
    170       by += SQ(points[i].Wy);
    171     }
    172     ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    173     ay /= 1.0 * Npoints;
    174     bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
    175     by /= 1.0 * (Npoints - data->Nterms);
    176  
    177     double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
    178     double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
    179  
    180     double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
    181     double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
    182 
    183     // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
    184     double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    185     double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    186 
    187     fit[0].dRo = sqrt(data->Cov[0][0]);
    188     fit[0].duR = sqrt(data->Cov[1][1]);
    189     fit[0].dDo = sqrt(data->Cov[2][2]);
    190     fit[0].duD = sqrt(data->Cov[3][3]);
    191     fit[0].dp  = sqrt(data->Cov[4][4]);
    192 
    193     fit[0].dRo *= sigma_final_x;
    194     fit[0].duR *= sigma_final_x;
    195     fit[0].dDo *= sigma_final_y;
    196     fit[0].duD *= sigma_final_y;
    197     fit[0].dp  *= sqrt(sigma_final_x * sigma_final_y);
     165    FitAstromResult fitErrors;
     166    FitAstromResultInit (&fitErrors);
     167    fitErrors.useWeight = FALSE;
     168
     169    FitPMandPar_MinChisq(&fitErrors, data, points, Npoints);
     170
     171    // we use the errors from a simple OLS, ignoring masked points
     172    fit[0].dRo = fitErrors.dRo;
     173    fit[0].duR = fitErrors.duR;
     174    fit[0].dDo = fitErrors.dDo;
     175    fit[0].duD = fitErrors.duD;
     176    fit[0].dp  = fitErrors.dp;
    198177  }
    199178
    … …  
    219198    Nfit ++;
    220199
    221     wx = points[i].qx;
    222     wy = points[i].qy;
     200    if (fit->useWeight) {
     201      wx = points[i].qx;
     202      wy = points[i].qy;
     203    } else {
     204      wx = points[i].Qx;
     205      wy = points[i].Qy;
     206    }
    223207
    224208    Wx += wx;

  • branches/czw_branch/20160809/Ohana/src/relastro/src/FitPosPMfixed.c

    r39612 r39731  
    5454  }
    5555
     56  fit->useWeight = FALSE; // Ordinary Least Squares
    5657  if (!FitPosPMfixed_MinChisq (fit, data, points, Npoints)) {
    5758    if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE;
    … …  
    8384 
    8485  // Solve OLS equation 
     86  fit->useWeight = FALSE; // Ordinary Least Squares
    8587  if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
    8688    return(FALSE);
    … …  
    126128
    127129    // Solve with the new weights
     130    fit->useWeight = TRUE; // Reweighted Least Squares
    128131    if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
    129132
    … …  
    193196  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
    194197  if (data->getError) {
    195     double ax = 0.0, ay = 0.0;
    196     double bx = 0.0, by = 0.0;
    197 
    198     for (i = 0; i < Npoints; i++) {
    199       ax += dpsi_cauchy(points[i].rx / points[i].dX);
    200       ay += dpsi_cauchy(points[i].ry / points[i].dY);
    201 
    202       bx += SQ(points[i].Wx);
    203       by += SQ(points[i].Wy);
    204     }
    205     ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    206     ay /= 1.0 * Npoints;
    207     bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
    208     by /= 1.0 * (Npoints - data->Nterms);
    209  
    210     double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
    211     double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
    212  
    213     double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
    214     double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
    215 
    216     // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
    217     double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    218     double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    219 
    220     fit[0].dRo = sqrt(data->Cov[0][0]);
    221     fit[0].dDo = sqrt(data->Cov[1][1]);
    222 
    223     fit[0].dRo *= sigma_final_x;
    224     fit[0].dDo *= sigma_final_y;
     198    FitAstromResult fitErrors;
     199    FitAstromResultInit (&fitErrors);
     200    fitErrors.useWeight = FALSE;
     201
     202    FitPosPMfixed_MinChisq(&fitErrors, data, points, Npoints);
     203
     204    // we use the errors from a simple OLS, ignoring masked points
     205    fit[0].dRo = fitErrors.dRo;
     206    fit[0].duR = fitErrors.duR;
     207    fit[0].dDo = fitErrors.dDo;
     208    fit[0].duD = fitErrors.duD;
     209    fit[0].dp  = fitErrors.dp;
    225210  }
    226211
    … …  
    242227    Nfit ++;
    243228
    244     wx = points[i].qx;
    245     wy = points[i].qy;
     229    if (fit->useWeight) {
     230      wx = points[i].qx;
     231      wy = points[i].qy;
     232    } else {
     233      wx = points[i].Qx;
     234      wy = points[i].Qy;
     235    }
    246236
    247237    Wx += wx;

  • branches/czw_branch/20160809/Ohana/src/relastro/src/ImageOps.c

    r39580 r39731  
    794794   
    795795    // if we are applying the galaxy model, move the reference position...
    796     if (USE_GALAXY_MODEL) {
     796    if (APPLY_PROPER_MOTION) {
    797797      // apply proper-motion from average position to measure epoch:
    798798      float dTime = (measure[0].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years
    … …  
    934934}
    935935
    936 # if (0)
    937 /** lifted from relphot/StarOps.clean_measures */
    938 void FlagOutliers2D(Catalog *catalog);
    939 
    940 // operates on Full values (not tiny)
    941 void FlagOutliers (Catalog *catalog) {
    942 
    943   // XXX FlagOutliers is now just using FlagOutliers2D
    944   FlagOutliers2D(catalog);
    945   return;
    946 
    947   int Ndel, Nave;
    948   off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
    949   double Ns, theta, x, y;
    950   double *R, *D, *dR, *dD;
    951   StatType statsR, statsD;
    952 
    953   Nsecfilt = GetPhotcodeNsecfilt();
    954   assert(catalog[0].Nsecfilt == Nsecfilt);
    955 
    956   if (VERBOSE2) fprintf (stderr, "marking poor measures\n");
    957   Nmax = 0;
    958   for (i = 0; i < catalog[0].Naverage; i++) {
    959     Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
    960   }
    961 
    962   ALLOCATE (R, double, Nmax);
    963   ALLOCATE (D, double, Nmax);
    964   ALLOCATE (dR, double, Nmax);
    965   ALLOCATE (dD, double, Nmax);
    966 
    967   /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    968   TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
    969 
    970   Ns = CLIP_THRESH;
    971   Ndel = Nave = 0;
    972      
    973   /* loop over each object in the catalog */
    974   for (j = 0; j < catalog[0].Naverage; j++) {
    975    
    976     // pointer to this set of measurements
    977     m = catalog[0].average[j].measureOffset;
    978     Measure *measure = &catalog[0].measure[m];
    979 
    980     /* accumulate list of valid measurements */
    981     N = 0;
    982     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    983       // skip measurements based on user selected criteria
    984       if (!MeasFilterTest(&measure[k], FALSE)) continue;
    985       R[N] = measure[k].R;
    986       D[N] = measure[k].D;
    987       dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA);
    988       dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC);
    989       if (isnan(R[N]) || isnan(D[N])) continue;
    990       N++;
    991     }
    992     if (N <= TOOFEW) continue;
    993    
    994     /* 3-sigma clip based on stats of inner 50% */
    995     initstats ("MEAN");
    996     liststats (R, dR, N, &statsR);
    997     liststats (D, dD, N, &statsD);
    998    
    999     statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    1000     statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    1001    
    1002     /* compare per-object distance to this standard deviation, and flag outliers*/
    1003     N = 0;
    1004     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    1005       // reset flag on each invocation
    1006       measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    1007 
    1008       // skip measurements based on user selected criteria
    1009       if (!MeasFilterTest(&measure[k], FALSE)) continue;
    1010      
    1011       x = measure[k].R - statsR.median;
    1012       y = measure[k].D - statsD.median;
    1013       theta = atan2(y,x);
    1014       if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) +
    1015                          SQR(statsD.sigma * Ns * sin(theta)))) {   
    1016         measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    1017         Ndel++;
    1018       }
    1019       N++;
    1020       Nave ++;
    1021     }
    1022 
    1023     // examine results
    1024     // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
    1025   }
    1026  
    1027   if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
    1028   free (R);
    1029   free (dR);
    1030   free (D);
    1031   free (dD);
    1032 }
    1033 
    1034 
    1035 /** an alternative outlier rejection scheme */
    1036 void FlagOutliers2D (Catalog *catalog) {
    1037 
    1038   int Ndel, Nave;
    1039   off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
    1040   double *index;
    1041   double Ns, theta, x, y;
    1042   double *R, *D, *dR, *dD, *d2;
    1043   StatType statsR, statsD;
    1044 
    1045   // XXX we are not going to use this for now
    1046   return;
    1047 
    1048   Nsecfilt = GetPhotcodeNsecfilt();
    1049   assert(catalog[0].Nsecfilt == Nsecfilt);
    1050 
    1051   if (VERBOSE2) fprintf (stderr, "marking poor measures\n");
    1052   Nmax = 0;
    1053   for (i = 0; i < catalog[0].Naverage; i++) {
    1054     Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
    1055   }
    1056 
    1057   ALLOCATE (R, double, Nmax);
    1058   ALLOCATE (D, double, Nmax);
    1059   ALLOCATE (dR, double, Nmax);
    1060   ALLOCATE (dD, double, Nmax);
    1061   ALLOCATE (d2, double, Nmax);
    1062   ALLOCATE (index, double, Nmax);
    1063 
    1064   /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    1065   TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
    1066 
    1067   Ns = CLIP_THRESH;
    1068   Ndel = Nave = 0;
    1069      
    1070   /* loop over each object in the catalog */
    1071   for (j = 0; j < catalog[0].Naverage; j++) {
    1072    
    1073     // pointer to this set of measurements
    1074     m = catalog[0].average[j].measureOffset;
    1075     Measure *measure = &catalog[0].measure[m];
    1076 
    1077     /* accumulate list of valid measurements */
    1078     N = 0;
    1079     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    1080 
    1081       // reset flag on each invocation
    1082       measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    1083      
    1084       // skip measurements based on user selected criteria
    1085       if (!MeasFilterTest(&measure[k], FALSE)) continue;
    1086       R[N] = measure[k].R;
    1087       D[N] = measure[k].D;
    1088       dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA);
    1089       dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC);
    1090       if (isnan(R[N]) || isnan(D[N])) continue;
    1091       N++;
    1092     }
    1093     if (N <= TOOFEW) continue;
    1094    
    1095     /* calculate mean of all points*/
    1096     initstats ("MEAN");
    1097     liststats (R, dR, N, &statsR);
    1098     liststats (D, dD, N, &statsD);
    1099     statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    1100     statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    1101    
    1102     /* calculate deviations of all points*/
    1103     N = 0;
    1104     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    1105       // skip bad measurements
    1106       if (!MeasFilterTest(&measure[k], FALSE)) continue; 
    1107       x = measure[k].R - statsR.median;
    1108       y = measure[k].D - statsD.median;
    1109       theta = atan2(y,x);
    1110       d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
    1111                              SQR(statsD.sigma * Ns * sin(theta)));     
    1112       index[N] = k;
    1113       N++;
    1114     }
    1115    
    1116     // sort d2
    1117     dsortpair(d2, index, N);
    1118     N = (N/2 > (N-1)) ? N/2 : N-1;
    1119 
    1120     // recalculate image center, sigma based on closest 50% of points
    1121     for (k = 0;  k < N; k++) {
    1122       off_t ind = (off_t) index[k];
    1123       R[k] = measure[ind].R;
    1124       D[k] = measure[ind].D;
    1125       dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA);
    1126       dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC);
    1127     }
    1128     liststats (R, dR, N, &statsR);
    1129     liststats (D, dD, N, &statsD);
    1130     statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    1131     statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    1132    
    1133     // use these new statistics to flag outliers
    1134     N = 0;
    1135     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    1136       //skip bad measurements
    1137       if (!MeasFilterTest(&measure[k], FALSE)) continue; 
    1138       x = measure[k].R - statsR.median;
    1139       y = measure[k].D - statsD.median;
    1140       theta = atan2(y,x);
    1141       d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
    1142                              SQR(statsD.sigma * Ns * sin(theta)));     
    1143       if ((d2[N]) > 1) {
    1144         measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    1145         Ndel ++;
    1146       }
    1147       N++;
    1148       Nave++;
    1149     }  // done rejecting outliers
    1150 
    1151     // examine results
    1152     // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
    1153    
    1154   } // done looping over objects
    1155  
    1156   if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
    1157   free (R);
    1158   free (dR);
    1159   free (D);
    1160   free (dD);
    1161   free (d2);
    1162   free (index);
    1163 }
    1164 
    1165 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */
    1166 // we only optionally apply the sigma limit: for object averages, this should not be used
    1167 int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) {
    1168   int found, k;
    1169   long mask;
    1170   PhotCode *code;
    1171   float mag;
    1172 
    1173   if (!finite(measure[0].R) || !finite(measure[0].D)) return FALSE;
    1174   if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
    1175   if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks?
    1176  
    1177   if ((MinBadQF > 0.0) && (isGPC1chip(measure[0].photcode) || isGPC1stack(measure[0].photcode))) {
    1178     if (!isfinite(measure[0].psfQF)) { return FALSE; };
    1179     if (measure[0].psfQF < MinBadQF) { return FALSE; };
    1180   }
    1181 
    1182   /* select measurements by photcode, or equiv photcode, if specified */
    1183   if (NphotcodesKeep > 0) {
    1184     found = FALSE;
    1185     for (k = 0; (k < NphotcodesKeep) && !found; k++) {
    1186       if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
    1187       if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    1188     }
    1189     if (!found) return FALSE;
    1190   }
    1191  
    1192   if (NphotcodesSkip > 0) {
    1193     found = FALSE;
    1194     for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    1195       if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
    1196       if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    1197     }
    1198     if (found) return FALSE;
    1199   } 
    1200  
    1201   /* select measurements by time */
    1202   if (TimeSelect) {
    1203     if (measure[0].t < TSTART) return FALSE;
    1204     if (measure[0].t > TSTOP) return FALSE;
    1205   }
    1206  
    1207   /* select measurements by quality */
    1208   if (PhotFlagSelect) {
    1209     if (PhotFlagBad) {
    1210       mask = PhotFlagBad;
    1211     } else {
    1212       code = GetPhotcodebyCode (measure[0].photcode);
    1213       if (!code) return FALSE;
    1214       mask = code[0].astromBadMask;
    1215     }
    1216     if (mask & measure[0].photFlags) return FALSE;
    1217   }
    1218 
    1219   /* select measurements by measurement error */
    1220   // this is a bit convoluted: applySigmaLim is only TRUE when this function is
    1221   // called by bcatalog.  for UpdateObjects, it is FALSE
    1222   if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
    1223     return FALSE;
    1224   }
    1225  
    1226   /* select measurements by mag limit */
    1227   if (ImagSelect) {
    1228     mag = PhotInstTiny (measure, MAG_CLASS_PSF);
    1229     if (mag < ImagMin || mag > ImagMax) return FALSE;
    1230   }
    1231  
    1232   return TRUE;
    1233 }
    1234 
    1235 # define SUPER_VERBOSE 0
    1236 
    1237 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */
    1238 // we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
    1239 int MeasFilterTest(Measure *measure, int applySigmaLim) {
    1240   int found, k;
    1241   long mask;
    1242   PhotCode *code;
    1243   float mag;
    1244 
    1245   if (!finite(measure[0].R) || !finite(measure[0].D)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 1\n"); return FALSE; };
    1246   if (!finite(measure[0].M)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 2\n"); return FALSE; }; //XXX is this necessary for all relastro tasks?
    1247   if (!finite(measure[0].dM)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 3\n"); return FALSE; }; //XXX is this necessary for all relastro tasks?
    1248  
    1249   /* select measurements by photcode, or equiv photcode, if specified */
    1250   if (NphotcodesKeep > 0) {
    1251     found = FALSE;
    1252     for (k = 0; (k < NphotcodesKeep) && !found; k++) {
    1253       if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
    1254       if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    1255     }
    1256     if (!found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 4\n"); return FALSE; };
    1257   }
    1258  
    1259   if (NphotcodesSkip > 0) {
    1260     found = FALSE;
    1261     for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    1262       if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
    1263       if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    1264     }
    1265     if (found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 5\n"); return FALSE; };
    1266   } 
    1267  
    1268   if ((MinBadQF > 0.0) && isGPC1chip(measure[0].photcode)) {
    1269     if (!isfinite(measure[0].psfQF)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 6\n"); return FALSE; };
    1270     if (measure[0].psfQF < MinBadQF) { if (SUPER_VERBOSE) fprintf (stderr, "filter 7\n"); return FALSE; };
    1271   }
    1272 
    1273   /* select measurements by time */
    1274   if (TimeSelect) {
    1275     if (measure[0].t < TSTART) { if (SUPER_VERBOSE) fprintf (stderr, "filter 8\n"); return FALSE; };
    1276     if (measure[0].t > TSTOP) { if (SUPER_VERBOSE) fprintf (stderr, "filter 9\n"); return FALSE; };
    1277   }
    1278  
    1279   /* select measurements by quality */
    1280   if (PhotFlagSelect) {
    1281     if (PhotFlagBad) {
    1282       mask = PhotFlagBad;
    1283     } else {
    1284       code = GetPhotcodebyCode (measure[0].photcode);
    1285       if (!code) { if (SUPER_VERBOSE) fprintf (stderr, "filter 10\n"); return FALSE; };
    1286       mask = code[0].astromBadMask;
    1287     }
    1288     if (mask & measure[0].photFlags) { if (SUPER_VERBOSE) fprintf (stderr, "filter 11\n"); return FALSE; };
    1289   }
    1290 
    1291   /* select measurements by measurement error */
    1292   if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
    1293     { if (SUPER_VERBOSE) fprintf (stderr, "filter 12\n"); return FALSE; };
    1294   }
    1295  
    1296   /* select measurements by mag limit */
    1297   if (ImagSelect) {
    1298     mag = PhotInst (measure, MAG_CLASS_PSF);
    1299     if (mag < ImagMin || mag > ImagMax) { if (SUPER_VERBOSE) fprintf (stderr, "filter 13\n"); return FALSE; };
    1300   }
    1301  
    1302   return TRUE;
    1303 }
    1304 # endif

  • branches/czw_branch/20160809/Ohana/src/relastro/src/StarMaps.c

    r39457 r39731  
    166166  dLmax = dMmax = 0.0;
    167167
     168  float plateScale;
     169  if (images[N].coords.mosaic) {
     170    // NOTE: for the full pixel to sky plate scale, use this:
     171    // float plateScaleX = 3600.0*images[N].coords.mosaic->cdelt1*images[N].coords.cdelt1;
     172    // float plateScaleY = 3600.0*images[N].coords.mosaic->cdelt2*images[N].coords.cdelt2;
     173
     174    // since we are compare L,M values, just need to compensate for focal plate to sky:
     175    float plateScaleX = 3600.0*fabs(images[N].coords.mosaic->cdelt1);
     176    float plateScaleY = 3600.0*fabs(images[N].coords.mosaic->cdelt2);
     177    plateScale = 0.5*(plateScaleX + plateScaleY);
     178  } else {
     179    // since we are compare L,M values, just need to compensate for arcsec vs degrees:
     180    plateScale = 3600.0;
     181  }
     182
    168183  for (i = 0; i < starmap[N].Npoints; i++) {
    169184
    … …  
    171186    XY_to_LM (&L, &M, starmap[N].points[i].X, starmap[N].points[i].Y, &images[N].coords);
    172187
    173     starmap[N].points[i].dL = starmap[N].points[i].L - L;
    174     starmap[N].points[i].dM = starmap[N].points[i].M - M;
     188    starmap[N].points[i].dL = plateScale*(starmap[N].points[i].L - L);
     189    starmap[N].points[i].dM = plateScale*(starmap[N].points[i].M - M);
    175190
    176191    dLmax = MAX(fabs(starmap[N].points[i].dL), dLmax);

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateChips.c

    r39457 r39731  
    6565    // helper functions that grab the next available chip....
    6666
     67    // allow certain cameras to stay static
     68    if (SKIP_PS1_CHIP  && isGPC1chip (image[i].photcode)) { Nskip ++;  mode[i] = 0; continue; }
     69    if (SKIP_PS1_STACK && isGPC1stack(image[i].photcode)) { Nskip ++;  mode[i] = 0; continue; }
     70    if (SKIP_HSC       && isHSCchip  (image[i].photcode)) { Nskip ++;  mode[i] = 0; continue; }
     71    if (SKIP_CFH       && isCFHchip  (image[i].photcode)) { Nskip ++;  mode[i] = 0; continue; }
     72   
    6773    /* skip all except WRP images */
    6874    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
    … …  
    160166    setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC);
    161167    if (USE_GALAXY_MODEL) {
    162       // XXX DEPRECATE?
     168      // the image calibration was calculated using a galaxy motion model
    163169      image[i].flags |= ID_IMAGE_ASTROM_GMM;
    164170    }
    … …  
    281287    }
    282288
     289    // allow certain cameras to stay static
     290    if (SKIP_PS1_CHIP  && isGPC1chip (image[i].photcode)) { threadinfo->Nskip ++;  threadinfo->mode[i] = 0; continue; }
     291    if (SKIP_PS1_STACK && isGPC1stack(image[i].photcode)) { threadinfo->Nskip ++;  threadinfo->mode[i] = 0; continue; }
     292    if (SKIP_HSC       && isHSCchip  (image[i].photcode)) { threadinfo->Nskip ++;  threadinfo->mode[i] = 0; continue; }
     293    if (SKIP_CFH       && isCFHchip  (image[i].photcode)) { threadinfo->Nskip ++;  threadinfo->mode[i] = 0; continue; }
     294   
    283295    /* skip all except WRP images */
    284296    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateMeasures.c

    r39624 r39731  
    115115
    116116      // only modify the chip coordinates
    117       if (isGPC1chip (measureT->photcode)) {
     117      if (UPDATE_PS1_CHIP_MEASURE && isGPC1chip (measureT->photcode)) {
     118        measureT->R = R;
     119        measureT->D = D;
     120        if (measureB) {
     121          measureB->R = R;
     122          measureB->D = D;
     123        }
     124      }
     125      // only modify the chip coordinates
     126      if (UPDATE_HSC_MEASURE && isGPC1chip (measureT->photcode)) {
     127        measureT->R = R;
     128        measureT->D = D;
     129        if (measureB) {
     130          measureB->R = R;
     131          measureB->D = D;
     132        }
     133      }
     134      // only modify the chip coordinates
     135      if (UPDATE_CFH_MEASURE && isGPC1chip (measureT->photcode)) {
    118136        measureT->R = R;
    119137        measureT->D = D;

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjectOffsets.c

    r39719 r39731  
    197197    if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures");
    198198
     199    if (UPDATE_PS1_CHIP_MEASURE) strextend (&command, "-update-ps1-chip");
     200    if (UPDATE_HSC_MEASURE)      strextend (&command, "-update-hsc");
     201    if (UPDATE_CFH_MEASURE)      strextend (&command, "-update-cfh");
     202
    199203    if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST);
    200204    if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST);
    … …  
    330334    if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures");
    331335
     336    if (UPDATE_PS1_CHIP_MEASURE) strextend (&command, "-update-ps1-chip");
     337    if (UPDATE_HSC_MEASURE)      strextend (&command, "-update-hsc");
     338    if (UPDATE_CFH_MEASURE)      strextend (&command, "-update-cfh");
     339
    332340    if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST);
    333341    if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST);

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjects.c

    r39617 r39731  
    852852  int i;
    853853
    854   if (USE_GALAXY_MODEL) {
     854  if (APPLY_PROPER_MOTION) {
    855855    for (i = 0; i < Nfit; i++) {
    856       fit[i].uR = average->uRgal;
    857       fit[i].uD = average->uDgal;
     856      if (USE_GALAXY_MODEL) {
     857        fit[i].uR = average->uRgal;
     858        fit[i].uD = average->uDgal;
     859      } else {
     860        fit[i].uR = average->uR;
     861        fit[i].uD = average->uD;
     862      }
    858863    }
    859864  } else {

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateStacks.c

    r39457 r39731  
    33// NOTE: we only measure the systematic floor of the astrometric scatter per stack, no change to the calibration
    44int UpdateStacks (Catalog *catalog, int Ncatalog) {
     5
     6  if (SKIP_PS1_STACK) return TRUE;
    57
    68  off_t Nimage;

  • branches/czw_branch/20160809/Ohana/src/relastro/src/args.c

    r39719 r39731  
    116116  }
    117117
     118  // used to decide if changes to the image parameters get applied to the measures
     119  APPLY_OFFSETS = FALSE;
     120  if ((N = get_argument (argc, argv, "-apply-offsets"))) {
     121    remove_argument (N, &argc, argv);
     122    APPLY_OFFSETS = TRUE;
     123  }
     124
     125  APPLY_PROPER_MOTION = FALSE;
     126  if ((N = get_argument (argc, argv, "-apply-proper-motion"))) {
     127    remove_argument (N, &argc, argv);
     128    APPLY_PROPER_MOTION = TRUE;
     129  }
     130
     131  SKIP_PS1_CHIP = FALSE;
     132  if ((N = get_argument (argc, argv, "-skip-ps1-chip"))) {
     133    remove_argument (N, &argc, argv);
     134    SKIP_PS1_CHIP = TRUE;
     135  }
     136  SKIP_PS1_STACK = FALSE;
     137  if ((N = get_argument (argc, argv, "-skip-ps1-stack"))) {
     138    remove_argument (N, &argc, argv);
     139    SKIP_PS1_STACK = TRUE;
     140  }
     141  SKIP_HSC = FALSE;
     142  if ((N = get_argument (argc, argv, "-skip-hsc"))) {
     143    remove_argument (N, &argc, argv);
     144    SKIP_HSC = TRUE;
     145  }
     146  SKIP_CFH = FALSE;
     147  if ((N = get_argument (argc, argv, "-skip-cfh"))) {
     148    remove_argument (N, &argc, argv);
     149    SKIP_CFH = TRUE;
     150  }
     151
     152  UPDATE_PS1_CHIP_MEASURE = FALSE;
     153  if ((N = get_argument (argc, argv, "-update-ps1-chip"))) {
     154    remove_argument (N, &argc, argv);
     155    UPDATE_PS1_CHIP_MEASURE = TRUE;
     156  }
     157  UPDATE_HSC_MEASURE = FALSE;
     158  if ((N = get_argument (argc, argv, "-update-hsc"))) {
     159    remove_argument (N, &argc, argv);
     160    UPDATE_HSC_MEASURE = TRUE;
     161  }
     162  UPDATE_CFH_MEASURE = FALSE;
     163  if ((N = get_argument (argc, argv, "-update-cfh"))) {
     164    remove_argument (N, &argc, argv);
     165    UPDATE_CFH_MEASURE = TRUE;
     166  }
     167  if (RELASTRO_OP == OP_UPDATE_OFFSETS) {
     168    if (!UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) {
     169      fprintf (stderr, "for -update-offsets, need to select at least one of -update-ps1-chip, -update-hsc, -update-cfh\n");
     170      exit (2);
     171    }
     172  }
     173  if ((RELASTRO_OP == OP_PARALLEL_IMAGES) || (RELASTRO_OP == OP_IMAGES)) {
     174    if (APPLY_OFFSETS && !UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) {
     175      fprintf (stderr, "for [-images or -parallel-images] with -apply-offsets, need to select at least one of -update-ps1-chip, -update-hsc, -update-cfh\n");
     176      exit (2);
     177    }
     178  }
     179
    118180  PARALLEL_REGIONS_MANUAL = FALSE;
    119181  if ((N = get_argument (argc, argv, "-parallel-regions"))) {
    … …  
    211273    remove_argument (N, &argc, argv);
    212274    FIT_TARGET = TARGET_MOSAICS;
     275  }
     276  if ((N = get_argument (argc, argv, "-set-chips"))) {
     277    remove_argument (N, &argc, argv);
     278    FIT_TARGET = SET_CHIPS;
    213279  }
    214280
    … …  
    315381    remove_argument (N, &argc, argv);
    316382    MaxDensityUse = TRUE;
    317   }
    318 
    319   APPLY_OFFSETS = FALSE;
    320   if ((N = get_argument (argc, argv, "-apply-offsets"))) {
    321     remove_argument (N, &argc, argv);
    322     APPLY_OFFSETS = TRUE;
    323383  }
    324384
    … …  
    724784    remove_argument (N, &argc, argv);
    725785    CHECK_MEASURE_TO_IMAGE = TRUE;
     786  }
     787
     788  SKIP_PS1_CHIP = FALSE;
     789  if ((N = get_argument (argc, argv, "-skip-ps1-chip"))) {
     790    remove_argument (N, &argc, argv);
     791    SKIP_PS1_CHIP = TRUE;
     792  }
     793  SKIP_PS1_STACK = FALSE;
     794  if ((N = get_argument (argc, argv, "-skip-ps1-stack"))) {
     795    remove_argument (N, &argc, argv);
     796    SKIP_PS1_STACK = TRUE;
     797  }
     798  SKIP_HSC = FALSE;
     799  if ((N = get_argument (argc, argv, "-skip-hsc"))) {
     800    remove_argument (N, &argc, argv);
     801    SKIP_HSC = TRUE;
     802  }
     803  SKIP_CFH = FALSE;
     804  if ((N = get_argument (argc, argv, "-skip-cfh"))) {
     805    remove_argument (N, &argc, argv);
     806    SKIP_CFH = TRUE;
     807  }
     808
     809  UPDATE_PS1_CHIP_MEASURE = FALSE;
     810  if ((N = get_argument (argc, argv, "-update-ps1-chip"))) {
     811    remove_argument (N, &argc, argv);
     812    UPDATE_PS1_CHIP_MEASURE = TRUE;
     813  }
     814  UPDATE_HSC_MEASURE = FALSE;
     815  if ((N = get_argument (argc, argv, "-update-hsc"))) {
     816    remove_argument (N, &argc, argv);
     817    UPDATE_HSC_MEASURE = TRUE;
     818  }
     819  UPDATE_CFH_MEASURE = FALSE;
     820  if ((N = get_argument (argc, argv, "-update-cfh"))) {
     821    remove_argument (N, &argc, argv);
     822    UPDATE_CFH_MEASURE = TRUE;
     823  }
     824  if (RELASTRO_OP == OP_UPDATE_OFFSETS) {
     825    if (!UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) {
     826      fprintf (stderr, "for -update-offsets, need to select at least one of -update-ps1-chip, -update-hsc, -update-cfh\n");
     827      exit (2);
     828    }
    726829  }
    727830

  • branches/czw_branch/20160809/Ohana/src/relastro/src/assign_images.c

    r39457 r39731  
    7272
    7373  for (j = 0; j < Nimage; j++) {
     74   
     75    // allow certain cameras to stay static
     76    if (SKIP_PS1_CHIP  && isGPC1chip (image[j].photcode)) continue;
     77    if (SKIP_PS1_STACK && isGPC1stack(image[j].photcode)) continue;
     78    if (SKIP_HSC       && isHSCchip  (image[j].photcode)) continue;
     79    if (SKIP_CFH       && isCFHchip  (image[j].photcode)) continue;
    7480   
    7581    /* select images by photcode, or equiv photcode, if specified */

  • branches/czw_branch/20160809/Ohana/src/relastro/src/bcatalog.c

    r39587 r39731  
    187187      if (isGPC1warp(catalog[0].measure[offset].photcode)) continue;
    188188
     189      // allow certain cameras to stay static
     190      if (SKIP_PS1_CHIP  && isGPC1chip(catalog[0].measure[offset].photcode)) continue;
     191      if (SKIP_PS1_STACK && isGPC1stack(catalog[0].measure[offset].photcode)) continue;
     192      if (SKIP_HSC       && isHSCchip(catalog[0].measure[offset].photcode)) continue;
     193      if (SKIP_CFH       && isCFHchip(catalog[0].measure[offset].photcode)) continue;
     194
    189195      // filter objects based on user supplied criteria, including SIGMA_LIM
    190196      if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {

  • branches/czw_branch/20160809/Ohana/src/relastro/src/extra.c

    r37261 r39731  
    3939  return FALSE;
    4040}
     41
     42// for now (20160925) I need to identify HSC chips explicitly.  generalize in the future
     43int isHSCchip (int photcode) {
     44
     45  if ((photcode >= 20000) && (photcode <= 20111)) return TRUE; // g-band
     46  if ((photcode >= 21000) && (photcode <= 21111)) return TRUE; // r-band
     47  if ((photcode >= 22000) && (photcode <= 22111)) return TRUE; // i-band
     48  if ((photcode >= 23000) && (photcode <= 23111)) return TRUE; // z-band
     49  if ((photcode >= 24000) && (photcode <= 24111)) return TRUE; // y-band
     50
     51  return FALSE;
     52}
     53
     54// for now (20160925) I need to identify CFH chips explicitly.  generalize in the future
     55int isCFHchip (int photcode) {
     56
     57  if ((photcode >= 100) && (photcode <= 152)) return TRUE; // g-band
     58  if ((photcode >= 200) && (photcode <= 252)) return TRUE; // r-band
     59  if ((photcode >= 300) && (photcode <= 352)) return TRUE; // i-band
     60  if ((photcode >= 400) && (photcode <= 452)) return TRUE; // z-band
     61  if ((photcode >= 500) && (photcode <= 552)) return TRUE; // y-band
     62
     63  return FALSE;
     64}

  • branches/czw_branch/20160809/Ohana/src/relastro/src/launch_region_hosts.c

    r39457 r39731  
    9898        strextend (&command, "-update-mosaics");
    9999        break;
     100      case SET_CHIPS:
     101        strextend (&command, "-set-chips");
     102        break;
    100103      case TARGET_NONE:
    101104        abort();
    … …  
    136139    if (LoopWeight2MASS) {  strextend (&command, "-loop-weights-2mass %s", LoopWeight2MASSstr); }
    137140    if (LoopWeightTycho) {  strextend (&command, "-loop-weights-tycho %s", LoopWeightTychostr); }
     141    if (APPLY_PROPER_MOTION) strextend (&command, "-apply-proper-motion");
    138142
    139143    strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW);

  • branches/czw_branch/20160809/Ohana/src/relastro/src/load_catalogs.c

    r39466 r39731  
    217217    if (PhotFlagBad)         strextend (&command, "+photflagbad %d", PhotFlagBad);
    218218    if (PhotFlagPoor)        strextend (&command, "+photflagpoor %d", PhotFlagPoor);
     219
     220    if (SKIP_PS1_CHIP)       strextend (&command, "-skip-ps1-chip");
     221    if (SKIP_PS1_STACK)      strextend (&command, "-skip-ps1-stack");
     222    if (SKIP_HSC)            strextend (&command, "-skip-hsc");
     223    if (SKIP_CFH)            strextend (&command, "-skip-cfh");
     224
    219225    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    220226

  • branches/czw_branch/20160809/Ohana/src/relastro/src/relastro_images.c

    r39580 r39731  
    9090      break;
    9191
     92    case SET_CHIPS:
     93      // we just want to fit the selected chips to the mean positions
     94      UpdateChips (catalog, Ncatalog, i);   // measure.X,Y -> R,D, fit image.coords
     95      MARKTIME("update chips: %f sec\n", dtime);
     96
     97      break;
     98
    9299    case TARGET_MOSAICS:
    93100      for (i = 0; i < NLOOP; i++) {

  • branches/czw_branch/20160809/Ohana/src/relastro/src/select_images.c

    r39474 r39731  
    121121    }
    122122
     123    // allow certain cameras to stay static
     124    if (SKIP_PS1_CHIP  && isGPC1chip (timage[i].photcode)) continue;
     125    if (SKIP_PS1_STACK && isGPC1stack(timage[i].photcode)) continue;
     126    if (SKIP_HSC       && isHSCchip  (timage[i].photcode)) continue;
     127    if (SKIP_CFH       && isCFHchip  (timage[i].photcode)) continue;
     128   
    123129    /* select images by photcode, or equiv photcode, if specified */
    124130    if (NphotcodesKeep > 0) {
Note: See TracChangeset for help on using the changeset viewer.