IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2015, 5:50:46 AM (11 years ago)
Author:
eugene
Message:

merged changes from trunk

Location:
branches/eam_branches/ipp-20151113
Files:
89 edited
9 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20151113

  • branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/loadgalphot_readstars.c

    r39125 r39266  
    204204  fit->ClipNsigma = 5.0;
    205205
     206  // use the list of index values from above to join entries with the same index
     207
    206208  int Nbad = 0;
    207209  for (i = 0; i < Nrow; i++) {

  • branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/mkcmf.c

    r39121 r39266  
    1313# define FLAGS 0x1101
    1414
    15 void gauss_init (int Nbin);
    16 double rnd_gauss (double mean, double sigma);
    1715void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1816void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    … …  
    282280  }
    283281   
    284   gauss_init (2048);
     282  ohana_gaussdev_init ();
    285283
    286284  // load test stars from a file:
    … …  
    564562
    565563
    566 static int Ngaussint = 0;
    567 static double *gaussint;
    568 
    569 extern double drand48();
    570 
    571 double gaussian (double x, double mean, double sigma) {
    572 
    573   double f;
    574 
    575   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    576 
    577   return (f);
    578 
    579 }
    580 
    581 /* integrate a gaussian from -5 sigma to +5 sigma */
    582 void gauss_init (int Nbin) {
    583  
    584   int i;
    585   double val, x, dx, dx1, dx2, dx3, df;
    586   double mean, sigma;
    587  
    588   /* no need to generate this if it already exists */
    589   if (Ngaussint == Nbin) return;
    590 
    591   // A = time(NULL);
    592   // // XXX this is expensive if called a lot (1 sec min)
    593   // // for (B = 0; A == time(NULL); B++);
    594   // B = A + 10000;
    595   // srand48(B);
    596  
    597   Ngaussint = Nbin;
    598   ALLOCATE (gaussint, double, Ngaussint + 1);
    599 
    600   val = 0;
    601   dx = 1.0 / Ngaussint;
    602   dx1 = dx / 3.0;
    603   dx2 = 2.0*dx/3.0;
    604   dx3 = dx;
    605   mean = 0.0;
    606   sigma = 1.0;
    607  
    608   for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
    609     df = (3.0*gaussian(x    , mean, sigma) +
    610           9.0*gaussian(x+dx1, mean, sigma) +
    611           9.0*gaussian(x+dx2, mean, sigma) +
    612           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    613     val += df;
    614     if (val > (i + 0.5) / (double) Ngaussint) {
    615       gaussint[i] = x + dx / 2.0;
    616       i++;
    617     }
    618   }
    619 }
    620 
    621 double rnd_gauss (double mean, double sigma) {
    622  
    623   int i;
    624   double y;
    625  
    626   y = drand48();
    627   i = Ngaussint*y;
    628   y = gaussint[i]*sigma + mean;
    629  
    630   return (y);
    631  
    632 }
    633  
    634 double int_gauss (int i) {
    635   double y;
    636   y = gaussint[i];
    637   return (y);
    638 }
    639  
    640564void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {
    641565
    … …  
    654578
    655579    if (ADDNOISE) {
    656       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    657       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    658       M[i] += fSN*rnd_gauss(0.0, 1.0);
     580      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     581      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     582      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    659583      flux = pow (10.0, -0.4*M[i]);
    660584      fSN = 1.0 / sqrt(flux);
    … …  
    700624
    701625    if (ADDNOISE) {
    702       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    703       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    704       M[i] += fSN*rnd_gauss(0.0, 1.0);
     626      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     627      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     628      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    705629      flux = pow (10.0, -0.4*M[i]);
    706630      fSN = 1.0 / sqrt(flux);
    … …  
    749673
    750674    if (ADDNOISE) {
    751       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    752       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    753       M[i] += fSN*rnd_gauss(0.0, 1.0);
     675      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     676      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     677      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    754678      flux = pow (10.0, -0.4*M[i]);
    755679      fSN = 1.0 / sqrt(flux);
    … …  
    800724
    801725    if (ADDNOISE) {
    802       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    803       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    804       M[i] += fSN*rnd_gauss(0.0, 1.0);
     726      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     727      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     728      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    805729      flux = pow (10.0, -0.4*M[i]);
    806730      fSN = 1.0 / sqrt(flux);
    … …  
    857781
    858782    if (ADDNOISE) {
    859       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    860       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    861       M[i] += fSN*rnd_gauss(0.0, 1.0);
     783      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     784      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     785      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    862786      flux = pow (10.0, -0.4*M[i]);
    863787      fSN = 1.0 / sqrt(flux);
    … …  
    918842
    919843    if (ADDNOISE) {
    920       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    921       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    922       M[i] += fSN*rnd_gauss(0.0, 1.0);
     844      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     845      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     846      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    923847      flux = pow (10.0, -0.4*M[i]);
    924848      fSN = 1.0 / sqrt(flux);
    … …  
    987911
    988912    if (ADDNOISE) {
    989       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    990       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    991       M[i] += fSN*rnd_gauss(0.0, 1.0);
     913      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     914      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     915      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    992916      flux = pow (10.0, -0.4*M[i]);
    993917      fSN = 1.0 / sqrt(flux);
    … …  
    10941018
    10951019    if (ADDNOISE) {
    1096       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1097       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1098       M[i] += fSN*rnd_gauss(0.0, 1.0);
     1020      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1021      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1022      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    10991023      flux = pow (10.0, -0.4*M[i]);
    11001024      fSN = 1.0 / sqrt(flux);
    … …  
    11981122
    11991123    if (ADDNOISE) {
    1200       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1201       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1202       M[i] += fSN * rnd_gauss(0.0, 1.0);
     1124      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1125      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1126      M[i] += fSN * ohana_gaussdev_rnd(0.0, 1.0);
    12031127      flux = pow(10.0, -0.4 * M[i]);
    12041128      fSN = 1.0 / sqrt(flux);
    … …  
    13241248
    13251249    if (ADDNOISE) {
    1326       X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
    1327       Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
    1328       M[i] += fSN*rnd_gauss(0.0, 1.0);
     1250      X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1251      Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0);
     1252      M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0);
    13291253      flux = pow (10.0, -0.4*M[i]);
    13301254      fSN = 1.0 / sqrt(flux);

  • branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r39220 r39266  
    8585    N1[Nstars] = i;
    8686    Nstars ++;
    87     input[0].found_t[i] = -1;
     87    input[0].found_t[i] = FALSE;
    8888  }
    8989  if (Nstars < 1) {
    … …  
    180180    off_t N = N1[i];
    181181
     182    input[0].found_t[N] = TRUE;
     183
    182184    /* make sure there is space for next Nmeasure entries */
    183185    if (Nmeas + input[0].average[N].Nmeasure >= NMEAS) {
    … …  
    209211      Nreplace = replace_tycho (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[Minp]);
    210212      if (Nreplace == 6) {
    211         input[0].found_t[N] = Nmeas;
    212213        output[0].found_t[n] = Nmeas;
    213214        i++;
    … …  
    224225        // XXX this does not support lensing, starpar, or galphot measurements
    225226        if (replace_match (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[offset])) {
    226           input[0].found_t[N] = output[0].average[n].measureOffset; 
    227227          continue;
    228228        }
    … …  
    274274        }
    275275      }
    276       input[0].found_t[N] = Nmeas;
    277276      output[0].average[n].Nmeasure ++;
    278277      Nmeas ++;
    … …  
    394393    }
    395394
    396     if (input[0].found_t[N] >= 0) continue;
     395    if (input[0].found_t[N]) continue;
    397396
    398397    // if we are using MATCHED_TABLES, we are going to leave the edge cases in their
    399398    // source catalog, even if they have leaked beyond the edge
    400399    if (!MATCHED_TABLES && !IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue;
     400
     401    input[0].found_t[N] = TRUE;
    401402
    402403    // XXX should we accept the input measurements for these fields?
    … …  
    450451
    451452      /* we set next[Nmeas] to -1 here, and update correctly below */
    452       input[0].found_t[N] = Nmeas;
    453453      next_meas[Nmeas] = -1;
    454454      Nmeas ++;
    … …  
    559559  }
    560560
    561   /* note stars which have been found in this catalog */
    562   for (i = 0; i < input[0].Naverage; i++) {
    563     if (input[0].found_t[i] > -1) {
    564       input[0].found_t[i] = -2;
    565     } else {
    566       input[0].found_t[i] = -3;
    567     }
    568   }
    569 
    570561  /* check if the catalog has changed?  if no change, no need to write */
    571562  output[0].objID    = objID; // new max value, save on catalog close

  • branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c

    r39126 r39266  
    229229
    230230    float sersic_index = galphot->index;
    231     if (galphot->modelType == 6) sersic_index = 1.0;
    232     if (galphot->modelType == 7) sersic_index = 4.0;
     231    if (galphot->modelType == 6) {
     232      sersic_index = 1.0;
     233    }
     234    else if (galphot->modelType == 7) {
     235      sersic_index = 4.0;
     236    }
     237    else {
     238      sersic_index = 1.0 / (2.0 * sersic_index);
     239    }
    233240
    234241    PRINT_FLOAT(gal_buffer, galphot->mag,              "%.6e, ");

  • branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FWobjects_dvopsps_catalog.c

    r39092 r39266  
    148148                 "FLUX_PSF_WRP_ERR    FLOAT,  "         
    149149                 "FLUX_PSF_WRP_STDEV  FLOAT,  "         
    150 
     150                 "MAG_PSF_WRP         FLOAT,  "
     151                 "MAG_PSF_WRP_ERR     FLOAT,  "
     152                 
    151153                 "FLUX_KRON_WRP       FLOAT,  "       
    152154                 "FLUX_KRON_WRP_ERR   FLOAT,  "   
    153155                 "FLUX_KRON_WRP_STDEV FLOAT,  "   
    154 
     156                 "MAG_KRON_WRP        FLOAT,  "
     157                 "MAG_KRON_WRP_ERR    FLOAT,  "
     158                 
    155159                 "FLUX_AP_WRP         FLOAT,  "             
    156160                 "FLUX_AP_WRP_ERR     FLOAT,  "             
    157                  "FLUX_AP_WRP_STDEV   FLOAT,  "             
     161                 "FLUX_AP_WRP_STDEV   FLOAT,  "
     162                 "MAG_AP_WRP          FLOAT,  "
     163                 "MAG_AP_WRP_ERR      FLOAT,  "
    158164
    159165                 "FLAGS               INT "     
    … …  
    188194                 "FLUX_STD_AP_R5     FLOAT, "
    189195                 "FLUX_FIL_AP_R5     FLOAT, "
     196                 "MAG_AP_R5          FLOAT, "
     197                 "MAG_ERR_AP_R5      FLOAT, "
    190198
    191199                 "FLUX_AP_R6         FLOAT, "
    … …  
    193201                 "FLUX_STD_AP_R6     FLOAT, "
    194202                 "FLUX_FIL_AP_R6     FLOAT, "
     203                 "MAG_AP_R6          FLOAT, "
     204                 "MAG_ERR_AP_R6      FLOAT, "
    195205
    196206                 "FLUX_AP_R7         FLOAT, "
    … …  
    198208                 "FLUX_STD_AP_R7     FLOAT, "
    199209                 "FLUX_FIL_AP_R7     FLOAT, "
     210                 "MAG_AP_R7          FLOAT, "
     211                 "MAG_ERR_AP_R7      FLOAT, "
    200212
    201213                 "X11_SM_OBJ         FLOAT, "
    … …  
    262274                 "FLUX_PSF_WRP, "             
    263275                 "FLUX_PSF_WRP_ERR, "         
    264                  "FLUX_PSF_WRP_STDEV, "         
     276                 "FLUX_PSF_WRP_STDEV, "
     277                 "MAG_PSF_WRP,  "
     278                 "MAG_PSF_WRP_ERR, "
    265279
    266280                 "FLUX_KRON_WRP, "       
    267281                 "FLUX_KRON_WRP_ERR, "   
    268                  "FLUX_KRON_WRP_STDEV, "   
     282                 "FLUX_KRON_WRP_STDEV, "
     283                 "MAG_KRON_WRP, "
     284                 "MAG_KRON_WRP_ERR, "
    269285
    270286                 "FLUX_AP_WRP, "             
    271287                 "FLUX_AP_WRP_ERR, "             
    272                  "FLUX_AP_WRP_STDEV, "             
     288                 "FLUX_AP_WRP_STDEV, "
     289                 "MAG_AP_WRP, "
     290                 "MAG_AP_WRP_ERR, "
    273291
    274292                 "FLAGS"     
    … …  
    287305                 "FLUX_STD_AP_R5, "
    288306                 "FLUX_FIL_AP_R5, "
     307                 "MAG_AP_R5, "
     308                 "MAG_ERR_AP_R5, "
    289309
    290310                 "FLUX_AP_R6, "
    … …  
    292312                 "FLUX_STD_AP_R6, "
    293313                 "FLUX_FIL_AP_R6, "
     314                 "MAG_AP_R6, "
     315                 "MAG_ERR_AP_R6, "
    294316
    295317                 "FLUX_AP_R7, "
    … …  
    297319                 "FLUX_STD_AP_R7, "
    298320                 "FLUX_FIL_AP_R7, "
     321                 "MAG_AP_R7, "
     322                 "MAG_ERR_AP_R7, "
    299323
    300324                 "X11_SM_OBJ, "
    … …  
    348372    float meanPSFFluxErr = NAN;
    349373    float meanPSFFluxStd = NAN;
     374    float meanPSFMag     = NAN;
     375    float meanPSFMagErr  = NAN;
    350376    if (isfinite(secfilt->dFpsfWrp) && isfinite(secfilt->FpsfWrp)) {
    351377      meanPSFFlux    = secfilt-> FpsfWrp;
    352378      meanPSFFluxErr = secfilt->dFpsfWrp;
    353379      meanPSFFluxStd = secfilt->sFpsfWrp;
     380      if (meanPSFFlux > 0.0) {
     381        meanPSFMag    = -2.5 * log10(meanPSFFlux / 3631.0);
     382        meanPSFMagErr = (2.5 * meanPSFFluxErr) / (meanPSFFlux * log(10));
     383      }
    354384    }
    355385
    … …  
    357387    float meanKronFluxErr = NAN;
    358388    float meanKronFluxStd = NAN;
     389    float meanKronMag     = NAN;
     390    float meanKronMagErr  = NAN;
    359391    if (isfinite(secfilt->dFkronWrp) && isfinite(secfilt->FkronWrp)) {
    360392      meanKronFlux    = secfilt-> FkronWrp;
    361393      meanKronFluxErr = secfilt->dFkronWrp;
    362394      meanKronFluxStd = secfilt->sFkronWrp;
     395      if (meanKronFlux > 0.0) {
     396        meanKronMag    = -2.5 * log10(meanKronFlux / 3631.0);
     397        meanKronMagErr = (2.5 * meanKronFluxErr) / (meanKronFlux * log(10));
     398      }     
    363399    }
    364400
    … …  
    366402    float meanApFluxErr = NAN;
    367403    float meanApFluxStd = NAN;
     404    float meanApMag     = NAN;
     405    float meanApMagErr  = NAN;
    368406    if (isfinite(secfilt->dFapWrp) && isfinite(secfilt->FapWrp)) {
    369407      meanApFlux    = secfilt-> FapWrp;
    370408      meanApFluxErr = secfilt->dFapWrp;
    371409      meanApFluxStd = secfilt->sFapWrp;
     410      if (meanApFlux > 0.0) {
     411        meanApMag    = -2.5 * log10(meanApFlux / 3631.0);
     412        meanApMagErr = (2.5 * meanApFluxErr) / (meanApFlux * log(10));
     413      }
    372414    }
    373415
    … …  
    382424    PRINT_FLOAT(sec_buffer, meanPSFFluxErr,   "%.6e, "); // uflux precision
    383425    PRINT_FLOAT(sec_buffer, meanPSFFluxStd,   "%.6e, "); // uflux precision
    384 
     426    PRINT_FLOAT(sec_buffer, meanPSFMag,       "%.6f, ");
     427    PRINT_FLOAT(sec_buffer, meanPSFMagErr,    "%.6f, ");   
     428   
    385429    PRINT_FLOAT(sec_buffer, meanKronFlux,     "%.6e, "); // uflux precision
    386430    PRINT_FLOAT(sec_buffer, meanKronFluxErr,  "%.6e, "); // uflux precision
    387431    PRINT_FLOAT(sec_buffer, meanKronFluxStd,  "%.6e, "); // uflux precision
     432    PRINT_FLOAT(sec_buffer, meanKronMag,      "%.6f, ");
     433    PRINT_FLOAT(sec_buffer, meanKronMagErr,   "%.6f, ");
    388434
    389435    PRINT_FLOAT(sec_buffer, meanApFlux,     "%.6e, "); // uflux precision
    390436    PRINT_FLOAT(sec_buffer, meanApFluxErr,  "%.6e, "); // uflux precision
    391437    PRINT_FLOAT(sec_buffer, meanApFluxStd,  "%.6e, "); // uflux precision
     438    PRINT_FLOAT(sec_buffer, meanApMag,      "%.6f, ");
     439    PRINT_FLOAT(sec_buffer, meanApMagErr,   "%.6f, ");
    392440
    393441    PrintIOBuffer (sec_buffer, "%u ", secfilt->flags);       
    … …  
    402450    PrintIOBuffer (cpy_buffer, "%hd, ", lensobj->photcode);         
    403451
     452    float magApR5    = NAN;
     453    float magErrApR5 = NAN;
     454    float magApR6    = NAN;
     455    float magErrApR6 = NAN;
     456    float magApR7    = NAN;
     457    float magErrApR7 = NAN;
     458
     459    if (lensobj->F_ApR5 > 0.0) {
     460      magApR5    = -2.5 * log10(lensobj->F_ApR5 / 3631.0);
     461      magErrApR5 = (2.5 * lensobj->dF_ApR5) / (lensobj->F_ApR5 * log(10));
     462    }
     463    if (lensobj->F_ApR6 > 0.0) {
     464      magApR6    = -2.5 * log10(lensobj->F_ApR6 / 3631.0);
     465      magErrApR6 = (2.5 * lensobj->dF_ApR6) / (lensobj->F_ApR6 * log(10));
     466    }
     467    if (lensobj->F_ApR7 > 0.0) {
     468      magApR7    = -2.5 * log10(lensobj->F_ApR7 / 3631.0);
     469      magErrApR7 = (2.5 * lensobj->dF_ApR7) / (lensobj->F_ApR7 * log(10));
     470    }
     471   
    404472    PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR5,     "%.6e, ");
    405473    PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR5,     "%.6e, ");
    406474    PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR5,     "%.6e, ");
    407475    PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR5,     "%.6e, ");
     476    PRINT_FLOAT(cpy_buffer, magApR5,              "%.6f, ");
     477    PRINT_FLOAT(cpy_buffer, magErrApR5,           "%.6f, ");
    408478    PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR6,     "%.6e, ");
    409479    PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR6,     "%.6e, ");
    410480    PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR6,     "%.6e, ");
    411481    PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR6,     "%.6e, ");
     482    PRINT_FLOAT(cpy_buffer, magApR6,              "%.6f, ");
     483    PRINT_FLOAT(cpy_buffer, magErrApR6,           "%.6f, ");
    412484    PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR7,     "%.6e, ");
    413485    PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR7,     "%.6e, ");
    414486    PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR7,     "%.6e, ");
    415487    PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR7,     "%.6e, ");
     488    PRINT_FLOAT(cpy_buffer, magApR7,              "%.6f, ");
     489    PRINT_FLOAT(cpy_buffer, magErrApR7,           "%.6f, ");
    416490    PRINT_FLOAT(cpy_buffer, lensobj->X11_sm_obj,  "%.6e, ");
    417491    PRINT_FLOAT(cpy_buffer, lensobj->X12_sm_obj,  "%.6e, ");

  • branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_diffobj_dvopsps_catalog.c

    r39092 r39266  
    9494                 "CAT_ID        INT,      "
    9595                 "RA_MEAN       DOUBLE,   "         
    96                  "DEC_MEAN      DOUBLE,   "         
     96                 "DEC_MEAN      DOUBLE,   "
    9797                 "RA_ERR        FLOAT,    "     
    98                  "DEC_ERR       FLOAT,    "     
     98                 "DEC_ERR       FLOAT,    "
     99                 "IAU_NAME      VARCHAR(32), "
     100                 "PSO_NAME      VARCHAR(32), "
    99101                 "FLAGS         INT       "       
    100102                 ")\n", basename);
    … …  
    167169                 "DEC_MEAN,      "         
    168170                 "RA_ERR,        "     
    169                  "DEC_ERR,       "     
     171                 "DEC_ERR,       "
     172                 "IAU_NAME,      "
     173                 "PSO_NAME,      "
    170174                 "FLAGS          "       
    171175                 ") VALUES \n", basename);
    … …  
    218222  PRINT_FLOAT(ave_buffer, average->dD,       "%.5f, ");  // 0.010 mas precision
    219223
     224  // Add names.
     225  int ra_hr,ra_min,ra_sec,ra_fracsec;
     226  int dec_deg,dec_min,dec_sec,dec_fracsec;
     227  double tmp_ra = average->R;
     228  double tmp_dec = average->D;
     229  double dec_sign = tmp_dec / fabs(tmp_dec);
     230  tmp_dec = fabs(tmp_dec);
     231
     232  tmp_ra /= 15.0;
     233  ra_hr   = (int) floor(tmp_ra);
     234
     235  tmp_ra -= ra_hr;
     236  tmp_ra *= 60.0;
     237  ra_min  = (int) floor(tmp_ra);
     238 
     239  tmp_ra -= ra_min;
     240  tmp_ra *= 60.0;
     241  ra_sec  = (int) floor(tmp_ra);
     242
     243  tmp_ra -= ra_sec;
     244  tmp_ra *= 1000.0;
     245  ra_fracsec = (int) floor(tmp_ra);
     246
     247  dec_deg = (int) floor(tmp_dec);
     248
     249  tmp_dec -= dec_deg;
     250  tmp_dec *= 60.0;
     251  dec_min = (int) floor(tmp_dec);
     252
     253  tmp_dec -= dec_min;
     254  tmp_dec *= 60.0;
     255  dec_sec = (int) floor(tmp_dec);
     256
     257  tmp_dec -= dec_sec;
     258  tmp_dec *= 1000.0;
     259  dec_fracsec = (int) floor(tmp_dec);
     260
     261  dec_deg *= (int) dec_sign;
     262
     263  // IAU NAME
     264  PrintIOBuffer(ave_buffer, "'PSO J%02d%02d%02d.%03d%+02d%02d%02d.%03d',  ",
     265                ra_hr,ra_min,ra_sec,ra_fracsec,
     266                dec_deg,dec_min,dec_sec,dec_fracsec);
     267
     268  // PSO NAME
     269  tmp_ra = average->R;
     270  tmp_dec = average->D;
     271  PrintIOBuffer(ave_buffer, "'PSO J%.5f%+.5f',  ",
     272                tmp_ra,tmp_dec);
     273
     274  // End names
     275 
    220276  PrintIOBuffer (ave_buffer, "%u ",  average->flags);         
    221277  PrintIOBuffer (ave_buffer, "),\n");

  • branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c

    r39098 r39266  
    101101                 "EPOCH_MEAN    DOUBLE,   "         
    102102                 "RA_ERR        FLOAT,    "     
    103                  "DEC_ERR       FLOAT,    "     
     103                 "DEC_ERR       FLOAT,    "
     104                 "IAU_NAME      VARCHAR(32), "
     105                 "PSO_NAME      VARCHAR(32), "
    104106                 "CHISQ_POS     FLOAT,    "   
    105107                 "CHISQ_PM      FLOAT,    "   
    … …  
    208210                 "EPOCH_MEAN,    "         
    209211                 "RA_ERR,        "     
    210                  "DEC_ERR,       "     
     212                 "DEC_ERR,       "
     213                 "IAU_NAME,      "
     214                 "PSO_NAME,      "
    211215                 "CHISQ_POS,     "   
    212216                 "CHISQ_PM,      "   
    … …  
    276280  PRINT_FLOAT(ave_buffer, average->dD,       "%.5f, ");  // 0.010 mas precision
    277281
     282  // Add names.
     283  int ra_hr,ra_min,ra_sec,ra_fracsec;
     284  int dec_deg,dec_min,dec_sec,dec_fracsec;
     285  double tmp_ra,tmp_dec;
     286
     287  if (secfilt->NstackDet != 0) {
     288    tmp_ra = average->Rstk;
     289    tmp_dec = average->Dstk;
     290  }
     291  else {
     292    tmp_ra = average->R;
     293    tmp_dec = average->D;
     294  }
     295 
     296  double dec_sign = tmp_dec / fabs(tmp_dec);
     297  tmp_dec = fabs(tmp_dec);
     298
     299
     300  tmp_ra /= 15.0;
     301  ra_hr   = (int) floor(tmp_ra);
     302
     303  tmp_ra -= ra_hr;
     304  tmp_ra *= 60.0;
     305  ra_min  = (int) floor(tmp_ra);
     306 
     307  tmp_ra -= ra_min;
     308  tmp_ra *= 60.0;
     309  ra_sec  = (int) floor(tmp_ra);
     310
     311  tmp_ra -= ra_sec;
     312  tmp_ra *= 1000.0;
     313  ra_fracsec = (int) floor(tmp_ra);
     314
     315  dec_deg = (int) floor(tmp_dec);
     316
     317  tmp_dec -= dec_deg;
     318  tmp_dec *= 60.0;
     319  dec_min = (int) floor(tmp_dec);
     320
     321  tmp_dec -= dec_min;
     322  tmp_dec *= 60.0;
     323  dec_sec = (int) floor(tmp_dec);
     324
     325  tmp_dec -= dec_sec;
     326  tmp_dec *= 1000.0;
     327  dec_fracsec = (int) floor(tmp_dec);
     328
     329  dec_deg *= (int) dec_sign;
     330
     331  // IAU NAME
     332  PrintIOBuffer(ave_buffer, "'PSO J%02d%02d%02d.%03d%+02d%02d%02d.%03d',  ",
     333                ra_hr,ra_min,ra_sec,ra_fracsec,
     334                dec_deg,dec_min,dec_sec,dec_fracsec);
     335
     336  // PSO NAME
     337  if (secfilt->NstackDet != 0) {
     338    tmp_ra = average->Rstk;
     339    tmp_dec = average->Dstk;
     340  }
     341  else {
     342    tmp_ra = average->R;
     343    tmp_dec = average->D;
     344  }
     345
     346  PrintIOBuffer(ave_buffer, "'PSO J%.5f%+.5f',  ",
     347                tmp_ra,tmp_dec);
     348
     349  // End names
     350 
    278351  PRINT_FLOAT(ave_buffer, average->ChiSqAve, "%.4f, ");
    279352  PRINT_FLOAT(ave_buffer, average->ChiSqPM,  "%.4f, ");

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/Makefile

    r37846 r39266  
    3838$(SRC)/fakestar_io.$(ARCH).o \
    3939$(SRC)/insert_fakestar.$(ARCH).o \
    40 $(SRC)/gaussian.$(ARCH).o \
    4140$(SRC)/fakeastro_images.$(ARCH).o \
    4241$(SRC)/fakeastro_images_region.$(ARCH).o \

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/fakeastro.c

    r37807 r39266  
    33int main (int argc, char **argv) {
    44
    5   gauss_init (50000);
     5  ohana_gaussdev_init ();
    66
    77  /* get configuration info, args */

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_2mass_measures.c

    r38986 r39266  
    8080
    8181    // uR,uD in linear arcsec
    82     double dRsee = rnd_gauss (0.0, 1.0 / SN);
    83     double dDsee = rnd_gauss (0.0, 1.0 / SN);
     82    double dRsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
     83    double dDsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
    8484
    8585    double dRoff = (dRpm + dRsee) / 3600.0;

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fake_stars_catalog.c

    r37807 r39266  
    113113
    114114    // uR,uD in linear arcsec
    115     double dRsee = rnd_gauss (0.0, 1.0 / SN);
    116     double dDsee = rnd_gauss (0.0, 1.0 / SN);
     115    double dRsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
     116    double dDsee = ohana_gaussdev_rnd (0.0, 1.0 / SN);
    117117
    118118    double dRoff = (dRpm + dRsee) / 3600.0;

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakeqsos.c

    r37807 r39266  
    5151
    5252    // Mr will be interpretted as m_r
    53     double Mr = rnd_gauss (18.0, 1.5);
     53    double Mr = ohana_gaussdev_rnd (18.0, 1.5);
    5454   
    5555    stars[i].R = R;

  • branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakestars.c

    r37807 r39266  
    4848      int inPatch = FALSE;
    4949      while (!inPatch) {
    50         z = rnd_gauss (0.0, FAKEASTRO_ZGAL);
     50        z = ohana_gaussdev_rnd (0.0, FAKEASTRO_ZGAL);
    5151        r = sqrt(drand48()) * FAKEASTRO_RGAL;
    5252        Lrad = drand48() * 2 * M_PI;
    … …  
    9898    double Mr;
    9999    if (bigPeak) {
    100       Mr = rnd_gauss (11.25, 1.0);
     100      Mr = ohana_gaussdev_rnd (11.25, 1.0);
    101101    } else {
    102       Mr = rnd_gauss (10.00, 2.5);
     102      Mr = ohana_gaussdev_rnd (10.00, 2.5);
    103103    }
    104104   

  • branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-sim.d

    r37807 r39266  
    2727FIELD XoffDCR,        X_OFF_DCR,     float,          X offset from correction,     pixels
    2828FIELD YoffDCR,        Y_OFF_DCR,     float,          Y offset from correction,     pixels
    29 FIELD RoffGAL,        R_OFF_GAL,     float,          RA offset from correction,    arcsec
    30 FIELD DoffGAL,        D_OFF_GAL,     float,          DEC offset from correction,   arcsec
     29
     30FIELD Mflat,          M_FLAT,        float,          Static Flat-field offset,     mag
     31FIELD dummy2,         PADDING,       int,            unused 4 bytes
    3132
    3233FIELD t,              TIME,         int,            time in seconds (UNIX)
    … …  
    4849FIELD dRsys,          POS_SYS_ERR,  short,          systematic error from astrom,   1/100 of pixels
    4950
    50 FIELD dummy,          PADDING,      short,          padding
     51FIELD dummy,          PADDING2,     short,          padding
    5152FIELD posangle,       POSANGLE,     short,          position angle sky to chip,     (0xffff/360) deg
    5253FIELD pltscale,       PLTSCALE,     float,          plate scale,                    arcsec/pixel

  • branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-v5.d

    r38986 r39266  
    4242FIELD XoffCAM,        X_OFF_CAM,     float,          X offset from correction,     pixels
    4343FIELD YoffCAM,        Y_OFF_CAM,     float,          Y offset from correction,     pixels
    44 FIELD RoffGAL,        R_OFF_GAL,     float,          RA offset from correction,    arcsec
    45 FIELD DoffGAL,        D_OFF_GAL,     float,          DEC offset from correction,   arcsec
     44
     45FIELD Mflat,          M_FLAT,        float,          Static Flat-field offset,     mag
     46FIELD dummy2,         PADDING,       int,            unused 4 bytes
    4647
    4748# could these be packed into fewer bits?

  • branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure.d

    r39205 r39266  
    4343FIELD YoffCAM,        Y_OFF_CAM,     float,          Y offset from correction,     pixels
    4444
    45 # XXX I can deprecate these as I am going to apply the correct uR,uD offset
    46 # XXX not sure how to use this yet...
    47 # XXX replace RoffGall with Moff -> static component of the zpt correction
    48 FIELD RoffGAL,        R_OFF_GAL,     float,          RA offset from correction,    arcsec
    49 FIELD DoffGAL,        D_OFF_GAL,     float,          DEC offset from correction,   arcsec
     45FIELD Mflat,          M_FLAT,        float,          Static Flat-field offset,     mag
     46FIELD dummy2,         PADDING,       int,            unused 4 bytes
    5047
    5148# could these be packed into fewer bits?

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/dvo.h

    r39125 r39266  
    476476  double         R;
    477477  double         D;
    478   float          RoffGAL;
    479   float          DoffGAL;
    480478  float          M;
    481479  float          Mcal;
     480  float          Mflat;
    482481  float          dM;
    483482  float          airmass;
    … …  
    620619  float            XoffDCR;              // X offset from correction (pixels)
    621620  float            YoffDCR;              // Y offset from correction (pixels)
    622   float            RoffGAL;              // RA offset from correction (arcsec)
    623   float            DoffGAL;              // DEC offset from correction (arcsec)
     621  float            Mflat;                // flat offset from correction (arcsec)
     622  int              padding2;             // dummy
    624623  float            Sky;                  // local estimate of sky flux (counts/sec)
    625624  float            dSky;                 // local estimate of sky flux (counts/sec)

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/libdvo_astro.h

    r38986 r39266  
    138138} Image;
    139139
     140CoordTransform *AllocTransform (double phi, double Xo, double xo);
    140141CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output);
    141142int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform);

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/coord_systems.c

    r37807 r39266  
    11# include "dvo.h"
    22   
     3CoordTransform *AllocTransform (double phi, double Xo, double xo) {
     4
     5  CoordTransform *transform;
     6  ALLOCATE (transform, CoordTransform, 1);
     7  transform->isIdentity = FALSE;
     8 
     9  transform->phi = phi*RAD_DEG;
     10  transform->Xo  = Xo*RAD_DEG;
     11  transform->xo  = xo;
     12
     13  // pre-calculated constants:
     14  transform->sin_phi_cos_Xo = sin(transform->phi)*cos(transform->Xo);
     15  transform->sin_phi_sin_Xo = sin(transform->phi)*sin(transform->Xo);
     16  transform->cos_phi        = cos(transform->phi);
     17 
     18  transform->cos_phi_cos_Xo = cos(transform->phi)*cos(transform->Xo);
     19  transform->cos_phi_sin_Xo = cos(transform->phi)*sin(transform->Xo);
     20  transform->sin_phi        = sin(transform->phi);
     21  transform->cos_Xo         = cos(transform->Xo);
     22  transform->sin_Xo         = sin(transform->Xo);
     23
     24  return transform;
     25}
     26
    327CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output) {
    428

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dbExtractMeasures.c

    r38986 r39266  
    597597    case MEAS_XOFF_CAM: { value.Flt = measure[0].XoffCAM; break; }
    598598    case MEAS_YOFF_CAM: { value.Flt = measure[0].YoffCAM; break; }
    599     case MEAS_ROFF_GAL: { value.Flt = measure[0].RoffGAL; break; }
    600     case MEAS_DOFF_GAL: { value.Flt = measure[0].DoffGAL; break; }
    601599
    602600    case MEAS_XFIELD: /* offset relative to exposure center in ra,dec space */
    … …  
    723721
    724722    case MEAS_MCAL_OFFSET: { value.Flt = measure[0].Mcal; break; }
    725 
    726     case MEAS_FLAT: /* OK */
    727       // flat = measure.Mcal - image.Mcal
    728       if (REMOTE_CLIENT) {
    729         ImageMetadata *image = MatchImageMetadataDVO (measure[0].imageID);
    730         if (image == NULL) break;
    731         value.Flt = measure[0].Mcal - image[0].Mcal;
    732       } else {
    733         Image *image = MatchImageDVO (measure[0].t, measure[0].photcode, measure[0].imageID);
    734         if (image == NULL) break;
    735         value.Flt = measure[0].Mcal - image[0].Mcal;
    736       }
    737       break;
     723    case MEAS_FLAT: { value.Flt = measure[0].Mflat; break; }
    738724
    739725      // we have measure[0].Xccd,Yccd and image[0].NX,NY.  Find the distance to the center

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_catalog.c

    r39125 r39266  
    332332 measure->XoffCAM   = NAN;
    333333 measure->YoffCAM   = NAN;
    334  measure->RoffGAL   = NAN;
    335  measure->DoffGAL   = NAN;
     334
     335 measure->Mflat     = 0.0;
     336 measure->dummy2    = 0;
    336337
    337338 measure->Sky       = NAN;

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c

    r38462 r39266  
    2929    out[i].XoffDCR    = in[i].XoffDCR;
    3030    out[i].YoffDCR    = in[i].YoffDCR;
    31     out[i].RoffGAL    = in[i].RoffGAL;
    32     out[i].DoffGAL    = in[i].DoffGAL;
     31    out[i].Mflat      = in[i].Mflat;
    3332    out[i].t          = in[i].t;
    3433    out[i].averef     = in[i].averef;
    … …  
    7574    out[i].XoffDCR    = in[i].XoffDCR;
    7675    out[i].YoffDCR    = in[i].YoffDCR;
    77     out[i].RoffGAL    = in[i].RoffGAL;
    78     out[i].DoffGAL    = in[i].DoffGAL;
     76    out[i].Mflat      = in[i].Mflat;
    7977    out[i].t          = in[i].t;
    8078    out[i].averef     = in[i].averef;

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c

    r39125 r39266  
    4242    out[i].XoffCAM    = in[i].XoffCAM;
    4343    out[i].YoffCAM    = in[i].YoffCAM;
    44     out[i].RoffGAL    = in[i].RoffGAL;
    45     out[i].DoffGAL    = in[i].DoffGAL;
     44    out[i].Mflat      = in[i].Mflat;
    4645    out[i].Sky        = in[i].Sky;
    4746    out[i].dSky       = in[i].dSky;
    … …  
    116115    out[i].XoffCAM    = in[i].XoffCAM;
    117116    out[i].YoffCAM    = in[i].YoffCAM;
    118     out[i].RoffGAL    = in[i].RoffGAL;
    119     out[i].DoffGAL    = in[i].DoffGAL;
     117    out[i].Mflat      = in[i].Mflat;
    120118    out[i].Sky        = in[i].Sky;
    121119    out[i].dSky       = in[i].dSky;
    … …  
    16781676    out[i].XoffDCR    = in[i].XoffDCR;
    16791677    out[i].YoffDCR    = in[i].YoffDCR;
    1680     out[i].RoffGAL    = in[i].RoffGAL;
    1681     out[i].DoffGAL    = in[i].DoffGAL;
     1678    out[i].Mflat      = in[i].Mflat;
    16821679    out[i].Sky        = in[i].Sky;
    16831680    out[i].dSky       = in[i].dSky;

  • branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_tiny_values.c

    r38062 r39266  
    2222  measureT[0].dM         = measure[0].dM;
    2323  measureT[0].airmass    = measure[0].airmass;
    24   measureT[0].RoffGAL    = measure[0].RoffGAL;
    25   measureT[0].DoffGAL    = measure[0].DoffGAL;
     24  measureT[0].Mflat      = measure[0].Mflat;
    2625  measureT[0].Xccd       = measure[0].Xccd;
    2726  measureT[0].Yccd       = measure[0].Yccd;

  • branches/eam_branches/ipp-20151113/Ohana/src/libohana/Makefile

    r39058 r39266  
    3535$(SRC)/sorts.$(ARCH).o           \
    3636$(SRC)/bisection.$(ARCH).o       \
     37$(SRC)/gaussdev.$(ARCH).o       \
    3738$(SRC)/string.$(ARCH).o          \
    3839$(SRC)/findexec.$(ARCH).o        \

  • branches/eam_branches/ipp-20151113/Ohana/src/libohana/include/ohana.h

    r39040 r39266  
    486486int sprintf_double (char *output, double value);
    487487
    488 # endif
     488/* in gaussdev.c */
     489double ohana_gaussian (double x, double mean, double sigma);
     490void ohana_gaussdev_init (void);
     491void ohana_gaussdev_free (void);
     492double ohana_gaussdev_rnd (double mean, double sigma);
     493
     494# endif

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakedist.c

    r15487 r39266  
    11# include <ohana.h>
    2 
    3 # define NGAUSS 2048
    4 double drand48();
    5 double gaussint[NGAUSS];
    6 double gaussian();
    7 double rnd_gauss();
    82
    93# define MMIN 1.0
    … …  
    4539  init_outmatrix ();
    4640 
    47   gauss_init ();
     41  ohana_gaussdev_init ();
    4842
    4943  read_datafiles ();
    … …  
    349343        if (noise < dMo) {
    350344          *Noise = noise;
    351           *mag = rnd_gauss (v, noise);
    352           *color = rnd_gauss (uv, 1.4*noise);
     345          *mag = ohana_gaussdev_rnd (v, noise);
     346          *color = ohana_gaussdev_rnd (uv, 1.4*noise);
    353347          return;
    354348        }
    … …  
    414408
    415409/*****************************************************************************/
    416 
    417 double
    418 rnd_gauss (mean, sigma)
    419 double mean, sigma;
    420 {
    421 
    422   int i;
    423   double y;
    424 
    425   y = drand48();
    426   i = NGAUSS*y;
    427   y = gaussint[i]*sigma + mean;
    428 
    429   return (y);
    430 
    431 }
    432 
    433 gauss_init ()
    434 {
    435 
    436   int i;
    437   long A, B;
    438   double val, x, dx, dx1, dx2, dx3, df;
    439   double mean, sigma;
    440 
    441   A = time(NULL);
    442   for (B = 0; A == time(NULL); B++);
    443   srand48(B);
    444 
    445   val = 0;
    446   dx = 0.001;
    447   dx1 = dx / 3.0;
    448   dx2 = 2.0*dx/3.0;
    449   dx3 = dx;
    450   mean = 0.0;
    451   sigma = 1.0;
    452 
    453   for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    454     df = (3.0*gaussian(x    , mean, sigma) +
    455           9.0*gaussian(x+dx1, mean, sigma) +
    456           9.0*gaussian(x+dx2, mean, sigma) +
    457           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    458     val += df;
    459     if (val > (i / (double) NGAUSS)) {
    460       gaussint[i] = x + dx / 2.0;
    461       i++;
    462     }
    463   }
    464   return (x + dx / 2.0);
    465 }
    466 
    467 double
    468 gaussian (x, mean, sigma)
    469 double x, mean, sigma;
    470 {
    471 
    472   double f, X;
    473 
    474   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    475 
    476   return (f);
    477 
    478 }
    479410
    480411/*****************************************************************************/
    … …  
    505436      Ngood = 1.0;
    506437      for (k = 0; k < NTRY; k++) {
    507         v = rnd_gauss (mag, noise);
    508         uv = rnd_gauss (color, 1.4*noise);
     438        v = ohana_gaussdev_rnd (mag, noise);
     439        uv = ohana_gaussdev_rnd (color, 1.4*noise);
    509440        x = (uv - UV0 - 0.7*A_V) / DUV;
    510441        y = (v - V0 - Dist - A_V) / DV;

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakepop.c

    r7080 r39266  
    22# define MMIN 1.0
    33# define MMAX 120.0
    4 extern double drand48();
    54extern double rnd_mass();
    65double rnd_mass();
    76double term;
    8 extern double rnd_gauss();
    9 extern double rnd_integrate ();
    10 extern double gaussian ();
    11 
    12 double gaussian();
    13 double rnd_gauss();
    14 double rnd_integrate();
    15 
    16 # define NGAUSS 2048
    17 double gaussint[NGAUSS];
    187
    198double AgeS, AgeE, alpha, d, Av, dMo, dVo, dVref;
    … …  
    4736  gfits_scan (&UV_h, "DEC_Y", "%lf", 1, &ldM);
    4837 
    49   gauss_init ();
     38  ohana_gaussdev_init ();
    5039
    5140  fprintf (stderr, "beginning main loop\n");
    … …  
    6857        v = gfits_get_matrix_value (&V_i, X, Y) + d + Av;
    6958        noise = dVo*sqrt(1.0 + pow (10.0, (0.4*(v - dVref))));
    70         V = rnd_gauss (v, noise);
    71         UV = rnd_gauss (uv, 1.4*noise);
     59        V = ohana_gaussdev_rnd (v, noise);
     60        UV = ohana_gaussdev_rnd (uv, 1.4*noise);
    7261        if (noise < dMo) {
    7362          fprintf (stdout, "%f %f   %f %f   %f  %f\n", V, noise, V+UV, noise, mass, age);
    … …  
    9685}
    9786
    98 
    99 double
    100 rnd_gauss (mean, sigma)
    101 double mean, sigma;
    102 {
    103 
    104   int i;
    105   double y;
    106 
    107   y = drand48();
    108   i = NGAUSS*y;
    109   y = gaussint[i]*sigma + mean;
    110 
    111   return (y);
    112 
    113 }
    114 
    115 
    116 gauss_init ()
    117 {
    118 
    119   int i;
    120   long A, B;
    121   double val, x, dx, dx1, dx2, dx3, df;
    122   double mean, sigma;
    123 
    124   A = time(NULL);
    125   for (B = 0; A == time(NULL); B++);
    126   srand48(B);
    127 
    128   val = 0;
    129   dx = 0.001;
    130   dx1 = dx / 3.0;
    131   dx2 = 2.0*dx/3.0;
    132   dx3 = dx;
    133   mean = 0.0;
    134   sigma = 1.0;
    135 
    136   for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    137     df = (3.0*gaussian(x    , mean, sigma) +
    138           9.0*gaussian(x+dx1, mean, sigma) +
    139           9.0*gaussian(x+dx2, mean, sigma) +
    140           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    141     val += df;
    142     if (val > (i / (double) NGAUSS)) {
    143       gaussint[i] = x + dx / 2.0;
    144       i++;
    145     }
    146   }
    147   return (x + dx / 2.0);
    148 }
    14987
    15088double
    … …  
    172110  return (x + dx / 2.0);
    173111}
    174 
    175 double
    176 gaussian (x, mean, sigma)
    177 double x, mean, sigma;
    178 {
    179 
    180   double f, X;
    181 
    182   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    183 
    184   return (f);
    185 
    186 }
    187 
    188112
    189113/*****************************************************************************/

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakestars.c

    r27435 r39266  
    1717  int found;
    1818} Stars;
    19 
    20 void   gauss_init ();
    21 double gaussian (double x, double mean, double sigma);
    22 double rnd_gauss (double mean, double sigma);
    2319
    2420  /*
    … …  
    4541  offset = atof (argv[4]);
    4642
    47   gauss_init ();
     43  ohana_gaussdev_init ();
    4844
    4945  /* load header, open file */
    … …  
    9894    dMr = 0.001 * stars[i].dM;
    9995    dMo = hypot (dMs, dMr);
    100     dM  = rnd_gauss (0.0, dMo);
     96    dM  = ohana_gaussdev_rnd (0.0, dMo);
    10197    stars[i].M += dM + offset;
    10298  }

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fitdist.c

    r8301 r39266  
    11# include <ohana.h>
    2 
    3 # define NGAUSS 2048
    4 double drand48();
    5 double gaussint[NGAUSS];
    6 double gaussian();
    7 double rnd_gauss();
    82
    93# define MMIN 1.0
    … …  
    4539  fprintf (stderr, "read in parameters\n");
    4640
    47   gauss_init ();
     41  ohana_gaussdev_init ();
    4842
    4943  read_datafiles ();
    … …  
    360354/*****************************************************************************/
    361355
    362 double
    363 rnd_gauss (mean, sigma)
    364 double mean, sigma;
    365 {
    366 
    367   int i;
    368   double y;
    369 
    370   y = drand48();
    371   i = NGAUSS*y;
    372   y = gaussint[i]*sigma + mean;
    373 
    374   return (y);
    375 
    376 }
    377 
    378 gauss_init ()
    379 {
    380 
    381   int i;
    382   long A, B;
    383   double val, x, dx, dx1, dx2, dx3, df;
    384   double mean, sigma;
    385 
    386   A = time(NULL);
    387   for (B = 0; A == time(NULL); B++);
    388   srand48(B);
    389 
    390   val = 0;
    391   dx = 0.001;
    392   dx1 = dx / 3.0;
    393   dx2 = 2.0*dx/3.0;
    394   dx3 = dx;
    395   mean = 0.0;
    396   sigma = 1.0;
    397 
    398   for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    399     df = (3.0*gaussian(x    , mean, sigma) +
    400           9.0*gaussian(x+dx1, mean, sigma) +
    401           9.0*gaussian(x+dx2, mean, sigma) +
    402           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    403     val += df;
    404     if (val > (i / (double) NGAUSS)) {
    405       gaussint[i] = x + dx / 2.0;
    406       i++;
    407     }
    408   }
    409   return (x + dx / 2.0);
    410 }
    411 
    412 double
    413 gaussian (x, mean, sigma)
    414 double x, mean, sigma;
    415 {
    416 
    417   double f, X;
    418 
    419   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    420 
    421   return (f);
    422 
    423 }
    424 
    425356/*****************************************************************************/
    426357
    … …  
    453384      Ngood = 1.0;
    454385      for (k = 0; k < NTRY; k++) {
    455         v = rnd_gauss (Mag, Noise);
    456         uv = rnd_gauss (Color, 1.4*Noise);
     386        v = ohana_gaussdev_rnd (Mag, Noise);
     387        uv = ohana_gaussdev_rnd (Color, 1.4*Noise);
    457388        x = (uv - UV0 - 0.7*A_V) / DUV;
    458389        y = (v - V0 - Dist - A_V) / DV;

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtoage.c

    r7080 r39266  
    22# define MMIN 1.0
    33# define MMAX 120.0
    4 extern double drand48();
    5 extern double rnd_gauss();
    6 extern double rnd_integrate ();
    7 extern double gaussian ();
    8 
    9 double gaussian();
    10 double rnd_gauss();
    11 double rnd_integrate();
    124
    135void main (argc, argv)
    … …  
    2921  long A, B;
    3022 
     23  ohana_gaussdev_init();
     24
    3125  lAo =   0.0;
    3226  ldA =   1.0;
    … …  
    118112      }
    119113      else {
    120         v = rnd_gauss (V, dV);
    121         uv = rnd_gauss ((U-V), dUV);
     114        v = ohana_gaussdev_rnd (V, dV);
     115        uv = ohana_gaussdev_rnd ((U-V), dUV);
    122116      }
    123117      x = (uv - UV0 - 0.7*Av) / DUV;
    … …  
    167161}
    168162
    169 double
    170 rnd_gauss (mean, sigma)
    171 double mean, sigma;
    172 {
    173 
    174   double range, x;
    175 
    176   range = drand48();
    177   x = rnd_integrate (*gaussian, range, mean, sigma);
    178 
    179   return (x);
    180 
    181 }
    182 
    183 
    184 double
    185 rnd_integrate (function, range, mean, sigma)
    186 double (*function) ();
    187 double range, mean, sigma;
    188 {
    189 
    190   double val, x, dx, dx1, dx2, dx3, df;
    191 
    192   range += 0.0001;
    193   val = 0;
    194   dx = sigma / 10.0;
    195   dx1 = dx / 3.0;
    196   dx2 = 2.0*dx/3.0;
    197   dx3 = dx;
    198 
    199   for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx)  {
    200     df = (3.0*function(x    , mean, sigma) +
    201           9.0*function(x+dx1, mean, sigma) +
    202           9.0*function(x+dx2, mean, sigma) +
    203           3.0*function(x+dx3, mean, sigma)) * (dx1/8.0);
    204     val += df;
    205   }
    206   return (x + dx / 2.0);
    207 }
    208 
    209 double
    210 gaussian (x, mean, sigma)
    211 double x, mean, sigma;
    212 {
    213 
    214   double f, X;
    215 
    216   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    217 
    218   return (f);
    219 
    220 }
    221 

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtomass.c

    r7080 r39266  
    11# include <ohana.h>
    2 extern double drand48();
    3 extern double rnd_gauss();
    4 extern double rnd_integrate ();
    5 extern double gaussian ();
    62# define NEWWAY 1
    7 
    8 double gaussian();
    9 double rnd_gauss();
    10 double rnd_integrate();
    11 
    12 # define NGAUSS 2048
    13 double gaussint[NGAUSS];
    143
    154void main (argc, argv)
    … …  
    5241  fprintf (stderr, "using mags in columns %d & %d\n", col1, col2);
    5342
    54   gauss_init ();
     43  ohana_gaussdev_init ();
    5544
    5645  lAo = 0.1;
    … …  
    147136        Ngood = 1.0;
    148137        for (k = 0; k < Ntry; k++) {
    149           v = rnd_gauss (V, dV);
    150           uv = rnd_gauss ((U-V), dUV);
     138          v = ohana_gaussdev_rnd (V, dV);
     139          uv = ohana_gaussdev_rnd ((U-V), dUV);
    151140          x = (uv - UV0 - 0.7*Av) / DUV;
    152141          y = (v - d - V0 - Av) / DV;
    … …  
    202191*/
    203192
    204 double
    205 rnd_gauss (mean, sigma)
    206 double mean, sigma;
    207 {
    208 
    209   int i;
    210   double y;
    211 
    212   y = drand48();
    213   i = NGAUSS*y;
    214   y = gaussint[i]*sigma + mean;
    215 
    216   return (y);
    217 
    218 }
    219 
    220 
    221 gauss_init ()
    222 {
    223 
    224   int i;
    225   long A, B;
    226   double val, x, dx, dx1, dx2, dx3, df;
    227   double mean, sigma;
    228 
    229   A = time(NULL);
    230   for (B = 0; A == time(NULL); B++);
    231   srand48(B);
    232 
    233   val = 0;
    234   dx = 0.001;
    235   dx1 = dx / 3.0;
    236   dx2 = 2.0*dx/3.0;
    237   dx3 = dx;
    238   mean = 0.0;
    239   sigma = 1.0;
    240 
    241   for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    242     df = (3.0*gaussian(x    , mean, sigma) +
    243           9.0*gaussian(x+dx1, mean, sigma) +
    244           9.0*gaussian(x+dx2, mean, sigma) +
    245           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    246     val += df;
    247     if (val > (i / (double) NGAUSS)) {
    248       gaussint[i] = x + dx / 2.0;
    249       i++;
    250     }
    251   }
    252   return (x + dx / 2.0);
    253 }
    254 
    255 double
    256 rnd_integrate (function, range, mean, sigma)
    257 double (*function) ();
    258 double range, mean, sigma;
    259 {
    260 
    261   double val, x, dx, dx1, dx2, dx3, df;
    262 
    263   range += 0.0001;
    264   val = 0;
    265   dx = sigma / 100.0;
    266   dx1 = dx / 3.0;
    267   dx2 = 2.0*dx/3.0;
    268   dx3 = dx;
    269 
    270   for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx)  {
    271     df = (3.0*function(x    , mean, sigma) +
    272           9.0*function(x+dx1, mean, sigma) +
    273           9.0*function(x+dx2, mean, sigma) +
    274           3.0*function(x+dx3, mean, sigma)) * (dx1/8.0);
    275     val += df;
    276   }
    277   return (x + dx / 2.0);
    278 }
    279 
    280 double
    281 gaussian (x, mean, sigma)
    282 double x, mean, sigma;
    283 {
    284 
    285   double f, X;
    286 
    287   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    288 
    289   return (f);
    290 
    291 }
    292 

  • branches/eam_branches/ipp-20151113/Ohana/src/misc/src/mkgauss.c

    r34088 r39266  
    11# include <ohana.h>
    2 
    3 void   gauss_init ();
    4 double gaussian (double x, double mean, double sigma);
    5 double rnd_gauss (double mean, double sigma);
    6 
    7 /* NGAUSS defines the resolution of the distribution */
    8 # define NGAUSS 2048
    9 static double gaussint[NGAUSS];
    102
    113int main (int argc, char **argv) {
    … …  
    2315  Npts = atoi (argv[3]);
    2416
    25   gauss_init ();
     17  ohana_gaussdev_init ();
    2618
    2719  for (i = 0; i < Npts; i++) {
    28     f = rnd_gauss (mean, sigma);
     20    f = ohana_gaussdev_rnd (mean, sigma);
    2921    fprintf (stdout, "%f\n", f);
    3022  }
    … …  
    3325}
    3426
    35 void gauss_init () {
    36 
    37   int i;
    38   long A, B;
    39   double val, x, dx, dx1, dx2, dx3, df;
    40 
    41   A = time(NULL);
    42   for (B = 0; A == time(NULL); B++);
    43   srand48(B);
    44 
    45   val = 0;
    46   dx = 0.001;
    47   dx1 = dx / 3.0;
    48   dx2 = 2.0*dx/3.0;
    49   dx3 = dx;
    50 
    51   for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    52     df = (3.0*gaussian(x    , 0.0, 1.0) +
    53           9.0*gaussian(x+dx1, 0.0, 1.0) +
    54           9.0*gaussian(x+dx2, 0.0, 1.0) +
    55           3.0*gaussian(x+dx3, 0.0, 1.0)) * (dx1/8.0);
    56     val += df;
    57     if (val > (i / (double) NGAUSS)) {
    58       gaussint[i] = x + dx / 2.0;
    59       i++;
    60     }
    61   }
    62   return;
    63 }
    64 
    65 double gaussian (double x, double mean, double sigma) {
    66   double f;
    67 
    68   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    69   return (f);
    70 }
    71 
    72 double rnd_gauss (double mean, double sigma) {
    73 
    74   int i;
    75   double y;
    76 
    77   i = 0;
    78   while ((i == 0) || (i == NGAUSS - 1)) {
    79     y = drand48();
    80     i = MIN (NGAUSS - 1, MAX (0, NGAUSS*y));
    81   }
    82   y = gaussint[i]*sigma + mean;
    83 
    84   return (y);
    85 }

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile

    r39223 r39266  
    2525$(SRC)/cneedles.$(ARCH).o                  \
    2626$(SRC)/cplot.$(ARCH).o             \
     27$(SRC)/crotation.$(ARCH).o         \
    2728$(SRC)/csystem.$(ARCH).o           \
    2829$(SRC)/ctimes.$(ARCH).o    \
    … …  
    3738$(SRC)/fitpm_irls.$(ARCH).o  \
    3839$(SRC)/fitplx_irls.$(ARCH).o  \
     40$(SRC)/astrom_ops.$(ARCH).o        \
    3941$(SRC)/fixwrap.$(ARCH).o           \
    4042$(SRC)/fixcols.$(ARCH).o           \

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdensify.c

    r38463 r39266  
    8484  gfits_free_matrix (&bf[0].matrix);
    8585  gfits_free_header (&bf[0].header);
    86   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     86  if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    8787  strcpy (bf[0].file, "(empty)");
    8888  PutCoords (&newcoords, &bf[0].header);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdhistogram.c

    r38062 r39266  
    101101  gfits_free_matrix (&bf[0].matrix);
    102102  gfits_free_header (&bf[0].header);
    103   CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0);
     103  if (!CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0)) return FALSE;
    104104  strcpy (bf[0].file, "(empty)");
    105105  PutCoords (&newcoords, &bf[0].header);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitplx.c

    r38986 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double *X;
    10   double *Y;
    11   double *t;
    12   double *pX;
    13   double *pY;
    14   double *dX;
    15   double *dY;
    16   int *index;
    17   int Npts;
    18 } PlxFitData;
    19 
    20 typedef struct {
    21   double Ro, dRo;
    22   double Do, dDo;
    23 
    24   double uR, duR;
    25   double uD, duD;
    26 
    27   double p, dp;
    28 
    29   double chisq;
    30   int Nfit;
    31   int getChisq;
    32 } PlxFit;
    33 
    34 int VectorRobustStats (Vector *vector, double *median, double *sigma);
    35 double VectorFractionInterpolate (double *values, float fraction, int Npts);
    36 
    37 int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
    38 int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
    39 int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
    40 
    41 int PlxFitDataAlloc (PlxFitData *data, int N);
    42 void PlxFitDataFree (PlxFitData *data);
    43 int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
    44 
    45 int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    46 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
    47 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time);
    482
    493int fitplx (int argc, char **argv) {
    … …  
    10054    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
    10155    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
     56    gprint (GP_ERR, "  -mask mask : excluded points are marked with a 0 mask value\n");
    10257    return (FALSE);
    10358  }
    … …  
    148103  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
    149104
    150   // determine dPsig for detections based on Noutlier attempts
     105  // determine dPsig for detections based on Noutlier attempts (mask is updated to mark the bad stars, mask == 0)
    151106  if (Noutlier) {
    152107    int clipRetry = TRUE;
    153108    for (i = 0; clipRetry && (i < 3); i++) {
     109      // XXX NOTE This will segfault is mask is not supplied:
    154110      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
    155111
    … …  
    163119  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
    164120    int n = fitdata.index[i];
    165     int maskValue = mask ? mask[n] : 0;
     121    int maskValue = mask ? mask[n] : 1;
    166122    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
    167123  }
    … …  
    401357  fit[0].Nfit = Npts;
    402358  return (TRUE);
    403 }
    404 
    405 /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
    406 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) {
    407 
    408   double n = mjd - J2000;                          // day number relative to standard epoch
    409   double L = 280.460 + 0.9856474 * n;              // mean solar longitute (corr. for aberration)
    410   double g = (357.528 + 0.9856003 * n)*RAD_DEG;    // Mean anomaly
    411 
    412   *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
    413   *beta = 0.0;                                     // approx latitude
    414   *epsilon = (23.439 - 0.0000004 * n);             // obliquity of ecliptic in degrees
    415   *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
    416   return TRUE;
    417 }
    418 
    419 /* given RA, DEC, Time, calculate the parallax factor */
    420 // RA,DEC are decimal degrees
    421 // Time is MJD
    422 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) {
    423 
    424   double lambda, beta, epsilon, Radius;
    425 
    426   // Time must be mjd
    427   sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius);
    428 
    429   double lambda_rad = lambda*RAD_DEG;
    430   double epsilon_rad = epsilon*RAD_DEG;
    431   double RA_rad = RA*RAD_DEG;
    432   double DEC_rad = DEC*RAD_DEG;
    433 
    434   double x = Radius*cos(lambda_rad);
    435   double y = Radius*cos(epsilon_rad)*sin(lambda_rad);
    436   double z = Radius*sin(epsilon_rad)*sin(lambda_rad);
    437 
    438   *pR = +(y*cos(RA_rad) - x*sin(RA_rad));
    439   *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad);
    440 
    441   return TRUE;
    442 }
    443 
    444 // allocate arrays but not the container
    445 int PlxFitDataAlloc (PlxFitData *data, int N) {
    446 
    447   data->Npts = N;
    448   ALLOCATE (data->X, double, N);
    449   ALLOCATE (data->Y, double, N);
    450   ALLOCATE (data->dX, double, N);
    451   ALLOCATE (data->dY, double, N);
    452   ALLOCATE (data->t, double, N);
    453   ALLOCATE (data->pX, double, N);
    454   ALLOCATE (data->pY, double, N);
    455   ALLOCATE (data->index, int, N);
    456   return TRUE;
    457 }
    458 
    459 void PlxFitDataFree (PlxFitData *data) {
    460   FREE (data->X);
    461   FREE (data->Y);
    462   FREE (data->dX);
    463   FREE (data->dY);
    464   FREE (data->t);
    465   FREE (data->pX);
    466   FREE (data->pY);
    467   FREE (data->index);
    468359}
    469360
    … …  
    481372    tgt->pX[i] = src->pX[N];
    482373    tgt->pY[i] = src->pY[N];
     374
     375    // *** make this optional?
     376    tgt->Wx[i] = src->Wx[N];
     377    tgt->Wy[i] = src->Wy[N];
    483378  }
    484379  return TRUE;
    … …  
    539434    pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
    540435    pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
     436
     437    fitdata->Wx[Nsubset] = 1.0;
     438    fitdata->Wy[Nsubset] = 1.0;   
    541439    fitdata->index[Nsubset] = i;
    542440    Nsubset++;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm.c

    r37807 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double Ro, dRo;
    10   double Do, dDo;
    11 
    12   double uR, duR;
    13   double uD, duD;
    14 
    15   double chisq;
    16   int Nfit;
    17 } PMFit;
    18 
    19 int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
    202
    213int fitpm (int argc, char **argv) {
    … …  
    122104  }
    123105
    124   PMFit fit;
     106  PlxFit fit;
    125107  if (!FitPMonly (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
    126108    return FALSE;
    … …  
    158140
    159141/* do we want an init function which does the alloc and a clear function to free? */
    160 int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
     142int FitPMonly (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
    161143
    162144  int i;
    … …  
    241223  fit[0].duD = sqrt(A[3][3]);
    242224 
     225  fit[0].p   = 0.0;
     226  fit[0].dp  = NAN;
     227
    243228  // add up the chi square for the fit
    244229  chisq = 0.0;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c

    r39223 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double Ro, dRo;
    10   double Do, dDo;
    11 
    12   double uR, duR;
    13   double uD, duD;
    14  
    15   double chisq;
    16   int Nfit;
    17 } PMFit_IRLS;
    18 
    19 int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
    20 int IRLS_converged (PMFit_IRLS *fit);
    21 int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
    22                  double **A, double **B, int VERBOSE);
    23 double weight_cauchy (double x);
    24 double dpsi_cauchy (double x);
    25 double MAD(double *in, int N);
    26      
    272
    283int fitpm_irls (int argc, char **argv) {
    … …  
    129104  }
    130105
    131   PMFit_IRLS fit;
     106  PlxFit fit;
    132107  if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
    133108    return FALSE;
    … …  
    165140
    166141/* do we want an init function which does the alloc and a clear function to free? */
    167 int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
     142int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
    168143
    169144  int i,j;
    … …  
    222197 
    223198  // Solve OLS equation 
    224   if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
     199  if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
    225200                   A,B,VERBOSE)) {
    226201    // Handle fail case
    … …  
    267242
    268243    // Solve
    269     if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
     244    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
    270245                     A,B,VERBOSE)) {
    271246      // Handle fail case
    … …  
    284259      u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
    285260    }
    286     sigma_hat = MAD(u,Npts) / 0.6745;
     261    sigma_hat = MedianAbsDeviation(u,Npts) / 0.6745;
    287262   
    288263    // Check convergence
    … …  
    310285  double Sum_Wx, Sum_Wy;
    311286 
     287  Sum_Wx = 0.0;
     288  Sum_Wy = 0.0;
    312289  ax = 0.0; ay = 0.0;
    313290  bx = 0.0; by = 0.0;
    … …  
    387364}
    388365
    389 
    390 double weight_cauchy (double x) {
    391   double r = x / 2.385;
    392   return (1.0 / (1.0 + SQ(r)));
    393 }
    394 
    395 // dpsi = (d/dx) (x * weight(x))
    396 double dpsi_cauchy (double x) {
    397   double r2 = SQ(x / 2.385);
    398   return ((1.0 - r2) / (SQ(1 + r2)));
    399 }
    400 
    401 
    402 // median absolute deviation
    403 // MAD = median(abs(x - median(x)))
    404 double MAD(double *in, int N) {
    405   double *x;
    406   double median = 0.0;
    407   int i;
    408  
    409   ALLOCATE(x,double,N);
    410   for (i = 0; i < N; i++) {
    411     x[i] = in[i];
    412   }
    413 
    414   dsort(x,N);
    415 
    416   if (N % 2) {
    417     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    418   } else {
    419     median = x[(int)(0.5*N)];
    420   }
    421 
    422   for (i = 0; i < N; i++ ) {
    423     x[i] = fabs(x[i] - median);
    424   }
    425 
    426   dsort(x,N);
    427 
    428   if (N % 2) {
    429     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    430   } else {
    431     median = x[(int)(0.5*N)];
    432   }
    433 
    434   return(median);
    435 }
    436    
    437  
    438  
    439 int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
    440                  double **A, double **B, int VERBOSE) {
     366int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE) {
    441367
    442368  int i,j;
    … …  
    497423  return TRUE;
    498424}
     425
     426double weight_cauchy (double x) {
     427  double r = x / 2.385;
     428  return (1.0 / (1.0 + SQ(r)));
     429}
     430
     431// dpsi = (d/dx) (x * weight(x))
     432double dpsi_cauchy (double x) {
     433  double r2 = SQ(x / 2.385);
     434  return ((1.0 - r2) / (SQ(1 + r2)));
     435}
     436
     437
     438// median absolute deviation
     439// MAD = median(abs(x - median(x)))
     440double MedianAbsDeviation(double *in, int N) {
     441  double *x;
     442  double median = 0.0;
     443  int i;
     444 
     445  ALLOCATE(x,double,N);
     446  for (i = 0; i < N; i++) {
     447    x[i] = in[i];
     448  }
     449
     450  dsort(x,N);
     451
     452  if (N % 2) {
     453    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
     454  } else {
     455    median = x[(int)(0.5*N)];
     456  }
     457
     458  for (i = 0; i < N; i++ ) {
     459    x[i] = fabs(x[i] - median);
     460  }
     461
     462  dsort(x,N);
     463
     464  if (N % 2) {
     465    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
     466  } else {
     467    median = x[(int)(0.5*N)];
     468  }
     469
     470  return(median);
     471}
     472 

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fiximage.c

    r34088 r39266  
    2626  gfits_free_matrix (&mask[0].matrix);
    2727  gfits_free_header (&mask[0].header);
    28   CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0);
     28  if (!CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    2929  strcpy (mask[0].file, "(empty)");
    3030  memset (mask[0].matrix.buffer, 0, Nx*Ny*sizeof(float));

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/getcoords.c

    r37807 r39266  
    7373  gfits_free_matrix (&outC1[0].matrix);
    7474  gfits_free_header (&outC1[0].header);
    75   CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0);
     75  if (!CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    7676
    7777  gfits_free_matrix (&outC2[0].matrix);
    7878  gfits_free_header (&outC2[0].header);
    79   CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0);
     79  if (!CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    8080
    8181  float *valC1 = (float *) outC1[0].matrix.buffer;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/imfit.c

    r34088 r39266  
    141141
    142142    strcpy (out[0].file, "(empty)");
    143     CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0);
     143    if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE;
    144144
    145145    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c

    r39223 r39266  
    99int cneedles                PROTO((int, char **));
    1010int cplot                   PROTO((int, char **));
     11int crotation               PROTO((int, char **));
    1112int csystem                 PROTO((int, char **));
    1213int ctimes                  PROTO((int, char **));
    … …  
    7879  {1, "cneedles",    cneedles,     "plot vectors in sky coordinates"},
    7980  {1, "cplot",       cplot,        "plot vectors in sky coordinates"},
     81  {1, "crotation",   crotation,    "rotate in 3D"},
    8082  {1, "csystem",     csystem,      "convert between coordinate systems"},
    8183  {1, "ctimes",      ctimes,       "convert between time formats"},

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/shimage.c

    r37807 r39266  
    3030  int Ny = src[0].header.Naxis[1];
    3131
    32   gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0);
    33   gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0);
     32  gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); if (!CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     33  gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); if (!CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3434
    3535  // l=0 allocates space for a single value

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/vshimage.c

    r37807 r39266  
    3232  int Ny = src[0].header.Naxis[1];
    3333
    34   gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0);
    35   gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0);
    36   gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0);
    37   gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0);
     34  gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); if (!CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     35  gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); if (!CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     36  gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); if (!CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     37  gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); if (!CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3838
    3939  // l=0 allocates space for a single value

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/densify.c

    r36489 r39266  
    9696  gfits_free_matrix (&bf[0].matrix);
    9797  gfits_free_header (&bf[0].header);
    98   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     98  if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    9999  strcpy (bf[0].file, "(empty)");
    100100 

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/fft2d.c

    r16107 r39266  
    5858 
    5959  /* fix up output headers (real) & allocate data buffer */
    60   CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0);
    61   CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0);
     60  if (!CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
     61  if (!CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    6262
    6363  gfits_copy_header (&Ire[0].header, &Ore[0].header);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gaussdeviate.c

    r34088 r39266  
    1717  ResetVector (vec, OPIHI_FLT, Npts);
    1818
    19   gauss_init (2048);
     19  ohana_gaussdev_init ();
    2020  for (i = 0; i < Npts; i++) {
    21     vec[0].elements.Flt[i] = rnd_gauss (mean, sigma);
     21    vec[0].elements.Flt[i] = ohana_gaussdev_rnd (mean, sigma);
    2222  }
    2323  return (TRUE);
    … …  
    2929}
    3030
    31 double int_gauss (int i);
    32 
    3331int gaussintegral (int argc, char **argv) {
    3432 
     33  gprint (GP_ERR, "fix gaussintegral\n");
     34  return (FALSE);
     35
     36/*
    3537  int i, Npts;
    3638  Vector *vec;
    … …  
    4648  ResetVector (vec, OPIHI_FLT, Npts);
    4749
    48   gauss_init (Npts);
     50  ohana_gaussdev_init ();
    4951  for (i = 0; i < Npts; i++) {
    50     vec[0].elements.Flt[i] = int_gauss (i);
     52    vec[0].elements.Flt[i] = gaussian_int (i);
    5153  }
    5254  return (TRUE);
    … …  
    5557  gprint (GP_ERR, "USAGE: gaussintegral Npts mean sigma\n");
    5658  return (FALSE);
     59*/
    5760   
    5861}

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gridify.c

    r33963 r39266  
    8585    gfits_free_matrix (&bf[0].matrix);
    8686    gfits_free_header (&bf[0].header);
    87     CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     87    if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    8888    strcpy (bf[0].file, "(empty)");
    8989  }
    … …  
    111111    gfits_free_matrix (&ct[0].matrix);
    112112    gfits_free_header (&ct[0].header);
    113     CreateBuffer (ct, Nx, Ny, -32, 0.0, 1.0);
     113    if (!CreateBuffer (ct, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    114114    strcpy (ct[0].file, "(empty)");
    115115

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mcreate.c

    r38062 r39266  
    2727
    2828  if (Nz) {
    29     CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0);
     29    if (!CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0)) return FALSE;
    3030  } else {
    31     CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0);
     31    if (!CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3232  }
    3333  return (TRUE);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/medimage_commands.c

    r36679 r39266  
    8282  gfits_free_matrix (&output->matrix);
    8383  gfits_free_header (&output->header);
    84   CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0);
     84  if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    8585
    8686  float *outvalue = (float *) output->matrix.buffer;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mslice.c

    r38234 r39266  
    3535  gfits_free_matrix (&out[0].matrix);
    3636  gfits_free_header (&out[0].header);
    37   CreateBuffer (out, Nx, Ny, -32, 1.0, 0.0);
     37  if (!CreateBuffer (out, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3838
    3939  float *inF  = (float *) in[0].matrix.buffer + plane*Nx*Ny;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/read_vectors.c

    r39223 r39266  
    348348}
    349349
     350# undef ESCAPE
    350351# define ESCAPE(...) {          \
    351352  gprint (GP_ERR, __VA_ARGS__); \

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/reindex.c

    r34461 r39266  
    11# include "data.h"
    2 
    3 # define ESCAPE(MSG,...){ gprint (GP_ERR, MSG, __VA_ARGS__); goto error; }
    42
    53int reindex (int argc, char **argv) {

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/shift.c

    r34088 r39266  
    4848  gfits_free_matrix (&out[0].matrix);
    4949  gfits_free_header (&out[0].header);
    50   CreateBuffer (out, nx, ny, -32, 0.0, 1.0);
     50  if (!CreateBuffer (out, nx, ny, -32, 0.0, 1.0)) return FALSE;
    5151
    5252  DXin = (dx < 0) ? -dx : 0;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/squash3d.c

    r38062 r39266  
    4646        int oNx = iNz;
    4747        int oNy = iNy;
    48         CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
     48        if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE;
    4949        float *oBuf  = (float *) tgt[0].matrix.buffer;
    5050
    … …  
    7070        int oNx = iNx;
    7171        int oNy = iNz;
    72         CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
     72        if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE;
    7373        float *oBuf  = (float *) tgt[0].matrix.buffer;
    7474
    … …  
    9494        int oNx = iNx;
    9595        int oNy = iNy;
    96         CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
     96        if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE;
    9797        float *oBuf  = (float *) tgt[0].matrix.buffer;
    9898

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/tdhistogram.c

    r38062 r39266  
    140140    gfits_free_matrix (&bf[0].matrix);
    141141    gfits_free_header (&bf[0].header);
    142     CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0);
     142    if (!CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0)) return FALSE;
    143143    strcpy (bf[0].file, "(empty)");
    144144

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/vgrid.c

    r20936 r39266  
    4343  gfits_free_matrix (&bf[0].matrix);
    4444  gfits_free_header (&bf[0].header);
    45   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     45  if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    4646  strcpy (bf[0].file, "(empty)");
    4747

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dimm/camera_cmds.c

    r23530 r39266  
    107107    gfits_free_matrix (&buf[0].matrix);
    108108    gfits_free_header (&buf[0].header);
    109     CreateBuffer (buf, dx, dy, -32, 0.0, 1.0);
     109    if (!CreateBuffer (buf, dx, dy, -32, 0.0, 1.0)) return FALSE;
    110110    strcpy (buf[0].file, "(empty)");
    111111

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordimage.c

    r36375 r39266  
    7979  gfits_free_matrix (&bufX[0].matrix);
    8080  gfits_free_header (&bufX[0].header);
    81   CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here
     81  if (!CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here
    8282  strcpy (bufX[0].file, "(empty)");
    8383
    8484  gfits_free_matrix (&bufY[0].matrix);
    8585  gfits_free_header (&bufY[0].header);
    86   CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here
     86  if (!CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here
    8787  strcpy (bufY[0].file, "(empty)");
    8888

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordmosaic.c

    r36833 r39266  
    8888  gfits_free_matrix (&bufX[0].matrix);
    8989  gfits_free_header (&bufX[0].header);
    90   CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here
     90  if (!CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here
    9191  strcpy (bufX[0].file, "(empty)");
    9292
    9393  gfits_free_matrix (&bufY[0].matrix);
    9494  gfits_free_header (&bufY[0].header);
    95   CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here
     95  if (!CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here
    9696  strcpy (bufY[0].file, "(empty)");
    9797
    … …  
    224224}
    225225
     226# undef ESCAPE
    226227# define ESCAPE(MSG) { \
    227228    if (src) fclose (src); \

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/fitcolors.c

    r38471 r39266  
    145145  gfits_free_matrix (&buf[0].matrix);
    146146  gfits_free_header (&buf[0].header);
    147   CreateBuffer (buf, NP, NP, -32, 0.0, 1.0);
     147  if (!CreateBuffer (buf, NP, NP, -32, 0.0, 1.0)) return FALSE;
    148148  strcpy (buf[0].file, "(empty)");
    149149

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/imphot.c

    r34088 r39266  
    3636  if (GreyScale) {
    3737    if ((buf = SelectBuffer (bufname, ANYBUFFER, TRUE)) == NULL) return (FALSE);
    38     CreateBuffer (buf, 100, 200, -32, 0.0, 1.0);
     38    if (!CreateBuffer (buf, 100, 200, -32, 0.0, 1.0)) return FALSE;
    3939  }
    4040

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/objectcoverage.c

    r38471 r39266  
    144144  gfits_free_matrix (&buf[0].matrix);
    145145  gfits_free_header (&buf[0].header);
    146   CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0);
     146  if (!CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    147147  strcpy (buf[0].file, "(empty)");
    148148

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/skycoverage.c

    r37807 r39266  
    201201  gfits_free_matrix (&buf[0].matrix);
    202202  gfits_free_header (&buf[0].header);
    203   CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0);
     203  if (!CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    204204  strcpy (buf[0].file, "(empty)");
    205205

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/astro.h

    r37807 r39266  
    1212void FreeAstro (void);
    1313
     14typedef struct {
     15  double *X;
     16  double *Y;
     17  double *t;
     18  double *pX;
     19  double *pY;
     20  double *dX;
     21  double *dY;
     22  double *Wx;
     23  double *Wy;
     24  int *index;
     25  int Npts;
     26} PlxFitData;
     27
     28typedef struct {
     29  double Ro, dRo;
     30  double Do, dDo;
     31
     32  double uR, duR;
     33  double uD, duD;
     34
     35  double p, dp;
     36
     37  double chisq;
     38  int Nfit;
     39  int getChisq;
     40} PlxFit;
     41
     42int VectorRobustStats (Vector *vector, double *median, double *sigma);
     43double VectorFractionInterpolate (double *values, float fraction, int Npts);
     44
     45int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
     46int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
     47int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
     48
     49int PlxFitDataAlloc (PlxFitData *data, int N);
     50void PlxFitDataFree (PlxFitData *data);
     51int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
     52
     53int FitPMonly (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
     54int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
     55int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
     56int ParFactor (double *pR, double *pD, double RA, double DEC, double Time);
     57
     58/***** */
     59
     60int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
     61int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, int VERBOSE);
     62
     63int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);
     64int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);
     65
     66double weight_cauchy (double x);
     67double dpsi_cauchy (double x);
     68double MedianAbsDeviation(double *in, int N);
     69
    1470# endif

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/shell.h

    r37049 r39266  
    44# ifndef SHELL_H
    55# define SHELL_H
     6
     7# define ESCAPE(MSG,...) { gprint (GP_ERR, MSG, __VA_ARGS__); return FALSE; }
    68
    79# define ISVAR(a) (isalnum (a) || (a == ':') || (a == '_'))

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/Makefile

    r36679 r39266  
    3030$(SDIR)/starfuncs.$(ARCH).o             \
    3131$(SDIR)/hermitian.$(ARCH).o             \
    32 $(SDIR)/gaussian.$(ARCH).o              \
    3332$(SDIR)/graphtools.$(ARCH).o            \
    3433$(SDIR)/queues.$(ARCH).o                \
    … …  
    3635$(SDIR)/open_kapa.$(ARCH).o             \
    3736$(SDIR)/style_args.$(ARCH).o
     37
     38# moved to libohana
     39# $(SDIR)/gaussian.$(ARCH).o            \
    3840
    3941# fix malloc

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/spline.c

    r39230 r39266  
    3232  float dx, a, b, value;
    3333 
     34  // saturate correction at high and low ends
     35  if (X < x[0]) return y[0];
     36  if (X > x[N-1]) return y[N-1];
     37
    3438  /* find correct element in array (x must be sorted) */
    3539  lo = 0;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/BufferOps.c

    r38062 r39266  
    9292  gfits_free_matrix (&buf[0].matrix);
    9393  gfits_free_header (&buf[0].header);
    94   CreateBuffer (buf, Nx, Ny, bitpix, bzero, bscale);
     94  if (!CreateBuffer (buf, Nx, Ny, bitpix, bzero, bscale)) return FALSE;
    9595  return TRUE;
    9696}
    9797
    9898int CreateBuffer (Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale) {
     99
     100  if (Nx < 0) { gprint (GP_ERR, "invalid negative matrix size Nx : %d\n", Nx);  return FALSE; }
     101  if (Ny < 0) { gprint (GP_ERR, "invalid negative matrix size Ny : %d\n", Ny);  return FALSE; }
    99102
    100103  /* store the default output values */
    … …  
    123126// buffer[x + Nx*y + Nx*Ny*z + ...]
    124127int CreateBuffer3D (Buffer *buf, int Nx, int Ny, int Nz, int bitpix, float bzero, float bscale) {
     128
     129  if (Nx < 0) { gprint (GP_ERR, "invalid negative matrix size Nx : %d\n", Nx);  return FALSE; }
     130  if (Ny < 0) { gprint (GP_ERR, "invalid negative matrix size Ny : %d\n", Ny);  return FALSE; }
     131  if (Nz < 0) { gprint (GP_ERR, "invalid negative matrix size Nz : %d\n", Nz);  return FALSE; }
    125132
    126133  /* store the default output values */

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/ListOps.c

    r36681 r39266  
    212212  if (temp != NULL) {
    213213      if (!strcmp (temp, "-x")) goto escape;
     214      if (!strcmp (temp, "-glob")) goto escape;
    214215      if (!strcmp (temp, "-split")) goto escape;
    215216      if (!strcmp (temp, "-splitbychar")) goto escape;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/mana/simsignal.c

    r7917 r39266  
    66  float *buf;
    77  double cvalue, dvalue, sigma, SN, period;
    8   double rnd_gauss ();
    98
    109  if (argc != 5) {
    … …  
    2827    ivalue = scale * 0.5 * (sin (i*2*M_PI/period) + 1) + 0.5;
    2928    /*
    30     dvalue = rnd_gauss (cvalue, sigma);
     29    dvalue = ohana_gaussdev_rnd (cvalue, sigma);
    3130    cvalue = (dvalue + range) / (2.0*range);
    3231    dvalue = MAX (0, MIN (0.99999, cvalue));
    … …  
    44438 bit = 2^8
    4544
    46   gauss_init (2*scale);
     45  ohana_gaussdev_init ();
    4746  sigma = 2.0 / SN;
    4847  range = 1 + 5*sigma;

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/Makefile

    r38986 r39266  
    1010MAN     =       $(HOME)/doc
    1111INC     =       $(HOME)/include
     12TESTDIR =       $(HOME)/test
     13TESTBIN =       $(HOME)/test
    1214include ../../Makefile.Common
    1315
    … …  
    2224install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client
    2325
    24 # $(DESTBIN)/testparallax
    25 
    2626RELASTRO = \
    2727$(SRC)/ConfigInit.$(ARCH).o          \
    … …  
    2929$(SRC)/FitMosaic.$(ARCH).o           \
    3030$(SRC)/FitSimple.$(ARCH).o           \
    31 $(SRC)/FitAstromOps.$(ARCH).o           \
     31$(SRC)/FitAstromOps.$(ARCH).o        \
    3232$(SRC)/FitPM.$(ARCH).o               \
     33$(SRC)/FitPM_IRLS.$(ARCH).o          \
    3334$(SRC)/FitPMandPar.$(ARCH).o         \
     35$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
    3436$(SRC)/FitPosPMfixed.$(ARCH).o       \
     37$(SRC)/BootstrapOps.$(ARCH).o        \
    3538$(SRC)/ImageOps.$(ARCH).o            \
    3639$(SRC)/MosaicOps.$(ARCH).o           \
    … …  
    103106$(SRC)/ConfigInit.$(ARCH).o          \
    104107$(SRC)/FitSimple.$(ARCH).o           \
    105 $(SRC)/FitAstromOps.$(ARCH).o           \
     108$(SRC)/FitAstromOps.$(ARCH).o        \
    106109$(SRC)/FitPM.$(ARCH).o               \
     110$(SRC)/FitPM_IRLS.$(ARCH).o          \
    107111$(SRC)/FitPMandPar.$(ARCH).o         \
     112$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
    108113$(SRC)/FitPosPMfixed.$(ARCH).o       \
     114$(SRC)/BootstrapOps.$(ARCH).o        \
    109115$(SRC)/ImageOps.$(ARCH).o            \
    110116$(SRC)/MosaicOps.$(ARCH).o           \
    … …  
    159165$(TESTPAR): $(INC)/relastro.h
    160166$(BIN)/testparallax.$(ARCH): $(TESTPAR)
     167
     168fitpm: $(BIN)/fitpm.$(ARCH)
     169
     170FITPM = \
     171$(SRC)/fitpm.$(ARCH).o        \
     172$(SRC)/ParFactor.$(ARCH).o           \
     173$(SRC)/FitAstromOps.$(ARCH).o        \
     174$(SRC)/FitPM.$(ARCH).o               \
     175$(SRC)/FitPM_IRLS.$(ARCH).o               \
     176$(SRC)/FitPMandPar.$(ARCH).o               \
     177$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
     178$(SRC)/BootstrapOps.$(ARCH).o        \
     179$(SRC)/mkpolyterm.$(ARCH).o            \
     180$(SRC)/fitpoly.$(ARCH).o
     181
     182$(FITPM): $(INC)/relastro.h
     183$(BIN)/fitpm.$(ARCH): $(FITPM)

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/include/relastro.h

    r39207 r39266  
    134134  double **A;
    135135  double **B;
     136  double **Cov;
     137  double *Beta;
     138  double *Beta_prev;
    136139  int Nterms;
     140  int getChisq;
     141  int getError;
    137142} FitAstromData;
    138143
     144// XXX do we need doubles for all of these?  I actually only have of order 100 of these
     145// allocated at a time, so size is not an issue.
    139146typedef struct {
    140147  double X, dX;
    … …  
    143150  double D, dD;
    144151  double T, dT;
     152  double Wx, Wy;
     153  double rx, ry;
     154  double u;
    145155  double pR;
    146156  double pD;
    … …  
    148158  double C_red;
    149159  int measure;
     160  int mask;
    150161} FitAstromPoint;
    151162
    … …  
    164175  FitAstromPoint *points;
    165176  FitAstromPoint *sample;
     177  FitAstromPoint *nomask;
    166178  int Npoints;
    167179  int NpointsAlloc;
    … …  
    723735int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode);
    724736int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints);
     737int BootstrapSaveUnmasked (FitAstromPoint *nomask, FitAstromPoint *points, int Npoints);
     738
    725739int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
    726740int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);
    … …  
    739753int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average);
    740754void AstromErrorSetLoop (int Nloop, int isImageMode);
     755
     756int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     757int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     758
     759double MedianAbsDeviation(FitAstromPoint *points, int Npoints);
     760
     761int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     762int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     763
     764double weight_cauchy (double x);
     765double dpsi_cauchy (double x);

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BrightCatalog.c

    r39212 r39266  
    6767    GET_COLUMN(R,         "RA",         double);
    6868    GET_COLUMN(D,         "DEC",        double);
    69     GET_COLUMN(RoffGAL,   "ROFF_GAL",   float);
    70     GET_COLUMN(DoffGAL,   "DOFF_GAL",   float);
    7169    GET_COLUMN(M,         "MAG_SYS",    float);
    7270    GET_COLUMN(Mcal,      "MAG_CAL",    float);
    … …  
    9492      measure[i].R         = R[i];
    9593      measure[i].D         = D[i];
    96       measure[i].RoffGAL   = RoffGAL[i];
    97       measure[i].DoffGAL   = DoffGAL[i];
    9894      measure[i].M         = M[i];
    9995      measure[i].Mcal      = Mcal[i];
    … …  
    121117    free (R       );
    122118    free (D       );
    123     free (RoffGAL );
    124     free (DoffGAL );
    125119    free (M       );
    126120    free (Mcal    );
    … …  
    362356    gfits_define_bintable_column (&theader, "D", "RA",       "ra",                         "degrees", 1.0, 0.0);
    363357    gfits_define_bintable_column (&theader, "D", "DEC",      "dec",                        "degrees", 1.0, 0.0);
    364     gfits_define_bintable_column (&theader, "E", "ROFF_GAL", "ra offset",                  "arcsec", 1.0, 0.0);
    365     gfits_define_bintable_column (&theader, "E", "DOFF_GAL", "dec offset",                 "arcsec", 1.0, 0.0);
    366358    gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
    367359    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
    … …  
    382374
    383375    // generate the output array that carries the data
    384     ohana_memcheck (1);
    385376    gfits_create_table (&theader, &ftable);
    386     ohana_memcheck (1);
    387377
    388378    // create intermediate storage arrays
    389379    double *R         ; ALLOCATE (R        ,  double, catalog->Nmeasure);
    390380    double *D         ; ALLOCATE (D        ,  double, catalog->Nmeasure);
    391     float  *RoffGAL   ; ALLOCATE (RoffGAL  ,  float,  catalog->Nmeasure);
    392     float  *DoffGAL   ; ALLOCATE (DoffGAL  ,  float,  catalog->Nmeasure);
    393381    float  *M         ; ALLOCATE (M        ,  float,  catalog->Nmeasure);
    394382    float  *Mcal      ; ALLOCATE (Mcal     ,  float,  catalog->Nmeasure);
    … …  
    407395    int    *catID     ; ALLOCATE (catID    ,  int  ,  catalog->Nmeasure);
    408396    short  *photcode  ; ALLOCATE (photcode ,  short,  catalog->Nmeasure);
    409     ohana_memcheck (1);
    410397
    411398    // assign the storage arrays
    … …  
    414401      R[i]        = measure[i].R        ;
    415402      D[i]        = measure[i].D        ;
    416       RoffGAL[i]  = measure[i].RoffGAL  ;
    417       DoffGAL[i]  = measure[i].DoffGAL  ;
    418403      M[i]        = measure[i].M        ;
    419404      Mcal[i]     = measure[i].Mcal     ;
    … …  
    433418      photcode[i] = measure[i].photcode ;
    434419    }
    435     ohana_memcheck (1);
    436420
    437421    // add the columns to the output array
    438422    gfits_set_bintable_column (&theader, &ftable, "RA",         R,         catalog->Nmeasure);
    439     ohana_memcheck (1);
    440423    gfits_set_bintable_column (&theader, &ftable, "DEC",        D,         catalog->Nmeasure);
    441     ohana_memcheck (1);
    442     gfits_set_bintable_column (&theader, &ftable, "ROFF_GAL",   RoffGAL,   catalog->Nmeasure);
    443     gfits_set_bintable_column (&theader, &ftable, "DOFF_GAL",   DoffGAL,   catalog->Nmeasure);
    444424    gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
    445425    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
    … …  
    461441    free (R       );
    462442    free (D       );
    463     free (RoffGAL );
    464     free (DoffGAL );
    465443    free (M       );
    466444    free (Mcal    );

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitAstromOps.c

    r38986 r39266  
    22
    33// I am modifying FitPM with an eye to (a) threaded operations and (b) bootstrap resampling tests.
     4
     5int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
     6
     7  int k;
     8
     9  int Npoints = fitStats->Npoints;
     10  FitAstromPoint *points = fitStats->points;
     11
     12  // find Tmin & Tmax from the list of accepted measurements
     13  double Tmin  = points[0].T;
     14  double Tmax  = points[0].T;
     15  double pRmin = +2.0;
     16  double pRmax = -2.0;
     17  double pDmin = +2.0;
     18  double pDmax = -2.0;
     19
     20  *Tmean = 0.0;
     21
     22  double Tsum = 0.0;
     23  double Wsum = 0.0;
     24  for (k = 0; k < Npoints; k++) {
     25    Tmin = MIN(Tmin, points[k].T);
     26    Tmax = MAX(Tmax, points[k].T);
     27
     28    float wx = 1.0 / SQ(points[k].dX);
     29
     30    Tsum += points[k].T * wx;
     31    Wsum += wx;
     32
     33    // at this point, T is in years since J2000
     34    ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
     35    pRmin = MIN (pRmin, points[k].pR);
     36    pRmax = MAX (pRmax, points[k].pR);
     37    pDmin = MIN (pDmin, points[k].pD);
     38    pDmax = MAX (pDmax, points[k].pD);
     39  }
     40  *Trange = Tmax - Tmin;
     41
     42  // mean epoch
     43  *Tmean = Tsum / Wsum;
     44
     45  // for HIGH_SPEED, just use the center of the range
     46  if (RELASTRO_OP == OP_HIGH_SPEED) {
     47    *Tmean = 0.5*(Tmax - Tmin);
     48  }
     49
     50  *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
     51
     52  /* we need to do the fit in a locally linear space; choose a ref coordinate */
     53  fitStats->coords.crval1 = points[0].R;
     54  fitStats->coords.crval2 = points[0].D;
     55
     56  // project all of the R,D coordinates to a plane centered on this coordinate. set
     57  // the times to be relative to Tmean
     58  for (k = 0; k < Npoints; k++) {
     59    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
     60    points[k].T -= *Tmean;
     61  }       
     62  return TRUE;
     63}
    464
    565FitStats *FitStatsInit (int Nmax, int Nboot) {
    … …  
    2787  fitStats->points = NULL;
    2888  fitStats->sample = NULL;
     89  fitStats->nomask = NULL;
    2990  fitStats->Npoints = 0;
    3091  fitStats->NpointsAlloc = Nmax;
    … …  
    3293    ALLOCATE (fitStats->points, FitAstromPoint, Nmax);
    3394    ALLOCATE (fitStats->sample, FitAstromPoint, Nmax);
     95    ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax);
    3496  }
    3597
    … …  
    79141  FREE (fitStats->points);
    80142  FREE (fitStats->sample);
     143  FREE (fitStats->nomask);
    81144
    82145  FitAstromDataFree (fitStats->fitdataPos);
    … …  
    93156
    94157  /* do I need to do this as 2 2x2 matrix equations? */
    95   fit->A = array_init (Nterms, Nterms);
    96   fit->B = array_init (Nterms, 1);
     158  fit->B   = array_init (Nterms, 1);
     159  fit->A   = array_init (Nterms, Nterms);
     160  fit->Cov = array_init (Nterms, Nterms);
     161
     162  ALLOCATE (fit->Beta, double, Nterms);
     163  ALLOCATE (fit->Beta_prev, double, Nterms);
    97164  fit->Nterms = Nterms;
     165
     166  fit->getChisq = TRUE;
     167  fit->getError = TRUE;
    98168
    99169  return fit;
    … …  
    106176  array_free (fit->A, fit->Nterms);
    107177  array_free (fit->B, fit->Nterms);
     178  array_free (fit->Cov, fit->Nterms);
     179
     180  free (fit->Beta);
     181  free (fit->Beta_prev);
     182
    108183  free (fit);
    109184  return;
    … …  
    126201  object->C_red  = 0.0;
    127202  object->measure= -1;
     203
     204  object->Wx     = 1.0;
     205  object->Wy     = 1.0;
     206
     207  object->rx     = 0.0;
     208  object->ry     = 0.0;
     209  object->u      = 0.0;
     210
     211  object->mask   = 0; // keep point if mask == 0
    128212  return;
    129213}
    … …  
    147231  return;
    148232}
     233
     234double weight_cauchy (double x) {
     235  double r = x / 2.385;
     236  return (1.0 / (1.0 + SQ(r)));
     237}
     238
     239// dpsi = (d/dx) (x * weight(x))
     240double dpsi_cauchy (double x) {
     241  double r2 = SQ(x / 2.385);
     242  return ((1.0 - r2) / (SQ(1 + r2)));
     243}
     244
     245
     246// median absolute deviation
     247// MAD = median(abs(x - median(x)))
     248double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
     249
     250  double *x;
     251  double median = 0.0;
     252  int i;
     253 
     254  ALLOCATE(x, double, Npoints);
     255  for (i = 0; i < Npoints; i++) {
     256    x[i] = points[i].u;
     257  }
     258  dsort(x, Npoints);
     259
     260  if (Npoints % 2) {
     261    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
     262  } else {
     263    median = x[(int)(0.5*Npoints)];
     264  }
     265
     266  for (i = 0; i < Npoints; i++ ) {
     267    x[i] = fabs(x[i] - median);
     268  }
     269  dsort(x, Npoints);
     270
     271  if (Npoints % 2) {
     272    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
     273  } else {
     274    median = x[(int)(0.5*Npoints)];
     275  }
     276
     277  return median;
     278}

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar.c

    r38986 r39266  
    1515
    1616  for (i = 0; i < Npoints; i++) {
     17    if (points[i].mask) continue; // respect the mask if set
     18    fit->Nfit ++;
     19
    1720    /* handle case where dX or dY = 0.0 */
    1821    wx = 1.0 / SQ(points[i].dX);
    … …  
    8992  fit->dp  = sqrt(data->A[4][4]);
    9093 
    91   fit->Nfit = Npoints;
    92 
    9394  return (TRUE);
    9495}

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/ImageOps.c

    r38986 r39266  
    769769      ref[i].R += dTime * catalog[c].average[n].uR / 3600.0 / cos(ref[i].D*RAD_DEG);
    770770      ref[i].D += dTime * catalog[c].average[n].uD / 3600.0;
    771     }
    772 
    773     // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position
    774     if (0) {
    775       myAssert (!isnan(measure[0].RoffGAL), "oops");
    776       myAssert (!isnan(measure[0].DoffGAL), "oops");
    777       ref[i].R += measure[0].RoffGAL / 3600.0;
    778       ref[i].D += measure[0].DoffGAL / 3600.0;
    779771    }
    780772

  • branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/UpdateObjects.c

    r38986 r39266  
    140140  // use a fit with fewer parameters.
    141141
     142  // fit the parallax + proper-motion model
     143  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
     144  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
     145  if (mode == FIT_PM_AND_PAR) {
     146    if (Trange < PM_DT_MIN) {
     147      mode = FIT_AVERAGE;
     148      goto justPosition;
     149    }
     150    if (parRange < PAR_FACTOR_MIN) {
     151      mode = FIT_PM_ONLY;
     152      goto skipParallax;
     153    }
     154    if (fitStats->Npoints <= PAR_TOOFEW) {
     155      mode = FIT_PM_ONLY;
     156      goto skipParallax;
     157    }
     158
     159    // we are going to use the IRLS analysis to calculate the mean solution and the masking
     160    // then run N_BOOTSTRAP_SAMPLES to measure the errors
     161    fitStats->fitdataPar->getError = (N_BOOTSTRAP_SAMPLES < 3);
     162    FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, VERBOSE);
     163
     164    if (N_BOOTSTRAP_SAMPLES >= 3) {
     165      fitStats->Nfit = 0;
     166      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     167      if (Nnomask < 5) {
     168        myAbort ("handle this case, please");
     169      }
     170      for (k = 0; k < fitStats->NfitAlloc; k++) {
     171        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     172        FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
     173        fitStats->Nfit ++;
     174      }
     175      // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     176      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     177      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     178      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     179      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     180      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
     181    }
     182
     183    // project Ro, Do back to RA,DEC
     184    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
     185    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     186
     187    average[0].flags |= ID_STAR_FIT_PAR;
     188    fitStats->Npar ++;
     189
     190    // XXX a hard-wired hack...
     191    // unless there is a clear problems (below) with the parallax fit, we will use it
     192    if ((fabs(fitPar.uR) > 4.0) || (fabs(fitPar.uD) > 4.0)) {
     193      mode = FIT_PM_ONLY;
     194    } else {
     195      average[0].flags |= ID_STAR_USE_PAR;
     196    }
     197  }       
     198
     199skipParallax:
     200
    142201  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
    143202  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
    … …  
    151210    }
    152211
    153     if (fitStats->NfitAlloc == 1) {
    154       // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
    155       FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
     212    if (average[0].flags & ID_STAR_USE_PAR) {
     213      if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
     214        myAbort ("oops");
     215      }
    156216    } else {
    157       fitStats->Nfit = 0;
    158       for (k = 0; k < fitStats->NfitAlloc; k++) {
    159         BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    160         if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
    161         fitStats->Nfit ++;
    162       }
    163       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    164       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    165       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    166       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    167     }
    168     FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
     217      fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3);
     218      FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
     219
     220      if (N_BOOTSTRAP_SAMPLES >= 3) {
     221        fitStats->Nfit = 0;
     222        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     223        if (Nnomask < 5) {
     224          myAbort ("handle this case, please");
     225        }
     226        for (k = 0; k < fitStats->NfitAlloc; k++) {
     227          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     228          if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
     229          fitStats->Nfit ++;
     230        }
     231        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     232        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     233        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     234        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     235      }
     236    }
    169237
    170238    // project Ro, Do back to RA,DEC
    … …  
    172240    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    173241
    174     fitPM.p  = fitPM.dp  = 0.0;
    175242    average[0].flags |= ID_STAR_FIT_PM;
    176243    fitStats->Npm ++;
    177244
    178245    // XXX a hard-wired hack...
    179     if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
     246    // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
     247    if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) {
    180248      mode = FIT_AVERAGE;
    181249      average[0].flags |= ID_STAR_BAD_PM;
     250    } else {
     251      if (!(average[0].flags & ID_STAR_USE_PAR)) {
     252        average[0].flags |= ID_STAR_USE_PM;
     253      }
    182254    }
    183255  }
    184256 
    185   // fit the parallax + proper-motion model
    186   // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
    187   // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
    188   if (mode == FIT_PM_AND_PAR) {
    189     if (parRange < PAR_FACTOR_MIN) {
    190       mode = FIT_PM_ONLY;
    191       goto justPosition;
    192     }
    193     if (fitStats->Npoints <= PAR_TOOFEW) {
    194       mode = FIT_PM_ONLY;
    195       goto justPosition;
    196     }
    197 
    198     if (fitStats->NfitAlloc == 1) {
    199       // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
    200       FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
    201     } else {
    202       fitStats->Nfit = 0;
    203       for (k = 0; k < fitStats->NfitAlloc; k++) {
    204         BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    205         FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
    206         fitStats->Nfit ++;
    207       }
    208       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    209       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    210       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    211       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    212       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
    213     }
    214     FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
    215 
    216     // project Ro, Do back to RA,DEC
    217     XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
    218     if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    219 
    220     average[0].flags |= ID_STAR_FIT_PAR;
    221     fitStats->Npar ++;
    222 
    223     // XXX a hard-wired hack...
    224     if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) {
    225       mode = FIT_PM_ONLY;
    226     }
    227   }       
    228 
    229257justPosition:
    230258  {
    … …  
    232260    // if we only have one point, this is silly...
    233261   
     262    // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway?
     263
    234264    if (fitStats->NfitAlloc == 1) {
    235265      FitAstromResultSetPM (&fitPos, 1, average);
    … …  
    618648}
    619649
    620 int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
    621 
    622   int k;
    623 
    624   int Npoints = fitStats->Npoints;
    625   FitAstromPoint *points = fitStats->points;
    626 
    627   // find Tmin & Tmax from the list of accepted measurements
    628   double Tmin  = points[0].T;
    629   double Tmax  = points[0].T;
    630   double pRmin = +2.0;
    631   double pRmax = -2.0;
    632   double pDmin = +2.0;
    633   double pDmax = -2.0;
    634 
    635   *Tmean = 0.0;
    636 
    637   double Tsum = 0.0;
    638   double Wsum = 0.0;
    639   for (k = 0; k < Npoints; k++) {
    640     Tmin = MIN(Tmin, points[k].T);
    641     Tmax = MAX(Tmax, points[k].T);
    642 
    643     float wx = 1.0 / SQ(points[k].dX);
    644 
    645     Tsum += points[k].T * wx;
    646     Wsum += wx;
    647 
    648     // at this point, T is in years since J2000
    649     ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
    650     pRmin = MIN (pRmin, points[k].pR);
    651     pRmax = MAX (pRmax, points[k].pR);
    652     pDmin = MIN (pDmin, points[k].pD);
    653     pDmax = MAX (pDmax, points[k].pD);
    654   }
    655   *Trange = Tmax - Tmin;
    656 
    657   // mean epoch
    658   *Tmean = Tsum / Wsum;
    659 
    660   // for HIGH_SPEED, just use the center of the range
    661   if (RELASTRO_OP == OP_HIGH_SPEED) {
    662     *Tmean = 0.5*(Tmax - Tmin);
    663   }
    664 
    665   *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
    666 
    667   /* we need to do the fit in a locally linear space; choose a ref coordinate */
    668   fitStats->coords.crval1 = points[0].R;
    669   fitStats->coords.crval2 = points[0].D;
    670 
    671   // project all of the R,D coordinates to a plane centered on this coordinate. set
    672   // the times to be relative to Tmean
    673   for (k = 0; k < Npoints; k++) {
    674     RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
    675     points[k].T -= *Tmean;
    676   }       
    677   return TRUE;
    678 }
    679 
    680650int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
    681651
    … …  
    689659  }
    690660  return Nmax;
    691 }
    692 
    693 int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) {
    694   int i;
    695 
    696   // I need to draw Npoints random entries from 'points' with replacement:
    697   for (i = 0; i < Npoints; i++) {
    698     int N = Npoints * drand48();
    699     sample[i] = points[N];
    700   }
    701   return TRUE;
    702 }
    703 
    704 // calculate mean and sigma points for the 5 fit parameter
    705 int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
    706 
    707   // generate a histogram for the selected element
    708   double *values = NULL;
    709   ALLOCATE (values, double, Nfit);
    710  
    711   int i;
    712 
    713   for (i = 0; i < Nfit; i++) {
    714     switch (mode) {
    715       case FIT_RESULT_RA:
    716         values[i] = fit[i].Ro;
    717         break;
    718       case FIT_RESULT_DEC:
    719         values[i] = fit[i].Do;
    720         break;
    721       case FIT_RESULT_uR:
    722         values[i] = fit[i].uR;
    723         break;
    724       case FIT_RESULT_uD:
    725         values[i] = fit[i].uD;
    726         break;
    727       case FIT_RESULT_PLX:
    728         values[i] = fit[i].p;
    729         break;
    730       default:
    731         myAbort ("invalid option");
    732     }
    733   }
    734 
    735   dsort (values, Nfit);
    736 
    737   double median;
    738   if (Nfit % 2) {
    739     int Ncenter = Nfit / 2;
    740     median = values[Ncenter];
    741   } else {
    742     int Ncenter = Nfit / 2 - 1;
    743     median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
    744   }
    745 
    746   double Slo = VectorFractionInterpolate (values, 0.158655, Nfit);
    747   double Shi = VectorFractionInterpolate (values, 0.841345, Nfit);
    748   double sigma = (Shi - Slo) / 2.0;
    749 
    750   switch (mode) {
    751     case FIT_RESULT_RA:
    752       result->Ro = median;
    753       result->dRo = sigma;
    754       break;
    755     case FIT_RESULT_DEC:
    756       result->Do = median;
    757       result->dDo = sigma;
    758       break;
    759     case FIT_RESULT_uR:
    760       result->uR = median;
    761       result->duR = sigma;
    762       break;
    763     case FIT_RESULT_uD:
    764       result->uD = median;
    765       result->duD = sigma;
    766       break;
    767     case FIT_RESULT_PLX:
    768       result->p = median;
    769       result->dp = sigma;
    770       break;
    771     default:
    772       myAbort ("invalid option");
    773   }
    774 
    775   return TRUE;
    776 }
    777 
    778 double VectorFractionInterpolate (double *values, float fraction, int Npts) {
    779 
    780   float F = fraction * Npts;
    781   int   N = fraction * Npts;
    782 
    783   if (N < 0        ) return NAN;
    784   if (N >= Npts - 2) return NAN;
    785 
    786   // interpolate between N,N+1
    787    
    788   double S = (F - N) * (values[N+1] - values[N]) + values[N];
    789   return S;
    790661}
    791662

  • branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/BrightCatalog.c

    r37037 r39266  
    6767    GET_COLUMN(M,         "MAG_SYS",    float);
    6868    GET_COLUMN(Mcal,      "MAG_CAL",    float);
     69    GET_COLUMN(Mflat,     "MAG_FLAT",   float);
    6970    GET_COLUMN(dM,        "MAG_ERR",    float);
    7071    GET_COLUMN(airmass,   "AIRMASS",    float);
    … …  
    8990      measure[i].M         = M[i];
    9091      measure[i].Mcal      = Mcal[i];
     92      measure[i].Mflat     = Mflat[i];
    9193      measure[i].dM        = dM[i];
    9294      measure[i].airmass   = airmass[i];
    … …  
    109111    free (M       );
    110112    free (Mcal    );
     113    free (Mflat   );
    111114    free (dM      );
    112115    free (airmass );
    … …  
    269272    gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
    270273    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
     274    gfits_define_bintable_column (&theader, "E", "MAG_FLAT", "magnitude (flat)",            NULL,    1.0, 0.0);
    271275    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "magnitude (err)",             NULL,    1.0, 0.0);
    272276    gfits_define_bintable_column (&theader, "E", "AIRMASS",  "airmass",                     NULL,    1.0, 0.0);
    … …  
    290294    float  *M         ; ALLOCATE (M        ,  float,  catalog->Nmeasure);
    291295    float  *Mcal      ; ALLOCATE (Mcal     ,  float,  catalog->Nmeasure);
     296    float  *Mflat     ; ALLOCATE (Mflat    ,  float,  catalog->Nmeasure);
    292297    float  *dM        ; ALLOCATE (dM       ,  float,  catalog->Nmeasure);
    293298    float  *airmass   ; ALLOCATE (airmass  ,  float,  catalog->Nmeasure);
    … …  
    310315      M[i]        = measure[i].M        ;
    311316      Mcal[i]     = measure[i].Mcal     ;
     317      Mflat[i]    = measure[i].Mflat    ;
    312318      dM[i]       = measure[i].dM       ;
    313319      airmass[i]  = measure[i].airmass  ;
    … …  
    329335    gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
    330336    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
     337    gfits_set_bintable_column (&theader, &ftable, "MAG_FLAT",   Mflat,     catalog->Nmeasure);
    331338    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",    dM,        catalog->Nmeasure);
    332339    gfits_set_bintable_column (&theader, &ftable, "AIRMASS",    airmass,   catalog->Nmeasure);
    … …  
    346353    free (M       );
    347354    free (Mcal    );
     355    free (Mflat   );
    348356    free (dM      );
    349357    free (airmass );

  • branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/setMrelCatalog.c

    r39107 r39266  
    11# include "relphot.h"
     2
     3/*
     4# define TEST_OBJ_ID 0x3ae2
     5# define TEST_CAT_ID 0x16c5f
     6*/
     7# define TEST_OBJ_ID 0
     8# define TEST_CAT_ID 0
     9
    210int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
    311
    … …  
    113121
    114122  // option for a test print
    115   if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) {
     123  if ((averageT[0].objID == TEST_OBJ_ID) && (averageT[0].catID == TEST_CAT_ID)) {
    116124    fprintf (stderr, "test obj\n");
    117125    print_measure_set_alt (average, secfilt, measure);
    … …  
    353361
    354362    if (isSetMrelFinal) {
    355       if (FALSE && (average->objID == 0x1) && (average->catID)) {
     363      if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
    356364        fprintf (stderr, "test obj\n");
    357365      }
    … …  
    723731
    724732  // option for a test print
    725   if (FALSE && (average[0].objID == 0x0000049c) && (average[0].catID == 0x00016c5e)) {
     733  if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
    726734    fprintf (stderr, "test obj\n");
    727735    print_measure_set_alt (average, secfilt, measure);

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/Makefile

    r38986 r39266  
    1 default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client
     1# default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client
     2default: uniphot setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client
     3
    24help:
    35@echo "make options: uniphot setphot setfwhm default help install default (uniphot setphot setfwhm)"
    … …  
    3436ckids_client: $(BIN)/ckids_client.$(ARCH)
    3537
    36 install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/setgalmodel $(DESTBIN)/setgalmodel_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client $(DESTBIN)/fixstkids $(DESTBIN)/fixstkids_client $(DESTBIN)/ckids $(DESTBIN)/ckids_client
     38# install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/setgalmodel $(DESTBIN)/setgalmodel_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client $(DESTBIN)/fixstkids $(DESTBIN)/fixstkids_client $(DESTBIN)/ckids $(DESTBIN)/ckids_client
     39install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/setgalmodel $(DESTBIN)/setgalmodel_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client $(DESTBIN)/fixstkids $(DESTBIN)/fixstkids_client $(DESTBIN)/ckids $(DESTBIN)/ckids_client
    3740
    3841UNIPHOT = \

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/include/setphot.h

    r39219 r39266  
    1717  unsigned int flags;
    1818} ImageSubset;
     19
     20// we have one correction (an image) for each filter and chip
     21typedef struct {
     22  int Nx;       // number of chips in x
     23  int Ny;       // number of chips in y
     24  int Nfilter;  // number of filters
     25  int Nseason;  // number of correction peridos
     26
     27  int Nchips;   // chip offset (Nx*Ny)
     28  int Nflats;   // season offset (Nx*Ny*Nfilters)
     29  int Nvalues;  // Nx*Ny*Nfilters*Nseason
     30
     31  int dX;       // superpixel size
     32  int dY;       // superpixel size
     33
     34  int NxCCD;    // number of pixels
     35  int NyCCD;    // number of pixels
     36
     37  Matrix **matrix; // allocate an array of pointers
     38  // index = ix + iy*Nx + filter*Nchips + dir*Ngroup
     39} CamCorrection;
    1940
    2041/* global variables set in parameter file */

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_correction.c

    r38062 r39266  
    2929
    3030  int Ndir, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD;
    31   if (!gfits_scan (&header, "NDIR",    "%d", 1, &Ndir))    return FALSE;
    32   if (!gfits_scan (&header, "NFILTER", "%d", 1, &Nfilter)) return FALSE;
     31  if (!gfits_scan (&header, "NDIR",     "%d", 1, &Ndir))    return FALSE;
     32  if (!gfits_scan (&header, "NFILTER",  "%d", 1, &Nfilter)) return FALSE;
    3333  if (!gfits_scan (&header, "NX",       "%d", 1, &Nx))      return FALSE;
    3434  if (!gfits_scan (&header, "NY",       "%d", 1, &Ny))      return FALSE;
    … …  
    6161  }
    6262
    63   int found = FALSE;
    6463  while (TRUE) {
    6564    Header theader;
    … …  
    9089    // assert that cam->matrix[index] is NULL?
    9190    cam->matrix[index] = matrix;
    92 
    93     if (!found) {
    94      
    95     }
    9691  }
    9792  return TRUE;

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot.c

    r39219 r39266  
    3333
    3434  int N;
     35
     36  CAM_FILE = NULL;
     37  if ((N = get_argument (argc, argv, "-CAM"))) {
     38    remove_argument (N, &argc, argv);
     39    char *tmpfile = strcreate (argv[N]);
     40    CAM_FILE = abspath (tmpfile, DVO_MAX_PATH);
     41    remove_argument (N, &argc, argv);
     42  }
     43  CAM_RESET = FALSE;
     44  if ((N = get_argument (argc, argv, "-CAM-reset"))) {
     45    remove_argument (N, &argc, argv);
     46    CAM_RESET = TRUE;
     47  }
    3548
    3649  VERBOSE = FALSE;

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot_client.c

    r38986 r39266  
    2626
    2727  int N;
     28
     29  CAM_FILE = NULL;
     30  if ((N = get_argument (argc, argv, "-CAM"))) {
     31    remove_argument (N, &argc, argv);
     32    char *tmpfile = strcreate (argv[N]);
     33    CAM_FILE = abspath (tmpfile, DVO_MAX_PATH);
     34    remove_argument (N, &argc, argv);
     35  }
     36  CAM_RESET = FALSE;
     37  if ((N = get_argument (argc, argv, "-CAM-reset"))) {
     38    remove_argument (N, &argc, argv);
     39    CAM_RESET = TRUE;
     40  }
    2841
    2942  VERBOSE = FALSE;

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/setphot.c

    r38986 r39266  
    88  ZptTable *zpts;
    99  Image *image;
    10   FlatCorrectionTable flatcorrTable;
     10
     11  // FlatCorrectionTable flatcorrTable;
    1112
    1213  /* get configuration info, args, lockfile */
    … …  
    3334  }
    3435
     36  // if supplied, read in the camera-static flat-field correction
     37  if (CAM_FILE) {
     38    if (!load_cam_correction (CAM_FILE)) {
     39      fprintf (stderr, "failed to load camera-static flat-field correction\n");
     40      exit (2);
     41    }
     42    merge_flatcorr_with_cam (flatcorrTable);
     43  } else {
     44    convert_flatcorr_to_cam (flatcorrTable);
     45  }
     46
     47  {
     48    char newflatfile[DVO_MAX_PATH];
     49    int size = snprintf (newflatfile, DVO_MAX_PATH, "%s/flatfield.fits", CATDIR);
     50    assert (size < DVO_MAX_PATH);
     51    save_cam_correction (newflatfile);
     52  }
     53
    3554  if (!zpts) Shutdown ("failed to load zero points, or empty table");
    3655
    … …  
    4261
    4362  if (UBERCAL) {
    44     match_flatcorr_to_images (image, Nimage, &flatcorrTable);
     63    // we are going to deprecate the flatcorr imaage lookup
     64    // match_flatcorr_to_images (image, Nimage, &flatcorrTable);
     65    match_camcorr_to_images (image, Nimage);
    4566  }
    4667
    4768  if (!IMAGES_ONLY) {
    48     status = update_dvo_setphot (image, Nimage, &flatcorrTable);
     69    status = update_dvo_setphot (image, Nimage);
    4970  }
    5071

  • branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/update_catalog_setphot.c

    r35425 r39266  
    11# include "setphot.h"
    22
    3 void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) {
     3void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage) {
    44
    55  off_t i, j, found;
    … …  
    77  // if we are resetting, reset all flags
    88  DVOMeasureFlags photomFlags = ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_PHOTOM_UBERCAL;
     9
     10  CamCorrection *camcorr = get_cam_correction_ptr (); 
    911
    1012  found = 0;   
    … …  
    2123      float Mcal = image[id].Mcal;
    2224      float dMcal = image[id].dMcal;
    23       float Mcal_offset = 0.0;
     25      float Mflat = 0.0;
    2426
    2527      // if we know about a flat-field correction, then we need to apply the sub-chip correction
    2628      int flat_id = image[id].photom_map_id;
    2729      if (flat_id > 0) {
    28           Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);
     30        Mflat = CamCorrectionOffset (camcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);
    2931      }
    3032
    31       catalog[0].measure[m].Mcal = Mcal - Mcal_offset;
     33      catalog[0].measure[m].Mcal = Mcal;
     34      catalog[0].measure[m].Mflat = Mflat;
    3235      catalog[0].measure[m].dMcal = dMcal;
    3336      myAssert(isfinite(catalog[0].measure[m].Mcal), "oops: ubercal made a nan");
Note: See TracChangeset for help on using the changeset viewer.