IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42821


Ignore:
Timestamp:
May 8, 2025, 4:29:52 PM (12 months ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20230313

Location:
trunk/Ohana
Files:
1 deleted
74 edited
11 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/addstar/Makefile

    r42389 r42821  
    9696$(SRC)/find_matches_refstars.$(ARCH).o \
    9797$(SRC)/find_matches_closest_refstars.$(ARCH).o \
     98$(SRC)/photcodes.$(ARCH).o \
    9899$(SRC)/getgsc.$(ARCH).o \
    99100$(SRC)/getusno.$(ARCH).o \
    … …  
    151152$(SRC)/find_matches_refstars.$(ARCH).o \
    152153$(SRC)/find_subset.$(ARCH).o \
     154$(SRC)/photcodes.$(ARCH).o \
    153155$(SRC)/getgsc.$(ARCH).o \
    154156$(SRC)/getusno.$(ARCH).o \
    … …  
    181183$(SRC)/find_matches_refstars.$(ARCH).o \
    182184$(SRC)/find_subset.$(ARCH).o \
     185$(SRC)/photcodes.$(ARCH).o \
    183186$(SRC)/getgsc.$(ARCH).o \
    184187$(SRC)/getusno.$(ARCH).o \

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

    r42712 r42821  
    142142int     OLD_RESORT;
    143143int     READ_XRAD_DATA;
     144int     CONVOLVED_APERTURES;
    144145int     DIFF_WITH_INV;
    145146
    … …  
    344345void GetConfig (char *config, char *field, char *format, int N, void *ptr);
    345346
     347int isGPC1warp (int photcode);
     348int isGPC1stack (int photcode);
     349
    346350/**
    347351    there is an inconsistency to be resolved: fixed structures (like Image)

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

    r42477 r42821  
    12901290    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
    12911291    catalog->measure[i].detID      = ps1data[i].detID;
    1292 
     1292   
    12931293    dvo_lensing_init (&catalog->lensing[i]);
    12941294

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

    r42077 r42821  
    9797  // first entry and compare to the rest
    9898
    99   float fwhmValues[3];
    100   fwhmValues[0] = PSFfwhm[0];
    101   fwhmValues[1] = PSFfwhm[1];
    102   fwhmValues[2] = PSFfwhm[2];
    103 
    104   int i;
    10599  int Nap = 0;
    106   for (i = 0; i < catalog->Nmeasure; i++) {
    107     dvo_lensing_init (&catalog->lensing[i]);
    108 
    109     if (catalog->measure[i].detID < RadID[Nap]) {
    110       continue;
     100
     101  if (CONVOLVED_APERTURES) {
     102    // XXX code to support the 3 convolved apertures
     103
     104    float fwhmValues[3];
     105    fwhmValues[0] = PSFfwhm[0];
     106    fwhmValues[1] = PSFfwhm[1];
     107    fwhmValues[2] = PSFfwhm[2];
     108
     109    for (int i = 0; i < catalog->Nmeasure; i++) {
     110      dvo_lensing_init (&catalog->lensing[i]);
     111
     112      if (catalog->measure[i].detID < RadID[Nap]) {
     113        continue;
     114        // this is a psf measurement which does not have a radial aperture
     115      }
     116      if (catalog->measure[i].detID > RadID[Nap]) {
     117        myAbort("radial apertures for source not in psf list? sources out of order?  seems like a bug\n");
     118        // this could be a radial aperture which does not have a PSF source, but that is not possible
     119      }
     120
     121      // confirm the 3 FWHM values:
     122      myAssert (fwhmValues[0] == PSFfwhm[Nap+0], "FWHM mismatch %f vs %f", fwhmValues[0], PSFfwhm[Nap+0]);
     123      myAssert (fwhmValues[1] == PSFfwhm[Nap+1], "FWHM mismatch %f vs %f", fwhmValues[1], PSFfwhm[Nap+1]);
     124      myAssert (fwhmValues[2] == PSFfwhm[Nap+2], "FWHM mismatch %f vs %f", fwhmValues[2], PSFfwhm[Nap+2]);
     125
     126      // EAM 2022.02.17 : here is the comment for the PV3 load:
     127      // XXX this is all hard-wired and should make use of the headers.
     128      // psphot cmfs have 5 radial apertures:
     129      // array 0, 1, 2, 3, 4
     130      // SDSS  3, 4, 5, 6, 7
     131
     132      // EAM 2022.02.17 : here is the situation for UNIONS DR3:
     133      // we have 3 convolutions (raw, 6", 8")
     134      // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec
     135      // I am going to save (4, 16, 32) which have index of (0, 2, 3)
     136
     137# define RAD_0 0
     138# define RAD_1 2
     139# define RAD_2 3
     140      catalog->lensing[i]. F_ApR5    = AperFlux   [(Nap + 0)*Ncol + RAD_0];
     141      catalog->lensing[i].dF_ApR5    = AperFluxErr[(Nap + 0)*Ncol + RAD_0];
     142      catalog->lensing[i].sF_ApR5    = AperFluxStd[(Nap + 0)*Ncol + RAD_0];
     143      catalog->lensing[i].fF_ApR5    = AperFill   [(Nap + 0)*Ncol + RAD_0];
     144                                   
     145      catalog->lensing[i]. F_ApR6    = AperFlux   [(Nap + 0)*Ncol + RAD_1];
     146      catalog->lensing[i].dF_ApR6    = AperFluxErr[(Nap + 0)*Ncol + RAD_1];
     147      catalog->lensing[i].sF_ApR6    = AperFluxStd[(Nap + 0)*Ncol + RAD_1];
     148      catalog->lensing[i].fF_ApR6    = AperFill   [(Nap + 0)*Ncol + RAD_1];
     149                                   
     150      catalog->lensing[i]. F_ApR7    = AperFlux   [(Nap + 0)*Ncol + RAD_2];
     151      catalog->lensing[i].dF_ApR7    = AperFluxErr[(Nap + 0)*Ncol + RAD_2];
     152      catalog->lensing[i].sF_ApR7    = AperFluxStd[(Nap + 0)*Ncol + RAD_2];
     153      catalog->lensing[i].fF_ApR7    = AperFill   [(Nap + 0)*Ncol + RAD_2];
     154
     155      catalog->lensing[i]. F_ApR5_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_0];
     156      catalog->lensing[i].dF_ApR5_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_0];
     157      catalog->lensing[i]. F_ApR6_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_1];
     158      catalog->lensing[i].dF_ApR6_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_1];
     159      catalog->lensing[i]. F_ApR7_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_2];
     160      catalog->lensing[i].dF_ApR7_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_2];
     161
     162      catalog->lensing[i]. F_ApR5_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_0];
     163      catalog->lensing[i].dF_ApR5_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_0];
     164      catalog->lensing[i]. F_ApR6_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_1];
     165      catalog->lensing[i].dF_ApR6_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_1];
     166      catalog->lensing[i]. F_ApR7_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_2];
     167      catalog->lensing[i].dF_ApR7_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_2];
     168
     169      catalog->lensing[i].detID = catalog->measure[i].detID;
     170
     171      // XXX set the measure, object, etc ID values here
     172      // catalog->lensing[i].objID : set in find_matches_closest.c
     173      // catalog->lensing[i].catID : set in find_matches_closest.c
     174      // catalog->lensing[i].averef : set in find_matches_closest.c
     175      // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c
     176
     177      Nap += 3;
     178    }
     179  } else {
     180    // PSF values for unconvolve data are undefined ?
     181
     182    for (int i = 0; i < catalog->Nmeasure; i++) {
     183      dvo_lensing_init (&catalog->lensing[i]);
     184
    111185      // this is a psf measurement which does not have a radial aperture
     186      if (catalog->measure[i].detID < RadID[Nap]) { continue; }
     187
     188      // this could be a radial aperture which does not have a PSF source, but that is not possible
     189      if (catalog->measure[i].detID > RadID[Nap]) {
     190        myAbort("radial apertures for source not in psf list? sources out of order?  seems like a bug\n");
     191      }
     192
     193      // EAM 2022.02.17 : here is the comment for the PV3 load:
     194      // XXX this is all hard-wired and should make use of the headers.
     195      // psphot cmfs have 5 radial apertures:
     196      // array 0, 1, 2, 3, 4
     197      // SDSS  3, 4, 5, 6, 7
     198
     199      // EAM 2022.02.17 : here is the situation for UNIONS DR4:
     200      // we have 1 convolution (raw)
     201      // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec
     202      // I am going to save (4, 16, 32) which have index of (0, 2, 3)
     203
     204# define RAD_0 0
     205# define RAD_1 2
     206# define RAD_2 3
     207      catalog->lensing[i]. F_ApR5    = AperFlux   [(Nap + 0)*Ncol + RAD_0];
     208      catalog->lensing[i].dF_ApR5    = AperFluxErr[(Nap + 0)*Ncol + RAD_0];
     209      catalog->lensing[i].sF_ApR5    = AperFluxStd[(Nap + 0)*Ncol + RAD_0];
     210      catalog->lensing[i].fF_ApR5    = AperFill   [(Nap + 0)*Ncol + RAD_0];
     211                                   
     212      catalog->lensing[i]. F_ApR6    = AperFlux   [(Nap + 0)*Ncol + RAD_1];
     213      catalog->lensing[i].dF_ApR6    = AperFluxErr[(Nap + 0)*Ncol + RAD_1];
     214      catalog->lensing[i].sF_ApR6    = AperFluxStd[(Nap + 0)*Ncol + RAD_1];
     215      catalog->lensing[i].fF_ApR6    = AperFill   [(Nap + 0)*Ncol + RAD_1];
     216                                   
     217      catalog->lensing[i]. F_ApR7    = AperFlux   [(Nap + 0)*Ncol + RAD_2];
     218      catalog->lensing[i].dF_ApR7    = AperFluxErr[(Nap + 0)*Ncol + RAD_2];
     219      catalog->lensing[i].sF_ApR7    = AperFluxStd[(Nap + 0)*Ncol + RAD_2];
     220      catalog->lensing[i].fF_ApR7    = AperFill   [(Nap + 0)*Ncol + RAD_2];
     221
     222      catalog->lensing[i].detID = catalog->measure[i].detID;
     223
     224      // XXX set the measure, object, etc ID values here
     225      // catalog->lensing[i].objID : set in find_matches_closest.c
     226      // catalog->lensing[i].catID : set in find_matches_closest.c
     227      // catalog->lensing[i].averef : set in find_matches_closest.c
     228      // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c
     229
     230      Nap += 1;
    112231    }
    113     if (catalog->measure[i].detID > RadID[Nap]) {
    114       myAbort("radial apertures for source not in psf list? sources out of order?  seems like a bug\n");
    115       // this could be a radial aperture which does not have a PSF source, but that is not possible
    116     }
    117 
    118     // confirm the 3 FWHM values:
    119     myAssert (fwhmValues[0] == PSFfwhm[Nap+0], "FWHM mismatch %f vs %f", fwhmValues[0], PSFfwhm[Nap+0]);
    120     myAssert (fwhmValues[1] == PSFfwhm[Nap+1], "FWHM mismatch %f vs %f", fwhmValues[1], PSFfwhm[Nap+1]);
    121     myAssert (fwhmValues[2] == PSFfwhm[Nap+2], "FWHM mismatch %f vs %f", fwhmValues[2], PSFfwhm[Nap+2]);
    122 
    123     // EAM 2022.02.17 : here is the comment for the PV3 load:
    124     // XXX this is all hard-wired and should make use of the headers.
    125     // psphot cmfs have 5 radial apertures:
    126     // array 0, 1, 2, 3, 4
    127     // SDSS  3, 4, 5, 6, 7
    128 
    129     // EAM 2022.02.17 : here is the situation for UNIONS DR3:
    130     // we have 3 convolutions (raw, 6", 8")
    131     // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec
    132     // I am going to save (4, 16, 32) which have index of (0, 2, 3)
    133 
    134     # define RAD_0 0
    135     # define RAD_1 2
    136     # define RAD_2 3
    137     catalog->lensing[i]. F_ApR5    = AperFlux   [(Nap + 0)*Ncol + RAD_0];
    138     catalog->lensing[i].dF_ApR5    = AperFluxErr[(Nap + 0)*Ncol + RAD_0];
    139     catalog->lensing[i].sF_ApR5    = AperFluxStd[(Nap + 0)*Ncol + RAD_0];
    140     catalog->lensing[i].fF_ApR5    = AperFill   [(Nap + 0)*Ncol + RAD_0];
    141                                    
    142     catalog->lensing[i]. F_ApR6    = AperFlux   [(Nap + 0)*Ncol + RAD_1];
    143     catalog->lensing[i].dF_ApR6    = AperFluxErr[(Nap + 0)*Ncol + RAD_1];
    144     catalog->lensing[i].sF_ApR6    = AperFluxStd[(Nap + 0)*Ncol + RAD_1];
    145     catalog->lensing[i].fF_ApR6    = AperFill   [(Nap + 0)*Ncol + RAD_1];
    146                                    
    147     catalog->lensing[i]. F_ApR7    = AperFlux   [(Nap + 0)*Ncol + RAD_2];
    148     catalog->lensing[i].dF_ApR7    = AperFluxErr[(Nap + 0)*Ncol + RAD_2];
    149     catalog->lensing[i].sF_ApR7    = AperFluxStd[(Nap + 0)*Ncol + RAD_2];
    150     catalog->lensing[i].fF_ApR7    = AperFill   [(Nap + 0)*Ncol + RAD_2];
    151 
    152     catalog->lensing[i]. F_ApR5_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_0];
    153     catalog->lensing[i].dF_ApR5_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_0];
    154     catalog->lensing[i]. F_ApR6_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_1];
    155     catalog->lensing[i].dF_ApR6_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_1];
    156     catalog->lensing[i]. F_ApR7_C1 = AperFlux   [(Nap + 1)*Ncol + RAD_2];
    157     catalog->lensing[i].dF_ApR7_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_2];
    158 
    159     catalog->lensing[i]. F_ApR5_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_0];
    160     catalog->lensing[i].dF_ApR5_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_0];
    161     catalog->lensing[i]. F_ApR6_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_1];
    162     catalog->lensing[i].dF_ApR6_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_1];
    163     catalog->lensing[i]. F_ApR7_C2 = AperFlux   [(Nap + 2)*Ncol + RAD_2];
    164     catalog->lensing[i].dF_ApR7_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_2];
    165 
    166     catalog->lensing[i].detID = catalog->measure[i].detID;
    167 
    168     // XXX set the measure, object, etc ID values here
    169     // catalog->lensing[i].objID : set in find_matches_closest.c
    170     // catalog->lensing[i].catID : set in find_matches_closest.c
    171     // catalog->lensing[i].averef : set in find_matches_closest.c
    172     // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c
    173 
    174     Nap += 3;
    175232  }
    176233  myAssert (Nap == Nrow, "did we go too far???");

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

    r42712 r42821  
    9191    remove_argument (N, &argc, argv);
    9292    READ_XRAD_DATA = TRUE;
     93  }
     94  CONVOLVED_APERTURES = FALSE;
     95  if ((N = get_argument (argc, argv, "-convolved-ap"))) {
     96    remove_argument (N, &argc, argv);
     97    CONVOLVED_APERTURES = TRUE;
    9398  }
    9499  DIFF_WITH_INV = FALSE;
    … …  
    304309  }
    305310  /* force using single time from PHU */
     311  FORCE_SINGLE_TIME = FALSE;
     312  if ((N = get_argument (argc, argv, "-force-single-time"))) {
     313    FORCE_SINGLE_TIME = TRUE;
     314    remove_argument (N, &argc, argv);
     315  }
     316  /* accept bad header astrometry */
    306317  FORCE_SINGLE_TIME = FALSE;
    307318  if ((N = get_argument (argc, argv, "-force-single-time"))) {

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

    r40291 r42821  
    248248      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    249249      if (Nsec > -1) {
     250        if (isGPC1warp(tgtcat[0].measure[Nmeas].photcode)) {
     251          if (FALSE && isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfWrp)) {
     252            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfWrp  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     253            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     254            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapWrp   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     255          }
     256          goto set_average;
     257        }
     258        if (isGPC1stack(tgtcat[0].measure[Nmeas].photcode)) {
     259          if (FALSE && isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfStk)) {
     260            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfStk  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     261            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     262            tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapStk   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     263          }
     264          goto set_average;
     265        }
    250266        if (isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp)) {
    251           tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
    252         }
    253       }
     267          tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     268          tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     269          tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapChp   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     270        }
     271      }
     272    set_average:
    254273
    255274      /* Nm is updated, but not written out in -update mode (for existing entries)
    … …  
    352371      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    353372      if (Nsec > -1) {
    354         tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
    355       }
     373        if (isGPC1warp(tgtcat[0].measure[Nmeas].photcode)) {
     374          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfWrp  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     375          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     376          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapWrp   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     377          goto set_average_2;
     378        }
     379        if (isGPC1stack(tgtcat[0].measure[Nmeas].photcode)) {
     380          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfStk  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     381          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     382          // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapStk   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     383          goto set_average_2;
     384        }
     385        if (isnan(tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp)) {
     386          tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp  = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     387          tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON);
     388          tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapChp   = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER);
     389        }
     390      }
     391    set_average_2:
    356392
    357393      Nmeas ++;

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

    r42133 r42821  
    22# define DEBUG 1
    33# define MACHINE_GROUPS 1
     4# define MEMCHECK_VERBOSE 0
    45
    56// in parallel mode, the target database may be distributed, while the input database is

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r41341 r42821  
    497497      int outputIndex;
    498498      if (secfiltMap) {
     499        if (secfiltMap[j] < 0) continue;  // skip secfilt entries from input that are not in the output
    499500        outputIndex = (Nave * NsecfiltOut) + secfiltMap[j];
    500501      } else {

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

    r27435 r42821  
    3535PhotCode *photcode;
    3636
     37float EpochMJD;
     38time_t EpochUnix;
     39
    3740int  args                    PROTO((int argc, char **argv));
    3841int  ConfigInit              PROTO((int *argc, char **argv));

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

    r34459 r42821  
    3737    VERBOSE = TRUE;
    3838    remove_argument (N, &argc, argv);
     39  }
     40
     41  EpochUnix = 0;
     42  EpochMJD = NAN;
     43  if ((N = get_argument (argc, argv, "-epoch-mjd"))) {
     44    remove_argument (N, &argc, argv);
     45    EpochMJD = atof (argv[N]);
     46    remove_argument (N, &argc, argv);
     47    EpochUnix = ohana_mjd_to_sec (EpochMJD);
    3948  }
    4049

  • trunk/Ohana/src/getstar/src/write_getstar_ps1_dev_0.c

    r40291 r42821  
    4545    output[i].R        = average[i].R;
    4646    output[i].D        = average[i].D;
     47    if (isfinite(EpochMJD) && isfinite(average[i].uR) && isfinite(average[i].uD)) {
     48      // RA(epoch) = RA(mean) + uR*(t - Tmean)/cos(dec)
     49      // XXX how does this behave near dec = 90?
     50      // double dT = (EpochUnix - average[i].Tmean) / (86400*365.25);
     51      // XXX refcat 20220324.v0 does not have Tmean correctly set to 2016.0 = MJD 57388.0
     52      // NOTE average[i].Tmean is in UNIX, but here I am using the MJD of the Gaia DR3 reference epoch
     53      double dT = (EpochMJD - 57388.0) / 365.25;
     54      double dR = average[i].uR * dT / cos(RAD_DEG*average[i].D) / 3600.0;
     55      double dD = average[i].uD * dT / 3600.0;
     56      output[i].R += dR;
     57      output[i].D += dD;
     58    }
    4759
    4860    output[i].code     = average[i].flags;

  • trunk/Ohana/src/kapa2/Makefile

    r41344 r42821  
    1010SRC     =       $(HOME)/src
    1111INC     =       $(HOME)/include
     12DATA    =       $(DESTDATA)/kapa
    1213include ../../Makefile.Common
     14
     15%.c : %.c.in
     16        sed "s|@DATADIR@|$(DATA)|" $< > $@
    1317
    1418# programs may add their own internal requirements here
    … …  
    2125
    2226kapa: $(BIN)/kapa.$(ARCH)
    23 install: $(DESTBIN)/kapa
     27install: $(DESTBIN)/kapa cmaps
    2428
    2529PDF = \
    … …  
    100104$(SRC)/CheckButtons.$(ARCH).o             $(SRC)/InvertButton.$(ARCH).o       \
    101105$(SRC)/UpdatePointer.$(ARCH).o            $(SRC)/JPEGit24.$(ARCH).o           \
    102 $(SRC)/ButtonFunctions.$(ARCH).o    \
    103 $(SRC)/SetChannel.$(ARCH).o         \
     106$(SRC)/ButtonFunctions.$(ARCH).o          $(SRC)/FindColormap.$(ARCH).o       \
    104107$(SRC)/SetColorScale.$(ARCH).o            $(SRC)/ColorCube.$(ARCH).o          \
    105 $(SRC)/ColorHistogram.$(ARCH).o           $(SRC)/CreateWide.$(ARCH).o
     108$(SRC)/ColorHistogram.$(ARCH).o           $(SRC)/CreateWide.$(ARCH).o         \
     109$(SRC)/SetChannel.$(ARCH).o
    106110
    107111OBJ  =  $(KAPA) $(PDF) $(BDRAW) $(PSFILES)
     112
     113cmaps:
     114        @echo "installing colormaps for kapa"
     115        @if [ ! -d $(DATA)/colormaps ]; then mkdir -p $(DATA)/colormaps; fi
     116        @if [   -d $(HOME)/colormaps ]; then rm -f $(HOME)/colormaps/*~; fi
     117        @if [   -d $(HOME)/colormaps ]; then rm -f $(HOME)/colormaps/#*; fi
     118        @if [   -d $(HOME)/colormaps ]; then for i in `find $(HOME)/colormaps -name CVS -prune -o -type f -print`; do cp -f $$i $(DATA)/colormaps; done; fi
    108119
    109120# dependancy rules for include files ########################
    … …  
    116127
    117128$(BIN)/kapa.$(ARCH): $(OBJ)
     129
     130.PRECIOUS: $(SRC)/FindColormap.c

  • trunk/Ohana/src/kapa2/include/constants.h

    r40501 r42821  
    5454} KapaChannels;
    5555
     56typedef enum {
     57  KAPA_CM_NONE,
     58  KAPA_CM_RAW,
     59  KAPA_CM_CSV,
     60  KAPA_CM_CET,
     61  KAPA_CM_CET_REV,
     62  KAPA_CM_LEGACY,
     63  KAPA_CM_STATIC,
     64} KapaColorMapMode;
     65
    5666// use an enum to identify the 3 dimensions:
    5767typedef enum {

  • trunk/Ohana/src/kapa2/include/prototypes.h

    r41344 r42821  
    328328int ResizeByImage (int sock);
    329329
     330char *ParseColormapName (char *name, KapaColorMapMode *mode);
     331int CheckFileAccess (char *name);
     332char *FindColormap (char *name);
     333char *ParseColormapModeFromEnum (KapaColorMapMode mode);
     334void ListColormaps (void);

  • trunk/Ohana/src/kapa2/src

    • Property svn:ignore set to
      FindColormap.c

  • trunk/Ohana/src/kapa2/src/SetColormap.c

    r41341 r42821  
    2121  float scale, blueRef, redRef, greenRef;
    2222  Graphic *graphic;
     23
     24  if (inName && !strcasecmp(inName, "list")) {
     25    ListColormaps ();
     26    return TRUE;
     27  }
    2328
    2429  graphic = GetGraphic();
    … …  
    214219
    215220  /* anuenue */
    216   if (!strncmp (graphic->colormapName, "file:", 5) ||
    217       !strncmp (graphic->colormapName, "lgcy:", 5) ||
    218       !strncmp (graphic->colormapName, "cetf:", 5) ||
    219       !strncmp (graphic->colormapName, "cetr:", 5) ||
    220       !strncmp (graphic->colormapName, "csvf:", 5)) {
    221 
    222     FILE *f = fopen (&graphic->colormapName[5], "r");
    223     if (!f) {
    224       fprintf (stderr, "failed to open colormap file %s\n", &graphic->colormapName[5]);
     221  KapaColorMapMode mode = KAPA_CM_STATIC;
     222
     223  // we should already have handled the static colormaps by this point
     224  char *basename = ParseColormapName (graphic->colormapName, &mode);
     225  if (mode == KAPA_CM_STATIC) {
     226    fprintf (stderr, "unknown static colormap name %s\n", graphic->colormapName);
     227    free (basename);
     228    return FALSE;
     229  }   
     230
     231  // look in local directory and installation directory
     232  char *truename = FindColormap (basename);
     233  FREE (basename);
     234
     235  if (!truename) {
     236      fprintf (stderr, "cannot find colormap file %s\n", &graphic->colormapName[5]);
    225237      return FALSE;
    226     }
    227     char line[1024];
    228 
    229     int lastIndex = 0;
    230     float lastRed = 0.0;
    231     float lastBlue = 0.0;
    232     float lastGreen = 0.0;
    233 
    234     int isCSV = !strncmp (graphic->colormapName, "csvf:", 5);
    235     int isCET = !strncmp (graphic->colormapName, "cetf:", 5);
    236     int isCETRev = !strncmp (graphic->colormapName, "cetr:", 5);
    237     int isLegacy = !strncmp (graphic->colormapName, "lgcy:", 5);
    238 
    239     float fracIndex, fracRed, fracBlue, fracGreen;
    240 
    241     if (isCET || isCETRev) fracIndex = 0.0;
    242 
    243     while (scan_line_maxlen (f, line, 1024) != EOF) {
    244       // file contains f R G B : f = 0.0 - 1.0, R,G,B = 0.0 - 1.0
    245 
    246       int Nscan;
    247       if (isCSV) {
    248         Nscan = sscanf (line, "%f,%f,%f,%f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
    249         if (Nscan != 4) continue;
    250       }
    251       if (isCET || isCETRev) {
    252         Nscan = sscanf (line, "%f,%f,%f", &fracRed, &fracGreen, &fracBlue);
    253         fracIndex += 1.0 / 256.0;
    254         if (Nscan != 3) continue;
    255       }
    256       if (isLegacy) {
    257         Nscan = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracBlue, &fracGreen);
    258         if (Nscan != 4) continue;
    259       }
    260       if (!isCSV && !isLegacy && !isCET && !isCETRev) {
    261         Nscan = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
    262         if (Nscan != 4) continue;
    263       }
    264 
    265       int nextIndex = fracIndex * MaxValue;
    266       if (nextIndex <= lastIndex) {
    267         lastRed = fracRed;
    268         lastBlue = fracBlue;
    269         lastGreen = fracGreen;
    270         continue;
    271       }
    272       float dRdI = (fracRed   - lastRed)   / (float) (nextIndex - lastIndex);
    273       float dBdI = (fracBlue  - lastBlue)  / (float) (nextIndex - lastIndex);
    274       float dGdI = (fracGreen - lastGreen) / (float) (nextIndex - lastIndex);
    275       for (i = lastIndex; i < nextIndex; i++) {
    276         float redValue   = 0xffff * (dRdI * (i - lastIndex) + lastRed);
    277         float blueValue  = 0xffff * (dBdI * (i - lastIndex) + lastBlue);
    278         float greenValue = 0xffff * (dGdI * (i - lastIndex) + lastGreen);
    279         SETVALUE (graphic[0].cmap[i].red,   redValue,   0, 0xffff);
    280         SETVALUE (graphic[0].cmap[i].blue,  blueValue,  0, 0xffff);
    281         SETVALUE (graphic[0].cmap[i].green, greenValue, 0, 0xffff);
    282       }
     238  }
     239 
     240  FILE *f = fopen (truename, "r");
     241  if (!f) {
     242    fprintf (stderr, "failed to open colormap file %s\n", truename);
     243    FREE (truename);
     244    return FALSE;
     245  }
     246  FREE (truename);
     247
     248  char line[1024];
     249 
     250  int lastIndex = 0;
     251  float lastRed = 0.0;
     252  float lastBlue = 0.0;
     253  float lastGreen = 0.0;
     254 
     255  float fracIndex, fracRed, fracBlue, fracGreen;
     256
     257  if (mode == KAPA_CM_CET || mode == KAPA_CM_CET_REV) fracIndex = 0.0;
     258
     259  int checkFormat = TRUE;
     260
     261  while (scan_line_maxlen (f, line, 1024) != EOF) {
     262    // file contains f R G B : f = 0.0 - 1.0, R,G,B = 0.0 - 1.0
     263    // skip any lines with comment character:
     264
     265    char *p = line;
     266    while ((*p != 0) && OHANA_WHITESPACE(*p)) { p++; }
     267    if (*p == 0) { continue; }
     268    if (*p == '#') { continue; }
     269
     270    // attempt to confirm the format (and warn if it does not seem to match)
     271    if (checkFormat) {
     272      int Nscan1 = sscanf (line, "%f,%f,%f,%f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
     273      int Nscan2 = sscanf (line, "%f,%f,%f", &fracRed, &fracGreen, &fracBlue);
     274      int Nscan3 = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracBlue, &fracGreen);
     275
     276      if (Nscan1 == 4) {
     277        if (mode != KAPA_CM_CSV) { fprintf (stderr, "warning: %s appears to be format csvf: not %s\n line: %s\n", graphic->colormapName, ParseColormapModeFromEnum(mode), line); }
     278      }
     279      if (Nscan2 == 3) {
     280        if ((mode != KAPA_CM_CET) && (mode != KAPA_CM_CET_REV)) { fprintf (stderr, "warning: %s appears to be format cetf: or cetf: not %s\n line: %s\n", graphic->colormapName, ParseColormapModeFromEnum(mode), line); }
     281      }
     282      if (Nscan3 == 4) {
     283        if (mode != KAPA_CM_RAW) { fprintf (stderr, "warning: %s appears to be format file: not %s\n line: %s\n", graphic->colormapName, ParseColormapModeFromEnum(mode), line); }
     284      }
     285      checkFormat = FALSE;
     286    }
     287
     288    int Nscan;
     289    if (mode == KAPA_CM_CSV) {
     290      Nscan = sscanf (line, "%f,%f,%f,%f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
     291      if (Nscan != 4) continue;
     292    }
     293    if ((mode == KAPA_CM_CET || mode == KAPA_CM_CET_REV)) {
     294      Nscan = sscanf (line, "%f,%f,%f", &fracRed, &fracGreen, &fracBlue);
     295      fracIndex += 1.0 / 256.0;
     296      if (Nscan != 3) continue;
     297    }
     298    if (mode == KAPA_CM_LEGACY) {
     299      Nscan = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracBlue, &fracGreen);
     300      if (Nscan != 4) continue;
     301    }
     302    if (mode == KAPA_CM_RAW) {
     303      Nscan = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
     304      if (Nscan != 4) continue;
     305    }
     306   
     307    int nextIndex = fracIndex * MaxValue;
     308    if (nextIndex <= lastIndex) {
    283309      lastRed = fracRed;
    284310      lastBlue = fracBlue;
    285311      lastGreen = fracGreen;
    286       lastIndex = nextIndex;
    287     }
    288     fclose (f);
    289 
    290     if ((int) (MaxValue*fracIndex) < MaxValue) {
    291       float dRdI = (1.0 - lastRed)   / (float) (MaxValue - lastIndex);
    292       float dBdI = (1.0 - lastBlue)  / (float) (MaxValue - lastIndex);
    293       float dGdI = (1.0 - lastGreen) / (float) (MaxValue - lastIndex);
    294       for (i = lastIndex; i < MaxValue; i++) {
    295         float redValue   = 0xffff * (dRdI * (i - lastIndex) + lastRed);
    296         float blueValue  = 0xffff * (dBdI * (i - lastIndex) + lastBlue);
    297         float greenValue = 0xffff * (dGdI * (i - lastIndex) + lastGreen);
    298         SETVALUE (graphic[0].cmap[i].red,   redValue,   0, 0xffff);
    299         SETVALUE (graphic[0].cmap[i].blue,  blueValue,  0, 0xffff);
    300         SETVALUE (graphic[0].cmap[i].green, greenValue, 0, 0xffff);
    301       }
    302     }
    303 
    304     // reverse the color sequence:
    305     if (isCETRev) {
    306       for (i = 0; i < MaxValue / 2; i++) { 
    307         unsigned short tmp;
    308         // fprintf (stderr, "swap: %d %d\n", graphic[0].cmap[i].red, graphic[0].cmap[MaxValue - 1 - i].red);
    309         tmp = graphic[0].cmap[i].red;   graphic[0].cmap[i].red   = graphic[0].cmap[MaxValue - 1 - i].red;   graphic[0].cmap[MaxValue - 1 - i].red   = tmp;
    310         tmp = graphic[0].cmap[i].blue;  graphic[0].cmap[i].blue  = graphic[0].cmap[MaxValue - 1 - i].blue;  graphic[0].cmap[MaxValue - 1 - i].blue  = tmp;
    311         tmp = graphic[0].cmap[i].green; graphic[0].cmap[i].green = graphic[0].cmap[MaxValue - 1 - i].green; graphic[0].cmap[MaxValue - 1 - i].green = tmp;
    312       }
    313     }
    314 
    315     goto store_colors;
    316   }
    317 
    318   return (FALSE);
     312      continue;
     313    }
     314    float dRdI = (fracRed   - lastRed)   / (float) (nextIndex - lastIndex);
     315    float dBdI = (fracBlue  - lastBlue)  / (float) (nextIndex - lastIndex);
     316    float dGdI = (fracGreen - lastGreen) / (float) (nextIndex - lastIndex);
     317    for (i = lastIndex; i < nextIndex; i++) {
     318      float redValue   = 0xffff * (dRdI * (i - lastIndex) + lastRed);
     319      float blueValue  = 0xffff * (dBdI * (i - lastIndex) + lastBlue);
     320      float greenValue = 0xffff * (dGdI * (i - lastIndex) + lastGreen);
     321      SETVALUE (graphic[0].cmap[i].red,   redValue,   0, 0xffff);
     322      SETVALUE (graphic[0].cmap[i].blue,  blueValue,  0, 0xffff);
     323      SETVALUE (graphic[0].cmap[i].green, greenValue, 0, 0xffff);
     324    }
     325    lastRed = fracRed;
     326    lastBlue = fracBlue;
     327    lastGreen = fracGreen;
     328    lastIndex = nextIndex;
     329  }
     330  fclose (f);
     331
     332  if ((int) (MaxValue*fracIndex) < MaxValue) {
     333    float dRdI = (1.0 - lastRed)   / (float) (MaxValue - lastIndex);
     334    float dBdI = (1.0 - lastBlue)  / (float) (MaxValue - lastIndex);
     335    float dGdI = (1.0 - lastGreen) / (float) (MaxValue - lastIndex);
     336    for (i = lastIndex; i < MaxValue; i++) {
     337      float redValue   = 0xffff * (dRdI * (i - lastIndex) + lastRed);
     338      float blueValue  = 0xffff * (dBdI * (i - lastIndex) + lastBlue);
     339      float greenValue = 0xffff * (dGdI * (i - lastIndex) + lastGreen);
     340      SETVALUE (graphic[0].cmap[i].red,   redValue,   0, 0xffff);
     341      SETVALUE (graphic[0].cmap[i].blue,  blueValue,  0, 0xffff);
     342      SETVALUE (graphic[0].cmap[i].green, greenValue, 0, 0xffff);
     343    }
     344  }
     345
     346  // reverse the color sequence:
     347  if (mode == KAPA_CM_CET_REV) {
     348    for (i = 0; i < MaxValue / 2; i++) { 
     349      unsigned short tmp;
     350      // fprintf (stderr, "swap: %d %d\n", graphic[0].cmap[i].red, graphic[0].cmap[MaxValue - 1 - i].red);
     351      tmp = graphic[0].cmap[i].red;   graphic[0].cmap[i].red   = graphic[0].cmap[MaxValue - 1 - i].red;   graphic[0].cmap[MaxValue - 1 - i].red   = tmp;
     352      tmp = graphic[0].cmap[i].blue;  graphic[0].cmap[i].blue  = graphic[0].cmap[MaxValue - 1 - i].blue;  graphic[0].cmap[MaxValue - 1 - i].blue  = tmp;
     353      tmp = graphic[0].cmap[i].green; graphic[0].cmap[i].green = graphic[0].cmap[MaxValue - 1 - i].green; graphic[0].cmap[MaxValue - 1 - i].green = tmp;
     354    }
     355  }
    319356
    320357 store_colors:
    … …  
    333370  return (TRUE);
    334371}
    335 

  • trunk/Ohana/src/libdvo/include/dvodb.h

    r42389 r42821  
    547547
    548548dbValue      dbExtractAverages      PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, GalPhot *galphot, dbField *field));
    549 dbValue      dbExtractMeasures      PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensing *lensing, StarPar *starpar, dbField *field));
     549dbValue      dbExtractMeasures      PROTO((Average *average, SecFilt *secfilt, Lensobj *lensobj, Measure *measure, Lensing *lensing, StarPar *starpar, GalPhot *galphot, dbField *field));
    550550dbValue      dbExtractImages        PROTO((Image *image, off_t Nimage, off_t N, dbField *field));
    551551

  • trunk/Ohana/src/libdvo/include/libdvo_astro.h

    r42075 r42821  
    2020  PROJ_TNX, // zenithal
    2121  PROJ_DIS, // zenithal (TAN + polyterms)
     22  PROJ_CAR, // cartesian (Plat-Carre)
    2223  PROJ_LIN, // cartesian
    2324  PROJ_PLY, // cartesian (allow polyterms)

  • trunk/Ohana/src/libdvo/src/cmf-ps1-v5-r0.c

    r42477 r42821  
    9696    SWAP_BYTE (236); // SRC_CHIP_Y     
    9797    SWAP_BYTE (238); // PADDING3       
    98     SWAP_WORD (242); // MOMENTS_R1
    99     SWAP_WORD (246); // MOMENTS_RH
    100     SWAP_WORD (250); // KRON_FLUX
    101     SWAP_WORD (254); // KRON_FLUX_ERR
    102     SWAP_WORD (258); // KRON_FLUX_INNER
    103     SWAP_WORD (262); // KRON_FLUX_OUTER
    104     SWAP_WORD (266); // SKY_LIMIT_RAD
    105     SWAP_WORD (270); // SKY_LIMIT_FLUX
    106     SWAP_WORD (274); // SKY_LIMIT_SLOPE
    107     SWAP_WORD (278); // FLAGS
     98    SWAP_WORD (240); // MOMENTS_R1
     99    SWAP_WORD (244); // MOMENTS_RH
     100    SWAP_WORD (248); // KRON_FLUX
     101    SWAP_WORD (252); // KRON_FLUX_ERR
     102    SWAP_WORD (256); // KRON_FLUX_INNER
     103    SWAP_WORD (260); // KRON_FLUX_OUTER
     104    SWAP_WORD (264); // SKY_LIMIT_RAD
     105    SWAP_WORD (268); // SKY_LIMIT_FLUX
     106    SWAP_WORD (272); // SKY_LIMIT_SLOPE
     107    SWAP_WORD (276); // FLAGS
    108108    SWAP_WORD (280); // FLAGS2
    109109    SWAP_BYTE (284); // N_FRAMES

  • trunk/Ohana/src/libdvo/src/coordops.c

    r42389 r42821  
    10721072  if (!strncmp(&ctype[4], "-TNX", 4)) return PROJ_TNX;
    10731073  if (!strncmp(&ctype[4], "-DIS", 4)) return PROJ_DIS;
     1074  if (!strncmp(&ctype[4], "-CAR", 4)) return PROJ_CAR;
    10741075  if (!strncmp(&ctype[4], "-LIN", 4)) return PROJ_LIN;
    10751076  if (!strncmp(&ctype[0], "GENE", 4)) return PROJ_LIN; // note ctype[0]
    … …  
    10941095    case PROJ_TNX: strcpy(&ctype[4], "-TNX"); return TRUE;
    10951096    case PROJ_DIS: strcpy(&ctype[4], "-DIS"); return TRUE;
     1097    case PROJ_CAR: strcpy(&ctype[4], "-CAR"); return TRUE;
    10961098    case PROJ_LIN: strcpy(&ctype[4], "-LIN"); return TRUE;
    10971099    case PROJ_PLY: strcpy(&ctype[4], "-PLY"); return TRUE;
    … …  
    11181120    case PROJ_DIS:
    11191121      return PROJ_MODE_ZENITHAL;
     1122    case PROJ_CAR:
    11201123    case PROJ_LIN:
    11211124    case PROJ_PLY:

  • trunk/Ohana/src/libdvo/src/dbExtractMeasures.c

    r42102 r42821  
    8585
    8686/* return measure.field based on the selection */
    87 dbValue dbExtractMeasures (Average *average, SecFilt *secfilt, Measure *measure, Lensing *lensing, StarPar *starpar, dbField *field) {
     87dbValue dbExtractMeasures (Average *average, SecFilt *secfilt, Lensobj *lensobj, Measure *measure, Lensing *lensing, StarPar *starpar, GalPhot *galphot, dbField *field) {
    8888
    8989  int Nsec;
    … …  
    264264         
    265265          // the following values come from the mean lensobj table
    266         case MAG_OPTION_X11_SM_OBJ:
    267         case MAG_OPTION_X12_SM_OBJ:
    268         case MAG_OPTION_X22_SM_OBJ:
    269         case MAG_OPTION_E1_SM_OBJ:
    270         case MAG_OPTION_E2_SM_OBJ:
    271         case MAG_OPTION_X11_SH_OBJ:
    272         case MAG_OPTION_X12_SH_OBJ:
    273         case MAG_OPTION_X22_SH_OBJ:
    274         case MAG_OPTION_E1_SH_OBJ:
    275         case MAG_OPTION_E2_SH_OBJ:
    276         case MAG_OPTION_X11_SM_PSF:
    277         case MAG_OPTION_X12_SM_PSF:
    278         case MAG_OPTION_X22_SM_PSF:
    279         case MAG_OPTION_E1_SM_PSF:
    280         case MAG_OPTION_E2_SM_PSF:
    281         case MAG_OPTION_X11_SH_PSF:
    282         case MAG_OPTION_X12_SH_PSF:
    283         case MAG_OPTION_X22_SH_PSF:
    284         case MAG_OPTION_E1_SH_PSF:
    285         case MAG_OPTION_E2_SH_PSF:
    286 
    287         case MAG_OPTION_E1_PSF:
    288         case MAG_OPTION_E2_PSF:
    289 
    290         case MAG_OPTION_F_AP_R5:
    291         case MAG_OPTION_F_ERR_AP_R5:
    292         case MAG_OPTION_F_STDEV_AP_R5:
    293         case MAG_OPTION_F_FILL_AP_R5:
    294         case MAG_OPTION_F_AP_R6:
    295         case MAG_OPTION_F_ERR_AP_R6:
    296         case MAG_OPTION_F_STDEV_AP_R6:
    297         case MAG_OPTION_F_FILL_AP_R6:
    298         case MAG_OPTION_F_AP_R7:
    299         case MAG_OPTION_F_ERR_AP_R7:
    300         case MAG_OPTION_F_STDEV_AP_R7:
    301         case MAG_OPTION_F_FILL_AP_R7:
    302         case MAG_OPTION_E1:
    303         case MAG_OPTION_E2:
    304 
    305         case MAG_OPTION_GAL_MAG:       
    306         case MAG_OPTION_GAL_MAG_ERR:   
    307         case MAG_OPTION_GAL_MAJ:       
    308         case MAG_OPTION_GAL_MAJ_ERR:   
    309         case MAG_OPTION_GAL_MIN:       
    310         case MAG_OPTION_GAL_MIN_ERR:   
    311         case MAG_OPTION_GAL_THETA:     
    312         case MAG_OPTION_GAL_THETA_ERR:
    313         case MAG_OPTION_GAL_INDEX:     
    314         case MAG_OPTION_GAL_CHISQ:     
    315         case MAG_OPTION_GAL_NPIX:     
    316         case MAG_OPTION_GAL_FLAGS:     
    317         case MAG_OPTION_GAL_TYPE:     
    318         case MAG_OPTION_GAL_OBJ_ID:
    319         case MAG_OPTION_GAL_CAT_ID:
    320         case MAG_OPTION_GAL_DET_ID:
    321         case MAG_OPTION_GAL_IMAGE_ID:
    322 
    323         case MAG_OPTION_NONE:
    324           break;
     266          // g:X11_SM_OBJ, etc, are only valid for average
     267        case MAG_OPTION_X11_SM_OBJ: { value.Flt = LensValue_X11_sm_obj (field->photcode, lensobj, average->Nlensobj); break; }
     268        case MAG_OPTION_X12_SM_OBJ: { value.Flt = LensValue_X12_sm_obj (field->photcode, lensobj, average->Nlensobj); break; }
     269        case MAG_OPTION_X22_SM_OBJ: { value.Flt = LensValue_X22_sm_obj (field->photcode, lensobj, average->Nlensobj); break; }
     270        case MAG_OPTION_E1_SM_OBJ:  { value.Flt = LensValue_E1_sm_obj  (field->photcode, lensobj, average->Nlensobj); break; }
     271        case MAG_OPTION_E2_SM_OBJ:  { value.Flt = LensValue_E2_sm_obj  (field->photcode, lensobj, average->Nlensobj); break; }
     272
     273        case MAG_OPTION_X11_SH_OBJ: { value.Flt = LensValue_X11_sh_obj (field->photcode, lensobj, average->Nlensobj); break; }
     274        case MAG_OPTION_X12_SH_OBJ: { value.Flt = LensValue_X12_sh_obj (field->photcode, lensobj, average->Nlensobj); break; }
     275        case MAG_OPTION_X22_SH_OBJ: { value.Flt = LensValue_X22_sh_obj (field->photcode, lensobj, average->Nlensobj); break; }
     276        case MAG_OPTION_E1_SH_OBJ:  { value.Flt = LensValue_E1_sh_obj  (field->photcode, lensobj, average->Nlensobj); break; }
     277        case MAG_OPTION_E2_SH_OBJ:  { value.Flt = LensValue_E2_sh_obj  (field->photcode, lensobj, average->Nlensobj); break; }
     278
     279        case MAG_OPTION_X11_SM_PSF: { value.Flt = LensValue_X11_sm_psf (field->photcode, lensobj, average->Nlensobj); break; }
     280        case MAG_OPTION_X12_SM_PSF: { value.Flt = LensValue_X12_sm_psf (field->photcode, lensobj, average->Nlensobj); break; }
     281        case MAG_OPTION_X22_SM_PSF: { value.Flt = LensValue_X22_sm_psf (field->photcode, lensobj, average->Nlensobj); break; }
     282        case MAG_OPTION_E1_SM_PSF:  { value.Flt = LensValue_E1_sm_psf  (field->photcode, lensobj, average->Nlensobj); break; }
     283        case MAG_OPTION_E2_SM_PSF:  { value.Flt = LensValue_E2_sm_psf  (field->photcode, lensobj, average->Nlensobj); break; }
     284
     285        case MAG_OPTION_X11_SH_PSF: { value.Flt = LensValue_X11_sh_psf (field->photcode, lensobj, average->Nlensobj); break; }
     286        case MAG_OPTION_X12_SH_PSF: { value.Flt = LensValue_X12_sh_psf (field->photcode, lensobj, average->Nlensobj); break; }
     287        case MAG_OPTION_X22_SH_PSF: { value.Flt = LensValue_X22_sh_psf (field->photcode, lensobj, average->Nlensobj); break; }
     288        case MAG_OPTION_E1_SH_PSF:  { value.Flt = LensValue_E1_sh_psf  (field->photcode, lensobj, average->Nlensobj); break; }
     289        case MAG_OPTION_E2_SH_PSF:  { value.Flt = LensValue_E2_sh_psf  (field->photcode, lensobj, average->Nlensobj); break; }
     290
     291        case MAG_OPTION_F_AP_R5:       { value.Flt = LensValue_F_ApR5  (field->photcode, lensobj, average->Nlensobj); break; }
     292        case MAG_OPTION_F_ERR_AP_R5:   { value.Flt = LensValue_dF_ApR5 (field->photcode, lensobj, average->Nlensobj); break; }
     293        case MAG_OPTION_F_STDEV_AP_R5: { value.Flt = LensValue_sF_ApR5 (field->photcode, lensobj, average->Nlensobj); break; }
     294        case MAG_OPTION_F_FILL_AP_R5:  { value.Flt = LensValue_fF_ApR5 (field->photcode, lensobj, average->Nlensobj); break; }
     295                                                                                               
     296        case MAG_OPTION_F_AP_R6:       { value.Flt = LensValue_F_ApR6  (field->photcode, lensobj, average->Nlensobj); break; }
     297        case MAG_OPTION_F_ERR_AP_R6:   { value.Flt = LensValue_dF_ApR6 (field->photcode, lensobj, average->Nlensobj); break; }
     298        case MAG_OPTION_F_STDEV_AP_R6: { value.Flt = LensValue_sF_ApR6 (field->photcode, lensobj, average->Nlensobj); break; }
     299        case MAG_OPTION_F_FILL_AP_R6:  { value.Flt = LensValue_fF_ApR6 (field->photcode, lensobj, average->Nlensobj); break; }
     300
     301        case MAG_OPTION_F_AP_R7:       { value.Flt = LensValue_F_ApR7  (field->photcode, lensobj, average->Nlensobj); break; }
     302        case MAG_OPTION_F_ERR_AP_R7:   { value.Flt = LensValue_dF_ApR7 (field->photcode, lensobj, average->Nlensobj); break; }
     303        case MAG_OPTION_F_STDEV_AP_R7: { value.Flt = LensValue_sF_ApR7 (field->photcode, lensobj, average->Nlensobj); break; }
     304        case MAG_OPTION_F_FILL_AP_R7:  { value.Flt = LensValue_fF_ApR7 (field->photcode, lensobj, average->Nlensobj); break; }
     305
     306        case MAG_OPTION_E1:            { value.Flt = LensValue_E1      (field->photcode, lensobj, average->Nlensobj); break; }
     307        case MAG_OPTION_E2:            { value.Flt = LensValue_E2      (field->photcode, lensobj, average->Nlensobj); break; }
     308
     309        case MAG_OPTION_GAL_MAG:      { value.Flt = GalphotValue_GAL_MAG      (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     310        case MAG_OPTION_GAL_MAG_ERR:  { value.Flt = GalphotValue_GAL_MAG_ERR  (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     311        case MAG_OPTION_GAL_MAJ:      { value.Flt = GalphotValue_GAL_MAJ      (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     312        case MAG_OPTION_GAL_MAJ_ERR:  { value.Flt = GalphotValue_GAL_MAJ_ERR  (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     313        case MAG_OPTION_GAL_MIN:      { value.Flt = GalphotValue_GAL_MIN      (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     314        case MAG_OPTION_GAL_MIN_ERR:  { value.Flt = GalphotValue_GAL_MIN_ERR  (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     315        case MAG_OPTION_GAL_THETA:    { value.Flt = GalphotValue_GAL_THETA    (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     316        case MAG_OPTION_GAL_THETA_ERR:{ value.Flt = GalphotValue_GAL_THETA_ERR(field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     317        case MAG_OPTION_GAL_INDEX:    { value.Flt = GalphotValue_GAL_INDEX    (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     318        case MAG_OPTION_GAL_CHISQ:    { value.Flt = GalphotValue_GAL_CHISQ    (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     319        case MAG_OPTION_GAL_NPIX:     { value.Flt = GalphotValue_GAL_NPIX     (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     320        case MAG_OPTION_GAL_TYPE:     { value.Flt = GalphotValue_GAL_TYPE     (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     321        case MAG_OPTION_GAL_FLAGS:    { value.Flt = GalphotValue_GAL_FLAGS    (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     322
     323        case MAG_OPTION_GAL_OBJ_ID:   { value.Flt = GalphotValue_GAL_OBJ_ID   (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     324        case MAG_OPTION_GAL_CAT_ID:   { value.Flt = GalphotValue_GAL_CAT_ID   (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     325        case MAG_OPTION_GAL_DET_ID:   { value.Flt = GalphotValue_GAL_DET_ID   (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     326        case MAG_OPTION_GAL_IMAGE_ID: { value.Flt = GalphotValue_GAL_IMAGE_ID (field->photcode, field->magClass, galphot, average->Ngalphot); break; }
     327
     328        case MAG_OPTION_E1_PSF: break;
     329        case MAG_OPTION_E2_PSF: break;
     330        case MAG_OPTION_NONE: break;
    325331      }
    326332      break;

  • trunk/Ohana/src/libdvo/src/dvo_photcode_utils.c

    r33649 r42821  
    11# include <dvo.h>
     2
     3/* the photcode table includes a number of SEC types, photcodes->Nsecfilt
     4   hashNsec[code] returns the number in the sequence 0 - Nsecfilt for a given SEC code
     5*/
     6
     7// these APIs originally just matched by number of the code, but they should
     8// match by the *name* of the photcode (e.g., g = g)
    29
    310// Create a map between the secfilt values from one photcode table to another
    … …  
    1118  // loop over entries in the input table
    1219  for (i = 0; i < NsecfiltIn; i++) {
    13     int code  = input[0].codeNsec[i];
     20    int inseq  = input[0].codeNsec[i];
     21    myAssert (inseq > -1, "invalid SEC");
     22    int incode   = input[0].hashcode[inseq];
    1423    int entry = -1;
    1524    // find the matching entry in the output table
    1625    for (j = 0; j < NsecfiltOut; j++) {
    17       int outcode = output[0].codeNsec[j];
    18       if (code == outcode) {
    19          entry = j;
    20          break;
    21       }
     26      int outseq = output[0].codeNsec[j];
     27      myAssert (outseq > -1, "invalid SEC");
     28      int outcode   = output[0].hashcode[outseq];
     29      if (strcmp (input[0].code[incode].name, output[0].code[outcode].name)) continue;
     30      entry = j;
     31      break;
    2232    }
    23     // if (entry == -1) {
    24     //   // entry is missing fail (no printfs in this file)
    25     //   free(map);
    26     //   return(FALSE);
    27     // }
    28 
    2933    // if entry is still -1, we will skip this one (not map into the output db)
    3034    map[i] = entry;

  • trunk/Ohana/src/libkapa/src/IOfuncs.c

    r36084 r42821  
    8989  if (Nbytes == 0) {
    9090    return TRUE;
     91  }
     92  if (Nbytes < 0) {
     93    fprintf (stderr, "error in kapa communication (KiiScanMessage) %s\n", buffer);
     94    return FALSE;
    9195  }
    9296

  • trunk/Ohana/src/libkapa/src/KapaWindow.c

    r41377 r42821  
    259259                  data[0].formatYm, data[0].formatYp);
    260260
     261  memset (&data[0].coords, 0, sizeof(data[0].coords));
     262
    261263  KiiScanMessage (fd, "%f %f %f %f",
    262264                  &data[0].coords.pc1_1, &data[0].coords.pc2_2,
    … …  
    301303  KiiSendCommand (fd, 4, "GSTY");
    302304 
    303   KapaScanGraphData (fd, data);
     305  if (!KapaScanGraphData (fd, data)) {
     306    fprintf (stderr, "error getting graph data KapaGetGraphData (sent GSTY)\n");
     307    return FALSE;
     308  }
    304309
    305310  KiiWaitAnswer (fd, "DONE");

  • trunk/Ohana/src/libohana/include/ohana.h

    r42389 r42821  
    170170
    171171/* if your build crashes on OFF_T_MODE, you probably need to add your 64bit hardware to this list */
    172 # ifdef _LARGEFILE_SOURCE
     172# if defined(_LARGEFILE_SOURCE) && !defined(OFF_T_FMT)
    173173#  define OFF_T_FMT "%jd"
    174174# endif
    175175
    176 # ifdef lin64
     176# if defined(lin64) && !defined(OFF_T_FMT)
    177177#  define OFF_T_FMT "%jd"
    178178# endif
    179179
    180180// mac is annoying
    181 # ifdef _DARWIN_C_SOURCE
     181# if defined(_DARWIN_C_SOURCE) && !defined(OFF_T_FMT)
    182182#  define OFF_T_FMT "%lld"
    183183# endif
    184 # ifdef darwin_x86
     184# if defined(darwin_x86) && !defined(OFF_T_FMT)
    185185#  define OFF_T_FMT "%lld"
    186186# endif
    187 # ifdef darwin
     187# if defined(darwin) && !defined(OFF_T_FMT)
    188188#  define OFF_T_FMT "%lld"
    189189# endif

  • trunk/Ohana/src/libohana/src/ohana_allocate.c

    r42074 r42821  
    448448  int status = TRUE;
    449449
     450  OhanaMemblock *prevBlock = NULL;
     451
    450452  while (thisBlock) {
    451453
    … …  
    475477    Nbytes += thisBlock->size;
    476478
     479    prevBlock = thisBlock; // save for debugging
    477480    thisBlock = thisBlock->prevBlock;
    478481  }
    … …  
    480483  if (Ntotal || VERBOSE) {
    481484    fprintf (stderr, "%zd memory blocks allocated (%zd bytes total), %zd good, %zd bad ", Ntotal, Nbytes, Ngood, Nbad);
    482     fprintf (stderr, "@ %s:%d (%s)\n", myFile, myLine, myFunc);
     485    fprintf (stderr, "@ %s:%d (%s) ", myFile, myLine, myFunc);
     486    fprintf (stderr, "(%zd penultimate block)\n", (size_t) prevBlock);
    483487  }
    484488

  • trunk/Ohana/src/opihi/cmd.astro/cdensify.c

    r41515 r42821  
    11# include "data.h"
     2# define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
    23
    34# define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
    … …  
    116117  PutCoords (&newcoords, &bf[0].header);
    117118 
     119  // use the mask to prevent double-counting
     120  ALLOCATE_PTR (mask, char, Nx*Ny);
     121
    118122  float scalescale = scale*scale;
    119123  float scale2 = (scale + 1.0) * (scale + 1.0);
    … …  
    142146    coords.crval1 = rn;
    143147    coords.crval2 = *d;
     148   
     149    // XXX do not oversample by more than a factor of 10
     150    float dXn = dX * MAX(fabs(dCOS(*d)), 0.1);
    144151
    145152    float F = 1.0;
    146153    if (vv) { F = isFloatScale ? Fs[i] : Is[i]; }
     154
     155    // reset the mask so we do not double-count
     156    memset (mask, 0, Nx*Ny);
    147157
    148158    switch (PSFTYPE) {
    … …  
    162172        break;
    163173      case IS_SQUARE:
    164         for (ix = -scale; ix <= scale; ix += dX) {
     174        for (ix = -scale; ix <= scale; ix += dXn) {
    165175          for (iy = -scale; iy <= scale; iy += dY) {
    166176            double rp, dp;
    … …  
    184194        break;
    185195      case IS_CIRCLE:
    186         for (ix = -scale; ix <= scale; ix += dX) {
    187           for (iy = -scale; iy <= scale; iy += dY) {
     196        for (ix = -scale - dXn; ix <= scale + dXn; ix += dXn) {
     197          for (iy = -scale - dY; iy <= scale + dY; iy += dY) {
    188198            float r2 = ix*ix + iy*iy;
    189199            double rp, dp;
    … …  
    199209            if (Xb < 0) continue;
    200210            if (Yb < 0) continue;
     211            if (mask[Xb + Yb*Nx]) continue;
     212            mask[Xb + Yb*Nx] = 1;
    201213            if (vv) {
    202214              val[Xb + Yb*Nx] += Normalize ? fCircle*F : F;
    … …  
    208220        break;
    209221      case IS_GAUSS:
    210         for (ix = -3.0*scale; ix <= 3.0*scale; ix += dX) {
     222        for (ix = -3.0*scale; ix <= 3.0*scale; ix += dXn) {
    211223          for (iy = -3.0*scale; iy <= 3.0*scale; iy += dY) {
    212224            float r2 = ix*ix + iy*iy;
    … …  
    232244    }
    233245  }
     246
     247  free (mask);
    234248  return (TRUE);
    235249}

  • trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss-pol.c

    r42078 r42821  
    5151  set_variable ("Zpk",  par[5]);
    5252  set_variable ("Sg",   par[6]);
     53
     54  /*
     55  set_variable ("dXg",   dpar[0]);
     56  set_variable ("dYg",   dpar[1]);
     57  set_variable ("dSXg",  NAN);
     58  set_variable ("dSYg",  NAN);
     59  set_variable ("dSXYg", NAN);
     60  set_variable ("dZpk",  dpar[5]);
     61  set_variable ("dSg",   dpar[6]);
     62  */
    5363}
    5464

  • trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c

    r42078 r42821  
    3939  set_variable ("Zpk",  par[5]);
    4040  set_variable ("Sg",   par[6]);
     41
     42/*
     43  set_variable ("dXg",   dpar[0]);
     44  set_variable ("dYg",   dpar[1]);
     45  set_variable ("dSXg",  2.35 / dpar[2]);
     46  set_variable ("dSYg",  2.35 / dpar[3]);
     47  set_variable ("dSXYg", dpar[4]);
     48  set_variable ("dZpk",  dpar[5]);
     49  set_variable ("dSg",   dpar[6]);
     50*/
    4151}
    4252

  • trunk/Ohana/src/opihi/cmd.astro/imfit-pgauss-psf.c

    r39457 r42821  
    3535  set_variable ("Zpk",  par[2]);
    3636  set_variable ("Sg",   par[3]);
     37
     38/*
     39  set_variable ("dXg",   dpar[0]);
     40  set_variable ("dYg",   dpar[1]);
     41  set_variable ("dZpk",  dpar[2]);
     42  set_variable ("dSg",   dpar[3]);
     43*/
    3744}
    3845

  • trunk/Ohana/src/opihi/cmd.astro/imfit-pgauss.c

    r39457 r42821  
    3737  set_variable ("Zpk",  par[5]);
    3838  set_variable ("Sg",   par[6]);
     39
     40/*
     41  set_variable ("dXg",   dpar[0]);
     42  set_variable ("dYg",   dpar[1]);
     43  set_variable ("dSXg",  2.35 / dpar[2]);
     44  set_variable ("dSYg",  2.35 / dpar[3]);
     45  set_variable ("dSXYg", dpar[4]);
     46  set_variable ("dZpk",  dpar[5]);
     47  set_variable ("dSg",   dpar[6]);
     48*/
    3949}
    4050

  • trunk/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c

    r42078 r42821  
    4141  set_variable ("Sg",   par[6]);
    4242  set_variable ("Sr", par[7]);
     43
     44/*
     45  set_variable ("dXg",   dpar[0]);
     46  set_variable ("dYg",   dpar[1]);
     47  set_variable ("dSXg",  2.35 / dpar[2]);
     48  set_variable ("dSYg",  2.35 / dpar[3]);
     49  set_variable ("dSXYg", dpar[4]);
     50  set_variable ("dZpk",  dpar[5]);
     51  set_variable ("dSg",   dpar[6]);
     52  set_variable ("dSr",   dpar[7]);
     53*/
    4354}
    4455

  • trunk/Ohana/src/opihi/cmd.astro/imfit-qfgauss.c

    r39457 r42821  
    4141  set_variable ("Zpk",  par[5]);
    4242  set_variable ("Sg",   par[6]);
     43
     44/*
     45  set_variable ("dXg",   dpar[0]);
     46  set_variable ("dYg",   dpar[1]);
     47  set_variable ("dSXg",  2.35 / dpar[2]);
     48  set_variable ("dSYg",  2.35 / dpar[3]);
     49  set_variable ("dSXYg", dpar[4]);
     50  set_variable ("dZpk",  dpar[5]);
     51  set_variable ("dSg",   dpar[6]);
     52*/
    4353}
    4454

  • trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss-psf.c

    r39457 r42821  
    3838  set_variable ("Zpk", par[2]);
    3939  set_variable ("Sg",  par[3]);
     40
     41  /*
     42  set_variable ("dXg",  dpar[0]);
     43  set_variable ("dYg",  dpar[1]);
     44  set_variable ("dZpk", dpar[2]);
     45  set_variable ("dSg",  dpar[3]);
     46  */
    4047}
    4148

  • trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c

    r42078 r42821  
    4141  set_variable ("Sg",   par[6]);
    4242  set_variable ("Sr",   par[7]);
     43
     44/*
     45  set_variable ("dXg",   dpar[0]);
     46  set_variable ("dYg",   dpar[1]);
     47  set_variable ("dSXg",  2.35 / dpar[2]);
     48  set_variable ("dSYg",  2.35 / dpar[3]);
     49  set_variable ("dSXYg", dpar[4]);
     50  set_variable ("dZpk",  dpar[5]);
     51  set_variable ("dSg",   dpar[6]);
     52  set_variable ("dSr",   dpar[7]);
     53*/
    4354}
    4455

  • trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c

    r42078 r42821  
    4242  set_variable ("Sg",   par[6]);
    4343  set_variable ("Npow", par[7]);
     44
     45/*
     46  set_variable ("dXg",   dpar[0]);
     47  set_variable ("dYg",   dpar[1]);
     48  set_variable ("dSXg",  2.35 / dpar[2]);
     49  set_variable ("dSYg",  2.35 / dpar[3]);
     50  set_variable ("dSXYg", dpar[4]);
     51  set_variable ("dZpk",  dpar[5]);
     52  set_variable ("dSg",   dpar[6]);
     53  set_variable ("dNpow",   dpar[7]);
     54*/
    4455}
    4556

  • trunk/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c

    r42078 r42821  
    8585  set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
    8686  set_variable ("SXYf", par[9]);
     87
     88/*
     89  set_variable ("dXg",   dpar[0]);
     90  set_variable ("dYg",   dpar[1]);
     91  set_variable ("dSXg",  2.35 * sqrt(2.0) / dpar[2]);
     92  set_variable ("dSYg",  2.35 * sqrt(2.0) / dpar[3]);
     93  set_variable ("dSXYg", dpar[4]);
     94  set_variable ("dZpk",  dpar[5]);
     95  set_variable ("dSg",   dpar[6]);
     96  set_variable ("dSXf", 2.35 * sqrt(2.0) / dpar[7]);
     97  set_variable ("dSYf", 2.35 * sqrt(2.0) / dpar[8]);
     98  set_variable ("dSXYf", dpar[9]);
     99*/
    87100}

  • trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c

    r42078 r42821  
    4242  set_variable ("Sg",   par[6]);
    4343  set_variable ("Sr",   par[7]);
     44
     45/*
     46  set_variable ("dXg",   dpar[0]);
     47  set_variable ("dYg",   dpar[1]);
     48  set_variable ("dSXg",  2.35 / dpar[2]);
     49  set_variable ("dSYg",  2.35 / dpar[3]);
     50  set_variable ("dSXYg", dpar[4]);
     51  set_variable ("dZpk",  dpar[5]);
     52  set_variable ("dSg",   dpar[6]);
     53  set_variable ("dSr",   dpar[7]);
     54*/
    4455}
    4556

  • trunk/Ohana/src/opihi/cmd.astro/imfit-sgauss.c

    r39457 r42821  
    4949  set_variable ("SYf", 2.35 / par[8]);
    5050  set_variable ("SXYf", par[9]);
     51
     52/*
     53  set_variable ("dXg",   dpar[0]);
     54  set_variable ("dYg",   dpar[1]);
     55  set_variable ("dSXg",  2.35 / dpar[2]);
     56  set_variable ("dSYg",  2.35 / dpar[3]);
     57  set_variable ("dSXYg", dpar[4]);
     58  set_variable ("dZpk",  dpar[5]);
     59  set_variable ("dSg",   dpar[6]);
     60  set_variable ("dSXf", 2.35 / dpar[7]);
     61  set_variable ("dSYf", 2.35 / dpar[8]);
     62  set_variable ("dSXYf", dpar[9]);
     63*/
    5164}
    5265

  • trunk/Ohana/src/opihi/cmd.astro/imfit-tgauss.c

    r39457 r42821  
    8080  set_variable ("SYf",  2.35 * sqrt(2.0) / par[8]);
    8181  set_variable ("SXYf", par[9]);
     82
     83/*
     84  set_variable ("dXg",   dpar[0]);
     85  set_variable ("dYg",   dpar[1]);
     86  set_variable ("dSXg",  2.35 * sqrt(2.0) / dpar[2]);
     87  set_variable ("dSYg",  2.35 * sqrt(2.0) / dpar[3]);
     88  set_variable ("dSXYg", dpar[4]);
     89  set_variable ("dZpk",  dpar[5]);
     90  set_variable ("dSg",   dpar[6]);
     91  set_variable ("dSXf",  2.35 * sqrt(2.0) / dpar[7]);
     92  set_variable ("dSYf",  2.35 * sqrt(2.0) / dpar[8]);
     93  set_variable ("dSXYf", dpar[9]);
     94*/
    8295}

  • trunk/Ohana/src/opihi/cmd.astro/imfit-vgauss.c

    r39457 r42821  
    6969  set_variable ("SXf", par[7]);
    7070  set_variable ("SYf", par[8]);
     71
     72  /*
     73  set_variable ("dXg",   dpar[0]);
     74  set_variable ("dYg",   dpar[1]);
     75  set_variable ("dSXg",  2.35 / dpar[2]);
     76  set_variable ("dSYg",  2.35 / dpar[3]);
     77  set_variable ("dSXYg", dpar[4]);
     78  set_variable ("dZpk",  dpar[5]);
     79  set_variable ("dSg",   dpar[6]);
     80  set_variable ("dSXf",  dpar[7]);
     81  set_variable ("dSYf",  dpar[8]);
     82  */
    7183}

  • trunk/Ohana/src/opihi/cmd.astro/imfit.c

    r42078 r42821  
    1313  }
    1414
     15  Buffer *var = NULL;
     16  if ((N = get_argument (argc, argv, "-var-image"))) {
     17    remove_argument (N, &argc, argv);
     18    var = SelectBuffer (argv[N], OLDBUFFER, TRUE);
     19    if (!var) {
     20      gprint (GP_ERR, "unknown buffer for variance %s\n", argv[N]);
     21      FREE (Save);
     22      return (FALSE);
     23    }
     24    remove_argument (N, &argc, argv);
     25  }
     26
    1527  int Insert = FALSE;
    1628  if ((N = get_argument (argc, argv, "-insert"))) {
    … …  
    1830    Insert = TRUE;
    1931    if (Save) { gprint (GP_ERR, "-save and -insert are mutually exclusive\n"); free (Save); return (FALSE); }
     32    if (var)  { gprint (GP_ERR, "-save and -var-image are mutually exclusive\n"); return (FALSE); }
    2033  }
    2134
    … …  
    4053    Gain = atof(argv[N]);
    4154    remove_argument (N, &argc, argv);
     55    if (var) { gprint (GP_ERR, "warning: -var-image selected, -gain will have no effect\n"); }
    4256  }
    4357
    … …  
    4862    RDnoise = atof(argv[N]);
    4963    remove_argument (N, &argc, argv);
     64    if (var) { gprint (GP_ERR, "warning: -var-image selected, -rdnoise will have no effect\n"); }
    5065  }
    5166
    … …  
    140155    if (j + sy >= Ny) continue;
    141156    float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
     157    float *dV = (var) ? ((float *)(var[0].matrix.buffer) + (j+sy)*var[0].matrix.Naxis[0] + sx) : NULL;
    142158    for (int i = 0; i < dX; i++) {
    143159      if (i + sx < 0) continue;
    … …  
    145161      if (*V > SatThreshold) goto next; // skip pixels above threshold
    146162      if (!isfinite(*V)) goto next; // skip nan pixels
    147       dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
     163      if (dV) {
     164        dz[N] = *dV;
     165      } else {
     166        dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
     167      }
    148168      if (dz[N] <= 0) goto next;
    149169      dz[N] = 1.0 / dz[N];
    … …  
    154174    next:
    155175      V++;
     176      if (dV) { dV++; }
    156177    }
    157178  }

  • trunk/Ohana/src/opihi/cmd.astro/mkgauss.c

    r41341 r42821  
    1111  Buffer *buf;
    1212
     13  // if TRUE, integrated flux is Flux, else Io is Flux
    1314  int Normalize = FALSE;
    1415  if ((N = get_argument (argc, argv, "-norm"))) {
    … …  
    1617    remove_argument (N, &argc, argv);
    1718  }   
     19
     20  // Io or integrated flux (depending on value of Normalize)
     21  float Flux = 1.0;
     22  if ((N = get_argument (argc, argv, "-flux"))) {
     23    remove_argument (N, &argc, argv);
     24    Flux = atof(argv[N]);
     25    remove_argument (N, &argc, argv);
     26  }   
     27
    1828
    1929  // this should be Nx/2, Ny/2 if not set
    … …  
    2939  if ((argc < 3) || (argc > 5)) {
    3040    gprint (GP_ERR, "USAGE: mkgauss (buffer) (sigma) [[sy/sx] angle]\n");
     41    gprint (GP_ERR, " -flux flux : integral or peak is flux (default is 1.0)\n");
     42    gprint (GP_ERR, " -norm : integral is flux (else peak)\n");
     43    gprint (GP_ERR, " -c X Y : place center at X,Y\n");
    3144    return (FALSE);
    3245  }
    … …  
    6679  /* f = exp (-r), r = (x^2 / 2Sx) + (y^2 / 2Sy) + Sxy*x*y */
    6780
    68   double Io = Normalize ? 1.0 / (2.0 * M_PI * Sig_x * Sig_y) : 1.0;
     81  double Io = Normalize ? Flux / (2.0 * M_PI * Sig_x * Sig_y) : Flux;
    6982
    7083  in = (float *) buf[0].matrix.buffer;

  • trunk/Ohana/src/opihi/cmd.astro/region.c

    r42078 r42821  
    160160    if (!strcasecmp (argv[CtypeArg], "PAR")) { strcpy (graphmode.coords.ctype, "DEC--PAR"); goto got_ctype; }
    161161    if (!strcasecmp (argv[CtypeArg], "MOL")) { strcpy (graphmode.coords.ctype, "DEC--MOL"); goto got_ctype; }
     162    if (!strcasecmp (argv[CtypeArg], "LIN")) { strcpy (graphmode.coords.ctype, "DEC--LIN"); goto got_ctype; }
     163    if (!strcasecmp (argv[CtypeArg], "CAR")) { strcpy (graphmode.coords.ctype, "DEC--CAT"); goto got_ctype; }
    162164    gprint (GP_ERR, "ERROR: invalid projection type %s\n", argv[CtypeArg]);
    163     gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR, MOL\n");
     165    gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR, MOL, LIN, CAR\n");
    164166    return FALSE;
    165167  }

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r42456 r42821  
    6767$(SRC)/ungridify.$(ARCH).o     \
    6868$(SRC)/histogram.$(ARCH).o      \
     69$(SRC)/histbins.$(ARCH).o       \
    6970$(SRC)/tdhistogram.$(ARCH).o    \
    7071$(SRC)/hermitian1d.$(ARCH).o    \
    … …  
    8081$(SRC)/imresample.$(ARCH).o     \
    8182$(SRC)/imcollapse.$(ARCH).o     \
     83$(SRC)/impoints.$(ARCH).o       \
    8284$(SRC)/integrate.$(ARCH).o      \
    8385$(SRC)/interpolate.$(ARCH).o    \

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r42456 r42821  
    5656int ungridify        PROTO((int, char **));
    5757int histogram        PROTO((int, char **));
     58int histbins         PROTO((int, char **));
    5859int tdhistogram      PROTO((int, char **));
    5960int hermitian1d      PROTO((int, char **));
    … …  
    7071int imresample       PROTO((int, char **));
    7172int imcollapse       PROTO((int, char **));
     73int impoints         PROTO((int, char **));
    7274int integrate        PROTO((int, char **));
    7375int interpolate      PROTO((int, char **));
    … …  
    265267  {1, "header",       header,           "print image header"},
    266268  {1, "histogram",    histogram,        "generate histogram from vector"},
     269  {1, "histbins",     histbins,         "generate histogram from vector, bins specified by vectors"},
    267270  {1, "tdhistogram",  tdhistogram,      "generate 2D histogram image from vector set"},
    268271  {1, "hermitian1d",  hermitian1d,      "generate 1-D Hermitian Polynomial"},
    … …  
    281284  {1, "imconvolve",   imconvolve,       "full 2D real-space convolution"},
    282285  {1, "imstats",      imstats,          "statistics on a portion of an image"},
     286  {1, "impoints",     impoints,         "insert points into an image by vector pair"},
    283287  {1, "integrate",    integrate,        "integrate a vector"},
    284288  {1, "interpolate_presort",  interpolate_presort,      "interpolate between vector pairs"},

  • trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c

    r41341 r42821  
    1717  /** check basic syntax **/
    1818  if (argc != 5) {
    19     gprint (GP_ERR, "USAGE: interpolate Xi Yi Xo Yo\n");
     19    gprint (GP_ERR, "USAGE: interpolate Xi Yi Xo Yo [-fill-ends]\n");
    2020    gprint (GP_ERR, "  Xi Yi - sorted reference vectors\n");
    2121    gprint (GP_ERR, "  Xo    - output positions (vector)\n");
    … …  
    2323    gprint (GP_ERR, "  (vectors must be pre-sorted)\n");
    2424    gprint (GP_ERR, "  (use 'threshold' to interpolate to a value)\n");
     25    gprint (GP_ERR, "  -fill-ends : values outside of range are set to end-point values\n");
    2526    return (FALSE);
    2627  }

  • trunk/Ohana/src/opihi/cmd.data/vgroup.c

    r33963 r42821  
    7070  ALLOCATE (values, double, xin[0].Nelements);
    7171
    72   for (i = 0; i < xout[0].Nelements - 1; i++) {
     72  // if we specify binsize, then we need to examine all bins
     73  int Nout = isnan(binsize) ? xout[0].Nelements - 1 : xout[0].Nelements;
     74  for (i = 0; i < Nout; i++) {
    7375    if (isnan(binsize)) {
    7476      xmin = xout[0].elements.Flt[i];
    … …  
    8183    N = 0;
    8284    for (j = 0; j < xin[0].Nelements; j++) {
    83       if (xin[0].elements.Flt[j] < xmin) continue;
    84       if (xin[0].elements.Flt[j] > xmax) continue;
     85      if (xin[0].elements.Flt[j] <  xmin) continue;
     86      if (xin[0].elements.Flt[j] >= xmax) continue;
    8587      if (yin) {
    8688        values[N] = yin[0].elements.Flt[j];

  • trunk/Ohana/src/opihi/dvo/avperiodogram.c

    r40574 r42821  
    274274        // The first 3 values[] / fields[] will be time, mag, dmag
    275275        for (off_t n = 0; n < Nmfields; n++) {
    276           mvalues[n] = dbExtractMeasures (average, secfilt, &measure[k], NULL, NULL, &mfields[n]);
     276          mvalues[n] = dbExtractMeasures (average, secfilt, NULL, &measure[k], NULL, NULL, NULL, &mfields[n]);
    277277        }
    278278        if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;

  • trunk/Ohana/src/opihi/dvo/avperiodomatch.c

    r41379 r42821  
    302302        // The first 3 values[] / fields[] will be time, mag, dmag
    303303        for (off_t n = 0; n < Nmfields; n++) {
    304           mvalues[n] = dbExtractMeasures (average, secfilt, &measure[k], NULL, NULL, &mfields[n]);
     304          mvalues[n] = dbExtractMeasures (average, secfilt, NULL, &measure[k], NULL, NULL, NULL, &mfields[n]);
    305305        }
    306306        if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;

  • trunk/Ohana/src/opihi/dvo/avselect.c

    r41341 r42821  
    1717  int VERBOSE;
    1818  char name[1024];
    19   float RADIUS;
     19  opihi_flt RADIUS;
    2020
    2121  Catalog catalog;
    2222
    2323  Vector **vec, **invec, *RAvec, *DECvec, *IDXvec, *RADvec;
     24  Vector *RINvec = NULL;
    2425  dbField *fields;
    2526  dbValue *values;
    … …  
    102103  }
    103104
    104   RADIUS = atof (argv[1]);
    105 
    106   /* load regions which contain all supplied RA,DEC coordinates */
     105  if (SelectScalar (argv[1], &RADIUS)) {
     106    remove_argument (1, &argc, argv);
     107  } else {
     108    gprint (GP_ERR, " RADIUS must be a numerical value\n");
     109    goto help;
     110    if ((RINvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
     111    RADIUS = 0.0; // find the max radius for region selection:
     112    for (i = 0; i < RINvec->Nelements; i++) {
     113      RADIUS = MAX(RINvec->elements.Flt[i], RADIUS);
     114    }
     115  }
     116
     117  /* load regions which contain all supplied RA,DEC coordinates (for RINvec, uses max radius) */
    107118  if ((skylist = SelectRegionsByCoordVectorsAndRadius (RAvec, DECvec, RADIUS/3600.0)) == NULL) goto escape;
    108119
    … …  
    167178    return status;
    168179  }
    169 
    170   RADIUS = atof (argv[1]);
    171   remove_argument (1, &argc, argv);
    172180
    173181  // parse the fields to be extracted and returned

  • trunk/Ohana/src/opihi/dvo/gstar.c

    r41427 r42821  
    401401          gprint (GP_LOG, "%20s ",     date);
    402402          gprint (GP_LOG, "%f   ",     catalog.average[k].Trange / 86400.0);
     403
     404          gprint (GP_LOG, "%ld ", catalog.average[k].extID);
     405          // gprint (GP_LOG, "0x%16x ", catalog.average[k].extID);
    403406        }
    404407

  • trunk/Ohana/src/opihi/dvo/mextract.c

    r42103 r42821  
    221221  }
    222222
     223  int needLensobj = dbFieldNeedLensobj (fields, Nfields);
    223224  int needLensing = dbFieldNeedLensing (fields, Nfields);
    224225  int needStarpar = dbFieldNeedStarpar (fields, Nfields, FALSE);
     226  int needGalphot = dbFieldNeedGalphot (fields, Nfields);
    225227
    226228  // the lensing table does not have a good index to/from the measure table.  if we need lensing
    … …  
    247249    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT;
    248250    catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
     251    catalog.catflags |= needLensobj ? DVO_LOAD_LENSOBJ : DVO_SKIP_LENSOBJ;
    249252    catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
     253    catalog.catflags |= needGalphot ? DVO_LOAD_GALPHOT : DVO_SKIP_GALPHOT;
    250254    catalog.Nsecfilt = Nsecfilt;
    251255
    … …  
    297301        dbExtractMeasuresInitMeas (); // reset counters for saved fields  (costs very little
    298302
    299         int Nstarpar = average->starparOffset;
    300         StarPar *starpar = needStarpar && average->Nstarpar ? &catalog.starpar[Nstarpar] : NULL;
     303        int nStar = average->starparOffset;
     304        StarPar *starpar = needStarpar && average->Nstarpar ? &catalog.starpar[nStar] : NULL;
     305
     306        int nLens = average->lensobjOffset;
     307        Lensobj *lensobj = needLensobj && average->Nlensobj ? &catalog.lensobj[nLens] : NULL;
     308
     309        int nPhot = average->galphotOffset;
     310        GalPhot *galphot = needGalphot && average->Ngalphot ? &catalog.galphot[nPhot] : NULL;
    301311
    302312        Lensing *lensing = NULL;
    … …  
    313323
    314324        for (n = 0; n < Nfields; n++) {
    315           values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], lensing, starpar, &fields[n]);
     325          values[n] = dbExtractMeasures (average, secfilt, lensobj, &catalog.measure[m], lensing, starpar, galphot, &fields[n]);
    316326        }
    317327        // fprintf (stderr, "object: ave: %f, cat: %f, averef %d\n", fields[n].name, values[2], values[3], catalog.measure[m].averef);

  • trunk/Ohana/src/opihi/dvo/mmatch.c

    r41341 r42821  
    232232  ALLOCATE (index, off_t, Nelem);
    233233
    234   // int needLensing = dbFieldNeedLensing (fields, Nfields);
     234  int needLensobj = dbFieldNeedLensobj (fields, Nfields);
     235//int needLensing = dbFieldNeedLensing (fields, Nfields);
    235236  int needStarpar = dbFieldNeedStarpar (fields, Nfields, FALSE);
     237  int needGalphot = dbFieldNeedGalphot (fields, Nfields);
    236238
    237239  // grab data from all selected sky regions
    … …  
    249251    catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
    250252    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT | DVO_LOAD_MEASURE;
     253    catalog.catflags |= needLensobj ? DVO_LOAD_LENSOBJ : DVO_SKIP_LENSOBJ;
    251254    catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
     255    catalog.catflags |= needGalphot ? DVO_LOAD_GALPHOT : DVO_SKIP_GALPHOT;
    252256    catalog.Nsecfilt = Nsecfilt;
    253257
    … …  
    289293        Average *average = &catalog.average[Ncat];
    290294
    291         int Nstarpar = average->starparOffset;
    292         StarPar *starpar = needStarpar ? &catalog.starpar[Nstarpar] : NULL;
     295        int nStar = average->starparOffset;
     296        StarPar *starpar = needStarpar && average->Nstarpar ? &catalog.starpar[nStar] : NULL;
     297
     298        int nLens = average->lensobjOffset;
     299        Lensobj *lensobj = needLensobj && average->Nlensobj ? &catalog.lensobj[nLens] : NULL;
     300
     301        int nPhot = average->galphotOffset;
     302        GalPhot *galphot = needGalphot && average->Ngalphot ? &catalog.galphot[nPhot] : NULL;
    293303
    294304        // int Nlensing = average->lensobjOffset;
    … …  
    299309
    300310        for (n = 0; n < Nfields; n++) {
    301           values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], NULL, starpar, &fields[n]);
     311          values[n] = dbExtractMeasures (average, secfilt, lensobj, &catalog.measure[m], NULL, starpar, galphot, &fields[n]);
    302312        }
    303313

  • trunk/Ohana/src/opihi/dvo/mmextract.c

    r39457 r42821  
    174174  }
    175175
    176   // int needLensing = dbFieldNeedLensing (fields, Nfields);
     176  int needLensobj = dbFieldNeedLensobj (fields, Nfields);
     177//int needLensing = dbFieldNeedLensing (fields, Nfields);
    177178  int needStarpar = dbFieldNeedStarpar (fields, Nfields, FALSE);
     179  int needGalphot = dbFieldNeedGalphot (fields, Nfields);
    178180
    179181  // grab data from all selected sky regions
    … …  
    185187    catalog.filename = skylist[0].filename[i];
    186188    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT;
     189//  catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
     190    catalog.catflags |= needLensobj ? DVO_LOAD_LENSOBJ : DVO_SKIP_LENSOBJ;
     191    catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
     192    catalog.catflags |= needGalphot ? DVO_LOAD_GALPHOT : DVO_SKIP_GALPHOT;
    187193    catalog.Nsecfilt = Nsecfilt;
    188194
    … …  
    220226        Average *average = &catalog.average[j];
    221227
    222         int Nstarpar = average->starparOffset;
    223         StarPar *starpar = needStarpar ? &catalog.starpar[Nstarpar] : NULL;
     228        int nStar = average->starparOffset;
     229        StarPar *starpar = needStarpar && average->Nstarpar ? &catalog.starpar[nStar] : NULL;
     230
     231        int nLens = average->lensobjOffset;
     232        Lensobj *lensobj = needLensobj && average->Nlensobj ? &catalog.lensobj[nLens] : NULL;
     233
     234        int nPhot = average->galphotOffset;
     235        GalPhot *galphot = needGalphot && average->Ngalphot ? &catalog.galphot[nPhot] : NULL;
    224236
    225237        // int Nlensing = average->lensobjOffset;
    … …  
    231243        for (n = 0; n < Nfields; n++) {
    232244          // values needs to be a pointer to a type with FLT and INT (with a union, we would save a bit of memory...)
    233           values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], NULL, starpar, &fields[n]);
     245          values[n] = dbExtractMeasures (average, secfilt, lensobj, &catalog.measure[m], NULL, starpar, galphot, &fields[n]);
    234246        }
    235247        // fprintf (stderr, "object: ave: %f, cat: %f, averef %d\n", fields[n].name, values[2], values[3], catalog.measure[m].averef);

  • trunk/Ohana/src/opihi/lib.data/mrq2dmin.c

    r41715 r42821  
    1717static opihi_flt *parmax = NULL;
    1818
     19/** only used locally **/
    1920opihi_flt mrq2dcof (opihi_flt *x, opihi_flt *t, opihi_flt *y, opihi_flt *dy, int Npts,
    2021              opihi_flt *par, int Npar, opihi_flt **ta, opihi_flt **tb,
    … …  
    6061}
    6162
     63/** only used in imfit-test.c ?? **/
    6264opihi_flt mrq2dchi (opihi_flt *x, opihi_flt *t, opihi_flt *y, opihi_flt *dy, int Npts,
    6365                opihi_flt *par, int Npar,
    … …  
    126128    chisq = ochisq;
    127129  }
     130
     131  /* note that the parameter errors are sqrt(alpha[j][j]) */
    128132
    129133  return (chisq);

  • trunk/Ohana/src/opihi/lib.data/spline.c

    r42332 r42821  
    113113  opihi_flt dx, a, b, value;
    114114 
    115   // saturate correction at high and low ends
    116   if (X < x[0]) return y[0];
    117   if (X > x[N-1]) return y[N-1];
     115  // linear extrapolation past endpoints
     116  if (X < x[0]) {
     117    hi = 1;
     118    lo = 0;
     119    goto evaluate;
     120    // alternative: saturate correction at high and low ends
     121    // return y[0];
     122  }
     123  if (X > x[N-1]) {
     124    hi = N - 1;
     125    lo = N - 2;
     126    goto evaluate;
     127    // alternative: saturate correction at high and low ends
     128    // return y[N-1];
     129  }
    118130
    119131  /* find correct element in array (x must be sorted) */
    … …  
    128140    }
    129141  }
     142
     143evaluate:
    130144
    131145  /* error condition: duplicate abssisca */

  • trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c

    r42080 r42821  
    7777    if (!strcmp (argv[i], "isflt"))  { type = ST_UNARY; goto gotit; }
    7878    if (!strcmp (argv[i], "length")) { type = ST_UNARY; goto gotit; }
     79    if (!strcmp (argv[i], "toupper")) { type = ST_UNARY; goto gotit; }
     80    if (!strcmp (argv[i], "tolower")) { type = ST_UNARY; goto gotit; }
    7981
    8082    /* binary operations */

  • trunk/Ohana/src/opihi/lib.shell/stack_math.c

    r42389 r42821  
    16151615  int i, Nx;
    16161616 
     1617  // handle string vectors in the L_unary function
     1618  if (V1->vector->type == OPIHI_STR) {
     1619    int status = L_unary (OUT, V1, op);
     1620    return status;
     1621  }
     1622
    16171623  Nx = V1[0].vector[0].Nelements;
    16181624
    16191625  OUT[0].vector = InitVector ();
    16201626  OUT[0].type = ST_VECTOR_TMP; /*** <<--- says this is a temporary matrix ***/
    1621 
    1622   if (V1->vector->type == OPIHI_STR) {
    1623     ResetVector (OUT->vector, V1->vector->type, V1->vector->Nelements);
    1624     for (i = 0; i < V1->vector->Nelements; i++) {
    1625       OUT->vector->elements.Str[i] = strcreate (V1->vector->elements.Str[i]);
    1626     }
    1627     goto escape;
    1628   }
    16291627
    16301628# define V_FUNC(OP,FTYPE) {                                             \
    … …  
    17441742    for (int i = 0; i < Nx; i++) {                     
    17451743      Ov[i] = strcreate (Iv[i]);
     1744    }                                                                   
     1745    goto escape;                                                       
     1746  }
     1747  if (!strcmp (op, "toupper")) {
     1748    for (int i = 0; i < Nx; i++) {                     
     1749      Ov[i] = strcreate (Iv[i]);
     1750      for (int j = 0; Ov[i][j]; j++) {
     1751        Ov[i][j] = toupper ((unsigned int) Ov[i][j]);
     1752      }
     1753    }                                                                   
     1754    goto escape;                                                       
     1755  }
     1756  if (!strcmp (op, "tolower")) {
     1757    for (int i = 0; i < Nx; i++) {                     
     1758      Ov[i] = strcreate (Iv[i]);
     1759      for (int j = 0; Ov[i][j]; j++) {
     1760        Ov[i][j] = tolower ((unsigned int) Ov[i][j]);
     1761      }
    17461762    }                                                                   
    17471763    goto escape;                                                       

  • trunk/Ohana/src/relphot

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

    r42788 r42821  
    417417int IS_DIFF_DB;
    418418
     419double MAG_MIN;
    419420double MAG_LIM;
    420421double SIGMA_LIM;
     422
     423double STAR_MAX_PSF_KRON;
     424double PSF_QF_MIN;
     425
    421426double IMAGE_SCATTER;
    422427double IMAGE_OFFSET;
    … …  
    471476int    RESET_FLATCORR;
    472477int    REPAIR_WARPS;
     478int    USE_REF_EQUIV;
    473479int    PRESERVE_PS1;
    474480int    REQUIRE_PSFFIT;
    … …  
    486492double IMAGE_GOOD_FRACTION;
    487493int    CALIBRATE_STACKS_AND_WARPS;
     494int    CALIBRATE_WARPS_LIKE_CHIPS;
     495int    CALIBRATE_STACKS_LIKE_CHIPS;
    488496
    489497int    KEEP_UBERCAL;

  • trunk/Ohana/src/relphot/src/ConfigInit.c

    r41647 r42821  
    2323  if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
    2424
    25   GetConfig (config, "MAG_LIM",                "%lf", 0, &MAG_LIM);
    26   GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
     25  // the following cuts are limit the objects used to calculate the image zero points.
     26  // These are used in bcatalog.c and apply to Mcat (nominal measurement magnitude)
     27  GetConfig (config, "SIGMA_LIM",              "%lf",    0, &SIGMA_LIM);
     28  GetConfig (config, "MAG_LIM",                "%lf",    0, &MAG_LIM);
     29
     30  DefConfig ("MAG_MIN",                        "%lf", 15.0, MAG_MIN);
     31  DefConfig ("STAR_MAX_PSF_KRON",              "%lf",  0.2, STAR_MAX_PSF_KRON);
     32  DefConfig ("PSF_QF_MIN",                     "%lf", 0.95, PSF_QF_MIN);
    2733
    2834  DefConfig ("STAR_SCATTER",                    "%lf",  0.1, STAR_SCATTER);

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

    r41664 r42821  
    531531
    532532  // until the analysis has converged a bit, do not use the IRLS analysis
    533   // default is MaxIterations = 10
    534   if (UseStandardOLS(ZPT_IMAGES)) {
     533  // default is MaxIterations = 10.
     534
     535  // However, for stack and warp calibration, we are not iteratively solving for the
     536  // zero points.  In this case, use IRLS in a single pass
     537  if (!CALIBRATE_STACKS_AND_WARPS && UseStandardOLS(ZPT_IMAGES)) {
    535538    brightStars.MaxIterations = 0;
    536539    kronStars.MaxIterations = 0;
    … …  
    672675
    673676      if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
    674         fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars.alldata->yVector[Nkron]);
     677        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f %6.3f %6.3f | %5.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, MsysKron, MrelKron, Mmos, Mgrid, Mflat, Moff, psfStars.alldata->dyVector[Nref]);
    675678      }
    676679
    … …  
    702705    if (mark) continue;
    703706
    704     // use liststats to find the 20-pct, median, 80-pct points
    705     StatType stats;
    706     liststats_setmode (&stats, "MEDIAN");
    707     liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &stats);
    708     double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
    709 
    710     // soften the individual errors with 10% of the scatter above
    711     for (j = 0; j < Nref; j++) {
    712       double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma);
    713       psfStars.alldata->dyVector[j] = newSigma;
    714     }
    715 
    716     // soften the errors based on the scatter
     707    // soften the errors based on the core scatter
    717708    FitDataSetSoften (&psfStars, Nref);
    718709

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

    r41647 r42821  
    3434} 
    3535
    36 //
     36// Nsecfilt is a property of the photcode table.  we only need to retrieve it once.
     37static int Nsecfilt = -1;
     38
     39// return the average photometry for this measurement
    3740float getMrel (Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source) {
    3841
    … …  
    4144
    4245  int ecode = GetPhotcodeEquivCodebyCode (photcode);
     46  if (USE_REF_EQUIV) {
     47    // tie this photometry to an alternative refernce
     48    // (e.g., tie UDR4.i to PV3.i)
     49    int refcode = GetPhotcodeEquivCodebyCode (ecode);
     50    if (refcode < 0) return NAN;
     51    ecode = refcode;
     52  }
     53
    4354  int Nsec = GetPhotcodeNsec(ecode);
    44   int Nsecfilt = GetPhotcodeNsecfilt ();
     55  if (Nsec < 0) return NAN;
     56
     57  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    4558
    4659  int entry = Nsecfilt*ave+Nsec;
    … …  
    312325  }
    313326
    314   int Nsecfilt = GetPhotcodeNsecfilt ();
     327  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    315328
    316329  SetMrelInfo summary, results;
    … …  
    334347  int i;
    335348
    336   int Nsecfilt = GetPhotcodeNsecfilt ();
     349  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    337350
    338351  SetMrelInfo summary, results;
    … …  
    424437  ThreadInfo *threadinfo = data;
    425438
    426   int Nsecfilt = GetPhotcodeNsecfilt ();
     439  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    427440
    428441  SetMrelInfo results;
    … …  
    451464  off_t k;
    452465
    453   int Nsecfilt = GetPhotcodeNsecfilt ();
     466  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    454467
    455468  off_t m = average[0].measureOffset;
    … …  
    559572  off_t j;
    560573
    561   int Nsecfilt = GetPhotcodeNsecfilt ();
     574  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    562575
    563576  for (i = 0; i < Ncatalog; i++) {
    … …  
    601614  ALLOCATE (slist, double, Ntot);
    602615
    603   int Nsecfilt = GetPhotcodeNsecfilt ();
     616  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    604617
    605618  // eliminate bad stars using the stats for a single secfilt at a time
    … …  
    666679  // N1 = N2 = N3 = N4 = N0 = 0;
    667680
    668   int Nsecfilt = GetPhotcodeNsecfilt ();
     681  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    669682
    670683  Ntot = 0;
    … …  
    719732  StatType stats;
    720733
    721   int Nsecfilt = GetPhotcodeNsecfilt ();
     734  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    722735
    723736  Ntot = 0;
    … …  
    760773  StatType stats;
    761774
    762   int Nsecfilt = GetPhotcodeNsecfilt ();
     775  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    763776
    764777  Ntot = 0;
    … …  
    805818  ALLOCATE (Mlist, double, NBIN);
    806819
    807   int Nsecfilt = GetPhotcodeNsecfilt ();
     820  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    808821
    809822  int Ns;
    … …  
    840853  Graphdata graphdata;
    841854
    842   int Nsecfilt = GetPhotcodeNsecfilt ();
     855  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    843856
    844857  Ntotal = 0;
    … …  
    911924  return;
    912925
    913   int Nsecfilt = GetPhotcodeNsecfilt ();
     926  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
    914927
    915928  FILE *fout = NULL;

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

    r42389 r42821  
    237237  }
    238238
     239  USE_REF_EQUIV = FALSE;
     240  if ((N = get_argument (argc, argv, "-use-ref-equiv"))) {
     241    remove_argument (N, &argc, argv);
     242    USE_REF_EQUIV = TRUE;
     243  }
    239244  PRESERVE_PS1 = FALSE;
    240245  if ((N = get_argument (argc, argv, "-preserve-ps1"))) {
    … …  
    521526    remove_argument (N, &argc, argv);
    522527    USE_REFERENCE_WEIGHT = TRUE;
     528  }
     529
     530  /* Normally, stacks and warps are treated differently from chips.  In this case,
     531     WHAT IS DIFFERENT?
     532   */
     533
     534  CALIBRATE_WARPS_LIKE_CHIPS = FALSE;
     535  if ((N = get_argument (argc, argv, "-calibrate-warps-like-chips"))) {
     536    remove_argument (N, &argc, argv);
     537    CALIBRATE_WARPS_LIKE_CHIPS = TRUE;
     538  }
     539  CALIBRATE_STACKS_LIKE_CHIPS = FALSE;
     540  if ((N = get_argument (argc, argv, "-calibrate-stacks-like-chips"))) {
     541    remove_argument (N, &argc, argv);
     542    CALIBRATE_STACKS_LIKE_CHIPS = TRUE;
    523543  }
    524544
    … …  
    836856  }
    837857
     858  USE_REF_EQUIV = FALSE;
     859  if ((N = get_argument (argc, argv, "-use-ref-equiv"))) {
     860    remove_argument (N, &argc, argv);
     861    USE_REF_EQUIV = TRUE;
     862  }
    838863  PRESERVE_PS1 = FALSE;
    839864  if ((N = get_argument (argc, argv, "-preserve-ps1"))) {
    … …  
    9781003  }
    9791004
     1005  CALIBRATE_WARPS_LIKE_CHIPS = FALSE;
     1006  if ((N = get_argument (argc, argv, "-calibrate-warps-like-chips"))) {
     1007    remove_argument (N, &argc, argv);
     1008    CALIBRATE_WARPS_LIKE_CHIPS = TRUE;
     1009  }
     1010  CALIBRATE_STACKS_LIKE_CHIPS = FALSE;
     1011  if ((N = get_argument (argc, argv, "-calibrate-stacks-like-chips"))) {
     1012    remove_argument (N, &argc, argv);
     1013    CALIBRATE_STACKS_LIKE_CHIPS = TRUE;
     1014  }
     1015
    9801016  STAGES = STAGE_CHIP | STAGE_WARP | STAGE_STACK;
    9811017  if ((N = get_argument (argc, argv, "-skip-chip"))) {

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

    r42389 r42821  
    11# include "relphot.h"
     2
     3# define BAD_COUNT(FIELD) \
     4  badCount.FIELD ++; \
     5  if (doTest1 && (catalog[0].measure[offset].imageID == TEST_IMAGE1)) { badTest1.FIELD ++; } \
     6  if (doTest2 && (catalog[0].measure[offset].imageID == TEST_IMAGE2)) { badTest2.FIELD ++; }
     7
     8typedef struct {
     9  int Ncode;
     10  int Ntime;
     11  int Ndophot;
     12  int Nmag;
     13  int Nsigma;
     14  int Nimag;
     15  int Nfew;
     16  int Ngalaxy;
     17  int Npsfqf;
     18  int Nnan;
     19  int Nbad;
     20  int Npoor;
     21} MyBadCountType;
    222
    323int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
    … …  
    929  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    1030  float mag;
    11   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
     31
     32  // int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
    1233
    1334  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    2647  int NmMin = catalog[0].Nmeasure;
    2748
    28   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
     49  MyBadCountType badCount = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     50  MyBadCountType badTest1 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     51  MyBadCountType badTest2 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     52
     53  int doTest1 = (TEST_IMAGE1 >= 0);
     54  int doTest2 = (TEST_IMAGE2 >= 0);
     55
     56  FILE *ftest1 = (doTest1) ? fopen ("test1.meas.dat", "a") : NULL;
     57  FILE *ftest2 = (doTest2) ? fopen ("test2.meas.dat", "a") : NULL;
     58
     59  // Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
    2960
    3061  // copy the following fields to the subcatalog:
    … …  
    5788
    5889      offset = catalog[0].average[i].measureOffset + j;
     90
     91      if (doTest1) {
     92        Measure *m = &catalog[0].measure[offset];
     93        if (m->imageID == TEST_IMAGE1) {
     94          fprintf (ftest1, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
     95        }
     96      }
     97      if (doTest2) {
     98        Measure *m = &catalog[0].measure[offset];
     99        if (m->imageID == TEST_IMAGE2) {
     100          fprintf (ftest2, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
     101        }
     102      }
    59103
    60104      /* select measurements by photcode */
    … …  
    69113      }
    70114      if (!found) {
    71         Ncode ++;
     115        BAD_COUNT(Ncode);
    72116        continue;
    73117      }
    … …  
    75119      /* select measurements by time */
    76120      if (TimeSelect) {
    77         if (catalog[0].measure[offset].t < TSTART) { Ntime ++; continue; }
    78         if (catalog[0].measure[offset].t > TSTOP)  { Ntime ++; continue; }
     121        if (catalog[0].measure[offset].t < TSTART) { BAD_COUNT(Ntime); continue; }
     122        if (catalog[0].measure[offset].t > TSTOP)  { BAD_COUNT(Ntime); continue; }
    79123      }
    80124
    81125      /* select measurements by quality */
    82       if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
     126      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { BAD_COUNT(Ndophot); continue; }
    83127
    84128      float Maper = USE_APER_FOR_STARGAL ? catalog[0].measure[offset].Map : catalog[0].measure[offset].Mkron;
    85129
    86130      // skip garbage measurements
    87       if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.95))     { Npsfqf ++; continue; }
    88       if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.95)) { Npsfqf ++; continue; }
    89       if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
    90       if (isnan(Maper)) { Nnan ++; continue; }
     131      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
     132      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
     133      if (isnan(catalog[0].measure[offset].M)) { BAD_COUNT(Nnan); continue; }
     134      if (isnan(Maper)) { BAD_COUNT(Nnan); continue; }
    91135
    92136      // require 0x01 in photFlags (fitted with a PSF) -- add to the test data
    93       if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
     137      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { BAD_COUNT(Nbad); continue; }
    94138
    95139      // very loose cut on PSF - APER mag (Map or Mkron) -- this lightly rejects galaxies & other oddities
    96140      float Mkp = catalog[0].measure[offset].M - Maper;
    97       if (fabs(Mkp) > 1.0) { Nbad ++; continue; }
     141      if (fabs(Mkp) > 1.0) { BAD_COUNT(Nbad); continue; }
    98142
    99143      if (catalog[0].measure[offset].photFlags & code->photomBadMask) {
    100         Nbad++;
     144        BAD_COUNT(Nbad);
    101145        continue;
    102146      }
    103       if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
     147      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { BAD_COUNT(Npoor); continue;}
    104148
    105149      // weak test for galaxies
    106       if (Mkp > 0.5) {
     150      if (Mkp > STAR_MAX_PSF_KRON) {
    107151        nEXT ++;
    108152      } else {
    … …  
    112156      /* select measurements by mag limit */
    113157      mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF);
    114       if (mag > MAG_LIM) { Nmag ++; continue; }
     158      if (mag > MAG_LIM) { BAD_COUNT(Nmag); continue; }
     159      if (mag < MAG_MIN) { BAD_COUNT(Nmag); continue; }
     160      // note that these two cuts are applied to all photcodes. these are also applied to
     161      // the nominal apparent mags : if the expected sensitivity is very different between
     162      // different filters, or if exposure times have a large range, these cuts may not be
     163      // effective (or may exclude too much).
    115164
    116165      /* select measurements by measurement error */
    117       if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { Nsigma ++; continue; }
     166      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { BAD_COUNT(Nsigma); continue; }
    118167
    119168      /* select measurements by mag limit */
    120169      if (ImagSelect) {
    121170        mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF);
    122         if (mag < ImagMin) { Nimag ++; continue; }
    123         if (mag > ImagMax) { Nimag ++; continue; }
    124       }
     171        if (mag < ImagMin) { BAD_COUNT(Nimag); continue; }
     172        if (mag > ImagMax) { BAD_COUNT(Nimag); continue; }
     173      }
     174      // note that instrumental mag cuts are reasonable for individual exposures using the same
     175      // detector, but are probably not appropriate for stacks with a large range of the number
     176      // of inputs. 
    125177
    126178      // count this measurement as valid for this secfilt entry
    … …  
    147199    if (nEXT >= nPSF) {
    148200      Nmeasure -= Nm;
    149       Ngalaxy ++;
     201      badCount.Ngalaxy++;
    150202      continue;
    151203    }
    … …  
    157209    if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
    158210      Nmeasure -= Nm;
    159       Nfew ++;
     211      badCount.Nfew++;
    160212      continue;
    161213    }
    … …  
    167219    if (!anySecfiltGood) {
    168220      Nmeasure -= Nm;
    169       Nfew ++;
     221      badCount.Nfew++;
    170222      continue;
    171223    }
    … …  
    204256             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
    205257    fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
    206              Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
     258             badCount.Nfew, badCount.Ncode, badCount.Ntime, badCount.Ndophot, badCount.Nmag, badCount.Nsigma, badCount.Nimag, badCount.Npsfqf, badCount.Nnan, badCount.Ngalaxy, badCount.Nbad, badCount.Npoor);
     259  }
     260
     261  if (doTest1) {
     262    fprintf (stderr, "TEST1 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     263             badTest1.Nfew, badTest1.Ncode, badTest1.Ntime, badTest1.Ndophot, badTest1.Nmag, badTest1.Nsigma, badTest1.Nimag, badTest1.Npsfqf, badTest1.Nnan, badTest1.Ngalaxy, badTest1.Nbad, badTest1.Npoor);
     264  }
     265  if (doTest2) {
     266    fprintf (stderr, "TEST2 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     267             badTest2.Nfew, badTest2.Ncode, badTest2.Ntime, badTest2.Ndophot, badTest2.Nmag, badTest2.Nsigma, badTest2.Nimag, badTest2.Npsfqf, badTest2.Nnan, badTest2.Ngalaxy, badTest2.Nbad, badTest2.Npoor);
    207268  }
    208269
    … …  
    212273  }
    213274
    214   free (Nvalid)
     275  free (Nvalid);
     276
     277  if (ftest1) { fclose (ftest1); }
     278  if (ftest2) { fclose (ftest2); }
    215279
    216280  return (TRUE);

  • trunk/Ohana/src/relphot/src/extra.c

    r41647 r42821  
    7070}
    7171
    72 // for now (20140710) I need to identify gpc1 stacks explicitly.  generalize in the future
     72// for now (20140710) I need to identify gpc1 & gpc2 stacks explicitly.  generalize in the future
    7373int isGPC1stack (int photcode) {
    7474
    75   if ((photcode == 11000) || (photcode == 14100)) return TRUE; // g-band
    76   if ((photcode == 11100) || (photcode == 14200)) return TRUE; // r-band
    77   if ((photcode == 11200) || (photcode == 14300)) return TRUE; // i-band
    78   if ((photcode == 11300) || (photcode == 14400)) return TRUE; // z-band
    79   if ((photcode == 11400) || (photcode == 14500)) return TRUE; // y-band
    80   if ((photcode == 11500) || (photcode == 14600)) return TRUE; // w-band
     75  // 11000 : GPC1.g.SkyChip
     76  // 11001 : VARIOUS.g.SkyChip (output of xcamera stacks)
     77  // 14100 : SIMTEST.g.SkyChip
     78  // 31000 : GPC2.g.SkyChip
     79 
     80  if (CALIBRATE_STACKS_LIKE_CHIPS) return FALSE;
     81
     82  if ((photcode == 11000) || (photcode == 11001) || (photcode == 14100) || (photcode == 31000)) return TRUE; // g-band
     83  if ((photcode == 11100) || (photcode == 11101) || (photcode == 14200) || (photcode == 31100)) return TRUE; // r-band
     84  if ((photcode == 11200) || (photcode == 11201) || (photcode == 14300) || (photcode == 31200)) return TRUE; // i-band
     85  if ((photcode == 11300) || (photcode == 11301) || (photcode == 14400) || (photcode == 31300)) return TRUE; // z-band
     86  if ((photcode == 11400) || (photcode == 11401) || (photcode == 14500) || (photcode == 31400)) return TRUE; // y-band
     87  if ((photcode == 11500) || (photcode == 11501) || (photcode == 14600) || (photcode == 31500)) return TRUE; // w-band
    8188
    8289  return FALSE;
    … …  
    8592// for now (20140710) I need to identify gpc1 stacks explicitly.  generalize in the future
    8693int isGPC1warp (int photcode) {
     94
     95  if (CALIBRATE_WARPS_LIKE_CHIPS) return FALSE;
     96
     97  // 11000 : GPC1.g.SkyChip
     98  // 14100 : SIMTEST.g.SkyChip
     99
     100  // 12000 : GPC1.g.ForcedWarp
     101  // 15100 : SIMTEST.g.ForcedWarp
     102  // 32000 : GPC2.g.ForcedWarp
    87103
    88104  // diff warps get stack-like photcodes (kind of lame)
    … …  
    97113  }
    98114
    99   if ((photcode == 12000) || (photcode == 15100)) return TRUE; // g-band
    100   if ((photcode == 12100) || (photcode == 15200)) return TRUE; // r-band
    101   if ((photcode == 12200) || (photcode == 15300)) return TRUE; // i-band
    102   if ((photcode == 12300) || (photcode == 15400)) return TRUE; // z-band
    103   if ((photcode == 12400) || (photcode == 15500)) return TRUE; // y-band
    104   if ((photcode == 12500) || (photcode == 15600)) return TRUE; // w-band
     115  if ((photcode == 12000) || (photcode == 15100) || (photcode == 32000)) return TRUE; // g-band
     116  if ((photcode == 12100) || (photcode == 15200) || (photcode == 32100)) return TRUE; // r-band
     117  if ((photcode == 12200) || (photcode == 15300) || (photcode == 32200)) return TRUE; // i-band
     118  if ((photcode == 12300) || (photcode == 15400) || (photcode == 32300)) return TRUE; // z-band
     119  if ((photcode == 12400) || (photcode == 15500) || (photcode == 32400)) return TRUE; // y-band
     120  if ((photcode == 12500) || (photcode == 15600) || (photcode == 32500)) return TRUE; // w-band
    105121
    106122  return FALSE;

  • trunk/Ohana/src/relphot/src/fit1d_irls.c

    r41557 r42821  
    374374  StatType stats;
    375375  liststats_setmode (&stats, "MEDIAN");
    376   liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats);
     376  liststats (dataset->alldata->yVector, dataset->alldata->dyVector, NULL, Nvalues, &stats);
     377  // liststats sorts yVector and dyVector together
    377378
    378379  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
    … …  
    380381  // soften the individual errors with 10% of the scatter above
    381382  for (int j = 0; j < Nvalues; j++) {
    382     double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma);
     383    double newSigma = hypot(dataset->alldata->dyVector[j], 0.5*altSigma);
    383384    dataset->alldata->dyVector[j] = newSigma;
    384385  }

  • trunk/Ohana/src/relphot/src/initialize.c

    r42134 r42821  
    5050      fprintf (stderr, "DophotSelect: FALSE\n");
    5151    }
    52     fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n");
     52    fprintf (stderr, "PSF_QF limit: %f\n", PSF_QF_MIN);
    5353
    5454    // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n");
    5555
    56     fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM);
     56    fprintf (stderr, "MAG_MIN: %f, MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_MIN, MAG_LIM, SIGMA_LIM);
    5757    fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax);
    5858
     59    fprintf (stderr, "STAR_MAX_PSF_KRON: %f\n", STAR_MAX_PSF_KRON);
    5960    fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW);
    6061

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

    r41647 r42821  
    151151    // TimeSelect -time
    152152    // DophotSelect
    153     // (note that psfQF is applied rigidly at 0.85, as is the galaxy test)
     153    // PSF_QF_MIN
     154    // STAR_MAX_PSF_KRON
     155    // MAG_MIN
    154156    // MAG_LIM
    155157    // SIGMA_LIM
    … …  
    158160
    159161    char *command = NULL;
    160     strextend (&command, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D CAMERA %s -D MAG_LIM %f -D SIGMA_LIM %f",
    161               PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA, MAG_LIM, SIGMA_LIM);
     162    strextend (&command, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s", PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname);
     163    strextend (&command, "-region %f %f %f %f -D CAMERA %s", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA);
     164    strextend (&command, "-D MAG_MIN %f -D MAG_LIM %f -D SIGMA_LIM %f", MAG_MIN, MAG_LIM, SIGMA_LIM);
     165    strextend (&command, "-D STAR_MAX_PSF_KRON %f -D PSF_QF_MIN %f", STAR_MAX_PSF_KRON, PSF_QF_MIN);
    162166
    163167    if (VERBOSE)             { strextend (&command, "-v"); }

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r42112 r42821  
    251251    if (!(STAGES & STAGE_WARP))  { strextend (&command, "-skip-warp"); }
    252252    if (!(STAGES & STAGE_STACK)) { strextend (&command, "-skip-stack"); }
     253    if (CALIBRATE_WARPS_LIKE_CHIPS) { strextend (&command, "-calibrate-warps-like-chips"); }
     254    if (CALIBRATE_STACKS_LIKE_CHIPS) { strextend (&command, "-calibrate-stacks-like-chips"); }
    253255
    254256    // deprecate

  • trunk/Ohana/src/relphot/src/relphot.c

    r42483 r42821  
    3333 
    3434  switch (mode) {
    35     case UPDATE_IMAGES:
     35    case UPDATE_IMAGES: // -images
    3636      // calculate zero points for images (may group by exposure [mosaic], night [tgroups]; may measure flat-correction)
    37       // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat
     37      // IF CALLED WITH NLOOP == 0 (and NOT -only-stacks-and-warps), this step DOES NOT LOAD bcatalog.  Instead it just loads
     38      // the image table and applies the Mcal values to the measurements.  In parallel context, it generates ImageSubset.dat.
    3839      // If called with -update, calls reload_catalogs() just like apply_offsets, UNLESS -only-stacks-and-warps is selected
    3940      relphot_images (skylist);

  • trunk/Ohana/src/tools/src/mpcorb_predict.c

    r42389 r42821  
    5151PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc);
    5252
    53 Planets *mpcorb_read_text (char *filename, int *nplanets);
     53Planets *mpcorb_read_text (char *filename, int *nplanets, int noHeader);
    5454Planets *mpcorb_read_fits (char *filename, int *nplanets);
    5555void     mpcorb_save_fits (char *filename, Planets *planets, int Nplanets);
    … …  
    7373
    7474int main (int argc, char **argv) {
     75
     76  if (get_argument (argc, argv, "-h")) goto usage;
     77  if (get_argument (argc, argv, "--help")) goto usage;
    7578
    7679  // XXX test
    … …  
    8487    // load the planets
    8588    int Nplanets = 0;
    86     Planets *planets = mpcorb_read_text (filename, &Nplanets);
     89    Planets *planets = mpcorb_read_text (filename, &Nplanets, FALSE);
    8790    for (int i = 0; i < 10; i++) {
    8891      // if it is in the region, use the more accurate calculation
    … …  
    9295  }
    9396
    94   if ((argc != 9) && (argc != 10)) goto usage;
    95 
     97  if (argc < 2) goto usage;
    9698  if (!strcasecmp (argv[1], "trange")) mpcorb_trange (argc, argv);
    9799  if (!strcasecmp (argv[1], "moment")) mpcorb_moment (argc, argv);
    … …  
    99101usage:
    100102  fprintf (stderr, "USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
    101   fprintf (stderr, "USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n", argv[0]);
     103  fprintf (stderr, "USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n\n", argv[0]);
    102104
    103105  fprintf (stderr, " trange : generate a table of asteroid orbital elements which fall within the region at Tmin or Tmax (in MJD)\n");
    104   fprintf (stderr, " moment : generate a table of asteroid positions which fall within the region at Tobs (in MJD)\n");
    105 
     106  fprintf (stderr, " moment : generate a table of asteroid positions which fall within the region at Tobs (in MJD)\n\n");
     107
     108  fprintf (stderr, " Tmin, Tmax, Tobs are in MJD; Rmin, Rmax, Dmin, Dmax are in decimal degrees\n");
     109  fprintf (stderr, " Orbital predictions are valid for the PS1 site\n");
     110
     111  fprintf (stderr, " trange options:\n");
     112  fprintf (stderr, " -no-header (for files other than MPCORB.DAT)\n");
     113  fprintf (stderr, " -min-obs (Nobs) : minimum number of observations to keep object\n");
     114  fprintf (stderr, " -min-opp (Nopp) : minimum number of oppositions to keep object\n");
    106115  exit (2);
    107116}
    … …  
    110119void mpcorb_trange (int argc, char **argv) {
    111120
    112   if (argc != 10) Shutdown ("USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
     121  int N;
     122 
     123  int noHeader = FALSE;
     124  if ((N = get_argument (argc, argv, "-no-header"))) {
     125    remove_argument (N, &argc, argv);
     126    noHeader = TRUE;
     127  }
     128  int MIN_OBS = 0;
     129  if ((N = get_argument (argc, argv, "-min-obs"))) {
     130    remove_argument (N, &argc, argv);
     131    MIN_OBS = atoi(argv[N]);
     132    remove_argument (N, &argc, argv);
     133  }
     134  int MIN_OPP = 0;
     135  if ((N = get_argument (argc, argv, "-min-opp"))) {
     136    remove_argument (N, &argc, argv);
     137    MIN_OPP = atoi(argv[N]);
     138    remove_argument (N, &argc, argv);
     139  }
     140
     141  if (argc != 10) Shutdown ("USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits) [-no-header] [-min-obs] [-min-opp]\n", argv[0]);
    113142
    114143  char  *filename = argv[2];
    … …  
    122151
    123152  int Nplanets = 0;
    124   Planets *planets = mpcorb_read_text (filename, &Nplanets);
     153  Planets *planets = mpcorb_read_text (filename, &Nplanets, noHeader);
    125154  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
    126155
    … …  
    133162    if (i % 1000 == 0) fprintf (stderr, ".");
    134163
    135     // XXX make these quality cuts user options
    136     if (SKIP_BAD && ((planets[i].Nobs < 20) || (planets[i].Nopp < 2))) continue;
     164    // skip poor-quality observations
     165    if (planets[i].Nobs < MIN_OBS) continue;
     166    if (planets[i].Nopp < MIN_OPP) continue;
    137167
    138168    // is the object in the region for the start of the period?
    … …  
    804834
    805835// read the full MPCORB file, storing desired information in Planets structure
    806 Planets *mpcorb_read_text (char *filename, int *nplanets) {
     836Planets *mpcorb_read_text (char *filename, int *nplanets, int noHeader) {
    807837
    808838  ALLOCATE_PTR (buffer, char, NBYTE);
    … …  
    838868
    839869      Nline ++;
    840       if ((Nline > 43) && (q - p > 5)) {
     870      int doParse = (noHeader || (Nline > 43));
     871      if (doParse && (q - p > 5)) {
    841872        mpcorb_parseline (p, &planets[Nplanets], Nbyte - offset);
     873        if ((Nplanets < 10) && DEBUG) fprintf (stderr, "object: %s", planets[Nplanets].ID);
    842874        Nplanets ++;
    843875        CHECK_REALLOCATE (planets, Planets, NPLANETS, Nplanets, 10000);
Note: See TracChangeset for help on using the changeset viewer.