Index: /trunk/Ohana/src/addstar/src/mkcmf.c =================================================================== --- /trunk/Ohana/src/addstar/src/mkcmf.c (revision 38937) +++ /trunk/Ohana/src/addstar/src/mkcmf.c (revision 38938) @@ -25,5 +25,5 @@ void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars); -int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); +int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime); int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); @@ -334,5 +334,5 @@ gfits_modify (&header, "UTC-OBS", "%s", 1, time); } - gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT); + gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT); // this is now not used gfits_modify (&header, "EXPTIME", "%lf", 1, exptime); gfits_modify (&header, "AIRMASS", "%lf", 1, airmass); @@ -370,4 +370,6 @@ gfits_modify (&header, "AST_CDY", "%f", 1, 0.04); + float zeroPt = 25.0 + 0.3*(drand48() - 0.5); + if (!isStack && !isForcedWarp) { gfits_modify (&header, "HIERARCH DETREND.MASK", "%s", 1, "detref615.XY33.fits"); @@ -377,15 +379,16 @@ gfits_modify (&header, "HIERARCH DETREND.NONLIN", "%s", 1, "linearity_data.XY33.fits"); gfits_modify (&header, "HIERARCH DETREND.VIDEODARK", "%s", 1, "GPC1.VIDEODARK.979.0.XY33.fits"); - gfits_modify (&header, "ZPT_OBS", "%f", 1, 25.0); + gfits_modify (&header, "ZPT_OBS", "%f", 1, zeroPt); } if (isStack) { - gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, 25.0); + gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt); gfits_modify (&header, "TESS_ID", "%s", 1, tess_id); gfits_modify (&header, "SKYCELL", "%s", 1, skycell); + gfits_modify (&header, "ZPT_OBS", "%f", 1, zeroPt); } if (isForcedWarp) { - gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, "r.00000"); - gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, 25.0); + gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, "r.00000"); // does this affect anything? I don't actually think so... + gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt); gfits_modify (&header, "TESS_ID", "%s", 1, tess_id); gfits_modify (&header, "SKYCELL", "%s", 1, skycell); @@ -472,5 +475,5 @@ gfits_fwrite_table (fits, &ftable); - if (WriteXSRC) WriteXSRCtable (fits, extroot, X, Y, M, Flag, Nstars); + if (WriteXSRC) WriteXSRCtable (fits, extroot, X, Y, M, Flag, Nstars, zeroPt, exptime); if (WriteXFIT) WriteXFITtable (fits, extroot, X, Y, M, Flag, Nstars); if (WriteXRAD) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars); @@ -572,13 +575,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -616,13 +621,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -663,13 +670,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -712,13 +721,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -767,14 +778,16 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; stars[i].Map = M[i] - 0.05; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -826,4 +839,12 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; @@ -837,10 +858,4 @@ stars[i].apFlux = pow(10.0, -0.4*stars[i].Map); stars[i].apFluxErr = stars[i].apFlux * fSN; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -893,4 +908,12 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].detID = i; stars[i].X = X[i]; @@ -972,14 +995,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -1001,4 +1014,12 @@ flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); + + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } stars[i].detID = i; @@ -1081,14 +1102,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -1118,4 +1129,12 @@ flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); + + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } stars[i].detID = i; @@ -1221,14 +1240,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -1237,5 +1246,5 @@ } -int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { +int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime) { int i, j; @@ -1306,4 +1315,6 @@ float *PROF_FILL ; ALLOCATE (PROF_FILL , float , 9*Nstars); + float magtime = 2.5*log10(exptime); + // assign the storage arrays for (i = 0; i < Nstars; i++) { @@ -1318,6 +1329,6 @@ F25_ARATIO [i] = 0.9; F25_THETA [i] = 10.0; - PETRO_MAG [i] = M[i]; - PETRO_MAG_ERR [i] = fSN; + PETRO_MAG [i] = M[i] + zeroPt + magtime; + PETRO_MAG_ERR [i] = 1.0 / fSN; PETRO_RADIUS [i] = 8.0; if (i == 2) { @@ -1528,5 +1539,5 @@ int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { - int i, j; + int i, j, k; Header header; FTable ftable; @@ -1553,12 +1564,14 @@ // create intermediate storage arrays - int *IPP_IDET ; ALLOCATE (IPP_IDET , int , Nstars); - float *X_APER ; ALLOCATE (X_APER , float , Nstars); - float *Y_APER ; ALLOCATE (Y_APER , float , Nstars); - float *PSF_FWHM ; ALLOCATE (PSF_FWHM , float , Nstars); - float *APER_FLUX ; ALLOCATE (APER_FLUX , float , Nstars*9); - float *APER_FLUX_ERR ; ALLOCATE (APER_FLUX_ERR , float , Nstars*9); - float *APER_FLUX_STDEV ; ALLOCATE (APER_FLUX_STDEV , float , Nstars*9); - float *APER_FILL ; ALLOCATE (APER_FILL , float , Nstars*9); + int *IPP_IDET ; ALLOCATE (IPP_IDET , int , 3*Nstars); + float *X_APER ; ALLOCATE (X_APER , float , 3*Nstars); + float *Y_APER ; ALLOCATE (Y_APER , float , 3*Nstars); + float *PSF_FWHM ; ALLOCATE (PSF_FWHM , float , 3*Nstars); + float *APER_FLUX ; ALLOCATE (APER_FLUX , float , 3*Nstars*9); + float *APER_FLUX_ERR ; ALLOCATE (APER_FLUX_ERR , float , 3*Nstars*9); + float *APER_FLUX_STDEV ; ALLOCATE (APER_FLUX_STDEV , float , 3*Nstars*9); + float *APER_FILL ; ALLOCATE (APER_FILL , float , 3*Nstars*9); + + float psfFWHM[3] = {4.8, 6.0, 8.0}; // assign the storage arrays @@ -1567,25 +1580,29 @@ // float fSN = 1.0 / sqrt(flux); - IPP_IDET [i] = i; - X_APER [i] = X[i]; - Y_APER [i] = Y[i]; - PSF_FWHM [i] = FX*2.8; - for (j = 0; j < 9; j ++) { - APER_FLUX [j + 9*i] = flux; - APER_FLUX_ERR [j + 9*i] = sqrt(flux); - APER_FLUX_STDEV[j + 9*i] = sqrt(flux); - APER_FILL [j + 9*i] = 0.9; + for (k = 0; k < 3; k++) { + int N = 3*i + k; + IPP_IDET [N] = i; + X_APER [N] = X[i]; + Y_APER [N] = Y[i]; + PSF_FWHM [N] = psfFWHM[k]; + for (j = 0; j < 9; j ++) { + float Flux = flux*(1.0 + 0.02*j - 0.2); + APER_FLUX [j + 9*N] = Flux; + APER_FLUX_ERR [j + 9*N] = sqrt(Flux); + APER_FLUX_STDEV[j + 9*N] = sqrt(Flux*1.1); + APER_FILL [j + 9*N] = 0.9; + } } } // set the table data - gfits_set_bintable_column (&header, &ftable, "IPP_IDET" , IPP_IDET , Nstars); - gfits_set_bintable_column (&header, &ftable, "X_APER" , X_APER , Nstars); - gfits_set_bintable_column (&header, &ftable, "Y_APER" , Y_APER , Nstars); - gfits_set_bintable_column (&header, &ftable, "PSF_FWHM" , PSF_FWHM , Nstars); - gfits_set_bintable_column (&header, &ftable, "APER_FLUX" , APER_FLUX , Nstars); - gfits_set_bintable_column (&header, &ftable, "APER_FLUX_ERR" , APER_FLUX_ERR , Nstars); - gfits_set_bintable_column (&header, &ftable, "APER_FLUX_STDEV" , APER_FLUX_STDEV , Nstars); - gfits_set_bintable_column (&header, &ftable, "APER_FILL" , APER_FILL , Nstars); + gfits_set_bintable_column (&header, &ftable, "IPP_IDET" , IPP_IDET , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "X_APER" , X_APER , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_APER" , Y_APER , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "PSF_FWHM" , PSF_FWHM , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX" , APER_FLUX , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX_ERR" , APER_FLUX_ERR , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX_STDEV" , APER_FLUX_STDEV , 3*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FILL" , APER_FILL , 3*Nstars); } /**/