IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38938


Ignore:
Timestamp:
Oct 22, 2015, 5:53:15 AM (11 years ago)
Author:
eugene
Message:

update mkcmf to make the ipptopsps test suite a bit clearer: add filter, randomize zero point, apply x,y,m random noise before setting output fields

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r38920 r38938  
    2525void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    2626
    27 int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
     27int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime);
    2828int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    2929int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    … …  
    334334    gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
    335335  }
    336   gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT);
     336  gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT); // this is now not used
    337337  gfits_modify (&header, "EXPTIME", "%lf", 1, exptime);
    338338  gfits_modify (&header, "AIRMASS", "%lf", 1, airmass);
    … …  
    370370  gfits_modify (&header, "AST_CDY",           "%f", 1, 0.04);
    371371
     372  float zeroPt = 25.0 + 0.3*(drand48() - 0.5);
     373
    372374  if (!isStack && !isForcedWarp) {
    373375    gfits_modify (&header, "HIERARCH DETREND.MASK",      "%s", 1, "detref615.XY33.fits");
    … …  
    377379    gfits_modify (&header, "HIERARCH DETREND.NONLIN",    "%s", 1, "linearity_data.XY33.fits");
    378380    gfits_modify (&header, "HIERARCH DETREND.VIDEODARK", "%s", 1, "GPC1.VIDEODARK.979.0.XY33.fits");
    379     gfits_modify (&header, "ZPT_OBS",           "%f", 1, 25.0);
     381    gfits_modify (&header, "ZPT_OBS",                    "%f", 1, zeroPt);
    380382  }
    381383
    382384  if (isStack) {
    383     gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, 25.0);
     385    gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt);
    384386    gfits_modify (&header, "TESS_ID",         "%s", 1, tess_id);
    385387    gfits_modify (&header, "SKYCELL",         "%s", 1, skycell);
     388    gfits_modify (&header, "ZPT_OBS",         "%f", 1, zeroPt);
    386389  }
    387390  if (isForcedWarp) {
    388     gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, "r.00000");
    389     gfits_modify (&header, "HIERARCH FPA.ZP",       "%f", 1, 25.0);
     391    gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, "r.00000"); // does this affect anything? I don't actually think so...
     392    gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt);
    390393    gfits_modify (&header, "TESS_ID",         "%s", 1, tess_id);
    391394    gfits_modify (&header, "SKYCELL",         "%s", 1, skycell);
    … …  
    472475  gfits_fwrite_table   (fits, &ftable);
    473476
    474   if (WriteXSRC) WriteXSRCtable (fits, extroot, X, Y, M, Flag, Nstars);
     477  if (WriteXSRC) WriteXSRCtable (fits, extroot, X, Y, M, Flag, Nstars, zeroPt, exptime);
    475478  if (WriteXFIT) WriteXFITtable (fits, extroot, X, Y, M, Flag, Nstars);
    476479  if (WriteXRAD) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars);
    … …  
    572575    fSN = 1.0 / sqrt(flux);
    573576
     577    if (ADDNOISE) {
     578      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     579      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     580      M[i] += fSN*rnd_gauss(0.0, 1.0);
     581      flux = pow (10.0, -0.4*M[i]);
     582      fSN = 1.0 / sqrt(flux);
     583    }
     584
    574585    stars[i].X = X[i];
    575586    stars[i].Y = Y[i];
    576587    stars[i].M = M[i];
    577 
    578     if (ADDNOISE) {
    579       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    580       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    581       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    582     }
    583588
    584589    stars[i].dX = FX * fSN;
    … …  
    616621    fSN = 1.0 / sqrt(flux);
    617622
     623    if (ADDNOISE) {
     624      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     625      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     626      M[i] += fSN*rnd_gauss(0.0, 1.0);
     627      flux = pow (10.0, -0.4*M[i]);
     628      fSN = 1.0 / sqrt(flux);
     629    }
     630
    618631    stars[i].X = X[i];
    619632    stars[i].Y = Y[i];
    620633    stars[i].M = M[i];
    621 
    622     if (ADDNOISE) {
    623       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    624       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    625       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    626     }
    627634
    628635    stars[i].dX = FX * fSN;
    … …  
    663670    fSN = 1.0 / sqrt(flux);
    664671
     672    if (ADDNOISE) {
     673      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     674      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     675      M[i] += fSN*rnd_gauss(0.0, 1.0);
     676      flux = pow (10.0, -0.4*M[i]);
     677      fSN = 1.0 / sqrt(flux);
     678    }
     679
    665680    stars[i].X = X[i];
    666681    stars[i].Y = Y[i];
    667682    stars[i].M = M[i];
    668 
    669     if (ADDNOISE) {
    670       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    671       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    672       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    673     }
    674683
    675684    stars[i].dX = FX * fSN;
    … …  
    712721    fSN = 1.0 / sqrt(flux);
    713722
     723    if (ADDNOISE) {
     724      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     725      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     726      M[i] += fSN*rnd_gauss(0.0, 1.0);
     727      flux = pow (10.0, -0.4*M[i]);
     728      fSN = 1.0 / sqrt(flux);
     729    }
     730
    714731    stars[i].X = X[i];
    715732    stars[i].Y = Y[i];
    716733    stars[i].M = M[i];
    717 
    718     if (ADDNOISE) {
    719       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    720       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    721       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    722     }
    723734
    724735    // randomly give poor PSFQF values
    … …  
    767778    fSN = 1.0 / sqrt(flux);
    768779
     780    if (ADDNOISE) {
     781      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     782      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     783      M[i] += fSN*rnd_gauss(0.0, 1.0);
     784      flux = pow (10.0, -0.4*M[i]);
     785      fSN = 1.0 / sqrt(flux);
     786    }
     787
    769788    stars[i].X = X[i];
    770789    stars[i].Y = Y[i];
    771790    stars[i].M = M[i];
    772791    stars[i].Map = M[i] - 0.05;
    773 
    774     if (ADDNOISE) {
    775       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    776       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    777       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    778     }
    779792
    780793    // randomly give poor PSFQF values
    … …  
    826839    fSN = 1.0 / sqrt(flux);
    827840
     841    if (ADDNOISE) {
     842      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     843      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     844      M[i] += fSN*rnd_gauss(0.0, 1.0);
     845      flux = pow (10.0, -0.4*M[i]);
     846      fSN = 1.0 / sqrt(flux);
     847    }
     848
    828849    stars[i].X = X[i];
    829850    stars[i].Y = Y[i];
    … …  
    837858    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
    838859    stars[i].apFluxErr = stars[i].apFlux * fSN;
    839 
    840     if (ADDNOISE) {
    841       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    842       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    843       stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    844     }
    845860
    846861    // randomly give poor PSFQF values
    … …  
    893908    fSN = 1.0 / sqrt(flux);
    894909
     910    if (ADDNOISE) {
     911      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     912      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     913      M[i] += fSN*rnd_gauss(0.0, 1.0);
     914      flux = pow (10.0, -0.4*M[i]);
     915      fSN = 1.0 / sqrt(flux);
     916    }
     917
    895918    stars[i].detID = i;
    896919    stars[i].X = X[i];
    … …  
    972995
    973996    stars[i].nFrames   = 1;
    974 
    975     if (ADDNOISE) {
    976       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    977       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    978       float Moff = fSN*rnd_gauss(0.0, 1.0);
    979       stars[i].M += Moff;
    980       stars[i].Map += Moff;
    981       stars[i].MapRaw += Moff;
    982       stars[i].Mcalib += Moff;
    983     }
    984997  }
    985998
    … …  
    10011014    flux = pow (10.0, -0.4*M[i]);
    10021015    fSN = 1.0 / sqrt(flux);
     1016
     1017    if (ADDNOISE) {
     1018      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     1019      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     1020      M[i] += fSN*rnd_gauss(0.0, 1.0);
     1021      flux = pow (10.0, -0.4*M[i]);
     1022      fSN = 1.0 / sqrt(flux);
     1023    }
    10031024
    10041025    stars[i].detID = i;
    … …  
    10811102
    10821103    stars[i].nFrames   = 1;
    1083 
    1084     if (ADDNOISE) {
    1085       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    1086       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    1087       float Moff = fSN*rnd_gauss(0.0, 1.0);
    1088       stars[i].M += Moff;
    1089       stars[i].Map += Moff;
    1090       stars[i].MapRaw += Moff;
    1091       stars[i].Mcalib += Moff;
    1092     }
    10931104  }
    10941105
    … …  
    11181129    flux = pow (10.0, -0.4*M[i]);
    11191130    fSN = 1.0 / sqrt(flux);
     1131
     1132    if (ADDNOISE) {
     1133      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
     1134      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
     1135      M[i] += fSN*rnd_gauss(0.0, 1.0);
     1136      flux = pow (10.0, -0.4*M[i]);
     1137      fSN = 1.0 / sqrt(flux);
     1138    }
    11201139
    11211140    stars[i].detID = i;
    … …  
    12211240
    12221241    stars[i].nFrames   = 1;
    1223 
    1224     if (ADDNOISE) {
    1225       stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
    1226       stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
    1227       float Moff = fSN*rnd_gauss(0.0, 1.0);
    1228       stars[i].M += Moff;
    1229       stars[i].Map += Moff;
    1230       stars[i].MapRaw += Moff;
    1231       stars[i].Mcalib += Moff;
    1232     }
    12331242  }
    12341243
    … …  
    12371246}
    12381247
    1239 int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
     1248int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime) {
    12401249
    12411250  int i, j;
    … …  
    13061315    float  *PROF_FILL           ; ALLOCATE (PROF_FILL           ,  float , 9*Nstars);
    13071316
     1317    float magtime = 2.5*log10(exptime);
     1318
    13081319    // assign the storage arrays
    13091320    for (i = 0; i < Nstars; i++) {
    … …  
    13181329      F25_ARATIO          [i] = 0.9;
    13191330      F25_THETA           [i] = 10.0;
    1320       PETRO_MAG           [i] = M[i];
    1321       PETRO_MAG_ERR       [i] = fSN;
     1331      PETRO_MAG           [i] = M[i] + zeroPt + magtime;
     1332      PETRO_MAG_ERR       [i] = 1.0 / fSN;
    13221333      PETRO_RADIUS        [i] = 8.0;
    13231334      if (i == 2) {
    … …  
    15281539int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
    15291540
    1530   int i, j;
     1541  int i, j, k;
    15311542  Header header;
    15321543  FTable ftable;
    … …  
    15531564
    15541565    // create intermediate storage arrays
    1555     int    *IPP_IDET           ; ALLOCATE (IPP_IDET         ,  int   , Nstars);
    1556     float  *X_APER             ; ALLOCATE (X_APER           ,  float , Nstars);
    1557     float  *Y_APER             ; ALLOCATE (Y_APER           ,  float , Nstars);
    1558     float  *PSF_FWHM           ; ALLOCATE (PSF_FWHM         ,  float , Nstars);
    1559     float  *APER_FLUX          ; ALLOCATE (APER_FLUX        ,  float , Nstars*9);
    1560     float  *APER_FLUX_ERR      ; ALLOCATE (APER_FLUX_ERR    ,  float , Nstars*9);
    1561     float  *APER_FLUX_STDEV    ; ALLOCATE (APER_FLUX_STDEV  ,  float , Nstars*9);
    1562     float  *APER_FILL          ; ALLOCATE (APER_FILL        ,  float , Nstars*9);
     1566    int    *IPP_IDET           ; ALLOCATE (IPP_IDET         ,  int   , 3*Nstars);
     1567    float  *X_APER             ; ALLOCATE (X_APER           ,  float , 3*Nstars);
     1568    float  *Y_APER             ; ALLOCATE (Y_APER           ,  float , 3*Nstars);
     1569    float  *PSF_FWHM           ; ALLOCATE (PSF_FWHM         ,  float , 3*Nstars);
     1570    float  *APER_FLUX          ; ALLOCATE (APER_FLUX        ,  float , 3*Nstars*9);
     1571    float  *APER_FLUX_ERR      ; ALLOCATE (APER_FLUX_ERR    ,  float , 3*Nstars*9);
     1572    float  *APER_FLUX_STDEV    ; ALLOCATE (APER_FLUX_STDEV  ,  float , 3*Nstars*9);
     1573    float  *APER_FILL          ; ALLOCATE (APER_FILL        ,  float , 3*Nstars*9);
     1574
     1575    float psfFWHM[3] = {4.8, 6.0, 8.0};
    15631576
    15641577    // assign the storage arrays
    … …  
    15671580      // float fSN = 1.0 / sqrt(flux);
    15681581
    1569       IPP_IDET         [i] = i;
    1570       X_APER           [i] = X[i];
    1571       Y_APER           [i] = Y[i];
    1572       PSF_FWHM         [i] = FX*2.8;
    1573       for (j = 0; j < 9; j ++) {
    1574         APER_FLUX      [j + 9*i] = flux;
    1575         APER_FLUX_ERR  [j + 9*i] = sqrt(flux);
    1576         APER_FLUX_STDEV[j + 9*i] = sqrt(flux);
    1577         APER_FILL      [j + 9*i] = 0.9;
     1582      for (k = 0; k < 3; k++) {
     1583        int N = 3*i + k;
     1584        IPP_IDET         [N] = i;
     1585        X_APER           [N] = X[i];
     1586        Y_APER           [N] = Y[i];
     1587        PSF_FWHM         [N] = psfFWHM[k];
     1588        for (j = 0; j < 9; j ++) {
     1589          float Flux = flux*(1.0 + 0.02*j - 0.2);
     1590          APER_FLUX      [j + 9*N] = Flux;
     1591          APER_FLUX_ERR  [j + 9*N] = sqrt(Flux);
     1592          APER_FLUX_STDEV[j + 9*N] = sqrt(Flux*1.1);
     1593          APER_FILL      [j + 9*N] = 0.9;
     1594        }
    15781595      }
    15791596    }
    15801597
    15811598    // set the table data
    1582     gfits_set_bintable_column (&header, &ftable, "IPP_IDET"         , IPP_IDET         , Nstars);
    1583     gfits_set_bintable_column (&header, &ftable, "X_APER"           , X_APER           , Nstars);
    1584     gfits_set_bintable_column (&header, &ftable, "Y_APER"           , Y_APER           , Nstars);
    1585     gfits_set_bintable_column (&header, &ftable, "PSF_FWHM"         , PSF_FWHM         , Nstars);
    1586     gfits_set_bintable_column (&header, &ftable, "APER_FLUX"        , APER_FLUX        , Nstars);
    1587     gfits_set_bintable_column (&header, &ftable, "APER_FLUX_ERR"    , APER_FLUX_ERR    , Nstars);
    1588     gfits_set_bintable_column (&header, &ftable, "APER_FLUX_STDEV"  , APER_FLUX_STDEV  , Nstars);
    1589     gfits_set_bintable_column (&header, &ftable, "APER_FILL"        , APER_FILL        , Nstars);
     1599    gfits_set_bintable_column (&header, &ftable, "IPP_IDET"         , IPP_IDET         , 3*Nstars);
     1600    gfits_set_bintable_column (&header, &ftable, "X_APER"           , X_APER           , 3*Nstars);
     1601    gfits_set_bintable_column (&header, &ftable, "Y_APER"           , Y_APER           , 3*Nstars);
     1602    gfits_set_bintable_column (&header, &ftable, "PSF_FWHM"         , PSF_FWHM         , 3*Nstars);
     1603    gfits_set_bintable_column (&header, &ftable, "APER_FLUX"        , APER_FLUX        , 3*Nstars);
     1604    gfits_set_bintable_column (&header, &ftable, "APER_FLUX_ERR"    , APER_FLUX_ERR    , 3*Nstars);
     1605    gfits_set_bintable_column (&header, &ftable, "APER_FLUX_STDEV"  , APER_FLUX_STDEV  , 3*Nstars);
     1606    gfits_set_bintable_column (&header, &ftable, "APER_FILL"        , APER_FILL        , 3*Nstars);
    15901607  }
    15911608  /**/
Note: See TracChangeset for help on using the changeset viewer.