IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37036


Ignore:
Timestamp:
Jul 17, 2014, 10:11:42 AM (12 years ago)
Author:
eugene
Message:

adding lensing parameters, new cmf formats, replace measure.dR,dD with measure.R,D

Location:
trunk/Ohana/src/addstar
Files:
24 edited

Legend:

Unmodified
Added
Removed

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

    r35579 r37036  
    4141
    4242typedef struct {
    43   Average average;
    44   Measure measure;
     43  Average  average;
     44  Measure  measure;
     45  Lensing *lensing; // optionally carry out the lensing measurements
    4546  int found;
    4647} Stars;
    … …  
    250251void       update_coords          PROTO((Average *average, Measure *measure, off_t *next));
    251252off_t     *init_measure_links     PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
     253off_t     *init_lensing_links     PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
    252254off_t     *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
    253255off_t      add_meas_link          PROTO((Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE));
    254256off_t      add_miss_link          PROTO((Average *average, off_t *next_miss, off_t Nmissing));
     257off_t      add_lens_link          PROTO((Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING));
    255258off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
     259off_t     *build_lensing_links    PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
    256260Measure   *sort_measure           PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas));
    257261Missing   *sort_missing           PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss));
     262Lensing   *sort_lensing           PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next_lens));
    258263int        ImageOptions           PROTO((AddstarClientOptions *options, Image *images, off_t Nimages));
    259264int        GetFileMode            PROTO((Header *header));

  • trunk/Ohana/src/addstar/src/FilterStars.c

    r34260 r37036  
    9999      stars[N].measure.dFluxKron /= image[0].exptime;
    100100    }
     101    if (!isnan(stars[N].measure.FluxAp)) {
     102      stars[N].measure.FluxAp /= image[0].exptime;
     103    }
     104    if (!isnan(stars[N].measure.dFluxAp)) {
     105      stars[N].measure.dFluxAp /= image[0].exptime;
     106    }
    101107   
    102108    // the external ID is supplied, but do we trust it?

  • trunk/Ohana/src/addstar/src/MatchHeaders.c

    r35416 r37036  
    5959    if (!strcmp (exttype, "PS1_V3")) goto keep;
    6060    if (!strcmp (exttype, "PS1_V4")) goto keep;
     61    if (!strcmp (exttype, "PS1_V5")) goto keep;
    6162    if (!strcmp (exttype, "PS1_SV1")) goto keep;
    6263    if (!strcmp (exttype, "PS1_SV2")) goto keep;
    63     if (!strcmp (exttype, "PS1_DV3")) {
     64    if (!strcmp (exttype, "PS1_SV3")) goto keep;
     65    if (!strcmp (exttype, "PS1_DV3")) goto keep;
     66    if (!strcmp (exttype, "PS1_DV4")) {
    6467      goto keep;
    6568    }

  • trunk/Ohana/src/addstar/src/ReadImageHeader.c

    r36833 r37036  
    4545  image[0].NY = Ny;
    4646
    47   image[0].refColor = NAN;
     47  image[0].refColorBlue = NAN;
     48  image[0].refColorRed = NAN;
    4849
    4950  if (!gfits_scan (header, "TZERO",   "%d",  1, &image[0].tzero) && !ACCEPT_TIME) {
    … …  
    135136
    136137  tmp = 0;
    137   gfits_scan (header, "FLIMIT",   "%lf", 1, &tmp);
     138  // gfits_scan (header, "FLIMIT",   "%lf", 1, &tmp);
     139  gfits_scan (header, "DETEFF.MAGREF", "%lf", 1, &tmp);
    138140  image[0].detection_limit = tmp * 10.0;
    139141

  • trunk/Ohana/src/addstar/src/ReadStarsFITS.c

    r35416 r37036  
    1212Stars     *Convert_PS1_V3         PROTO((FTable *table, unsigned int *nstars));
    1313Stars     *Convert_PS1_V4         PROTO((FTable *table, unsigned int *nstars));
     14Stars     *Convert_PS1_V5         PROTO((FTable *table, unsigned int *nstars));
     15Stars     *Convert_PS1_V5_Lensing PROTO((FTable *table, unsigned int *nstars));
    1416Stars     *Convert_PS1_SV1        PROTO((FTable *table, unsigned int *nstars));
    1517Stars     *Convert_PS1_SV1_Alt    PROTO((FTable *table, unsigned int *nstars));
    1618Stars     *Convert_PS1_SV2        PROTO((FTable *table, unsigned int *nstars));
     19Stars     *Convert_PS1_SV3        PROTO((FTable *table, unsigned int *nstars));
    1720Stars     *Convert_PS1_DV3        PROTO((FTable *table, unsigned int *nstars));
     21Stars     *Convert_PS1_DV4        PROTO((FTable *table, unsigned int *nstars));
    1822
    1923// given a file with the pointer at the start of the table block and the
    … …  
    7377    stars = Convert_PS1_V4 (&table, &Nstars);
    7478  }
     79  if (!strcmp (type, "PS1_V5")) {
     80    if (table.header[0].Naxis[0] == 312) {
     81      stars = Convert_PS1_V5_Lensing (&table, &Nstars);
     82    } else {
     83      stars = Convert_PS1_V5 (&table, &Nstars);
     84    }
     85  }
    7586  if (!strcmp (type, "PS1_SV1")) {
    7687    stars = Convert_PS1_SV1 (&table, &Nstars);
    … …  
    7990    stars = Convert_PS1_SV2 (&table, &Nstars);
    8091  }
     92  if (!strcmp (type, "PS1_SV3")) {
     93    stars = Convert_PS1_SV3 (&table, &Nstars);
     94  }
    8195  if (!strcmp (type, "PS1_DV3")) {
    8296    stars = Convert_PS1_DV3 (&table, &Nstars);
    8397  }
     98  if (!strcmp (type, "PS1_DV4")) {
     99    stars = Convert_PS1_DV4 (&table, &Nstars);
     100  }
    84101  if (stars == NULL) {
    85102    fprintf (stderr, "invalid table type %s\n", type);
    … …  
    91108
    92109  return stars;
     110}
     111
     112float GetFluxFromFluxOrMag (float flux, float mag) {
     113
     114  if (isnan(mag) && isnan(flux)) return NAN;
     115
     116  if (isnan(flux)) return pow(10.0,-0.4*mag);
     117  return flux;
     118}
     119
     120float GetFluxErrFromFluxOrMag (float dFlux, float flux, float dMag) {
     121
     122  if (isnan(dMag) && isnan(dFlux)) return NAN;
     123  if (isnan(dFlux) && isnan(flux)) return NAN;
     124
     125  if (isnan(dFlux)) return (fabs(dMag * flux));
     126  return dFlux;
    93127}
    94128
    … …  
    360394    stars[i].measure.dMcal      = ps1data[i].dMcal;
    361395    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     396    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
    362397                       
    363398    stars[i].measure.Mkron      = NAN; // not provided by PS1_V1:
    … …  
    445480    stars[i].measure.dMcal      = ps1data[i].dMcal;
    446481    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     482    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
    447483                       
    448484    stars[i].measure.Mkron      = NAN; // not provided by PS1_V1_Alt:
    … …  
    522558    stars[i].measure.dMcal      = ps1data[i].dMcal;
    523559    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     560    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
    524561                       
    525562    stars[i].measure.Mkron      = NAN; // not provided by PS1_V2:
    … …  
    595632    stars[i].measure.dMcal      = ps1data[i].dMcal;
    596633    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     634    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
    597635                       
    598636    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    … …  
    600638                       
    601639    // these fluxes are converted from counts to counts/sec in FilterStars.c
    602     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    603     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     640    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     641    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    604642    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    605643    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     644    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (NAN, ps1data[i].Map);
     645    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (NAN, stars[i].measure.FluxAp, stars[i].measure.dMap);
    606646
    607647    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    676716    stars[i].measure.dMcal      = ps1data[i].dMcal;
    677717    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     718    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
    678719                       
    679720    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    … …  
    681722                       
    682723    // these fluxes are converted from counts to counts/sec in FilterStars.c
    683     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    684     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     724    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     725    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    685726    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    686727    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     728    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     729    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
    687730
    688731    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    724767}
    725768
    726 Stars *Convert_PS1_SV1 (FTable *table, unsigned int *nstars) {
     769Stars *Convert_PS1_V5 (FTable *table, unsigned int *nstars) {
    727770
    728771  off_t Nstars;
    … …  
    730773  double ZeroPt;
    731774  Stars *stars;
    732   CMF_PS1_SV1 *ps1data;
    733 
    734   if (table[0].header[0].Naxis[0] == 196) {
    735     stars = Convert_PS1_SV1_Alt (table, nstars);
    736     return (stars);
    737   }
    738 
    739   ps1data = gfits_table_get_CMF_PS1_SV1 (table, &Nstars, NULL);
     775  CMF_PS1_V5 *ps1data;
     776
     777  ps1data = gfits_table_get_CMF_PS1_V5 (table, &Nstars, NULL);
    740778  if (!ps1data) {
    741779    fprintf (stderr, "skipping inconsistent entry\n");
    … …  
    763801    stars[i].measure.dMcal      = ps1data[i].dMcal;
    764802    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     803    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
    765804                       
    766805    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    767806    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
    768 
     807                       
    769808    // these fluxes are converted from counts to counts/sec in FilterStars.c
    770     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    771     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     809    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     810    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    772811    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    773812    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     813    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     814    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
    774815
    775816    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    799840
    800841    // the Average fields and the following Measure fields are set in FilterStars after
    801     // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID,
    802     // averef is set in find_matches, dbFlags is zero on ingest.
     842    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID.
     843
     844    // averef is set in find_matches
     845
     846    // dbFlags is zero on ingest.
    803847
    804848    // the following fields are currently not being set anywhere: t_msec
    … …  
    808852}
    809853
    810 Stars *Convert_PS1_SV1_Alt (FTable *table, unsigned int *nstars) {
     854Stars *Convert_PS1_V5_Lensing (FTable *table, unsigned int *nstars) {
    811855
    812856  off_t Nstars;
    … …  
    814858  double ZeroPt;
    815859  Stars *stars;
    816   CMF_PS1_SV1 *ps1data;
    817 
    818   // some test output files were produced called CMF_PS1_SV1 but with mismatch byte boundaries
    819 
    820   ps1data = gfits_table_get_CMF_PS1_SV1_Alt (table, &Nstars, NULL);
     860  CMF_PS1_V5_Lensing *ps1data;
     861
     862  ps1data = gfits_table_get_CMF_PS1_V5_Lensing (table, &Nstars, NULL);
    821863  if (!ps1data) {
    822864    fprintf (stderr, "skipping inconsistent entry\n");
    … …  
    844886    stars[i].measure.dMcal      = ps1data[i].dMcal;
    845887    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     888    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
    846889                       
    847890    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    848891    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
    849 
     892                       
    850893    // these fluxes are converted from counts to counts/sec in FilterStars.c
    851     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    852     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     894    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     895    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    853896    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    854897    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     898    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     899    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
     900
     901    stars[i].measure.Sky        = ps1data[i].sky;
     902    stars[i].measure.dSky       = ps1data[i].dSky;
     903                       
     904    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
     905    stars[i].measure.psfQF      = ps1data[i].psfQF;
     906    stars[i].measure.psfQFperf  = ps1data[i].psfQFperf;
     907
     908    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
     909    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
     910    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
     911    stars[i].measure.extNsigma  = ps1data[i].extNsigma;
     912
     913    stars[i].measure.FWx        = ToShortPixels(ps1data[i].fx);
     914    stars[i].measure.FWy        = ToShortPixels(ps1data[i].fy);
     915    stars[i].measure.theta      = ToShortDegrees(ps1data[i].df);
     916
     917    stars[i].measure.Mxx        = ToShortPixels(ps1data[i].Mxx);
     918    stars[i].measure.Mxy        = ToShortPixels(ps1data[i].Mxy);
     919    stars[i].measure.Myy        = ToShortPixels(ps1data[i].Myy);
     920                       
     921    stars[i].measure.photFlags  = ps1data[i].flags;
     922
     923    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
     924    stars[i].measure.detID      = ps1data[i].detID;
     925
     926    ALLOCATE (stars[i].lensing, Lensing, 1);
     927    dvo_lensing_init (stars[i].lensing);
     928
     929    stars[i].lensing->X11_sm_obj  = ps1data[i].X11_sm_obj;
     930    stars[i].lensing->X12_sm_obj  = ps1data[i].X12_sm_obj;
     931    stars[i].lensing->X22_sm_obj  = ps1data[i].X22_sm_obj;
     932    stars[i].lensing->E1_sm_obj   = ps1data[i].E1_sm_obj;
     933    stars[i].lensing->E2_sm_obj   = ps1data[i].E2_sm_obj;
     934
     935    stars[i].lensing->X11_sh_obj  = ps1data[i].X11_sh_obj;
     936    stars[i].lensing->X12_sh_obj  = ps1data[i].X12_sh_obj;
     937    stars[i].lensing->X22_sh_obj  = ps1data[i].X22_sh_obj;
     938    stars[i].lensing->E1_sh_obj   = ps1data[i].E1_sh_obj;
     939    stars[i].lensing->E2_sh_obj   = ps1data[i].E2_sh_obj;
     940
     941    stars[i].lensing->X11_sm_psf  = ps1data[i].X11_sm_psf;
     942    stars[i].lensing->X12_sm_psf  = ps1data[i].X12_sm_psf;
     943    stars[i].lensing->X22_sm_psf  = ps1data[i].X22_sm_psf;
     944    stars[i].lensing->E1_sm_psf   = ps1data[i].E1_sm_psf;
     945    stars[i].lensing->E2_sm_psf   = ps1data[i].E2_sm_psf;
     946
     947    stars[i].lensing->X11_sh_psf  = ps1data[i].X11_sh_psf;
     948    stars[i].lensing->X12_sh_psf  = ps1data[i].X12_sh_psf;
     949    stars[i].lensing->X22_sh_psf  = ps1data[i].X22_sh_psf;
     950    stars[i].lensing->E1_sh_psf   = ps1data[i].E1_sh_psf;
     951    stars[i].lensing->E2_sh_psf   = ps1data[i].E2_sh_psf;
     952
     953    // stars[i].lensing->F_ApR5    = ps1data[i].F_ApR5;
     954    // stars[i].lensing->dF_ApR5   = ps1data[i].dF_ApR5;
     955    // stars[i].lensing->sF_ApR5   = ps1data[i].sF_ApR5;
     956    // stars[i].lensing->fF_ApR5   = ps1data[i].fF_ApR5;
     957    //
     958    // stars[i].lensing->F_ApR6    = ps1data[i].F_ApR6;
     959    // stars[i].lensing->dF_ApR6   = ps1data[i].dF_ApR6;
     960    // stars[i].lensing->sF_ApR6   = ps1data[i].sF_ApR6;
     961    // stars[i].lensing->fF_ApR6   = ps1data[i].fF_ApR6;
     962
     963    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
     964    stars[i].lensing->detID      = ps1data[i].detID;
     965
     966    // the Average fields and the following Measure fields are set in FilterStars after
     967    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID.
     968
     969    // averef is set in find_matches
     970
     971    // dbFlags is zero on ingest.
     972
     973    // the following fields are currently not being set anywhere: t_msec
     974  }   
     975  *nstars = Nstars;
     976  return (stars);
     977}
     978
     979Stars *Convert_PS1_SV1 (FTable *table, unsigned int *nstars) {
     980
     981  off_t Nstars;
     982  unsigned int i;
     983  double ZeroPt;
     984  Stars *stars;
     985  CMF_PS1_SV1 *ps1data;
     986
     987  if (table[0].header[0].Naxis[0] == 196) {
     988    stars = Convert_PS1_SV1_Alt (table, nstars);
     989    return (stars);
     990  }
     991
     992  ps1data = gfits_table_get_CMF_PS1_SV1 (table, &Nstars, NULL);
     993  if (!ps1data) {
     994    fprintf (stderr, "skipping inconsistent entry\n");
     995    return (NULL);
     996  }
     997  ZeroPt = GetZeroPoint();
     998
     999  ALLOCATE (stars, Stars, Nstars);
     1000  for (i = 0; i < Nstars; i++) {
     1001    InitStar (&stars[i]);
     1002    stars[i].measure.Xccd       = ps1data[i].X;
     1003    stars[i].measure.Yccd       = ps1data[i].Y;
     1004    stars[i].measure.dXccd      = ToShortPixels(ps1data[i].dX);
     1005    stars[i].measure.dYccd      = ToShortPixels(ps1data[i].dY);
     1006
     1007    stars[i].measure.posangle   = ToShortDegrees(ps1data[i].posangle);
     1008    stars[i].measure.pltscale   = ps1data[i].pltscale;
     1009
     1010    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
     1011      stars[i].measure.M      = NAN;
     1012    } else {
     1013      stars[i].measure.M      = ps1data[i].M + ZeroPt;
     1014    }
     1015    stars[i].measure.dM         = ps1data[i].dM;
     1016    stars[i].measure.dMcal      = ps1data[i].dMcal;
     1017    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1018    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
     1019                       
     1020    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
     1021    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
     1022
     1023    // these fluxes are converted from counts to counts/sec in FilterStars.c
     1024    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1025    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
     1026    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
     1027    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1028    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (NAN, ps1data[i].Map);
     1029    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (NAN, stars[i].measure.FluxAp, stars[i].measure.dMap);
     1030
     1031    stars[i].measure.Sky        = ps1data[i].sky;
     1032    stars[i].measure.dSky       = ps1data[i].dSky;
     1033                       
     1034    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
     1035    stars[i].measure.psfQF      = ps1data[i].psfQF;
     1036    stars[i].measure.psfQFperf  = ps1data[i].psfQFperf;
     1037
     1038    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
     1039    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
     1040    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
     1041    stars[i].measure.extNsigma  = ps1data[i].extNsigma;
     1042
     1043    stars[i].measure.FWx        = ToShortPixels(ps1data[i].fx);
     1044    stars[i].measure.FWy        = ToShortPixels(ps1data[i].fy);
     1045    stars[i].measure.theta      = ToShortDegrees(ps1data[i].df);
     1046
     1047    stars[i].measure.Mxx        = ToShortPixels(ps1data[i].Mxx);
     1048    stars[i].measure.Mxy        = ToShortPixels(ps1data[i].Mxy);
     1049    stars[i].measure.Myy        = ToShortPixels(ps1data[i].Myy);
     1050                       
     1051    stars[i].measure.photFlags  = ps1data[i].flags;
     1052
     1053    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
     1054    stars[i].measure.detID      = ps1data[i].detID;
     1055
     1056    // the Average fields and the following Measure fields are set in FilterStars after
     1057    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID,
     1058    // averef is set in find_matches, dbFlags is zero on ingest.
     1059
     1060    // the following fields are currently not being set anywhere: t_msec
     1061  }   
     1062  *nstars = Nstars;
     1063  return (stars);
     1064}
     1065
     1066Stars *Convert_PS1_SV1_Alt (FTable *table, unsigned int *nstars) {
     1067
     1068  off_t Nstars;
     1069  unsigned int i;
     1070  double ZeroPt;
     1071  Stars *stars;
     1072  CMF_PS1_SV1 *ps1data;
     1073
     1074  // some test output files were produced called CMF_PS1_SV1 but with mismatch byte boundaries
     1075
     1076  ps1data = gfits_table_get_CMF_PS1_SV1_Alt (table, &Nstars, NULL);
     1077  if (!ps1data) {
     1078    fprintf (stderr, "skipping inconsistent entry\n");
     1079    return (NULL);
     1080  }
     1081  ZeroPt = GetZeroPoint();
     1082
     1083  ALLOCATE (stars, Stars, Nstars);
     1084  for (i = 0; i < Nstars; i++) {
     1085    InitStar (&stars[i]);
     1086    stars[i].measure.Xccd       = ps1data[i].X;
     1087    stars[i].measure.Yccd       = ps1data[i].Y;
     1088    stars[i].measure.dXccd      = ToShortPixels(ps1data[i].dX);
     1089    stars[i].measure.dYccd      = ToShortPixels(ps1data[i].dY);
     1090
     1091    stars[i].measure.posangle   = ToShortDegrees(ps1data[i].posangle);
     1092    stars[i].measure.pltscale   = ps1data[i].pltscale;
     1093
     1094    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
     1095      stars[i].measure.M      = NAN;
     1096    } else {
     1097      stars[i].measure.M      = ps1data[i].M + ZeroPt;
     1098    }
     1099    stars[i].measure.dM         = ps1data[i].dM;
     1100    stars[i].measure.dMcal      = ps1data[i].dMcal;
     1101    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1102    stars[i].measure.dMap       = ps1data[i].dM; // a proxy measure
     1103                       
     1104    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
     1105    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
     1106
     1107    // these fluxes are converted from counts to counts/sec in FilterStars.c
     1108    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1109    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
     1110    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
     1111    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1112    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (NAN, ps1data[i].Map);
     1113    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (NAN, stars[i].measure.FluxAp, stars[i].measure.dMap);
    8551114
    8561115    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    9211180    stars[i].measure.dMcal      = ps1data[i].dMcal;
    9221181    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1182    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
    9231183                       
    9241184    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    … …  
    9261186
    9271187    // these fluxes are converted from counts to counts/sec in FilterStars.c
    928     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    929     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     1188    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1189    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    9301190    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    9311191    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1192    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     1193    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
    9321194
    9331195    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    9661228}
    9671229
    968 Stars *Convert_PS1_DV3 (FTable *table, unsigned int *nstars) {
     1230Stars *Convert_PS1_SV3 (FTable *table, unsigned int *nstars) {
    9691231
    9701232  off_t Nstars;
    … …  
    9721234  double ZeroPt;
    9731235  Stars *stars;
    974   CMF_PS1_DV3 *ps1data;
    975 
    976   ps1data = gfits_table_get_CMF_PS1_DV3 (table, &Nstars, NULL);
     1236  CMF_PS1_SV3 *ps1data;
     1237
     1238  ps1data = gfits_table_get_CMF_PS1_SV3 (table, &Nstars, NULL);
    9771239  if (!ps1data) {
    9781240    fprintf (stderr, "skipping inconsistent entry\n");
    … …  
    10001262    stars[i].measure.dMcal      = ps1data[i].dMcal;
    10011263    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1264    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
    10021265                       
    10031266    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
    10041267    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
    1005                        
     1268
    10061269    // these fluxes are converted from counts to counts/sec in FilterStars.c
    1007     stars[i].measure.FluxPSF    = ps1data[i].Flux;
    1008     stars[i].measure.dFluxPSF   = ps1data[i].dFlux;
     1270    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1271    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
    10091272    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
    10101273    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1274    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     1275    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
     1276
     1277    stars[i].measure.Sky        = ps1data[i].sky;
     1278    stars[i].measure.dSky       = ps1data[i].dSky;
     1279                       
     1280    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
     1281    stars[i].measure.psfQF      = ps1data[i].psfQF;
     1282    stars[i].measure.psfQFperf  = ps1data[i].psfQFperf;
     1283
     1284    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
     1285    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
     1286    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
     1287    stars[i].measure.extNsigma  = ps1data[i].extNsigma;
     1288
     1289    stars[i].measure.FWx        = ToShortPixels(ps1data[i].fx);
     1290    stars[i].measure.FWy        = ToShortPixels(ps1data[i].fy);
     1291    stars[i].measure.theta      = ToShortDegrees(ps1data[i].df);
     1292
     1293    stars[i].measure.Mxx        = ToShortPixels(ps1data[i].Mxx);
     1294    stars[i].measure.Mxy        = ToShortPixels(ps1data[i].Mxy);
     1295    stars[i].measure.Myy        = ToShortPixels(ps1data[i].Myy);
     1296                       
     1297    stars[i].measure.photFlags  = ps1data[i].flags;
     1298
     1299    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
     1300    stars[i].measure.detID      = ps1data[i].detID;
     1301
     1302    // the Average fields and the following Measure fields are set in FilterStars after
     1303    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID,
     1304    // averef is set in find_matches, dbFlags is zero on ingest.
     1305
     1306    // the following fields are currently not being set anywhere: t_msec
     1307  }   
     1308  *nstars = Nstars;
     1309  return (stars);
     1310}
     1311
     1312Stars *Convert_PS1_DV3 (FTable *table, unsigned int *nstars) {
     1313
     1314  off_t Nstars;
     1315  unsigned int i;
     1316  double ZeroPt;
     1317  Stars *stars;
     1318  CMF_PS1_DV3 *ps1data;
     1319
     1320  ps1data = gfits_table_get_CMF_PS1_DV3 (table, &Nstars, NULL);
     1321  if (!ps1data) {
     1322    fprintf (stderr, "skipping inconsistent entry\n");
     1323    return (NULL);
     1324  }
     1325  ZeroPt = GetZeroPoint();
     1326
     1327  ALLOCATE (stars, Stars, Nstars);
     1328  for (i = 0; i < Nstars; i++) {
     1329    InitStar (&stars[i]);
     1330    stars[i].measure.Xccd       = ps1data[i].X;
     1331    stars[i].measure.Yccd       = ps1data[i].Y;
     1332    stars[i].measure.dXccd      = ToShortPixels(ps1data[i].dX);
     1333    stars[i].measure.dYccd      = ToShortPixels(ps1data[i].dY);
     1334
     1335    stars[i].measure.posangle   = ToShortDegrees(ps1data[i].posangle);
     1336    stars[i].measure.pltscale   = ps1data[i].pltscale;
     1337
     1338    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
     1339      stars[i].measure.M      = NAN;
     1340    } else {
     1341      stars[i].measure.M      = ps1data[i].M + ZeroPt;
     1342    }
     1343    stars[i].measure.dM         = ps1data[i].dM;
     1344    stars[i].measure.dMcal      = ps1data[i].dMcal;
     1345    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1346    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
     1347                       
     1348    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
     1349    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
     1350                       
     1351    // these fluxes are converted from counts to counts/sec in FilterStars.c
     1352    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1353    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
     1354    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
     1355    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1356    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     1357    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
    10111358
    10121359    stars[i].measure.Sky        = ps1data[i].sky;
    … …  
    10471394}
    10481395
    1049 
    1050 
     1396Stars *Convert_PS1_DV4 (FTable *table, unsigned int *nstars) {
     1397
     1398  off_t Nstars;
     1399  unsigned int i;
     1400  double ZeroPt;
     1401  Stars *stars;
     1402  CMF_PS1_DV4 *ps1data;
     1403
     1404  ps1data = gfits_table_get_CMF_PS1_DV4 (table, &Nstars, NULL);
     1405  if (!ps1data) {
     1406    fprintf (stderr, "skipping inconsistent entry\n");
     1407    return (NULL);
     1408  }
     1409  ZeroPt = GetZeroPoint();
     1410
     1411  fprintf (stderr, "WARNING: Convert_PS1_DV4 not yet updated to match real format\n");
     1412
     1413  ALLOCATE (stars, Stars, Nstars);
     1414  for (i = 0; i < Nstars; i++) {
     1415    InitStar (&stars[i]);
     1416    stars[i].measure.Xccd       = ps1data[i].X;
     1417    stars[i].measure.Yccd       = ps1data[i].Y;
     1418    stars[i].measure.dXccd      = ToShortPixels(ps1data[i].dX);
     1419    stars[i].measure.dYccd      = ToShortPixels(ps1data[i].dY);
     1420
     1421    stars[i].measure.posangle   = ToShortDegrees(ps1data[i].posangle);
     1422    stars[i].measure.pltscale   = ps1data[i].pltscale;
     1423
     1424    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
     1425      stars[i].measure.M      = NAN;
     1426    } else {
     1427      stars[i].measure.M      = ps1data[i].M + ZeroPt;
     1428    }
     1429    stars[i].measure.dM         = ps1data[i].dM;
     1430    stars[i].measure.dMcal      = ps1data[i].dMcal;
     1431    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
     1432    stars[i].measure.dMap       = (ps1data[i].apFlux > 0.0) ? fabs(ps1data[i].apFluxErr / ps1data[i].apFlux) : NAN;
     1433                       
     1434    stars[i].measure.Mkron      = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN;
     1435    stars[i].measure.dMkron     = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN;
     1436                       
     1437    // these fluxes are converted from counts to counts/sec in FilterStars.c
     1438    stars[i].measure.FluxPSF    = GetFluxFromFluxOrMag (ps1data[i].Flux, ps1data[i].M);
     1439    stars[i].measure.dFluxPSF   = GetFluxErrFromFluxOrMag (ps1data[i].dFlux, stars[i].measure.FluxPSF, ps1data[i].dM);
     1440    stars[i].measure.FluxKron   = ps1data[i].kronFlux;
     1441    stars[i].measure.dFluxKron  = ps1data[i].kronFluxErr;
     1442    stars[i].measure.FluxAp     = GetFluxFromFluxOrMag (ps1data[i].apFlux, ps1data[i].Map);
     1443    stars[i].measure.dFluxAp    = GetFluxErrFromFluxOrMag (ps1data[i].apFluxErr, stars[i].measure.FluxAp, stars[i].measure.dMap);
     1444
     1445    stars[i].measure.Sky        = ps1data[i].sky;
     1446    stars[i].measure.dSky       = ps1data[i].dSky;
     1447                       
     1448    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
     1449    stars[i].measure.psfQF      = ps1data[i].psfQF;
     1450    stars[i].measure.psfQFperf  = ps1data[i].psfQFperf;
     1451    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
     1452    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
     1453    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
     1454    stars[i].measure.extNsigma  = ps1data[i].extNsigma;
     1455
     1456    stars[i].measure.FWx        = ToShortPixels(ps1data[i].fx);
     1457    stars[i].measure.FWy        = ToShortPixels(ps1data[i].fy);
     1458    stars[i].measure.theta      = ToShortDegrees(ps1data[i].df);
     1459
     1460    stars[i].measure.Mxx        = ToShortPixels(ps1data[i].Mxx);
     1461    stars[i].measure.Mxy        = ToShortPixels(ps1data[i].Mxy);
     1462    stars[i].measure.Myy        = ToShortPixels(ps1data[i].Myy);
     1463                       
     1464    stars[i].measure.photFlags  = ps1data[i].flags;
     1465
     1466    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
     1467    stars[i].measure.detID      = ps1data[i].detID;
     1468
     1469    // the Average fields and the following Measure fields are set in FilterStars after
     1470    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID.
     1471
     1472    // averef is set in find_matches
     1473
     1474    // dbFlags is zero on ingest.
     1475
     1476    // the following fields are currently not being set anywhere: t_msec
     1477  }   
     1478  *nstars = Nstars;
     1479  return (stars);
     1480}
     1481
     1482
     1483

  • trunk/Ohana/src/addstar/src/SEDfit.c

    r34405 r37036  
    162162    }
    163163
     164    double R = incat[0].average[i].R;
     165    double D = incat[0].average[i].D;
     166
    164167    for (j = 0; valid && (j < Nmodel); j++) {
    165168      n = modelRow[j];
    166169      dvo_measure_init (&outcat[0].measure[Nmeas]);
    167       outcat[0].measure[Nmeas].dR        = 0.0;
    168       outcat[0].measure[Nmeas].dD        = 0.0;
     170      outcat[0].measure[Nmeas].R         = R;
     171      outcat[0].measure[Nmeas].D         = D;
    169172      outcat[0].measure[Nmeas].M         = table[0].row[minFit.row][0].mags[n] + minFit.Md;
    170173      outcat[0].measure[Nmeas].dM        = 0.0;

  • trunk/Ohana/src/addstar/src/StarOps.c

    r34260 r37036  
    88    star[0].found = -1; // found == -1 -> not yet found (use enums?)
    99
     10    star[0].lensing = NULL; // we only populate this if needed
    1011    return TRUE;
    1112}

  • trunk/Ohana/src/addstar/src/addstar.c

    r34405 r37036  
    77  int Nmatch, status, loadObjects;
    88  off_t i, Nimages;
    9   off_t Naverage, Nmeasure;
     9  off_t Naverage, Nmeasure, Nlensing;
    1010  Stars *stars, **subset;
    1111  Image *images;
    … …  
    100100
    101101  /* match stars to existing catalog data (or otherwise manipulate catalog data) */
    102   Nmatch = Naverage = Nmeasure = 0;
     102  Nmatch = Naverage = Nmeasure = Nlensing = 0;
    103103  for (i = 0; loadObjects && (i < skylist[0].Nregions); i++) {
    104104
    … …  
    107107    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    108108    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    109     catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     109    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING;
    110110    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    111111
    … …  
    182182    Naverage += catalog.Naverage;
    183183    Nmeasure += catalog.Nmeasure;
     184    Nlensing += catalog.Nlensing;
    184185
    185186    // write out catalog, if appropriate
    … …  
    187188      SetProtect (TRUE);
    188189      if (options.update) {
    189         catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     190        catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING;
    190191        dvo_catalog_update (&catalog, VERBOSE);
    191192      } else {
    … …  
    231232  gettimeofday (&stop, NULL);
    232233  dtime = DTIME (stop, start);
    233   fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for %5d stars (%5d matches), "OFF_T_FMT" average, "OFF_T_FMT" measure\n", dtime, Nstars, Nmatch,  Naverage,  Nmeasure);
     234  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for %5d stars (%5d matches), "OFF_T_FMT" average, "OFF_T_FMT" measure, "OFF_T_FMT" lensing\n", dtime, Nstars, Nmatch,  Naverage,  Nmeasure, Nlensing);
    234235
    235236  exit (0);

  • trunk/Ohana/src/addstar/src/build_links.c

    r27435 r37036  
    225225}
    226226
     227/*******************************************************************************************/
     228
     229/* build the initial links assuming the table is sorted,
     230   not partial, and has a correct set of average[].lensingOffset,Nlensing values */
     231off_t *init_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
     232
     233  off_t i, j, N;
     234  off_t *next_lens;
     235
     236  N = 0;
     237
     238  ALLOCATE (next_lens, off_t, Nlensing);
     239  for (i = 0; i < Naverage; i++) {
     240    if (!average[i].Nlensing) continue;
     241    for (j = 0; j < average[i].Nlensing - 1; j++, N++) {
     242      next_lens[N] = N + 1;
     243      if (N >= Nlensing) {
     244        fprintf (stderr, "WARNING: N out of bounds (1)\n");
     245      }
     246    }
     247    next_lens[N] = -1;
     248    if (N >= Nlensing) {
     249      fprintf (stderr, "WARNING: N out of bounds (2)\n");
     250    }
     251
     252    if (N >= Nlensing) {
     253      fprintf (stderr, "overflow in init_lensing_links\n");
     254      abort ();
     255    }
     256    N++;
     257  }
     258  return (next_lens);
     259}
     260
     261/* construct lensing links which are valid FOR THIS LOAD
     262 * - if we have a full load, we will get links which can
     263 *   be used by other programs (eg, relphot, etc)
     264 * - if we have a partial load, the links are only valid
     265 *   for that partial load
     266 */
     267
     268off_t *build_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
     269
     270  off_t i, m, k, Nm, averef;
     271  off_t *next_lens;
     272
     273  ALLOCATE (next_lens, off_t, Nlensing);
     274
     275  /* reset the Nm, offset values for average */
     276  for (i = 0; i < Naverage; i++) {
     277    average[i].lensingOffset = -1;
     278    average[i].Nlensing     =  0;
     279  }
     280
     281  for (Nm = 0; Nm < Nlensing; Nm++) {
     282    averef = lensing[Nm].averef;
     283    m = average[averef].lensingOffset; 
     284    next_lens[Nm] = -1;
     285
     286    if (m == -1) { /* no links yet for source */
     287      average[averef].lensingOffset = Nm;
     288      average[averef].Nlensing     = 1;
     289      continue;
     290    }
     291
     292    for (k = 0; next_lens[m] != -1; k++) {
     293      m = next_lens[m];
     294      if (m >= Nlensing) {
     295        fprintf (stderr, "WARNING: m out of bounds (1)\n");
     296      }
     297    }
     298
     299    average[averef].Nlensing = k + 2;
     300    next_lens[m] = Nm;
     301    if (m >= Nlensing) {
     302      fprintf (stderr, "WARNING: m out of bounds (2)\n");
     303    }
     304  }
     305  return (next_lens);
     306}
     307
     308/* average[].lensingOffset, average[].Nlensing are valid within an addstar run */
     309off_t add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) {
     310
     311  off_t k, m;
     312
     313  /* if we have trouble, check validity of next_lens[m] : m < Nlensing */
     314  m = average[0].lensingOffset; 
     315
     316  for (k = 0; k < average[0].Nlensing - 1; k++)  {
     317    m = next_lens[m];
     318    if (m >= NLENSING) {
     319      fprintf (stderr, "WARNING: m out of bounds (3)\n");
     320    }
     321  }
     322
     323  /* set up references */
     324  next_lens[Nlensing] = -1;
     325  if (Nlensing >= NLENSING) {
     326    fprintf (stderr, "WARNING: Nlensing out of bounds (1)\n");
     327  }
     328
     329  if (m == -1) {
     330    average[0].lensingOffset = Nlensing;
     331  } else {
     332    next_lens[m] = Nlensing;
     333    if (m >= NLENSING) {
     334      fprintf (stderr, "WARNING: m out of bounds (4)\n");
     335    }
     336  }
     337
     338  return (TRUE);
     339}
     340
     341Lensing *sort_lensing (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next_lens) {
     342
     343  off_t i, k, n, N;
     344  Lensing *tmplensing;
     345
     346  /* fix order of Lensing (memory intensive, but fast) */
     347  N = 0;
     348  ALLOCATE (tmplensing, Lensing, Nlensing);
     349  for (i = 0; i < Naverage; i++) {
     350    if (!average[i].Nlensing) continue;
     351    n = average[i].lensingOffset;
     352    average[i].lensingOffset = N;
     353    for (k = 0; k < average[i].Nlensing; k++, N++) {
     354      if (n == -1) abort();
     355      tmplensing[N] = lensing[n];
     356      if (lensing[n].averef != i) abort();
     357      tmplensing[N].averef = i;
     358      n = next_lens[n];
     359    }
     360  }
     361  free (lensing);
     362  return (tmplensing);
     363}
     364
     365

  • trunk/Ohana/src/addstar/src/fakeimage.c

    r36485 r37036  
    66  double pltscale, pixscale;
    77  double Rmin, Rmax, Dmin, Dmax;
    8   // double Xmin, Xmax, Ymin, Ymax;
    98  double dX, dY, r, d;
    109  char chipname[80], chipdata[256], name[80];
    … …  
    168167  image[0].cerror = 0.0;
    169168   
    170   // RD_to_XY (&Xmax, &Ymax, Rmax, Dmax, MOSAIC);
    171   // RD_to_XY (&Xmin, &Ymin, Rmin, Dmin, MOSAIC);
    172169  image[0].NX = Rmax - Rmin;
    173170  image[0].NY = Dmax - Dmin;

  • trunk/Ohana/src/addstar/src/find_matches.c

    r36485 r37036  
    1414
    1515  if (NSTAR_GROUP <= 0) {
    16       fprintf (stderr, "ERROR: NSTAR_GROUP NOT SET!\n");
    17       exit (1);
     16    fprintf (stderr, "ERROR: NSTAR_GROUP NOT SET!\n");
     17    exit (1);
    1818  }
    1919
    … …  
    3636  ALLOCATE (Y2, double, NAVE);
    3737  ALLOCATE (N2, off_t,  NAVE);
    38   ALLOCATE (catalog[0].found, off_t, NAVE);
     38  ALLOCATE (catalog[0].found_t, off_t, NAVE);
    3939  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    4040
    … …  
    7777  if (Nstars < 1) {
    7878    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
    79     free (catalog[0].found);
    8079    free (X1);
    8180    free (Y1);
    … …  
    9291    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    9392    N2[i] = i;
    94     catalog[0].found[N2[i]] = -1;
     93    catalog[0].found_t[N2[i]] = -1;
    9594  }
    9695  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    … …  
    160159
    161160      // the following measure elements cannot be set until here:
    162       catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    163       catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N].average.D);
     161      catalog[0].measure[Nmeas].R        = stars[N].average.R;
     162      catalog[0].measure[Nmeas].D        = stars[N].average.D;
    164163      catalog[0].measure[Nmeas].dbFlags  = 0;
    165164      catalog[0].measure[Nmeas].averef   = n; // this must be an absolute sequence number, if partial average is loaded
    … …  
    167166      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    168167
    169       // rationalize dR:
    170       if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
    171           // average on high end of boundary, move star up
    172           stars[N].average.R += 360.0;
    173           catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    174       }
    175       if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
    176           // average on low end of boundary, move star down
    177           stars[N].average.R -= 360.0;
    178           catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    179       }
    180       if (fabs(catalog[0].measure[Nmeas].dR) > 10*RADIUS) {
     168      float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
     169
     170      // rationalize R:
     171      if (dRoff > +180.0*3600.0) {
     172        // average on high end of boundary, move star up
     173        catalog[0].measure[Nmeas].R += 360.0;
     174        dRoff -= 360.0*3600.0;
     175      }
     176      if (dRoff < -180.0*3600.0) {
     177        // average on low end of boundary, move star down
     178        catalog[0].measure[Nmeas].R -= 360.0;
     179        dRoff += 360.0*3600.0;
     180      }
     181      if (fabs(dRoff) > 10*RADIUS) {
     182        // take declination into account and check again.
     183        double cosD = cos(RAD_DEG*catalog[0].average[n].D);
     184        if (fabs(dRoff*cosD) > 10*RADIUS) {
    181185          fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
    182186                   catalog[0].average[n].R, catalog[0].average[n].D,
    … …  
    184188                   X1[i], X2[J],
    185189                   Y1[i], Y2[J]);
     190        }
    186191      }
    187192
    … …  
    196201      if (Nsec > -1) {
    197202        if (isnan(catalog[0].secfilt[n*Nsecfilt+Nsec].M)) {
    198           catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
     203          catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
    199204        }
    200205      }
    … …  
    213218      }
    214219      /* this catalog star matches more than one image star */
    215       if (catalog[0].found[n] > -1) {
    216         catalog[0].measure[catalog[0].found[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
     220      if (catalog[0].found_t[n] > -1) {
     221        catalog[0].measure[catalog[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
    217222        catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
    218223      } else {
    219         catalog[0].found[n] = Nmeas;
     224        catalog[0].found_t[n] = Nmeas;
    220225      }
    221226      /* Nm is updated, but not written out in -update mode (for existing entries)
    … …  
    232237  }
    233238
    234   /* incorporate unmatched image stars, if this star is in field of this catalog */
    235   /* these new entries are all written out in UPDATE mode */
     239/* incorporate unmatched image stars, if this star is in field of this catalog */
     240/* these new entries are all written out in UPDATE mode */
    236241  for (i = 0; (i < Nstars) && !options.only_match; i += NSTAR_GROUP) {
    237242    /* make sure there is space for next entry */
    … …  
    260265
    261266    if (PSPS_ID) {
    262         catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
     267      catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
    263268    }
    264269
    … …  
    274279
    275280      // the following measure elements cannot be set until here:
    276       catalog[0].measure[Nmeas].dR              = 0.0;
    277       catalog[0].measure[Nmeas].dD              = 0.0;
     281      catalog[0].measure[Nmeas].R               = stars[i].average.R;
     282      catalog[0].measure[Nmeas].D               = stars[i].average.D;
    278283      catalog[0].measure[Nmeas].dbFlags         = 0;
    279284      catalog[0].measure[Nmeas].averef          = Nave; // XXX EAM : must be absolute Nave if partial read
    … …  
    284289      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    285290      if (Nsec > -1) {
    286           catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
     291        catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
    287292      }
    288293
    289294      /* next[Nmeas] should always be -1 in this context (it is always the only
    290         measurement for the star) */
     295        measurement for the star) */
    291296      stars[i].found = Nmeas;
    292297      next_meas[Nmeas] = -1;  // initial value here update below
    … …  
    309314  }
    310315
    311   /* note stars which have been found in this catalog */
     316/* note stars which have been found in this catalog */
    312317  for (i = 0; i < NstarsIn; i++) {
    313318    if (stars[i].found > -1) {
    … …  
    316321  }
    317322
    318   /* check if the catalog has changed?  if no change, no need to write */
     323/* check if the catalog has changed?  if no change, no need to write */
    319324  catalog[0].objID    = objID; // new max value, save on catalog close
    320325  catalog[0].Naverage = Nave;
    … …  
    323328  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    324329
    325   free (catalog[0].found);
    326330  free (X1);
    327331  free (Y1);

  • trunk/Ohana/src/addstar/src/find_matches_closest.c

    r36485 r37036  
    77  double *X1, *Y1, *X2, *Y2;
    88  double dX, dY, dR;
    9   off_t *N1, *N2, *next_meas;
    10   off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
     9  off_t *N1, *N2, *next_meas, *next_lens;
     10  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS;
    1111  int Nsecfilt, Nsec;
    1212  unsigned int objID, catID;
    … …  
    3636  ALLOCATE (Y2, double, NAVE);
    3737  ALLOCATE (N2, off_t,  NAVE);
    38   ALLOCATE (catalog[0].found, off_t, NAVE);
     38  ALLOCATE (catalog[0].found_t, off_t, NAVE);
    3939  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    4040
    … …  
    4242  Nmatch = 0;
    4343  NMEAS = Nmeas = catalog[0].Nmeasure;
     44  NLENS = Nlens = catalog[0].Nlensing;
    4445
    4546  // current max obj ID for this catalog
    … …  
    7677  if (Nstars < 1) {
    7778    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
    78     free (catalog[0].found);
    7979    free (X1);
    8080    free (Y1);
    … …  
    9191    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    9292    N2[i] = i;
    93     catalog[0].found[N2[i]] = -1;
     93    catalog[0].found_t[N2[i]] = -1;
    9494  }
    9595  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    … …  
    100100    // is sorted while processed
    101101    next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     102    next_lens = init_lensing_links (catalog[0].average, Nave, catalog[0].lensing, Nlens);
    102103  } else {
    103104    next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     105    next_lens = build_lensing_links (catalog[0].average, Nave, catalog[0].lensing, Nlens);
    104106  }   
    105107
    … …  
    177179      REALLOCATE (catalog[0].measure, Measure, NMEAS);
    178180    }
     181    if (Nlens >= NLENS) {
     182      NLENS = Nlens + 1000;
     183      REALLOCATE (next_lens, off_t, NLENS);
     184      REALLOCATE (catalog[0].lensing, Lensing, NLENS);
     185    }
    179186
    180187    /* add to end of measurement list */
    … …  
    186193    catalog[0].measure[Nmeas]          = stars[N].measure;
    187194
    188     /** dR,dD now represent arcsec **/
    189     catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    190     catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N].average.D);
     195    // measure now carries R,D (not dR,dD)
     196    // note that ReadStarsFITS does not set measure.R,D and average.R,D
     197    catalog[0].measure[Nmeas].R        = stars[N].average.R;
     198    catalog[0].measure[Nmeas].D        = stars[N].average.D;
    191199    catalog[0].measure[Nmeas].dbFlags  = 0;
    192200    catalog[0].measure[Nmeas].averef   = n;
    … …  
    194202    catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    195203
     204    float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
     205
    196206    // rationalize dR
    197     if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
     207    if (dRoff > +180.0*3600.0) {
    198208      // average on high end of boundary, move star up
    199       stars[N].average.R += 360.0;
    200       catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    201     }
    202     if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
     209      catalog[0].measure[Nmeas].R += 360.0;
     210      dRoff -= 360.0*3600.0;
     211    }
     212    if (dRoff < -180.0*3600.0) {
    203213      // average on low end of boundary, move star down
    204       stars[N].average.R -= 360.0;
    205       catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].average.R);
    206     }
    207     if (fabs(catalog[0].measure[Nmeas].dR) > 10*RADIUS) {
    208       fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
    209                catalog[0].average[n].R, catalog[0].average[n].D,
    210                stars[N].average.R, stars[N].average.D,
    211                X1[i], X2[Jmin],
    212                Y1[i], Y2[Jmin]);
    213         // XXX abort on this? -- this is a bad failure...
     214      catalog[0].measure[Nmeas].R -= 360.0;
     215      dRoff += 360.0*3600.0;
     216    }
     217    if (fabs(dRoff) > 10*RADIUS) {
     218        // take declination into account and check again.
     219        double cosD = cos(RAD_DEG*catalog[0].average[n].D);
     220        if (fabs(dRoff*cosD) > 10*RADIUS) {
     221            fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
     222                     catalog[0].average[n].R, catalog[0].average[n].D,
     223                     stars[N].average.R, stars[N].average.D,
     224                     X1[i], X2[Jmin],
     225                     Y1[i], Y2[Jmin]);
     226            // XXX abort on this? -- this is a bad failure...
     227        }
     228    }
     229
     230    // add the lensing values if they exist
     231    if (stars[N].lensing) {
     232      add_lens_link (&catalog[0].average[n], next_lens, Nlens, NLENS); // ?
     233      catalog[0].lensing[Nlens] = stars[N].lensing[0];
     234     
     235      catalog[0].lensing[Nlens].averef = n;
     236      catalog[0].lensing[Nlens].objID = catalog[0].average[n].objID;
     237      catalog[0].lensing[Nlens].catID = catalog[0].catID;
     238      catalog[0].average[n].Nlensing ++;
     239      Nlens ++;
    214240    }
    215241
    … …  
    224250    if (Nsec > -1) {
    225251      if (isnan(catalog[0].secfilt[n*Nsecfilt+Nsec].M)) {
    226         catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
     252        catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
    227253      }
    228254    }
    … …  
    234260       Nm is recalculated in build_meas_links if loaded table is not sorted */
    235261    stars[N].found = Nmeas;
    236     catalog[0].found[n] = Nmeas;
     262    catalog[0].found_t[n] = Nmeas;
    237263    catalog[0].average[n].Nmeasure ++;
    238264    Nmeas ++;
    … …  
    249275      REALLOCATE (catalog[0].measure, Measure, NMEAS);
    250276    }
     277    if (Nlens >= NLENS - NSTAR_GROUP) {
     278      NLENS = Nlens + 1000;
     279      REALLOCATE (next_lens, off_t, NLENS);
     280      REALLOCATE (catalog[0].lensing, Lensing, NLENS);
     281    }
    251282    if (Nave >= NAVE) {
    252283      NAVE = Nave + 1000;
    … …  
    278309
    279310    for (j = 0; j < NSTAR_GROUP; j++) {
    280         // supply the measurments from this detection
    281         catalog[0].measure[Nmeas]           = stars[i + j].measure;
    282 
    283         // the following measure elements cannot be set until here:
    284         catalog[0].measure[Nmeas].dR       = 0.0; // astrometric offset, not error
    285         catalog[0].measure[Nmeas].dD       = 0.0; // astrometric offset, not error
    286         catalog[0].measure[Nmeas].dbFlags  = 0;
    287         catalog[0].measure[Nmeas].averef   = Nave;
    288         catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
    289         catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    290 
    291         /* set the average magnitude if not already set and the photcode.equiv is not 0 */
    292         /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    293         if (Nsec > -1) {
    294             catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
    295         }
    296 
    297         /* next[Nmeas] should always be -1 in this context (it is always the only
    298            measurement for the star) */
    299         stars[i+j].found = Nmeas;
    300         next_meas[Nmeas] = -1;  // inital value here update below
    301         Nmeas ++;
     311      // supply the measurments from this detection
     312      dvo_measure_init (&catalog[0].measure[Nmeas]);
     313      catalog[0].measure[Nmeas]           = stars[i + j].measure;
     314
     315      // the following measure elements cannot be set until here:
     316      catalog[0].measure[Nmeas].R        = stars[i].average.R;
     317      catalog[0].measure[Nmeas].D        = stars[i].average.D;
     318      catalog[0].measure[Nmeas].dbFlags  = 0;
     319      catalog[0].measure[Nmeas].averef   = Nave;
     320      catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
     321      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
     322
     323      /* set the average magnitude if not already set and the photcode.equiv is not 0 */
     324      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
     325      if (Nsec > -1) {
     326        catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
     327      }
     328
     329      /* next[Nmeas] should always be -1 in this context (it is always the only
     330         measurement for the star) */
     331      stars[i+j].found = Nmeas;
     332      next_meas[Nmeas] = -1;  // inital value here update below
     333      Nmeas ++;
    302334    }
    303335    for (j = 0; j < NSTAR_GROUP - 1; j++) {
    304         next_meas[Nmeas - NSTAR_GROUP + j] = Nmeas - NSTAR_GROUP + j + 1;
    305     }
     336      next_meas[Nmeas - NSTAR_GROUP + j] = Nmeas - NSTAR_GROUP + j + 1;
     337    }
     338
     339    // if we have lensing data, insert that as well
     340    if (stars[i].lensing) {
     341      catalog[0].average[Nave].Nlensing = NSTAR_GROUP;
     342      catalog[0].average[Nave].lensingOffset = Nlens;
     343      for (j = 0; j < NSTAR_GROUP; j++) {
     344        // add the lensing values if they exist
     345        if (stars[i + j].lensing) {
     346          dvo_lensing_init (&catalog[0].lensing[Nlens]);
     347          catalog[0].lensing[Nlens] = stars[i + j].lensing[0];
     348         
     349          catalog[0].lensing[Nlens].averef = Nave;
     350          catalog[0].lensing[Nlens].objID = catalog[0].average[Nave].objID;
     351          catalog[0].lensing[Nlens].catID = catalog[0].catID;
     352          next_lens[Nlens] = -1;
     353          Nlens ++;
     354        }
     355      }
     356      for (j = 0; j < NSTAR_GROUP - 1; j++) {
     357        next_lens[Nlens - NSTAR_GROUP + j] = Nlens - NSTAR_GROUP + j + 1;
     358      }
     359    }
     360
    306361    Nave ++;
    307362  }
    … …  
    309364  REALLOCATE (catalog[0].average, Average, Nave);
    310365  REALLOCATE (catalog[0].measure, Measure, Nmeas);
     366  REALLOCATE (catalog[0].lensing, Lensing, Nlens);
    311367 
    312368  if (options.nosort) {
    … …  
    315371    catalog[0].sorted = TRUE;
    316372    catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
     373    catalog[0].lensing = sort_lensing (catalog[0].average, Nave, catalog[0].lensing, Nlens, next_lens);
    317374  }
    318375
    … …  
    328385  catalog[0].Naverage = Nave;
    329386  catalog[0].Nmeasure = Nmeas;
     387  catalog[0].Nlensing = Nlens;
    330388  catalog[0].Nsecf_mem = Nave*Nsecfilt;
    331   if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    332 
    333   free (catalog[0].found);
     389  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Nmatch);
     390
    334391  free (X1);
    335392  free (Y1);
    … …  
    339396  free (Y2);
    340397  free (next_meas);
     398  free (next_lens);
    341399
    342400  return (Nmatch);

  • trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c

    r36485 r37036  
    3131  ALLOCATE (Y2, double, NAVE);
    3232  ALLOCATE (N2, off_t,  NAVE);
    33   ALLOCATE (catalog[0].found, off_t, NAVE);
     33  ALLOCATE (catalog[0].found_t, off_t, NAVE);
    3434  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    3535
    … …  
    7171  if (Nstars < 1) {
    7272    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
    73     free (catalog[0].found);
    7473    free (X1);
    7574    free (Y1);
    … …  
    8685    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    8786    N2[i] = i;
    88     catalog[0].found[N2[i]] = -1;
     87    catalog[0].found_t[N2[i]] = -1;
    8988  }
    9089  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    … …  
    180179
    181180    /** *** dR,dD now in arcsec *** **/
    182     catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
    183     catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].average.D);
     181    catalog[0].measure[Nmeas].R        = stars[N][0].average.R;
     182    catalog[0].measure[Nmeas].D        = stars[N][0].average.D;
    184183    catalog[0].measure[Nmeas].dbFlags  = 0;
    185184    catalog[0].measure[Nmeas].averef   = n;
    … …  
    187186    catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    188187
     188    float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
     189
    189190    // rationalize dR:
    190     if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
     191    if (dRoff > +180.0*3600.0) {
    191192      // average on high end of boundary, move star up
    192       stars[N][0].average.R += 360.0;
    193       catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
    194     }
    195     if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
     193      catalog[0].measure[Nmeas].R += 360.0;
     194      dRoff -= 360.0*3600.0;
     195    }
     196    if (dRoff < -180.0*3600.0) {
    196197      // average on low end of boundary, move star down
    197       stars[N][0].average.R -= 360.0;
    198       catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
     198      catalog[0].measure[Nmeas].R -= 360.0;
     199      dRoff += 360.0*3600.0;
    199200    }
    200201
    … …  
    218219    /* this catalog star matches more than one image star */
    219220    // XXX should this be an average flag?
    220     if (catalog[0].found[n] > -1) {
    221       catalog[0].measure[catalog[0].found[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
     221    if (catalog[0].found_t[n] > -1) {
     222      catalog[0].measure[catalog[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
    222223      catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
    223224    } else {
    224       catalog[0].found[n] = Nmeas;
     225      catalog[0].found_t[n] = Nmeas;
    225226    }
    226227
    … …  
    290291      catalog[0].measure[Nmeas]          = stars[N][0].measure;
    291292
    292       catalog[0].measure[Nmeas].dR       = 0.0; // astrometric offset, not error
    293       catalog[0].measure[Nmeas].dD       = 0.0; // astrometric offset, not error
     293      catalog[0].measure[Nmeas].R       = catalog[0].average[Nave].R;
     294      catalog[0].measure[Nmeas].D       = catalog[0].average[Nave].D;
    294295      catalog[0].measure[Nmeas].dbFlags  = 0;
    295296      catalog[0].measure[Nmeas].averef   = Nave;
    … …  
    329330  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    330331
    331   free (catalog[0].found);
    332332  free (X1);
    333333  free (Y1);

  • trunk/Ohana/src/addstar/src/find_matches_refstars.c

    r36485 r37036  
    3737  ALLOCATE (Y2, double, NAVE);
    3838  ALLOCATE (N2, off_t,  NAVE);
    39   ALLOCATE (catalog[0].found, off_t, NAVE);
     39  ALLOCATE (catalog[0].found_t, off_t, NAVE);
    4040  REALLOCATE (catalog[0].average, Average, NAVE);
    4141  REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    … …  
    6969    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    7070    N2[i] = i;
    71     catalog[0].found[N2[i]] = -1;
     71    catalog[0].found_t[N2[i]] = -1;
    7272  }
    7373  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    … …  
    141141      catalog[0].measure[Nmeas]          = stars[N][0].measure;
    142142
    143       /** *** dR,dD now in arcsec *** **/
    144       catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
    145       catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].average.D);
     143      /** measure now stores R,D **/
     144      catalog[0].measure[Nmeas].R        = stars[N][0].average.R;
     145      catalog[0].measure[Nmeas].D        = stars[N][0].average.D;
    146146      catalog[0].measure[Nmeas].dbFlags  = 0;
    147147      catalog[0].measure[Nmeas].averef   = n;
    … …  
    149149      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    150150
     151      float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
     152
    151153      // rationalize dR:
    152       if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
     154      if (dRoff > +180.0*3600.0) {
    153155          // average on high end of boundary, move star up
    154           stars[N][0].average.R += 360.0;
    155           catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
    156       }
    157       if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
     156          catalog[0].measure[Nmeas].R += 360.0;
     157          dRoff -= 360.0*3600.0;
     158      }
     159      if (dRoff < -180.0*3600.0) {
    158160          // average on low end of boundary, move star down
    159           stars[N][0].average.R -= 360.0;
    160           catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
     161          catalog[0].measure[Nmeas].R -= 360.0;
     162          dRoff += 360.0*3600.0;
    161163      }
    162164
    … …  
    188190      }
    189191      /* this catalog star matches more than one image star */
    190       if (catalog[0].found[n] > -1) {
    191         catalog[0].measure[catalog[0].found[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
     192      if (catalog[0].found_t[n] > -1) {
     193        catalog[0].measure[catalog[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
    192194        catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
    193195      } else {
    194         catalog[0].found[n] = Nmeas;
     196        catalog[0].found_t[n] = Nmeas;
    195197      }
    196198
    … …  
    262264      catalog[0].measure[Nmeas]          = stars[N][0].measure;
    263265
    264       catalog[0].measure[Nmeas].dR       = 0.0;
    265       catalog[0].measure[Nmeas].dD       = 0.0;
     266      catalog[0].measure[Nmeas].R        = stars[N][0].average.R;
     267      catalog[0].measure[Nmeas].D        = stars[N][0].average.D;
    266268
    267269      catalog[0].measure[Nmeas].t        = (stars[N][0].measure.t == 0) ? TIMEREF : stars[N][0].measure.t;      /** careful : time_t vs e_time **/
    … …  
    303305  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    304306
    305   free (catalog[0].found);
    306307  free (X1);
    307308  free (Y1);

  • trunk/Ohana/src/addstar/src/findskycell.c

    r36497 r37036  
    161161  exit (0);
    162162}
    163 
    164 # define MARKTIME(MSG,...) {                    \
    165     float dtime;                                \
    166     gettimeofday (&stop, (void *) NULL);        \
    167     dtime = DTIME (stop, start);                \
    168     fprintf (stderr, MSG, __VA_ARGS__); }
    169163
    170164int mktree (char *treefile, char *catdir) {
    … …  
    402396  }
    403397
    404   struct timeval start, stop;
    405   gettimeofday (&start, (void *) NULL);
     398  INITTIME;
    406399
    407400  int Npts = 10000000;

  • trunk/Ohana/src/addstar/src/load2mass_catalog.c

    r35416 r37036  
    1717  for (i = 0; i < Nstars; i+=3) {
    1818
     19    double R = stars[i].average.R;
     20    double D = stars[i].average.D;
     21
    1922    // construct an average object for this object
    2023    // XXX for now, the output objects will have limited astrometric interpretation...
    2124    // XXX every 3 stars represents 3 measurements and 1 average
    2225    dvo_average_init (&catalog[0].average[Nave]);
    23     catalog[0].average[Nave].R     = stars[i].average.R;
    24     catalog[0].average[Nave].D     = stars[i].average.D;
     26    catalog[0].average[Nave].R     = R;
     27    catalog[0].average[Nave].D     = D;
    2528    catalog[0].average[Nave].measureOffset = Nmeas;
    2629
    … …  
    3336      catalog[0].measure[Nmeas]           = stars[i+j].measure;
    3437
    35       catalog[0].measure[Nmeas].dR        = 0.0;
    36       catalog[0].measure[Nmeas].dD        = 0.0;
     38      catalog[0].measure[Nmeas].R         = R;
     39      catalog[0].measure[Nmeas].D         = D;
    3740      catalog[0].measure[Nmeas].dt        = NAN_S_SHORT;
    3841

  • trunk/Ohana/src/addstar/src/loadsupercos_ops.c

    r33653 r37036  
    11# include "addstar.h"
    22# include "supercos.h"
    3 
    4 // this is defined in libohana/src/string.c but not generally exposed
    5 char *_parse_nextword_csv (char *string);
    63
    74int loadsupercos_getFilterInfo (char *line, char *emulsion, char *filterID) {
    … …  
    1310  // emulsion is field 11, filterID is field 12
    1411  for (i = 1; i < 11; i++) {
    15     p1 = _parse_nextword_csv (p1);
     12    p1 = parse_nextword_csv (p1);
    1613    if (!p1) {
    1714      fprintf (stderr, "error parsing filter info for line %s\n", line);
    … …  
    2421  }
    2522
    26   char *p2 = _parse_nextword_csv (p1);
     23  char *p2 = parse_nextword_csv (p1);
    2724  if (!p2) {
    2825    fprintf (stderr, "error parsing filter info for line %s\n", line);
    … …  
    3431  }
    3532
    36   char *p3 = _parse_nextword_csv (p2);
     33  char *p3 = parse_nextword_csv (p2);
    3734  if (!p3) {
    3835    fprintf (stderr, "error parsing filter info for line %s\n", line);
    … …  
    6057  // lstObs is field 20
    6158  for (i = 1; i < 20; i++) {
    62     p1 = _parse_nextword_csv (p1);
     59    p1 = parse_nextword_csv (p1);
    6360    if (!p1) {
    6461      fprintf (stderr, "error parsing filter info for line %s\n", line);

  • trunk/Ohana/src/addstar/src/mkcmf.c

    r35760 r37036  
    11# include "mkcmf.h"
    22
     3# define ZERO_POINT 25.0
    34# define SKY 100.0
    45# define DSKY 2.0
    … …  
    1415void gauss_init (int Nbin);
    1516double rnd_gauss (double mean, double sigma);
     17void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
     18void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1619void writeStars_PS1_V4 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1720void writeStars_PS1_V3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1821void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1922void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
     23void writeStars_PS1_SV3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    2024void writeStars_PS1_DEV_1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    2125void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    … …  
    2327int ADDNOISE = TRUE;
    2428float BAD_PSFQF_FRAC = 0.0;
     29
     30static float exptime = 1.0;
     31static Coords coords;
    2532
    2633int main (int argc, char **argv) {
    … …  
    3845  Matrix matrix;
    3946  FTable ftable;
    40   Coords coords;
    4147
    4248  int APPEND = FALSE;
    … …  
    105111  }
    106112
    107   float exptime = 1.0;
    108113  if ((N = get_argument (argc, argv, "-exptime"))) {
    109114    remove_argument (N, &argc, argv);
    … …  
    288293    gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
    289294  }
    290   gfits_modify (&header, "ZERO_PT", "%lf", 1, 25.0);
     295  gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT);
    291296  gfits_modify (&header, "EXPTIME", "%lf", 1, exptime);
    292297  gfits_modify (&header, "AIRMASS", "%lf", 1, airmass);
    … …  
    327332  if (!strcmp(type, "PS1_V4")) {
    328333    writeStars_PS1_V4 (&ftable, X, Y, M, Flag, Nstars);
     334    found = TRUE;
     335  }
     336  if (!strcmp(type, "PS1_V5")) {
     337    writeStars_PS1_V5 (&ftable, X, Y, M, Flag, Nstars);
     338    found = TRUE;
     339  }
     340  if (!strcmp(type, "PS1_V5_Lensing")) {
     341    writeStars_PS1_V5_Lensing (&ftable, X, Y, M, Flag, Nstars);
     342    found = TRUE;
     343  }
     344  if (!strcmp(type, "PS1_SV3")) {
     345    writeStars_PS1_SV3 (&ftable, X, Y, M, Flag, Nstars);
    329346    found = TRUE;
    330347  }
    … …  
    764781}
    765782
     783void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
     784
     785  int i;
     786  CMF_PS1_V5 *stars;
     787  float flux, fSN;
     788
     789  // XXX add gaussian-distributed noise based on counts
     790  // this needs to make different output 'stars' entries depending on the desired type
     791  ALLOCATE (stars, CMF_PS1_V5, Nstars);
     792  for (i = 0; i < Nstars; i++) {
     793
     794    flux = pow (10.0, -0.4*M[i]);
     795    fSN = 1.0 / sqrt(flux);
     796
     797    stars[i].detID = i;
     798    stars[i].X = X[i];
     799    stars[i].Y = Y[i];
     800    stars[i].dX = FX * fSN;
     801    stars[i].dY = FY * fSN;
     802
     803    stars[i].posangle = 10.0;
     804    stars[i].pltscale = 0.25;
     805
     806    stars[i].M = M[i];
     807    stars[i].dM = fSN;
     808
     809    stars[i].Flux = flux;
     810    stars[i].dFlux = flux * fSN;
     811
     812    stars[i].Map = M[i] - 0.05;
     813    stars[i].MapRaw = M[i] - 0.10;
     814
     815    stars[i].apRadius = 8.0;
     816
     817    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
     818    stars[i].apFluxErr = stars[i].apFlux * fSN;
     819
     820    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
     821    stars[i].dMcal = 0.05;
     822
     823    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
     824    stars[i].apNpix = 3.14*8.0*8.0;
     825
     826    stars[i].Mpeak     = M[i] + 1.0;
     827    stars[i].sky       = SKY;
     828    stars[i].dSky      = DSKY;
     829    stars[i].psfChisq  = PSFCHI;
     830    stars[i].crNsigma  = CRN;
     831    stars[i].extNsigma = EXTN;
     832    stars[i].fx        = FX;
     833    stars[i].fy        = FY;
     834    stars[i].df        = DF;
     835
     836    stars[i].k         = 1.0;
     837    stars[i].fwhmMaj   = FX*2.8;
     838    stars[i].fwhmMin   = FY*2.8;
     839
     840    // randomly give poor PSFQF values
     841    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
     842      stars[i].psfQF     = 0.25;
     843      stars[i].psfQFperf = 0.24;
     844    } else {
     845      stars[i].psfQF     = PSFQUAL;
     846      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
     847    }
     848
     849    stars[i].psfNdof   = 1;
     850    stars[i].psfNpix   = 2;
     851
     852    stars[i].Mxx       = FX;
     853    stars[i].Mxy       = 0.01;
     854    stars[i].Myy       = FX;
     855    stars[i].M3c       = FX;
     856    stars[i].M3s       = FX;
     857    stars[i].M4c       = FX;
     858    stars[i].M4s       = FX;
     859    stars[i].Mr1       = FX;
     860    stars[i].Mrh       = FX;
     861
     862    stars[i].kronFlux  = flux * 1.25;
     863    stars[i].kronFluxErr = fSN * flux * 1.25;
     864
     865    stars[i].kronInner = fSN * flux * 0.9;
     866    stars[i].kronOuter = fSN * flux * 1.5;
     867
     868    stars[i].skyLimitRad = 1;
     869    stars[i].skyLimitFlux = 2;
     870    stars[i].skyLimitSlope = 0.2;
     871
     872    stars[i].flags     = Flag[i];
     873    stars[i].flags2    = 0x80;
     874
     875    stars[i].nFrames   = 1;
     876
     877    if (ADDNOISE) {
     878      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     879      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     880      float Moff = fSN*rnd_gauss(0.0, 1.0);
     881      stars[i].M += Moff;
     882      stars[i].Map += Moff;
     883      stars[i].MapRaw += Moff;
     884      stars[i].Mcalib += Moff;
     885    }
     886  }
     887
     888  gfits_table_set_CMF_PS1_V5 (ftable, stars, Nstars);
     889  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V5");
     890}
     891
     892void writeStars_PS1_SV3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
     893
     894  int i;
     895  CMF_PS1_SV3 *stars;
     896  float flux, fSN;
     897
     898  // XXX add gaussian-distributed noise based on counts
     899  // this needs to make different output 'stars' entries depending on the desired type
     900  ALLOCATE (stars, CMF_PS1_SV3, Nstars);
     901  for (i = 0; i < Nstars; i++) {
     902
     903    flux = pow (10.0, -0.4*M[i]);
     904    fSN = 1.0 / sqrt(flux);
     905
     906    stars[i].detID = i;
     907    stars[i].X = X[i];
     908    stars[i].Y = Y[i];
     909    stars[i].dX = FX * fSN;
     910    stars[i].dY = FY * fSN;
     911
     912    stars[i].posangle = 10.0;
     913    stars[i].pltscale = 0.25;
     914
     915    stars[i].M = M[i];
     916    stars[i].dM = fSN;
     917
     918    stars[i].Flux = flux;
     919    stars[i].dFlux = flux * fSN;
     920
     921    stars[i].Map = M[i] - 0.05;
     922    stars[i].MapRaw = M[i] - 0.10;
     923
     924    stars[i].apRadius = 8.0;
     925
     926    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
     927    stars[i].apFluxErr = stars[i].apFlux * fSN;
     928
     929    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
     930    stars[i].dMcal = 0.05;
     931
     932    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
     933    stars[i].apNpix = 3.14*8.0*8.0;
     934
     935    stars[i].Mpeak     = M[i] + 1.0;
     936    stars[i].sky       = SKY;
     937    stars[i].dSky      = DSKY;
     938    stars[i].psfChisq  = PSFCHI;
     939    stars[i].crNsigma  = CRN;
     940    stars[i].extNsigma = EXTN;
     941    stars[i].fx        = FX;
     942    stars[i].fy        = FY;
     943    stars[i].df        = DF;
     944
     945    stars[i].k         = 1.0;
     946    stars[i].fwhmMaj   = FX*2.8;
     947    stars[i].fwhmMin   = FY*2.8;
     948
     949    // randomly give poor PSFQF values
     950    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
     951      stars[i].psfQF     = 0.25;
     952      stars[i].psfQFperf = 0.24;
     953    } else {
     954      stars[i].psfQF     = PSFQUAL;
     955      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
     956    }
     957
     958    stars[i].psfNdof   = 1;
     959    stars[i].psfNpix   = 2;
     960
     961    stars[i].Mxx       = FX;
     962    stars[i].Mxy       = 0.01;
     963    stars[i].Myy       = FX;
     964    stars[i].M3c       = FX;
     965    stars[i].M3s       = FX;
     966    stars[i].M4c       = FX;
     967    stars[i].M4s       = FX;
     968    stars[i].Mr1       = FX;
     969    stars[i].Mrh       = FX;
     970
     971    stars[i].kronFlux  = flux * 1.25;
     972    stars[i].kronFluxErr = fSN * flux * 1.25;
     973
     974    stars[i].kronInner = fSN * flux * 0.9;
     975    stars[i].kronOuter = fSN * flux * 1.5;
     976
     977    // stars[i].skyLimitRad = 1;
     978    // stars[i].skyLimitFlux = 2;
     979    // stars[i].skyLimitSlope = 0.2;
     980
     981    stars[i].flags     = Flag[i];
     982    stars[i].flags2    = 0x80;
     983
     984    stars[i].nFrames   = 1;
     985
     986    if (ADDNOISE) {
     987      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     988      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     989      float Moff = fSN*rnd_gauss(0.0, 1.0);
     990      stars[i].M += Moff;
     991      stars[i].Map += Moff;
     992      stars[i].MapRaw += Moff;
     993      stars[i].Mcalib += Moff;
     994    }
     995  }
     996
     997  gfits_table_set_CMF_PS1_SV3 (ftable, stars, Nstars);
     998  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_SV3");
     999}
     1000
     1001void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
     1002
     1003  int i;
     1004  CMF_PS1_V5_Lensing *stars;
     1005  float flux, fSN;
     1006
     1007  // XXX add gaussian-distributed noise based on counts
     1008  // this needs to make different output 'stars' entries depending on the desired type
     1009  ALLOCATE (stars, CMF_PS1_V5_Lensing, Nstars);
     1010  for (i = 0; i < Nstars; i++) {
     1011
     1012    flux = pow (10.0, -0.4*M[i]);
     1013    fSN = 1.0 / sqrt(flux);
     1014
     1015    stars[i].detID = i;
     1016    stars[i].X = X[i];
     1017    stars[i].Y = Y[i];
     1018    stars[i].dX = FX * fSN;
     1019    stars[i].dY = FY * fSN;
     1020
     1021    stars[i].posangle = 10.0;
     1022    stars[i].pltscale = 0.25;
     1023
     1024    stars[i].M = M[i];
     1025    stars[i].dM = fSN;
     1026
     1027    stars[i].Flux = flux;
     1028    stars[i].dFlux = flux * fSN;
     1029
     1030    stars[i].Map = M[i] - 0.05;
     1031    stars[i].MapRaw = M[i] - 0.10;
     1032
     1033    stars[i].apRadius = 8.0;
     1034
     1035    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
     1036    stars[i].apFluxErr = stars[i].apFlux * fSN;
     1037
     1038    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
     1039    stars[i].dMcal = 0.05;
     1040
     1041    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
     1042    stars[i].apNpix = 3.14*8.0*8.0;
     1043
     1044    stars[i].Mpeak     = M[i] + 1.0;
     1045    stars[i].sky       = SKY;
     1046    stars[i].dSky      = DSKY;
     1047    stars[i].psfChisq  = PSFCHI;
     1048    stars[i].crNsigma  = CRN;
     1049    stars[i].extNsigma = EXTN;
     1050    stars[i].fx        = FX;
     1051    stars[i].fy        = FY;
     1052    stars[i].df        = DF;
     1053
     1054    stars[i].k         = 1.0;
     1055    stars[i].fwhmMaj   = FX*2.8;
     1056    stars[i].fwhmMin   = FY*2.8;
     1057
     1058    // randomly give poor PSFQF values
     1059    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
     1060      stars[i].psfQF     = 0.25;
     1061      stars[i].psfQFperf = 0.24;
     1062    } else {
     1063      stars[i].psfQF     = PSFQUAL;
     1064      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
     1065    }
     1066
     1067    stars[i].psfNdof   = 1;
     1068    stars[i].psfNpix   = 2;
     1069
     1070    stars[i].Mxx       = FX;
     1071    stars[i].Mxy       = 0.01;
     1072    stars[i].Myy       = FX;
     1073    stars[i].M3c       = FX;
     1074    stars[i].M3s       = FX;
     1075    stars[i].M4c       = FX;
     1076    stars[i].M4s       = FX;
     1077    stars[i].Mr1       = FX;
     1078    stars[i].Mrh       = FX;
     1079
     1080    stars[i].X11_sm_obj = FX;
     1081    stars[i].X12_sm_obj = FX;
     1082    stars[i].X22_sm_obj = FY;
     1083    stars[i].E1_sm_obj  = FX;
     1084    stars[i].E2_sm_obj  = FX;
     1085
     1086    stars[i].X11_sh_obj = FX;
     1087    stars[i].X12_sh_obj = 0.2;
     1088    stars[i].X22_sh_obj = FY;
     1089    stars[i].E1_sh_obj  = FX*0.9;
     1090    stars[i].E2_sh_obj  = FX*0.8;
     1091
     1092    stars[i].X11_sm_psf = FX;
     1093    stars[i].X12_sm_psf = FX;
     1094    stars[i].X22_sm_psf = FX;
     1095    stars[i].E1_sm_psf  = FX;
     1096    stars[i].E2_sm_psf  = FX;
     1097
     1098    stars[i].X11_sh_psf = FX;
     1099    stars[i].X12_sh_psf = FX;
     1100    stars[i].X22_sh_psf = FX;
     1101    stars[i].E1_sh_psf  = FX;
     1102    stars[i].E2_sh_psf  = FX;
     1103
     1104    stars[i].kronFlux  = flux * 1.25;
     1105    stars[i].kronFluxErr = fSN * flux * 1.25;
     1106
     1107    stars[i].kronInner = fSN * flux * 0.9;
     1108    stars[i].kronOuter = fSN * flux * 1.5;
     1109
     1110    stars[i].skyLimitRad = 1;
     1111    stars[i].skyLimitFlux = 2;
     1112    stars[i].skyLimitSlope = 0.2;
     1113
     1114    stars[i].flags     = Flag[i];
     1115    stars[i].flags2    = 0x80;
     1116
     1117    stars[i].nFrames   = 1;
     1118
     1119    if (ADDNOISE) {
     1120      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     1121      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     1122      float Moff = fSN*rnd_gauss(0.0, 1.0);
     1123      stars[i].M += Moff;
     1124      stars[i].Map += Moff;
     1125      stars[i].MapRaw += Moff;
     1126      stars[i].Mcalib += Moff;
     1127    }
     1128  }
     1129
     1130  gfits_table_set_CMF_PS1_V5_Lensing (ftable, stars, Nstars);
     1131  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V5");
     1132}
     1133

  • trunk/Ohana/src/addstar/src/replace_match.c

    r26278 r37036  
    1212    j = i + m;
    1313    if (measure[j].photcode != star[0].measure.photcode) continue;
    14     measure[j].dR = 3600.0*(average[0].R - star[0].average.R);
    15     measure[j].dD = 3600.0*(average[0].D - star[0].average.D);
     14    measure[j].R = star[0].average.R;
     15    measure[j].D = star[0].average.D;
    1616    measure[j].M  = star[0].measure.M;
    1717    measure[j].dM = star[0].measure.dM;

  • trunk/Ohana/src/addstar/src/resort_catalog.c

    r35760 r37036  
    11# include "addstar.h"
     2
     3void SortAveMatch (off_t *MEAS, off_t *AVE, off_t N);
     4void resort_catalog_measure (Catalog *catalog);
     5void resort_catalog_lensing (Catalog *catalog);
    26
    37void resort_catalog_old (Catalog *catalog) {
    … …  
    59  off_t *next_meas;
    610  off_t Naves, Nmeas;
    7   double dtime;
    8   struct timeval start, stop;
    911
    1012  if (catalog[0].sorted == TRUE) return;
    1113
    12   gettimeofday (&start, NULL);
     14  INITTIME;
    1315
    1416  /* internal counters */
    … …  
    2224  catalog[0].measure = sort_measure (catalog[0].average, Naves, catalog[0].measure, Nmeas, next_meas);
    2325
    24   gettimeofday (&stop, NULL);
    25   dtime = DTIME (stop, start);
    26   fprintf (stderr, "  match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves);
     26  MARKTIME ("  match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves);
    2727
    2828  return;
    2929}
    3030
    31 // sort the measure Sequence based on the average Sequence entries
    32 void SortAveMeasMatch (off_t *MEAS, off_t *AVE, off_t N) {
    33 
    34 # define SWAPFUNC(A,B){ off_t tmp_meas; off_t tmp_ave;          \
    35     tmp_meas = MEAS[A]; MEAS[A] = MEAS[B]; MEAS[B] = tmp_meas;  \
    36     tmp_ave  = AVE[A];  AVE[A]  = AVE[B];  AVE[B]  = tmp_ave;   \
    37   }
    38 # define COMPARE(A,B)(AVE[A] < AVE[B])
    39   OHANA_SORT (N, COMPARE, SWAPFUNC);
    40 # undef SWAPFUNC
    41 # undef COMPARE
    42 }
    43 
    44 # define MARKTIME(MSG,...) {                    \
    45     float dtime;                                \
    46     gettimeofday (&stop, (void *) NULL);        \
    47     dtime = DTIME (stop, start); start = stop;  \
    48     fprintf (stderr, MSG, __VA_ARGS__); }
    49 
    50 // XXX : where is the time going?  perhaps the ALLOCATE?
    51 // XXX : I don't thnk his is getting the right answer yet.
    52 
    5331void resort_catalog (Catalog *catalog) {
     32
     33  if (catalog[0].sorted == TRUE) return;
     34
     35  resort_catalog_measure (catalog);
     36  resort_catalog_lensing (catalog);
     37}
     38
     39void resort_catalog_measure (Catalog *catalog) {
    5440
    5541  off_t Naverage, Nmeasure;
    … …  
    6248  Measure *measureTMP = NULL;
    6349
    64   if (catalog[0].sorted == TRUE) return;
    65 
    6650  // struct timeval start, stop;
    6751  // gettimeofday (&start, NULL);
    … …  
    7054  Nmeasure = catalog[0].Nmeasure;
    7155  Naverage = catalog[0].Naverage;
     56
     57  if (!Nmeasure) return;
    7258
    7359  measure = catalog[0].measure;
    … …  
    11399  // fprintf (stderr, "\n");
    114100
    115   SortAveMeasMatch(measureSeq, averageSeq, Nmeasure);
     101  SortAveMatch(measureSeq, averageSeq, Nmeasure);
    116102  // MARKTIME("sort : %f sec\n", dtime);
    117103
    … …  
    188174}
    189175
     176void resort_catalog_lensing (Catalog *catalog) {
     177
     178  off_t Naverage, Nlensing;
     179  Lensing *lensing;
     180  Average *average;
     181  off_t i, j, N, currentAve;
     182
     183  off_t *lensingSeq = NULL;
     184  off_t *averageSeq = NULL;
     185  Lensing *lensingTMP = NULL;
     186
     187  // struct timeval start, stop;
     188  // gettimeofday (&start, NULL);
     189
     190  /* internal counters */
     191  Nlensing = catalog[0].Nlensing;
     192  Naverage = catalog[0].Naverage;
     193
     194  if (!Nlensing) return;
     195
     196  lensing = catalog[0].lensing;
     197  average = catalog[0].average;
     198 
     199  // we have a table of average objects and an unsorted table of measurements.  each measurement
     200  // has a reference to the average object sequence (as well as an ID)
     201  // lensing[i].averef -> average[averef]
     202  // lensing[i].objID = average[averef].objID
     203  // lensing[i].catID = average[averef].catID
     204
     205  // we want a sorted lensing array with all averef entries in sequence
     206
     207  ALLOCATE (lensingSeq, off_t,   Nlensing);
     208  ALLOCATE (averageSeq, off_t,   Nlensing);
     209
     210  for (i = 0; i < Nlensing; i++) {
     211    lensingSeq[i] = i;
     212    averageSeq[i] = lensing[i].averef;
     213   
     214    if (catalog[0].catformat >= DVO_FORMAT_PS1_V1) {
     215      // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't)
     216      myAssert(average[averageSeq[i]].catID == lensing[lensingSeq[i]].catID, "object / detection mismatch");
     217# if (1)
     218      myAssert(average[averageSeq[i]].objID == lensing[lensingSeq[i]].objID, "object / detection mismatch");
     219# else
     220      // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID. 
     221      // this happened for 28 detections in the dbs on /data/stsci1?.0/eugene/dvo3pi.20130616, but not at all (as far as I know) in the rest of LAP DVO
     222      if (average[averageSeq[i]].objID != lensing[lensingSeq[i]].objID) {
     223        fprintf (stderr, "R");
     224        lensing[lensingSeq[i]].objID = average[averageSeq[i]].objID; // XXX I don't really like this...
     225      }
     226# endif
     227    }
     228  }
     229 
     230  // check that averageSeq is now in order
     231  // for (i = 1; i < Nlensing; i++) {
     232  //   if (averageSeq[i] < averageSeq[i-1]) {
     233  //     fprintf (stderr, "%d ", (int) i);
     234  //   }
     235  // }
     236  // fprintf (stderr, "\n");
     237
     238  SortAveMatch(lensingSeq, averageSeq, Nlensing);
     239  // MARKTIME("sort : %f sec\n", dtime);
     240
     241  // check that averageSeq is now in order
     242  // for (i = 1; i < Nlensing; i++) {
     243  //   if (averageSeq[i] < averageSeq[i-1]) {
     244  //     fprintf (stderr, "%d ", (int) i);
     245  //   }
     246  // }
     247  // fprintf (stderr, "\n");
     248
     249  // copy the lensing entries in the sorted order
     250  ALLOCATE (lensingTMP, Lensing, Nlensing);
     251  for (i = 0; i < Nlensing; i++) {
     252    j = lensingSeq[i];
     253    lensingTMP[i] = lensing[j];
     254  }
     255  // MARKTIME("assign lensing : %f sec\n", dtime);
     256
     257  // update the values of average.lensingOffset and average.Nlensing
     258  FREE(lensing);
     259  catalog[0].lensing = lensingTMP;
     260
     261  N = 0;
     262  currentAve = averageSeq[0];
     263  average[currentAve].lensingOffset = 0;
     264  for (i = 0; i < Nlensing; i++) {
     265    if (averageSeq[i] != currentAve) {
     266      // we have hit the next entry in the list
     267      average[currentAve].Nlensing = N;
     268      N = 0;
     269      currentAve = averageSeq[i];
     270      average[currentAve].lensingOffset = i;
     271    }
     272    N++;
     273  }
     274  // N++;
     275  average[currentAve].Nlensing = N;
     276  // MARKTIME("update Nlensing : %f sec\n", dtime);
     277
     278  int NlensingTotal = 0;
     279  int lensingOffsetOK = TRUE;
     280  for (i = 0; i < Naverage; i++) {
     281    NlensingTotal += catalog[0].average[i].Nlensing;
     282    if (VERBOSE && !(NlensingTotal <= catalog[0].Nlensing)) {
     283      fprintf (stderr, "too few lensing: %d %d %d\n", (int) i, NlensingTotal, (int) catalog[0].Nlensing);
     284    }
     285    lensingOffsetOK &= (catalog[0].average[i].lensingOffset < catalog[0].Nlensing);
     286    if (VERBOSE && !(catalog[0].average[i].lensingOffset < catalog[0].Nlensing)) {
     287      fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nlensing, (int) catalog[0].Nlensing);
     288    }
     289    lensingOffsetOK &= (catalog[0].average[i].lensingOffset + catalog[0].average[i].Nlensing <= catalog[0].Nlensing);
     290    if (VERBOSE && !(catalog[0].average[i].lensingOffset + catalog[0].average[i].Nlensing <= catalog[0].Nlensing)) {
     291      fprintf (stderr, "orrset + Nlensing too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].lensingOffset, catalog[0].average[i].Nlensing, (int) catalog[0].Nlensing);
     292    }
     293  }
     294
     295  if (!lensingOffsetOK) {
     296    fprintf (stderr, "ERROR: catalog %s has an invalid lensingOffset\n", catalog[0].filename);
     297  }
     298
     299  if (NlensingTotal != catalog[0].Nlensing) {
     300    fprintf (stderr, "ERROR: catalog %s has an invalid Nlensing\n", catalog[0].filename);
     301  }
     302
     303  // MARKTIME("  match time %9.4f sec for %7lld lensings, %6lld average\n", dtime, (long long) Nlensing, (long long) Naverage);
     304
     305  catalog[0].sorted = TRUE;
     306
     307  FREE (lensingSeq);
     308  FREE (averageSeq);
     309
     310  return;
     311}
     312
     313// sort the measure or lensing Sequence based on the average Sequence entries
     314void SortAveMatch (off_t *MEAS, off_t *AVE, off_t N) {
     315
     316# define SWAPFUNC(A,B){ off_t tmp_meas; off_t tmp_ave;          \
     317    tmp_meas = MEAS[A]; MEAS[A] = MEAS[B]; MEAS[B] = tmp_meas;  \
     318    tmp_ave  = AVE[A];  AVE[A]  = AVE[B];  AVE[B]  = tmp_ave;   \
     319  }
     320# define COMPARE(A,B)(AVE[A] < AVE[B])
     321  OHANA_SORT (N, COMPARE, SWAPFUNC);
     322# undef SWAPFUNC
     323# undef COMPARE
     324}
     325

  • trunk/Ohana/src/addstar/src/resort_threaded.c

    r29938 r37036  
    4141    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    4242    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    43     catalog.catflags  = LOAD_AVES | LOAD_MEAS;
     43    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_LENSING;
    4444    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    4545 

  • trunk/Ohana/src/addstar/src/resort_unthreaded_catalogs.c

    r33963 r37036  
    3030    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    3131    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    32     catalog.catflags  = LOAD_AVES | LOAD_MEAS;
     32    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_LENSING;
    3333    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    3434

  • trunk/Ohana/src/addstar/src/update_coords.c

    r34405 r37036  
    1919      continue;
    2020    }
    21     R = measure[m].dR;
    22     D = measure[m].dD;
     21    R = measure[m].R;
     22    D = measure[m].D;
    2323    r += R;
    2424    d += D;
    … …  
    3333  r = r / Npt;  /* these are corrections in 1/100 arcsec to RA and DEC */
    3434  d = d / Npt;
    35   average[0].R -= r / 3600.0;
    36   average[0].D -= d / 3600.0;
     35  average[0].R = r;
     36  average[0].D = d;
    3737  m = average[0].measureOffset;  /* first measurement of this star */
    38   for (i = 0; i < average[0].Nmeasure; i++) {
    39     measure[m].dR -= r;
    40     measure[m].dD -= d;
    41     m = next[m];
    42   }
    43  
     38
    4439  /* measure scatter, if possible */
    4540  if (Npt < 2) return;
    … …  
    4742  dR2 = r2 / Npt - r*r;
    4843  dD2 = d2 / Npt - d*d;
    49   average[0].ChiSqAve = sqrt (dD2 + dR2 / SQ(cos(d*RAD_DEG)));
     44  average[0].ChiSqAve = 3600.0*sqrt (dD2 + dR2 / SQ(cos(d*RAD_DEG)));
    5045  /* ChiSqAve is supposed to be a chisq */
    5146

  • trunk/Ohana/src/addstar/test/simple.dvo

    r35263 r37036  
    77  test.fields PS1_DEV_0 PS1_V1
    88  test.fields PS1_DEV_1 PS1_V1
    9   test.fields PS1_V1    PS1_V1
    10   test.fields PS1_V2    PS1_V1
    11   test.fields PS1_V3    PS1_V1
    12   test.fields PS1_V4    PS1_V1
     9  test.fields PS1_V1    PS1_V1
     10  test.fields PS1_V2    PS1_V1
     11  test.fields PS1_V3    PS1_V1
     12  test.fields PS1_V4    PS1_V1
    1313
    1414  test.fields PS1_DEV_0 PS1_V2
    1515  test.fields PS1_DEV_1 PS1_V2
    16   test.fields PS1_V1    PS1_V2
    17   test.fields PS1_V2    PS1_V2
    18   test.fields PS1_V3    PS1_V2
    19   test.fields PS1_V4    PS1_V4
     16  test.fields PS1_V1    PS1_V2
     17  test.fields PS1_V2    PS1_V2
     18  test.fields PS1_V3    PS1_V2
     19  test.fields PS1_V4    PS1_V4
    2020
    2121  test.fields PS1_DEV_0 PS1_V3
    2222  test.fields PS1_DEV_1 PS1_V3
    23   test.fields PS1_V1    PS1_V3
    24   test.fields PS1_V2    PS1_V3
    25   test.fields PS1_V3    PS1_V3
    26   test.fields PS1_V4    PS1_V4
     23  test.fields PS1_V1    PS1_V3
     24  test.fields PS1_V2    PS1_V3
     25  test.fields PS1_V3    PS1_V3
     26  test.fields PS1_V4    PS1_V4
    2727
    2828  test.fields PS1_DEV_0 PS1_V4
    2929  test.fields PS1_DEV_1 PS1_V4
    30   test.fields PS1_V1    PS1_V4
    31   test.fields PS1_V2    PS1_V4
    32   test.fields PS1_V3    PS1_V4
    33   test.fields PS1_V4    PS1_V4
     30  test.fields PS1_V1    PS1_V4
     31  test.fields PS1_V2    PS1_V4
     32  test.fields PS1_V3    PS1_V4
     33  test.fields PS1_V4    PS1_V4
    3434end 
     35
     36if (not($?NO_NOISE))
     37  $NO_NOISE = ""
     38end
    3539
    3640# create a populated catdir with a single cmf -- test each field
    … …  
    5054
    5155  mkinput
    52   exec mkcmf test.in.txt test.cmf -date 2008/1/1 -time 01:00:00 -radec $RA $DEC -type $1
     56  echo mkcmf test.in.txt test.cmf -date 2008/1/1 -time 01:00:00 -radec $RA $DEC -type $1 $NO_NOISE
     57  exec mkcmf test.in.txt test.cmf -date 2008/1/1 -time 01:00:00 -radec $RA $DEC -type $1 $NO_NOISE
    5358  if ($TAP_VERBOSE)
    5459    echo exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass
    5560    exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass
    5661  else
     62    echo addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass
    5763    exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass >& tmp.log
    5864  end
    … …  
    112118
    113119    # some fields require arithmetic manipulations
    114     if ("$name:0" == "KRON_FLUX")
    115      set v1 = -2.5*log(v1)
    116     end
    117     if ("$name:0" == "KRON_FLUX_ERR")
    118      set v1 = KRON_FLUX_ERR / KRON_FLUX
     120    ## if ("$name:0" == "KRON_FLUX")
     121    ##  set v1 = -2.5*log(v1)
     122    ## end
     123    ## if ("$name:0" == "KRON_FLUX_ERR")
     124    ##  set v1 = KRON_FLUX_ERR / KRON_FLUX
     125    ## end
     126    if (("$name:0" == "X_PSF_SIG") || ("$name:0" == "Y_PSF_SIG"))
     127     set v1 = int(100*(v1))
     128     set v2 = int(100*(v2 + 0.0001))
     129    end
     130    if ("$name:0" == "POSANGLE")
     131     set v1 = int((0xffff/360.0)*(v1))
     132     set v2 = int((0xffff/360.0)*(v2 + 0.0028))
    119133    end
    120134
    … …  
    206220# this list is good for PS1_V1
    207221list testfields_PS1_V1
    208   IPP_IDET          : detid     
    209   X_PSF             : xccd     
    210   Y_PSF             : yccd     
     222  IPP_IDET          : detid     
     223  X_PSF             : xccd     
     224  Y_PSF             : yccd     
    211225  X_PSF_SIG         : xccd:err #  FAIL
    212226  Y_PSF_SIG         : yccd:err # FAIL
    … …  
    215229  POSANGLE          : SKIP # astrometry is not calibrated in the cmf
    216230  PLTSCALE          : SKIP # astrometry is not calibrated in the cmf
    217   PSF_INST_MAG      : mag:inst 
    218   PSF_INST_MAG_SIG  : mag:err   
     231  PSF_INST_MAG      : mag:inst 
     232  PSF_INST_MAG_SIG  : mag:err   
    219233  AP_MAG_STANDARD   : mag:ap # FAIL
    220234  AP_MAG_RADIUS     : SKIP # no accessor
    … …  
    222236  CAL_PSF_MAG       : SKIP # photometry is not calibrated in the cmf
    223237  CAL_PSF_MAG_SIG   : SKIP # photometry is not calibrated in the cmf
    224   SKY               : sky       
    225   SKY_SIG           : sky_err   
    226   PSF_CHISQ         : psf_chisq 
    227   CR_NSIGMA         : cr_nsigma 
     238  SKY               : sky       
     239  SKY_SIG           : sky_err   
     240  PSF_CHISQ         : psf_chisq 
     241  CR_NSIGMA         : cr_nsigma 
    228242  EXT_NSIGMA        : ext_nsigma
    229   PSF_MAJOR         : FWHM_MAJ 
    230   PSF_MINOR         : FWHM_MIN 
     243  PSF_MAJOR         : FWHM_MAJ 
     244  PSF_MINOR         : FWHM_MIN 
    231245  PSF_THETA         : THETA # FAIL
    232   PSF_QF            : PSF_QF   
     246  PSF_QF            : PSF_QF   
    233247  PSF_NDOF          : SKIP # no accessor
    234248  PSF_NPIX          : SKIP # no accessor
    … …  
    237251  MOMENTS_YY        : SKIP # no accessor
    238252  FLAGS             : phot_flags
    239   N_FRAMES          : SKIP # no accessor       
     253  N_FRAMES          : SKIP # no accessor       
    240254end
    241255
    242256# this list is good for PS1_V2
    243257list testfields_PS1_V2
    244   IPP_IDET          : detid     
    245   X_PSF             : xccd     
    246   Y_PSF             : yccd     
     258  IPP_IDET          : detid     
     259  X_PSF             : xccd     
     260  Y_PSF             : yccd     
    247261  X_PSF_SIG         : xccd:err # FAIL
    248262  Y_PSF_SIG         : yccd:err # FAIL
    … …  
    251265  POSANGLE          : SKIP # astrometry is not calibrated in the cmf
    252266  PLTSCALE          : SKIP # astrometry is not calibrated in the cmf
    253   PSF_INST_MAG      : mag:inst 
    254   PSF_INST_MAG_SIG  : mag:err   
     267  PSF_INST_MAG      : mag:inst 
     268  PSF_INST_MAG_SIG  : mag:err   
    255269  AP_MAG_STANDARD   : mag:ap # FAIL
    256270  AP_MAG_RADIUS     : SKIP # no accessor
    … …  
    258272  CAL_PSF_MAG       : SKIP # photometry is not calibrated in the cmf
    259273  CAL_PSF_MAG_SIG   : SKIP # photometry is not calibrated in the cmf
    260   SKY               : sky       
    261   SKY_SIG           : sky_err   
    262   PSF_CHISQ         : psf_chisq 
    263   CR_NSIGMA         : cr_nsigma 
     274  SKY               : sky       
     275  SKY_SIG           : sky_err   
     276  PSF_CHISQ         : psf_chisq 
     277  CR_NSIGMA         : cr_nsigma 
    264278  EXT_NSIGMA        : ext_nsigma
    265   PSF_MAJOR         : FWHM_MAJ 
    266   PSF_MINOR         : FWHM_MIN 
     279  PSF_MAJOR         : FWHM_MAJ 
     280  PSF_MINOR         : FWHM_MIN 
    267281  PSF_THETA         : THETA # FAIL
    268   PSF_QF            : PSF_QF   
     282  PSF_QF            : PSF_QF   
    269283  PSF_NDOF          : SKIP # no accessor
    270284  PSF_NPIX          : SKIP # no accessor
    … …  
    273287  MOMENTS_YY        : SKIP # no accessor
    274288  FLAGS             : phot_flags
    275   N_FRAMES          : SKIP # no accessor       
     289  N_FRAMES          : SKIP # no accessor       
    276290end
    277291
    278292# this list is good for PS1_V3
    279293list testfields_PS1_V3
    280   IPP_IDET          : detid     
    281   X_PSF             : xccd     
    282   Y_PSF             : yccd     
     294  IPP_IDET          : detid     
     295  X_PSF             : xccd     
     296  Y_PSF             : yccd     
    283297  X_PSF_SIG         : xccd:err # FAIL
    284298  Y_PSF_SIG         : yccd:err # FAIL
    285299  POSANGLE          : SKIP # astrometry is not calibrated in the cmf
    286300  PLTSCALE          : SKIP # astrometry is not calibrated in the cmf
    287   PSF_INST_MAG      : mag:inst 
    288   PSF_INST_MAG_SIG  : mag:err   
     301  PSF_INST_MAG      : mag:inst 
     302  PSF_INST_MAG_SIG  : mag:err   
    289303  PSF_INST_FLUX     : SKIP # not ingested into DVO
    290304  PSF_INST_FLUX_SIG : SKIP # not ingested into DVO
    … …  
    297311  DEC_PSF           : SKIP # astrometry is not calibrated in the cmf
    298312  PEAK_FLUX_AS_MAG  : SKIP # not ingested into DVO
    299   SKY               : sky       
    300   SKY_SIG           : sky_err   
    301   PSF_CHISQ         : psf_chisq 
    302   CR_NSIGMA         : cr_nsigma 
     313  SKY               : sky       
     314  SKY_SIG           : sky_err   
     315  PSF_CHISQ         : psf_chisq 
     316  CR_NSIGMA         : cr_nsigma 
    303317  EXT_NSIGMA        : ext_nsigma
    304   PSF_MAJOR         : FWHM_MAJ 
    305   PSF_MINOR         : FWHM_MIN 
     318  PSF_MAJOR         : FWHM_MAJ 
     319  PSF_MINOR         : FWHM_MIN 
    306320  PSF_THETA         : THETA # FAIL
    307   PSF_QF            : PSF_QF   
     321  PSF_QF            : PSF_QF   
    308322  PSF_QF_PERFECT    : SKIP # not ingested into DVO
    309323  PSF_NDOF          : PSF_NDOF
    … …  
    326340end
    327341
    328 # this list is good for PS1_V4
    329 list testfields_PS1_V4
    330   IPP_IDET          : detid     
    331   X_PSF             : xccd     
    332   Y_PSF             : yccd     
    333   X_PSF_SIG         : xccd:err # FAIL
    334   Y_PSF_SIG         : yccd:err # FAIL
    335   POSANGLE          : SKIP # astrometry is not calibrated in the cmf
    336   PLTSCALE          : SKIP # astrometry is not calibrated in the cmf
    337   PSF_INST_MAG      : mag:inst 
    338   PSF_INST_MAG_SIG  : mag:err   
    339   PSF_INST_FLUX     : SKIP # not ingested into DVO
    340   PSF_INST_FLUX_SIG : SKIP # not ingested into DVO
    341   AP_MAG            : mag:aperinst # FAIL
     342# this list is good for PS1_V5
     343list testfields_PS1_V5
     344  IPP_IDET          : detid     
     345  X_PSF             : xccd     
     346  Y_PSF             : yccd     
     347  X_PSF_SIG         : xccd:err
     348  Y_PSF_SIG         : yccd:err
     349  POSANGLE          : posangle
     350  PLTSCALE          : platescale
     351  PSF_INST_MAG      : mag:inst 
     352  PSF_INST_MAG_SIG  : mag:err   
     353  PSF_INST_FLUX     : flux:inst
     354  PSF_INST_FLUX_SIG : flux:inst:err
     355  AP_MAG            : mag:aper:inst
    342356  AP_MAG_RAW        : SKIP # not ingested into DVO
    343357  AP_MAG_RADIUS     : SKIP # not ingested into DVO
    344   CAL_PSF_MAG       : SKIP # photometry is not calibrated in the cmf
    345   CAL_PSF_MAG_SIG   : SKIP # photometry is not calibrated in the cmf
    346   RA_PSF            : SKIP # astrometry is not calibrated in the cmf
    347   DEC_PSF           : SKIP # astrometry is not calibrated in the cmf
    348   PEAK_FLUX_AS_MAG  : SKIP # not ingested into DVO
    349   SKY               : sky       
    350   SKY_SIGMA         : sky_err   
    351   PSF_CHISQ         : psf_chisq
    352   CR_NSIGMA         : cr_nsigma
    353   EXT_NSIGMA        : ext_nsigma
    354   PSF_MAJOR         : FWHM_MAJ 
    355   PSF_MINOR         : FWHM_MIN 
    356   PSF_THETA         : THETA # FAIL
    357   PSF_QF            : PSF_QF   
    358   PSF_QF_PERFECT    : SKIP # not ingested into DVO
    359   PSF_NDOF          : PSF_NDOF
    360   PSF_NPIX          : PSF_NPIX
    361   MOMENTS_XX        : MXX
    362   MOMENTS_XY        : MXY
    363   MOMENTS_YY        : MYY
    364   MOMENTS_M3C       : SKIP # not ingested into DVO
    365   MOMENTS_M3S       : SKIP # not ingested into DVO
    366   MOMENTS_M4C       : SKIP # not ingested into DVO
    367   MOMENTS_M4S       : SKIP # not ingested into DVO
    368   MOMENTS_R1        : SKIP # not ingested into DVO
    369   MOMENTS_RH        : SKIP # not ingested into DVO
    370   KRON_FLUX         : mag:kroninst
    371   KRON_FLUX_ERR     : mag:kronerr
    372   KRON_FLUX_INNER   : SKIP # not ingested into DVO
    373   KRON_FLUX_OUTER   : SKIP # not ingested into DVO
    374   FLAGS             : phot_flags
    375   N_FRAMES          : SKIP # not ingested into DVO
    376 end
     358  AP_FLUX           : flux:aper:inst 
     359  AP_FLUX_SIG       : flux:aper:inst:err
     360  CAL_PSF_MAG       : mag:cal
     361  CAL_PSF_MAG_SIG   : mag:cal:err
     362  RA_PSF            : ra
     363  DEC_PSF           : dec
     364  AP_NPIX           : SKIP # not ingested into DVO                   
     365  PEAK_FLUX_AS_MAG  : SKIP # not ingested into DVO                   
     366  SKY               : sky                                             
     367  SKY_SIGMA         : sky_err                                         
     368  PSF_CHISQ         : psf_chisq                                       
     369  CR_NSIGMA         : cr_nsigma                                       
     370  EXT_NSIGMA        : ext_nsigma                                     
     371  PSF_MAJOR         : fwhm_maj # XXX technically, this is not FWHM but Sigma
     372  PSF_MINOR         : fwhm_min # XXX technically, this is not FWHM but Sigma
     373  PSF_THETA         : theta # FAIL                                   
     374  PSF_CORE          : SKIP # not ingested into DVO
     375  PSF_FWHM_MAJ      : SKIP # not ingested into DVO
     376  PSF_FWHM_MIN      : SKIP # not ingested into DVO
     377  PSF_QF            : psf_qf                                         
     378  PSF_QF_PERFECT    : SKIP # not ingested into DVO                   
     379  PSF_NDOF          : psf_ndof                               
     380  PSF_NPIX          : psf_npix                               
     381  MOMENTS_XX        : Mxx                                             
     382  MOMENTS_XY        : Mxy                                             
     383  MOMENTS_YY        : Myy                                             
     384  MOMENTS_M3C       : SKIP # not ingested into DVO                   
     385  MOMENTS_M3S       : SKIP # not ingested into DVO                   
     386  MOMENTS_M4C       : SKIP # not ingested into DVO                   
     387  MOMENTS_M4S       : SKIP # not ingested into DVO                   
     388  MOMENTS_R1        : SKIP # not ingested into DVO                   
     389  MOMENTS_RH        : SKIP # not ingested into DVO                   
     390  KRON_FLUX         : flux:kron:inst                                 
     391  KRON_FLUX_ERR     : flux:kron:inst:err                                     
     392  KRON_FLUX_INNER   : SKIP # not ingested into DVO                       
     393  KRON_FLUX_OUTER   : SKIP # not ingested into DVO                   
     394  SKY_LIMIT_RAD     : SKIP # not ingested into DVO                       
     395  SKY_LIMIT_FLUX    : SKIP # not ingested into DVO                       
     396  SKY_LIMIT_SLOPE   : SKIP # not ingested into DVO                       
     397  FLAGS             : phot_flags                                     
     398  FLAGS2            : SKIP # not ingested into DVO                       
     399  N_FRAMES          : SKIP # not ingested into DVO                   
     400end                                                                   
     401
     402# this list is good for PS1_SV3
     403list testfields_PS1_SV3
     404  IPP_IDET          : detid     
     405  X_PSF             : xccd     
     406  Y_PSF             : yccd     
     407  X_PSF_SIG         : xccd:err
     408  Y_PSF_SIG         : yccd:err
     409  POSANGLE          : posangle
     410  PLTSCALE          : platescale
     411  PSF_INST_MAG      : mag:inst 
     412  PSF_INST_MAG_SIG  : mag:err   
     413  PSF_INST_FLUX     : flux:inst
     414  PSF_INST_FLUX_SIG : flux:inst:err
     415  AP_MAG            : mag:aper:inst
     416  AP_MAG_RAW        : SKIP # not ingested into DVO
     417  AP_MAG_RADIUS     : SKIP # not ingested into DVO
     418  AP_FLUX           : flux:aper:inst 
     419  AP_FLUX_SIG       : flux:aper:inst:err
     420  CAL_PSF_MAG       : mag:cal
     421  CAL_PSF_MAG_SIG   : mag:cal:err
     422  RA_PSF            : ra
     423  DEC_PSF           : dec
     424  AP_NPIX           : SKIP # not ingested into DVO                   
     425  PEAK_FLUX_AS_MAG  : SKIP # not ingested into DVO                   
     426  SKY               : sky                                             
     427  SKY_SIGMA         : sky_err                                         
     428  PSF_CHISQ         : psf_chisq                                       
     429  CR_NSIGMA         : cr_nsigma                                       
     430  EXT_NSIGMA        : ext_nsigma                                     
     431  PSF_MAJOR         : fwhm_maj # XXX technically, this is not FWHM but Sigma
     432  PSF_MINOR         : fwhm_min # XXX technically, this is not FWHM but Sigma
     433  PSF_THETA         : theta # FAIL                                   
     434  PSF_CORE          : SKIP # not ingested into DVO
     435  PSF_FWHM_MAJ      : SKIP # not ingested into DVO
     436  PSF_FWHM_MIN      : SKIP # not ingested into DVO
     437  PSF_QF            : psf_qf                                         
     438  PSF_QF_PERFECT    : SKIP # not ingested into DVO                   
     439  PSF_NDOF          : psf_ndof                               
     440  PSF_NPIX          : psf_npix                               
     441  MOMENTS_XX        : Mxx                                             
     442  MOMENTS_XY        : Mxy                                             
     443  MOMENTS_YY        : Myy                                             
     444  MOMENTS_M3C       : SKIP # not ingested into DVO                   
     445  MOMENTS_M3S       : SKIP # not ingested into DVO                   
     446  MOMENTS_M4C       : SKIP # not ingested into DVO                   
     447  MOMENTS_M4S       : SKIP # not ingested into DVO                   
     448  MOMENTS_R1        : SKIP # not ingested into DVO                   
     449  MOMENTS_RH        : SKIP # not ingested into DVO                   
     450  KRON_FLUX         : flux:kron:inst                                 
     451  KRON_FLUX_ERR     : flux:kron:inst:err                                     
     452  KRON_FLUX_INNER   : SKIP # not ingested into DVO                       
     453  KRON_FLUX_OUTER   : SKIP # not ingested into DVO                   
     454  FLAGS             : phot_flags                                     
     455  FLAGS2            : SKIP # not ingested into DVO                       
     456  N_FRAMES          : SKIP # not ingested into DVO                   
     457end                                                                   
Note: See TracChangeset for help on using the changeset viewer.