Index: branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/loadgalphot_readstars.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 39266) @@ -204,4 +204,6 @@ fit->ClipNsigma = 5.0; + // use the list of index values from above to join entries with the same index + int Nbad = 0; for (i = 0; i < Nrow; i++) { Index: branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/mkcmf.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/mkcmf.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/addstar/src/mkcmf.c (revision 39266) @@ -13,6 +13,4 @@ # define FLAGS 0x1101 -void gauss_init (int Nbin); -double rnd_gauss (double mean, double sigma); void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); @@ -282,5 +280,5 @@ } - gauss_init (2048); + ohana_gaussdev_init (); // load test stars from a file: @@ -564,78 +562,4 @@ -static int Ngaussint = 0; -static double *gaussint; - -extern double drand48(); - -double gaussian (double x, double mean, double sigma) { - - double f; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} - -/* integrate a gaussian from -5 sigma to +5 sigma */ -void gauss_init (int Nbin) { - - int i; - double val, x, dx, dx1, dx2, dx3, df; - double mean, sigma; - - /* no need to generate this if it already exists */ - if (Ngaussint == Nbin) return; - - // A = time(NULL); - // // XXX this is expensive if called a lot (1 sec min) - // // for (B = 0; A == time(NULL); B++); - // B = A + 10000; - // srand48(B); - - Ngaussint = Nbin; - ALLOCATE (gaussint, double, Ngaussint + 1); - - val = 0; - dx = 1.0 / Ngaussint; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - mean = 0.0; - sigma = 1.0; - - for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , mean, sigma) + - 9.0*gaussian(x+dx1, mean, sigma) + - 9.0*gaussian(x+dx2, mean, sigma) + - 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - if (val > (i + 0.5) / (double) Ngaussint) { - gaussint[i] = x + dx / 2.0; - i++; - } - } -} - -double rnd_gauss (double mean, double sigma) { - - int i; - double y; - - y = drand48(); - i = Ngaussint*y; - y = gaussint[i]*sigma + mean; - - return (y); - -} - -double int_gauss (int i) { - double y; - y = gaussint[i]; - return (y); -} - void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars) { @@ -654,7 +578,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -700,7 +624,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -749,7 +673,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -800,7 +724,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -857,7 +781,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -918,7 +842,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -987,7 +911,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -1094,7 +1018,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); @@ -1198,7 +1122,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN * ohana_gaussdev_rnd(0.0, 1.0); flux = pow(10.0, -0.4 * M[i]); fSN = 1.0 / sqrt(flux); @@ -1324,7 +1248,7 @@ 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); + X[i] += FX * fSN * ohana_gaussdev_rnd(0.0, 1.0); + Y[i] += FY * fSN * ohana_gaussdev_rnd(0.0, 1.0); + M[i] += fSN*ohana_gaussdev_rnd(0.0, 1.0); flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); Index: branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 39266) @@ -85,5 +85,5 @@ N1[Nstars] = i; Nstars ++; - input[0].found_t[i] = -1; + input[0].found_t[i] = FALSE; } if (Nstars < 1) { @@ -180,4 +180,6 @@ off_t N = N1[i]; + input[0].found_t[N] = TRUE; + /* make sure there is space for next Nmeasure entries */ if (Nmeas + input[0].average[N].Nmeasure >= NMEAS) { @@ -209,5 +211,4 @@ Nreplace = replace_tycho (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[Minp]); if (Nreplace == 6) { - input[0].found_t[N] = Nmeas; output[0].found_t[n] = Nmeas; i++; @@ -224,5 +225,4 @@ // XXX this does not support lensing, starpar, or galphot measurements if (replace_match (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[offset])) { - input[0].found_t[N] = output[0].average[n].measureOffset; continue; } @@ -274,5 +274,4 @@ } } - input[0].found_t[N] = Nmeas; output[0].average[n].Nmeasure ++; Nmeas ++; @@ -394,9 +393,11 @@ } - if (input[0].found_t[N] >= 0) continue; + if (input[0].found_t[N]) continue; // if we are using MATCHED_TABLES, we are going to leave the edge cases in their // source catalog, even if they have leaked beyond the edge if (!MATCHED_TABLES && !IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue; + + input[0].found_t[N] = TRUE; // XXX should we accept the input measurements for these fields? @@ -450,5 +451,4 @@ /* we set next[Nmeas] to -1 here, and update correctly below */ - input[0].found_t[N] = Nmeas; next_meas[Nmeas] = -1; Nmeas ++; @@ -559,13 +559,4 @@ } - /* note stars which have been found in this catalog */ - for (i = 0; i < input[0].Naverage; i++) { - if (input[0].found_t[i] > -1) { - input[0].found_t[i] = -2; - } else { - input[0].found_t[i] = -3; - } - } - /* check if the catalog has changed? if no change, no need to write */ output[0].objID = objID; // new max value, save on catalog close Index: branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c (revision 39266) @@ -229,6 +229,13 @@ float sersic_index = galphot->index; - if (galphot->modelType == 6) sersic_index = 1.0; - if (galphot->modelType == 7) sersic_index = 4.0; + if (galphot->modelType == 6) { + sersic_index = 1.0; + } + else if (galphot->modelType == 7) { + sersic_index = 4.0; + } + else { + sersic_index = 1.0 / (2.0 * sersic_index); + } PRINT_FLOAT(gal_buffer, galphot->mag, "%.6e, "); Index: branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FWobjects_dvopsps_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FWobjects_dvopsps_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_FWobjects_dvopsps_catalog.c (revision 39266) @@ -148,12 +148,18 @@ "FLUX_PSF_WRP_ERR FLOAT, " "FLUX_PSF_WRP_STDEV FLOAT, " - + "MAG_PSF_WRP FLOAT, " + "MAG_PSF_WRP_ERR FLOAT, " + "FLUX_KRON_WRP FLOAT, " "FLUX_KRON_WRP_ERR FLOAT, " "FLUX_KRON_WRP_STDEV FLOAT, " - + "MAG_KRON_WRP FLOAT, " + "MAG_KRON_WRP_ERR FLOAT, " + "FLUX_AP_WRP FLOAT, " "FLUX_AP_WRP_ERR FLOAT, " - "FLUX_AP_WRP_STDEV FLOAT, " + "FLUX_AP_WRP_STDEV FLOAT, " + "MAG_AP_WRP FLOAT, " + "MAG_AP_WRP_ERR FLOAT, " "FLAGS INT " @@ -188,4 +194,6 @@ "FLUX_STD_AP_R5 FLOAT, " "FLUX_FIL_AP_R5 FLOAT, " + "MAG_AP_R5 FLOAT, " + "MAG_ERR_AP_R5 FLOAT, " "FLUX_AP_R6 FLOAT, " @@ -193,4 +201,6 @@ "FLUX_STD_AP_R6 FLOAT, " "FLUX_FIL_AP_R6 FLOAT, " + "MAG_AP_R6 FLOAT, " + "MAG_ERR_AP_R6 FLOAT, " "FLUX_AP_R7 FLOAT, " @@ -198,4 +208,6 @@ "FLUX_STD_AP_R7 FLOAT, " "FLUX_FIL_AP_R7 FLOAT, " + "MAG_AP_R7 FLOAT, " + "MAG_ERR_AP_R7 FLOAT, " "X11_SM_OBJ FLOAT, " @@ -262,13 +274,19 @@ "FLUX_PSF_WRP, " "FLUX_PSF_WRP_ERR, " - "FLUX_PSF_WRP_STDEV, " + "FLUX_PSF_WRP_STDEV, " + "MAG_PSF_WRP, " + "MAG_PSF_WRP_ERR, " "FLUX_KRON_WRP, " "FLUX_KRON_WRP_ERR, " - "FLUX_KRON_WRP_STDEV, " + "FLUX_KRON_WRP_STDEV, " + "MAG_KRON_WRP, " + "MAG_KRON_WRP_ERR, " "FLUX_AP_WRP, " "FLUX_AP_WRP_ERR, " - "FLUX_AP_WRP_STDEV, " + "FLUX_AP_WRP_STDEV, " + "MAG_AP_WRP, " + "MAG_AP_WRP_ERR, " "FLAGS" @@ -287,4 +305,6 @@ "FLUX_STD_AP_R5, " "FLUX_FIL_AP_R5, " + "MAG_AP_R5, " + "MAG_ERR_AP_R5, " "FLUX_AP_R6, " @@ -292,4 +312,6 @@ "FLUX_STD_AP_R6, " "FLUX_FIL_AP_R6, " + "MAG_AP_R6, " + "MAG_ERR_AP_R6, " "FLUX_AP_R7, " @@ -297,4 +319,6 @@ "FLUX_STD_AP_R7, " "FLUX_FIL_AP_R7, " + "MAG_AP_R7, " + "MAG_ERR_AP_R7, " "X11_SM_OBJ, " @@ -348,8 +372,14 @@ float meanPSFFluxErr = NAN; float meanPSFFluxStd = NAN; + float meanPSFMag = NAN; + float meanPSFMagErr = NAN; if (isfinite(secfilt->dFpsfWrp) && isfinite(secfilt->FpsfWrp)) { meanPSFFlux = secfilt-> FpsfWrp; meanPSFFluxErr = secfilt->dFpsfWrp; meanPSFFluxStd = secfilt->sFpsfWrp; + if (meanPSFFlux > 0.0) { + meanPSFMag = -2.5 * log10(meanPSFFlux / 3631.0); + meanPSFMagErr = (2.5 * meanPSFFluxErr) / (meanPSFFlux * log(10)); + } } @@ -357,8 +387,14 @@ float meanKronFluxErr = NAN; float meanKronFluxStd = NAN; + float meanKronMag = NAN; + float meanKronMagErr = NAN; if (isfinite(secfilt->dFkronWrp) && isfinite(secfilt->FkronWrp)) { meanKronFlux = secfilt-> FkronWrp; meanKronFluxErr = secfilt->dFkronWrp; meanKronFluxStd = secfilt->sFkronWrp; + if (meanKronFlux > 0.0) { + meanKronMag = -2.5 * log10(meanKronFlux / 3631.0); + meanKronMagErr = (2.5 * meanKronFluxErr) / (meanKronFlux * log(10)); + } } @@ -366,8 +402,14 @@ float meanApFluxErr = NAN; float meanApFluxStd = NAN; + float meanApMag = NAN; + float meanApMagErr = NAN; if (isfinite(secfilt->dFapWrp) && isfinite(secfilt->FapWrp)) { meanApFlux = secfilt-> FapWrp; meanApFluxErr = secfilt->dFapWrp; meanApFluxStd = secfilt->sFapWrp; + if (meanApFlux > 0.0) { + meanApMag = -2.5 * log10(meanApFlux / 3631.0); + meanApMagErr = (2.5 * meanApFluxErr) / (meanApFlux * log(10)); + } } @@ -382,12 +424,18 @@ PRINT_FLOAT(sec_buffer, meanPSFFluxErr, "%.6e, "); // uflux precision PRINT_FLOAT(sec_buffer, meanPSFFluxStd, "%.6e, "); // uflux precision - + PRINT_FLOAT(sec_buffer, meanPSFMag, "%.6f, "); + PRINT_FLOAT(sec_buffer, meanPSFMagErr, "%.6f, "); + PRINT_FLOAT(sec_buffer, meanKronFlux, "%.6e, "); // uflux precision PRINT_FLOAT(sec_buffer, meanKronFluxErr, "%.6e, "); // uflux precision PRINT_FLOAT(sec_buffer, meanKronFluxStd, "%.6e, "); // uflux precision + PRINT_FLOAT(sec_buffer, meanKronMag, "%.6f, "); + PRINT_FLOAT(sec_buffer, meanKronMagErr, "%.6f, "); PRINT_FLOAT(sec_buffer, meanApFlux, "%.6e, "); // uflux precision PRINT_FLOAT(sec_buffer, meanApFluxErr, "%.6e, "); // uflux precision PRINT_FLOAT(sec_buffer, meanApFluxStd, "%.6e, "); // uflux precision + PRINT_FLOAT(sec_buffer, meanApMag, "%.6f, "); + PRINT_FLOAT(sec_buffer, meanApMagErr, "%.6f, "); PrintIOBuffer (sec_buffer, "%u ", secfilt->flags); @@ -402,16 +450,42 @@ PrintIOBuffer (cpy_buffer, "%hd, ", lensobj->photcode); + float magApR5 = NAN; + float magErrApR5 = NAN; + float magApR6 = NAN; + float magErrApR6 = NAN; + float magApR7 = NAN; + float magErrApR7 = NAN; + + if (lensobj->F_ApR5 > 0.0) { + magApR5 = -2.5 * log10(lensobj->F_ApR5 / 3631.0); + magErrApR5 = (2.5 * lensobj->dF_ApR5) / (lensobj->F_ApR5 * log(10)); + } + if (lensobj->F_ApR6 > 0.0) { + magApR6 = -2.5 * log10(lensobj->F_ApR6 / 3631.0); + magErrApR6 = (2.5 * lensobj->dF_ApR6) / (lensobj->F_ApR6 * log(10)); + } + if (lensobj->F_ApR7 > 0.0) { + magApR7 = -2.5 * log10(lensobj->F_ApR7 / 3631.0); + magErrApR7 = (2.5 * lensobj->dF_ApR7) / (lensobj->F_ApR7 * log(10)); + } + PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR5, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR5, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR5, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR5, "%.6e, "); + PRINT_FLOAT(cpy_buffer, magApR5, "%.6f, "); + PRINT_FLOAT(cpy_buffer, magErrApR5, "%.6f, "); PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR6, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR6, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR6, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR6, "%.6e, "); + PRINT_FLOAT(cpy_buffer, magApR6, "%.6f, "); + PRINT_FLOAT(cpy_buffer, magErrApR6, "%.6f, "); PRINT_FLOAT(cpy_buffer, lensobj-> F_ApR7, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->dF_ApR7, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->sF_ApR7, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->fF_ApR7, "%.6e, "); + PRINT_FLOAT(cpy_buffer, magApR7, "%.6f, "); + PRINT_FLOAT(cpy_buffer, magErrApR7, "%.6f, "); PRINT_FLOAT(cpy_buffer, lensobj->X11_sm_obj, "%.6e, "); PRINT_FLOAT(cpy_buffer, lensobj->X12_sm_obj, "%.6e, "); Index: branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_diffobj_dvopsps_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_diffobj_dvopsps_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_diffobj_dvopsps_catalog.c (revision 39266) @@ -94,7 +94,9 @@ "CAT_ID INT, " "RA_MEAN DOUBLE, " - "DEC_MEAN DOUBLE, " + "DEC_MEAN DOUBLE, " "RA_ERR FLOAT, " - "DEC_ERR FLOAT, " + "DEC_ERR FLOAT, " + "IAU_NAME VARCHAR(32), " + "PSO_NAME VARCHAR(32), " "FLAGS INT " ")\n", basename); @@ -167,5 +169,7 @@ "DEC_MEAN, " "RA_ERR, " - "DEC_ERR, " + "DEC_ERR, " + "IAU_NAME, " + "PSO_NAME, " "FLAGS " ") VALUES \n", basename); @@ -218,4 +222,56 @@ PRINT_FLOAT(ave_buffer, average->dD, "%.5f, "); // 0.010 mas precision + // Add names. + int ra_hr,ra_min,ra_sec,ra_fracsec; + int dec_deg,dec_min,dec_sec,dec_fracsec; + double tmp_ra = average->R; + double tmp_dec = average->D; + double dec_sign = tmp_dec / fabs(tmp_dec); + tmp_dec = fabs(tmp_dec); + + tmp_ra /= 15.0; + ra_hr = (int) floor(tmp_ra); + + tmp_ra -= ra_hr; + tmp_ra *= 60.0; + ra_min = (int) floor(tmp_ra); + + tmp_ra -= ra_min; + tmp_ra *= 60.0; + ra_sec = (int) floor(tmp_ra); + + tmp_ra -= ra_sec; + tmp_ra *= 1000.0; + ra_fracsec = (int) floor(tmp_ra); + + dec_deg = (int) floor(tmp_dec); + + tmp_dec -= dec_deg; + tmp_dec *= 60.0; + dec_min = (int) floor(tmp_dec); + + tmp_dec -= dec_min; + tmp_dec *= 60.0; + dec_sec = (int) floor(tmp_dec); + + tmp_dec -= dec_sec; + tmp_dec *= 1000.0; + dec_fracsec = (int) floor(tmp_dec); + + dec_deg *= (int) dec_sign; + + // IAU NAME + PrintIOBuffer(ave_buffer, "'PSO J%02d%02d%02d.%03d%+02d%02d%02d.%03d', ", + ra_hr,ra_min,ra_sec,ra_fracsec, + dec_deg,dec_min,dec_sec,dec_fracsec); + + // PSO NAME + tmp_ra = average->R; + tmp_dec = average->D; + PrintIOBuffer(ave_buffer, "'PSO J%.5f%+.5f', ", + tmp_ra,tmp_dec); + + // End names + PrintIOBuffer (ave_buffer, "%u ", average->flags); PrintIOBuffer (ave_buffer, "),\n"); Index: branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 39266) @@ -101,5 +101,7 @@ "EPOCH_MEAN DOUBLE, " "RA_ERR FLOAT, " - "DEC_ERR FLOAT, " + "DEC_ERR FLOAT, " + "IAU_NAME VARCHAR(32), " + "PSO_NAME VARCHAR(32), " "CHISQ_POS FLOAT, " "CHISQ_PM FLOAT, " @@ -208,5 +210,7 @@ "EPOCH_MEAN, " "RA_ERR, " - "DEC_ERR, " + "DEC_ERR, " + "IAU_NAME, " + "PSO_NAME, " "CHISQ_POS, " "CHISQ_PM, " @@ -276,4 +280,73 @@ PRINT_FLOAT(ave_buffer, average->dD, "%.5f, "); // 0.010 mas precision + // Add names. + int ra_hr,ra_min,ra_sec,ra_fracsec; + int dec_deg,dec_min,dec_sec,dec_fracsec; + double tmp_ra,tmp_dec; + + if (secfilt->NstackDet != 0) { + tmp_ra = average->Rstk; + tmp_dec = average->Dstk; + } + else { + tmp_ra = average->R; + tmp_dec = average->D; + } + + double dec_sign = tmp_dec / fabs(tmp_dec); + tmp_dec = fabs(tmp_dec); + + + tmp_ra /= 15.0; + ra_hr = (int) floor(tmp_ra); + + tmp_ra -= ra_hr; + tmp_ra *= 60.0; + ra_min = (int) floor(tmp_ra); + + tmp_ra -= ra_min; + tmp_ra *= 60.0; + ra_sec = (int) floor(tmp_ra); + + tmp_ra -= ra_sec; + tmp_ra *= 1000.0; + ra_fracsec = (int) floor(tmp_ra); + + dec_deg = (int) floor(tmp_dec); + + tmp_dec -= dec_deg; + tmp_dec *= 60.0; + dec_min = (int) floor(tmp_dec); + + tmp_dec -= dec_min; + tmp_dec *= 60.0; + dec_sec = (int) floor(tmp_dec); + + tmp_dec -= dec_sec; + tmp_dec *= 1000.0; + dec_fracsec = (int) floor(tmp_dec); + + dec_deg *= (int) dec_sign; + + // IAU NAME + PrintIOBuffer(ave_buffer, "'PSO J%02d%02d%02d.%03d%+02d%02d%02d.%03d', ", + ra_hr,ra_min,ra_sec,ra_fracsec, + dec_deg,dec_min,dec_sec,dec_fracsec); + + // PSO NAME + if (secfilt->NstackDet != 0) { + tmp_ra = average->Rstk; + tmp_dec = average->Dstk; + } + else { + tmp_ra = average->R; + tmp_dec = average->D; + } + + PrintIOBuffer(ave_buffer, "'PSO J%.5f%+.5f', ", + tmp_ra,tmp_dec); + + // End names + PRINT_FLOAT(ave_buffer, average->ChiSqAve, "%.4f, "); PRINT_FLOAT(ave_buffer, average->ChiSqPM, "%.4f, "); Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/Makefile (revision 39266) @@ -38,5 +38,4 @@ $(SRC)/fakestar_io.$(ARCH).o \ $(SRC)/insert_fakestar.$(ARCH).o \ -$(SRC)/gaussian.$(ARCH).o \ $(SRC)/fakeastro_images.$(ARCH).o \ $(SRC)/fakeastro_images_region.$(ARCH).o \ Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/fakeastro.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/fakeastro.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/fakeastro.c (revision 39266) @@ -3,5 +3,5 @@ int main (int argc, char **argv) { - gauss_init (50000); + ohana_gaussdev_init (); /* get configuration info, args */ Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_2mass_measures.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_2mass_measures.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_2mass_measures.c (revision 39266) @@ -80,6 +80,6 @@ // uR,uD in linear arcsec - double dRsee = rnd_gauss (0.0, 1.0 / SN); - double dDsee = rnd_gauss (0.0, 1.0 / SN); + double dRsee = ohana_gaussdev_rnd (0.0, 1.0 / SN); + double dDsee = ohana_gaussdev_rnd (0.0, 1.0 / SN); double dRoff = (dRpm + dRsee) / 3600.0; Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fake_stars_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 39266) @@ -113,6 +113,6 @@ // uR,uD in linear arcsec - double dRsee = rnd_gauss (0.0, 1.0 / SN); - double dDsee = rnd_gauss (0.0, 1.0 / SN); + double dRsee = ohana_gaussdev_rnd (0.0, 1.0 / SN); + double dDsee = ohana_gaussdev_rnd (0.0, 1.0 / SN); double dRoff = (dRpm + dRsee) / 3600.0; Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakeqsos.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakeqsos.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakeqsos.c (revision 39266) @@ -51,5 +51,5 @@ // Mr will be interpretted as m_r - double Mr = rnd_gauss (18.0, 1.5); + double Mr = ohana_gaussdev_rnd (18.0, 1.5); stars[i].R = R; Index: branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakestars.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakestars.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/fakeastro/src/make_fakestars.c (revision 39266) @@ -48,5 +48,5 @@ int inPatch = FALSE; while (!inPatch) { - z = rnd_gauss (0.0, FAKEASTRO_ZGAL); + z = ohana_gaussdev_rnd (0.0, FAKEASTRO_ZGAL); r = sqrt(drand48()) * FAKEASTRO_RGAL; Lrad = drand48() * 2 * M_PI; @@ -98,7 +98,7 @@ double Mr; if (bigPeak) { - Mr = rnd_gauss (11.25, 1.0); + Mr = ohana_gaussdev_rnd (11.25, 1.0); } else { - Mr = rnd_gauss (10.00, 2.5); + Mr = ohana_gaussdev_rnd (10.00, 2.5); } Index: branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-sim.d =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-sim.d (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-sim.d (revision 39266) @@ -27,6 +27,7 @@ FIELD XoffDCR, X_OFF_DCR, float, X offset from correction, pixels FIELD YoffDCR, Y_OFF_DCR, float, Y offset from correction, pixels -FIELD RoffGAL, R_OFF_GAL, float, RA offset from correction, arcsec -FIELD DoffGAL, D_OFF_GAL, float, DEC offset from correction, arcsec + +FIELD Mflat, M_FLAT, float, Static Flat-field offset, mag +FIELD dummy2, PADDING, int, unused 4 bytes FIELD t, TIME, int, time in seconds (UNIX) @@ -48,5 +49,5 @@ FIELD dRsys, POS_SYS_ERR, short, systematic error from astrom, 1/100 of pixels -FIELD dummy, PADDING, short, padding +FIELD dummy, PADDING2, short, padding FIELD posangle, POSANGLE, short, position angle sky to chip, (0xffff/360) deg FIELD pltscale, PLTSCALE, float, plate scale, arcsec/pixel Index: branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-v5.d =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-v5.d (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure-ps1-v5.d (revision 39266) @@ -42,6 +42,7 @@ FIELD XoffCAM, X_OFF_CAM, float, X offset from correction, pixels FIELD YoffCAM, Y_OFF_CAM, float, Y offset from correction, pixels -FIELD RoffGAL, R_OFF_GAL, float, RA offset from correction, arcsec -FIELD DoffGAL, D_OFF_GAL, float, DEC offset from correction, arcsec + +FIELD Mflat, M_FLAT, float, Static Flat-field offset, mag +FIELD dummy2, PADDING, int, unused 4 bytes # could these be packed into fewer bits? Index: branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure.d =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure.d (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libautocode/def/measure.d (revision 39266) @@ -43,9 +43,6 @@ FIELD YoffCAM, Y_OFF_CAM, float, Y offset from correction, pixels -# XXX I can deprecate these as I am going to apply the correct uR,uD offset -# XXX not sure how to use this yet... -# XXX replace RoffGall with Moff -> static component of the zpt correction -FIELD RoffGAL, R_OFF_GAL, float, RA offset from correction, arcsec -FIELD DoffGAL, D_OFF_GAL, float, DEC offset from correction, arcsec +FIELD Mflat, M_FLAT, float, Static Flat-field offset, mag +FIELD dummy2, PADDING, int, unused 4 bytes # could these be packed into fewer bits? Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/dvo.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/dvo.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/dvo.h (revision 39266) @@ -476,8 +476,7 @@ double R; double D; - float RoffGAL; - float DoffGAL; float M; float Mcal; + float Mflat; float dM; float airmass; @@ -620,6 +619,6 @@ float XoffDCR; // X offset from correction (pixels) float YoffDCR; // Y offset from correction (pixels) - float RoffGAL; // RA offset from correction (arcsec) - float DoffGAL; // DEC offset from correction (arcsec) + float Mflat; // flat offset from correction (arcsec) + int padding2; // dummy float Sky; // local estimate of sky flux (counts/sec) float dSky; // local estimate of sky flux (counts/sec) Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/libdvo_astro.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/libdvo_astro.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/include/libdvo_astro.h (revision 39266) @@ -138,4 +138,5 @@ } Image; +CoordTransform *AllocTransform (double phi, double Xo, double xo); CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output); int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform); Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/coord_systems.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/coord_systems.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/coord_systems.c (revision 39266) @@ -1,4 +1,28 @@ # include "dvo.h" +CoordTransform *AllocTransform (double phi, double Xo, double xo) { + + CoordTransform *transform; + ALLOCATE (transform, CoordTransform, 1); + transform->isIdentity = FALSE; + + transform->phi = phi*RAD_DEG; + transform->Xo = Xo*RAD_DEG; + transform->xo = xo; + + // pre-calculated constants: + transform->sin_phi_cos_Xo = sin(transform->phi)*cos(transform->Xo); + transform->sin_phi_sin_Xo = sin(transform->phi)*sin(transform->Xo); + transform->cos_phi = cos(transform->phi); + + transform->cos_phi_cos_Xo = cos(transform->phi)*cos(transform->Xo); + transform->cos_phi_sin_Xo = cos(transform->phi)*sin(transform->Xo); + transform->sin_phi = sin(transform->phi); + transform->cos_Xo = cos(transform->Xo); + transform->sin_Xo = sin(transform->Xo); + + return transform; +} + CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output) { Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 39266) @@ -597,6 +597,4 @@ case MEAS_XOFF_CAM: { value.Flt = measure[0].XoffCAM; break; } case MEAS_YOFF_CAM: { value.Flt = measure[0].YoffCAM; break; } - case MEAS_ROFF_GAL: { value.Flt = measure[0].RoffGAL; break; } - case MEAS_DOFF_GAL: { value.Flt = measure[0].DoffGAL; break; } case MEAS_XFIELD: /* offset relative to exposure center in ra,dec space */ @@ -723,17 +721,5 @@ case MEAS_MCAL_OFFSET: { value.Flt = measure[0].Mcal; break; } - - case MEAS_FLAT: /* OK */ - // flat = measure.Mcal - image.Mcal - if (REMOTE_CLIENT) { - ImageMetadata *image = MatchImageMetadataDVO (measure[0].imageID); - if (image == NULL) break; - value.Flt = measure[0].Mcal - image[0].Mcal; - } else { - Image *image = MatchImageDVO (measure[0].t, measure[0].photcode, measure[0].imageID); - if (image == NULL) break; - value.Flt = measure[0].Mcal - image[0].Mcal; - } - break; + case MEAS_FLAT: { value.Flt = measure[0].Mflat; break; } // we have measure[0].Xccd,Yccd and image[0].NX,NY. Find the distance to the center Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_catalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_catalog.c (revision 39266) @@ -332,6 +332,7 @@ measure->XoffCAM = NAN; measure->YoffCAM = NAN; - measure->RoffGAL = NAN; - measure->DoffGAL = NAN; + + measure->Mflat = 0.0; + measure->dummy2 = 0; measure->Sky = NAN; Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c (revision 39266) @@ -29,6 +29,5 @@ out[i].XoffDCR = in[i].XoffDCR; out[i].YoffDCR = in[i].YoffDCR; - out[i].RoffGAL = in[i].RoffGAL; - out[i].DoffGAL = in[i].DoffGAL; + out[i].Mflat = in[i].Mflat; out[i].t = in[i].t; out[i].averef = in[i].averef; @@ -75,6 +74,5 @@ out[i].XoffDCR = in[i].XoffDCR; out[i].YoffDCR = in[i].YoffDCR; - out[i].RoffGAL = in[i].RoffGAL; - out[i].DoffGAL = in[i].DoffGAL; + out[i].Mflat = in[i].Mflat; out[i].t = in[i].t; out[i].averef = in[i].averef; Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 39266) @@ -42,6 +42,5 @@ out[i].XoffCAM = in[i].XoffCAM; out[i].YoffCAM = in[i].YoffCAM; - out[i].RoffGAL = in[i].RoffGAL; - out[i].DoffGAL = in[i].DoffGAL; + out[i].Mflat = in[i].Mflat; out[i].Sky = in[i].Sky; out[i].dSky = in[i].dSky; @@ -116,6 +115,5 @@ out[i].XoffCAM = in[i].XoffCAM; out[i].YoffCAM = in[i].YoffCAM; - out[i].RoffGAL = in[i].RoffGAL; - out[i].DoffGAL = in[i].DoffGAL; + out[i].Mflat = in[i].Mflat; out[i].Sky = in[i].Sky; out[i].dSky = in[i].dSky; @@ -1678,6 +1676,5 @@ out[i].XoffDCR = in[i].XoffDCR; out[i].YoffDCR = in[i].YoffDCR; - out[i].RoffGAL = in[i].RoffGAL; - out[i].DoffGAL = in[i].DoffGAL; + out[i].Mflat = in[i].Mflat; out[i].Sky = in[i].Sky; out[i].dSky = in[i].dSky; Index: branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_tiny_values.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 39266) @@ -22,6 +22,5 @@ measureT[0].dM = measure[0].dM; measureT[0].airmass = measure[0].airmass; - measureT[0].RoffGAL = measure[0].RoffGAL; - measureT[0].DoffGAL = measure[0].DoffGAL; + measureT[0].Mflat = measure[0].Mflat; measureT[0].Xccd = measure[0].Xccd; measureT[0].Yccd = measure[0].Yccd; Index: branches/eam_branches/ipp-20151113/Ohana/src/libohana/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libohana/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libohana/Makefile (revision 39266) @@ -35,4 +35,5 @@ $(SRC)/sorts.$(ARCH).o \ $(SRC)/bisection.$(ARCH).o \ +$(SRC)/gaussdev.$(ARCH).o \ $(SRC)/string.$(ARCH).o \ $(SRC)/findexec.$(ARCH).o \ Index: branches/eam_branches/ipp-20151113/Ohana/src/libohana/include/ohana.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libohana/include/ohana.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/libohana/include/ohana.h (revision 39266) @@ -486,3 +486,9 @@ int sprintf_double (char *output, double value); -# endif +/* in gaussdev.c */ +double ohana_gaussian (double x, double mean, double sigma); +void ohana_gaussdev_init (void); +void ohana_gaussdev_free (void); +double ohana_gaussdev_rnd (double mean, double sigma); + +# endif Index: branches/eam_branches/ipp-20151113/Ohana/src/libohana/src/gaussdev.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/libohana/src/gaussdev.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/libohana/src/gaussdev.c (revision 39266) @@ -0,0 +1,83 @@ +# include + +// NOTE: this code uses a local static array which must be freed on exit +static int Ngaussint = 0; +static double *gaussint = NULL; + +// extern double drand48(); + +// return the value of a gaussian at position x with mean,sigma +double ohana_gaussian (double x, double mean, double sigma) { + + double f; + + double S2 = sigma*sigma; + + f = exp (-0.5 * SQ(x - mean) / S2) / sqrt(2 * M_PI * S2); + + return (f); +} + +/* integrate a gaussian from -10 sigma to +10 sigma */ +// how long does this take with a reasonably high resolution? +void ohana_gaussdev_init (void) { + + int i; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + if (gaussint) return; + + // 16k bins takes ~0.1 sec to generate on a ~3GHz core and gives good resolution + Ngaussint = 0x4000; + ALLOCATE (gaussint, double, Ngaussint + 1); + + val = 0; + dx = 1.0 / Ngaussint; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -10.0; (i < Ngaussint) && (x <= 10.0); x += dx) { + df = (3.0*ohana_gaussian(x , mean, sigma) + + 9.0*ohana_gaussian(x+dx1, mean, sigma) + + 9.0*ohana_gaussian(x+dx2, mean, sigma) + + 3.0*ohana_gaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i + 0.5) / (double) Ngaussint) { + gaussint[i] = x + dx / 2.0; + myAssert (i < Ngaussint + 1, "oops"); + i++; + } + } +} + +void ohana_gaussdev_free (void) { + free (gaussint); +} + +// return a number drawn from a gaussian deviate with mean, sigma +// must first call gaussdev_init() +double ohana_gaussdev_rnd (double mean, double sigma) { + + int i; + double y; + + myAssert (gaussint, "need to call gaussdev_init before calling gaussdev"); + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]*sigma + mean; + + return (y); +} + +// return the value of \int gaussian(x,mean,sigma) +// double gaussian_integral (double x, double mean, double sigma) { +double ohana_gaussian_integral (int i) { + + // int i = (x - mean) / (sigma * Ngaussint); + return gaussint[i]; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakedist.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakedist.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakedist.c (revision 39266) @@ -1,9 +1,3 @@ # include - -# define NGAUSS 2048 -double drand48(); -double gaussint[NGAUSS]; -double gaussian(); -double rnd_gauss(); # define MMIN 1.0 @@ -45,5 +39,5 @@ init_outmatrix (); - gauss_init (); + ohana_gaussdev_init (); read_datafiles (); @@ -349,6 +343,6 @@ if (noise < dMo) { *Noise = noise; - *mag = rnd_gauss (v, noise); - *color = rnd_gauss (uv, 1.4*noise); + *mag = ohana_gaussdev_rnd (v, noise); + *color = ohana_gaussdev_rnd (uv, 1.4*noise); return; } @@ -414,67 +408,4 @@ /*****************************************************************************/ - -double -rnd_gauss (mean, sigma) -double mean, sigma; -{ - - int i; - double y; - - y = drand48(); - i = NGAUSS*y; - y = gaussint[i]*sigma + mean; - - return (y); - -} - -gauss_init () -{ - - int i; - long A, B; - double val, x, dx, dx1, dx2, dx3, df; - double mean, sigma; - - A = time(NULL); - for (B = 0; A == time(NULL); B++); - srand48(B); - - val = 0; - dx = 0.001; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - mean = 0.0; - sigma = 1.0; - - for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , mean, sigma) + - 9.0*gaussian(x+dx1, mean, sigma) + - 9.0*gaussian(x+dx2, mean, sigma) + - 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - if (val > (i / (double) NGAUSS)) { - gaussint[i] = x + dx / 2.0; - i++; - } - } - return (x + dx / 2.0); -} - -double -gaussian (x, mean, sigma) -double x, mean, sigma; -{ - - double f, X; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} /*****************************************************************************/ @@ -505,6 +436,6 @@ Ngood = 1.0; for (k = 0; k < NTRY; k++) { - v = rnd_gauss (mag, noise); - uv = rnd_gauss (color, 1.4*noise); + v = ohana_gaussdev_rnd (mag, noise); + uv = ohana_gaussdev_rnd (color, 1.4*noise); x = (uv - UV0 - 0.7*A_V) / DUV; y = (v - V0 - Dist - A_V) / DV; Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakepop.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakepop.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakepop.c (revision 39266) @@ -2,18 +2,7 @@ # define MMIN 1.0 # define MMAX 120.0 -extern double drand48(); extern double rnd_mass(); double rnd_mass(); double term; -extern double rnd_gauss(); -extern double rnd_integrate (); -extern double gaussian (); - -double gaussian(); -double rnd_gauss(); -double rnd_integrate(); - -# define NGAUSS 2048 -double gaussint[NGAUSS]; double AgeS, AgeE, alpha, d, Av, dMo, dVo, dVref; @@ -47,5 +36,5 @@ gfits_scan (&UV_h, "DEC_Y", "%lf", 1, &ldM); - gauss_init (); + ohana_gaussdev_init (); fprintf (stderr, "beginning main loop\n"); @@ -68,6 +57,6 @@ v = gfits_get_matrix_value (&V_i, X, Y) + d + Av; noise = dVo*sqrt(1.0 + pow (10.0, (0.4*(v - dVref)))); - V = rnd_gauss (v, noise); - UV = rnd_gauss (uv, 1.4*noise); + V = ohana_gaussdev_rnd (v, noise); + UV = ohana_gaussdev_rnd (uv, 1.4*noise); if (noise < dMo) { fprintf (stdout, "%f %f %f %f %f %f\n", V, noise, V+UV, noise, mass, age); @@ -96,55 +85,4 @@ } - -double -rnd_gauss (mean, sigma) -double mean, sigma; -{ - - int i; - double y; - - y = drand48(); - i = NGAUSS*y; - y = gaussint[i]*sigma + mean; - - return (y); - -} - - -gauss_init () -{ - - int i; - long A, B; - double val, x, dx, dx1, dx2, dx3, df; - double mean, sigma; - - A = time(NULL); - for (B = 0; A == time(NULL); B++); - srand48(B); - - val = 0; - dx = 0.001; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - mean = 0.0; - sigma = 1.0; - - for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , mean, sigma) + - 9.0*gaussian(x+dx1, mean, sigma) + - 9.0*gaussian(x+dx2, mean, sigma) + - 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - if (val > (i / (double) NGAUSS)) { - gaussint[i] = x + dx / 2.0; - i++; - } - } - return (x + dx / 2.0); -} double @@ -172,18 +110,4 @@ return (x + dx / 2.0); } - -double -gaussian (x, mean, sigma) -double x, mean, sigma; -{ - - double f, X; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} - /*****************************************************************************/ Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakestars.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakestars.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fakestars.c (revision 39266) @@ -17,8 +17,4 @@ int found; } Stars; - -void gauss_init (); -double gaussian (double x, double mean, double sigma); -double rnd_gauss (double mean, double sigma); /* @@ -45,5 +41,5 @@ offset = atof (argv[4]); - gauss_init (); + ohana_gaussdev_init (); /* load header, open file */ @@ -98,5 +94,5 @@ dMr = 0.001 * stars[i].dM; dMo = hypot (dMs, dMr); - dM = rnd_gauss (0.0, dMo); + dM = ohana_gaussdev_rnd (0.0, dMo); stars[i].M += dM + offset; } Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fitdist.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fitdist.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/fitdist.c (revision 39266) @@ -1,9 +1,3 @@ # include - -# define NGAUSS 2048 -double drand48(); -double gaussint[NGAUSS]; -double gaussian(); -double rnd_gauss(); # define MMIN 1.0 @@ -45,5 +39,5 @@ fprintf (stderr, "read in parameters\n"); - gauss_init (); + ohana_gaussdev_init (); read_datafiles (); @@ -360,67 +354,4 @@ /*****************************************************************************/ -double -rnd_gauss (mean, sigma) -double mean, sigma; -{ - - int i; - double y; - - y = drand48(); - i = NGAUSS*y; - y = gaussint[i]*sigma + mean; - - return (y); - -} - -gauss_init () -{ - - int i; - long A, B; - double val, x, dx, dx1, dx2, dx3, df; - double mean, sigma; - - A = time(NULL); - for (B = 0; A == time(NULL); B++); - srand48(B); - - val = 0; - dx = 0.001; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - mean = 0.0; - sigma = 1.0; - - for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , mean, sigma) + - 9.0*gaussian(x+dx1, mean, sigma) + - 9.0*gaussian(x+dx2, mean, sigma) + - 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - if (val > (i / (double) NGAUSS)) { - gaussint[i] = x + dx / 2.0; - i++; - } - } - return (x + dx / 2.0); -} - -double -gaussian (x, mean, sigma) -double x, mean, sigma; -{ - - double f, X; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} - /*****************************************************************************/ @@ -453,6 +384,6 @@ Ngood = 1.0; for (k = 0; k < NTRY; k++) { - v = rnd_gauss (Mag, Noise); - uv = rnd_gauss (Color, 1.4*Noise); + v = ohana_gaussdev_rnd (Mag, Noise); + uv = ohana_gaussdev_rnd (Color, 1.4*Noise); x = (uv - UV0 - 0.7*A_V) / DUV; y = (v - V0 - Dist - A_V) / DV; Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtoage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtoage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtoage.c (revision 39266) @@ -2,12 +2,4 @@ # define MMIN 1.0 # define MMAX 120.0 -extern double drand48(); -extern double rnd_gauss(); -extern double rnd_integrate (); -extern double gaussian (); - -double gaussian(); -double rnd_gauss(); -double rnd_integrate(); void main (argc, argv) @@ -29,4 +21,6 @@ long A, B; + ohana_gaussdev_init(); + lAo = 0.0; ldA = 1.0; @@ -118,6 +112,6 @@ } else { - v = rnd_gauss (V, dV); - uv = rnd_gauss ((U-V), dUV); + v = ohana_gaussdev_rnd (V, dV); + uv = ohana_gaussdev_rnd ((U-V), dUV); } x = (uv - UV0 - 0.7*Av) / DUV; @@ -167,55 +161,2 @@ } -double -rnd_gauss (mean, sigma) -double mean, sigma; -{ - - double range, x; - - range = drand48(); - x = rnd_integrate (*gaussian, range, mean, sigma); - - return (x); - -} - - -double -rnd_integrate (function, range, mean, sigma) -double (*function) (); -double range, mean, sigma; -{ - - double val, x, dx, dx1, dx2, dx3, df; - - range += 0.0001; - val = 0; - dx = sigma / 10.0; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - - for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx) { - df = (3.0*function(x , mean, sigma) + - 9.0*function(x+dx1, mean, sigma) + - 9.0*function(x+dx2, mean, sigma) + - 3.0*function(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - } - return (x + dx / 2.0); -} - -double -gaussian (x, mean, sigma) -double x, mean, sigma; -{ - - double f, X; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} - Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtomass.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtomass.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/magtomass.c (revision 39266) @@ -1,15 +1,4 @@ # include -extern double drand48(); -extern double rnd_gauss(); -extern double rnd_integrate (); -extern double gaussian (); # define NEWWAY 1 - -double gaussian(); -double rnd_gauss(); -double rnd_integrate(); - -# define NGAUSS 2048 -double gaussint[NGAUSS]; void main (argc, argv) @@ -52,5 +41,5 @@ fprintf (stderr, "using mags in columns %d & %d\n", col1, col2); - gauss_init (); + ohana_gaussdev_init (); lAo = 0.1; @@ -147,6 +136,6 @@ Ngood = 1.0; for (k = 0; k < Ntry; k++) { - v = rnd_gauss (V, dV); - uv = rnd_gauss ((U-V), dUV); + v = ohana_gaussdev_rnd (V, dV); + uv = ohana_gaussdev_rnd ((U-V), dUV); x = (uv - UV0 - 0.7*Av) / DUV; y = (v - d - V0 - Av) / DV; @@ -202,91 +191,2 @@ */ -double -rnd_gauss (mean, sigma) -double mean, sigma; -{ - - int i; - double y; - - y = drand48(); - i = NGAUSS*y; - y = gaussint[i]*sigma + mean; - - return (y); - -} - - -gauss_init () -{ - - int i; - long A, B; - double val, x, dx, dx1, dx2, dx3, df; - double mean, sigma; - - A = time(NULL); - for (B = 0; A == time(NULL); B++); - srand48(B); - - val = 0; - dx = 0.001; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - mean = 0.0; - sigma = 1.0; - - for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , mean, sigma) + - 9.0*gaussian(x+dx1, mean, sigma) + - 9.0*gaussian(x+dx2, mean, sigma) + - 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - if (val > (i / (double) NGAUSS)) { - gaussint[i] = x + dx / 2.0; - i++; - } - } - return (x + dx / 2.0); -} - -double -rnd_integrate (function, range, mean, sigma) -double (*function) (); -double range, mean, sigma; -{ - - double val, x, dx, dx1, dx2, dx3, df; - - range += 0.0001; - val = 0; - dx = sigma / 100.0; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - - for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx) { - df = (3.0*function(x , mean, sigma) + - 9.0*function(x+dx1, mean, sigma) + - 9.0*function(x+dx2, mean, sigma) + - 3.0*function(x+dx3, mean, sigma)) * (dx1/8.0); - val += df; - } - return (x + dx / 2.0); -} - -double -gaussian (x, mean, sigma) -double x, mean, sigma; -{ - - double f, X; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - - return (f); - -} - Index: branches/eam_branches/ipp-20151113/Ohana/src/misc/src/mkgauss.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/misc/src/mkgauss.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/misc/src/mkgauss.c (revision 39266) @@ -1,11 +1,3 @@ # include - -void gauss_init (); -double gaussian (double x, double mean, double sigma); -double rnd_gauss (double mean, double sigma); - -/* NGAUSS defines the resolution of the distribution */ -# define NGAUSS 2048 -static double gaussint[NGAUSS]; int main (int argc, char **argv) { @@ -23,8 +15,8 @@ Npts = atoi (argv[3]); - gauss_init (); + ohana_gaussdev_init (); for (i = 0; i < Npts; i++) { - f = rnd_gauss (mean, sigma); + f = ohana_gaussdev_rnd (mean, sigma); fprintf (stdout, "%f\n", f); } @@ -33,53 +25,2 @@ } -void gauss_init () { - - int i; - long A, B; - double val, x, dx, dx1, dx2, dx3, df; - - A = time(NULL); - for (B = 0; A == time(NULL); B++); - srand48(B); - - val = 0; - dx = 0.001; - dx1 = dx / 3.0; - dx2 = 2.0*dx/3.0; - dx3 = dx; - - for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { - df = (3.0*gaussian(x , 0.0, 1.0) + - 9.0*gaussian(x+dx1, 0.0, 1.0) + - 9.0*gaussian(x+dx2, 0.0, 1.0) + - 3.0*gaussian(x+dx3, 0.0, 1.0)) * (dx1/8.0); - val += df; - if (val > (i / (double) NGAUSS)) { - gaussint[i] = x + dx / 2.0; - i++; - } - } - return; -} - -double gaussian (double x, double mean, double sigma) { - double f; - - f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); - return (f); -} - -double rnd_gauss (double mean, double sigma) { - - int i; - double y; - - i = 0; - while ((i == 0) || (i == NGAUSS - 1)) { - y = drand48(); - i = MIN (NGAUSS - 1, MAX (0, NGAUSS*y)); - } - y = gaussint[i]*sigma + mean; - - return (y); -} Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile (revision 39266) @@ -25,4 +25,5 @@ $(SRC)/cneedles.$(ARCH).o \ $(SRC)/cplot.$(ARCH).o \ +$(SRC)/crotation.$(ARCH).o \ $(SRC)/csystem.$(ARCH).o \ $(SRC)/ctimes.$(ARCH).o \ @@ -37,4 +38,5 @@ $(SRC)/fitpm_irls.$(ARCH).o \ $(SRC)/fitplx_irls.$(ARCH).o \ +$(SRC)/astrom_ops.$(ARCH).o \ $(SRC)/fixwrap.$(ARCH).o \ $(SRC)/fixcols.$(ARCH).o \ Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/astrom_ops.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/astrom_ops.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/astrom_ops.c (revision 39266) @@ -0,0 +1,72 @@ +# include "astro.h" +# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ + +/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { + + double n = mjd - J2000; // day number relative to standard epoch + double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) + double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly + + *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees + *beta = 0.0; // approx latitude + *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees + *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU + return TRUE; +} + +/* given RA, DEC, Time, calculate the parallax factor */ +// RA,DEC are decimal degrees +// Time is MJD +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { + + double lambda, beta, epsilon, Radius; + + // Time must be mjd + sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius); + + double lambda_rad = lambda*RAD_DEG; + double epsilon_rad = epsilon*RAD_DEG; + double RA_rad = RA*RAD_DEG; + double DEC_rad = DEC*RAD_DEG; + + double x = Radius*cos(lambda_rad); + double y = Radius*cos(epsilon_rad)*sin(lambda_rad); + double z = Radius*sin(epsilon_rad)*sin(lambda_rad); + + *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); + *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); + + return TRUE; +} + +// allocate arrays but not the container +int PlxFitDataAlloc (PlxFitData *data, int N) { + + data->Npts = N; + ALLOCATE (data->X, double, N); + ALLOCATE (data->Y, double, N); + ALLOCATE (data->dX, double, N); + ALLOCATE (data->dY, double, N); + ALLOCATE (data->t, double, N); + ALLOCATE (data->pX, double, N); + ALLOCATE (data->pY, double, N); + ALLOCATE (data->Wx, double, N); + ALLOCATE (data->Wy, double, N); + ALLOCATE (data->index, int, N); + return TRUE; +} + +void PlxFitDataFree (PlxFitData *data) { + FREE (data->X); + FREE (data->Y); + FREE (data->dX); + FREE (data->dY); + FREE (data->t); + FREE (data->pX); + FREE (data->pY); + FREE (data->Wx); + FREE (data->Wy); + FREE (data->index); +} + Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdensify.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdensify.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdensify.c (revision 39266) @@ -84,5 +84,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); PutCoords (&newcoords, &bf[0].header); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdhistogram.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdhistogram.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdhistogram.c (revision 39266) @@ -101,5 +101,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0); + if (!CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); PutCoords (&newcoords, &bf[0].header); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/crotation.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/crotation.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/crotation.c (revision 39266) @@ -0,0 +1,52 @@ +# include "astro.h" + +int crotation (int argc, char **argv) { + + int i; + double X, Y, x, y; + opihi_flt *xptr, *yptr; + Vector *xvec, *yvec; + CoordTransform *transform; + + if (argc != 6) { + gprint (GP_ERR, "USAGE: crotation (phi) (Xo) (xo) X Y\n"); + return (FALSE); + } + + double phi = atof(argv[1]); + double Xo = atof(argv[2]); + double xo = atof(argv[3]); + char *Xname = argv[4]; + char *Yname = argv[5]; + + transform = AllocTransform (phi, Xo, xo); + if (SelectScalar (Xname, &X)) { + if (!SelectScalar (Yname, &Y)) return (FALSE); + + ApplyTransform (&x, &y, X, Y, transform); + + gprint (GP_LOG, "%10.6f %10.6f\n", x, y); + return (TRUE); + } + + /* find vectors */ + if ((xvec = SelectVector (Xname, OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (Yname, OLDVECTOR, TRUE)) == NULL) return (FALSE); + + if (xvec[0].Nelements != yvec[0].Nelements) { + gprint (GP_ERR, "vectors %s and %s not the same length\n", Xname, Yname); + return (FALSE); + } + + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + + xptr = xvec[0].elements.Flt; + yptr = yvec[0].elements.Flt; + + for (i = 0; i < xvec[0].Nelements; i++, xptr++, yptr++) { + ApplyTransform (xptr, yptr, *xptr, *yptr, transform); + } + + return (TRUE); +} Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitplx.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitplx.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitplx.c (revision 39266) @@ -1,49 +1,3 @@ # include "astro.h" -# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ - -# define ESCAPE(MSG,...) { \ - gprint (GP_ERR, MSG, __VA_ARGS__); \ - return FALSE; } - -typedef struct { - double *X; - double *Y; - double *t; - double *pX; - double *pY; - double *dX; - double *dY; - int *index; - int Npts; -} PlxFitData; - -typedef struct { - double Ro, dRo; - double Do, dDo; - - double uR, duR; - double uD, duD; - - double p, dp; - - double chisq; - int Nfit; - int getChisq; -} PlxFit; - -int VectorRobustStats (Vector *vector, double *median, double *sigma); -double VectorFractionInterpolate (double *values, float fraction, int Npts); - -int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal); -int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean); -int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE); - -int PlxFitDataAlloc (PlxFitData *data, int N); -void PlxFitDataFree (PlxFitData *data); -int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt); - -int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); -int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); -int ParFactor (double *pR, double *pD, double RA, double DEC, double Time); int fitplx (int argc, char **argv) { @@ -100,4 +54,5 @@ gprint (GP_ERR, " -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n"); gprint (GP_ERR, " -dPsig vec : save path deviations in vec\n"); + gprint (GP_ERR, " -mask mask : excluded points are marked with a 0 mask value\n"); return (FALSE); } @@ -148,8 +103,9 @@ PlxFit fit; memset (&fit, 0, sizeof(PlxFit)); - // determine dPsig for detections based on Noutlier attempts + // determine dPsig for detections based on Noutlier attempts (mask is updated to mark the bad stars, mask == 0) if (Noutlier) { int clipRetry = TRUE; for (i = 0; clipRetry && (i < 3); i++) { + // XXX NOTE This will segfault is mask is not supplied: clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE); @@ -163,5 +119,5 @@ for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) { int n = fitdata.index[i]; - int maskValue = mask ? mask[n] : 0; + int maskValue = mask ? mask[n] : 1; 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]); } @@ -401,69 +357,4 @@ fit[0].Nfit = Npts; return (TRUE); -} - -/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ -int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { - - double n = mjd - J2000; // day number relative to standard epoch - double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) - double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly - - *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees - *beta = 0.0; // approx latitude - *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees - *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU - return TRUE; -} - -/* given RA, DEC, Time, calculate the parallax factor */ -// RA,DEC are decimal degrees -// Time is MJD -int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { - - double lambda, beta, epsilon, Radius; - - // Time must be mjd - sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius); - - double lambda_rad = lambda*RAD_DEG; - double epsilon_rad = epsilon*RAD_DEG; - double RA_rad = RA*RAD_DEG; - double DEC_rad = DEC*RAD_DEG; - - double x = Radius*cos(lambda_rad); - double y = Radius*cos(epsilon_rad)*sin(lambda_rad); - double z = Radius*sin(epsilon_rad)*sin(lambda_rad); - - *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); - *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); - - return TRUE; -} - -// allocate arrays but not the container -int PlxFitDataAlloc (PlxFitData *data, int N) { - - data->Npts = N; - ALLOCATE (data->X, double, N); - ALLOCATE (data->Y, double, N); - ALLOCATE (data->dX, double, N); - ALLOCATE (data->dY, double, N); - ALLOCATE (data->t, double, N); - ALLOCATE (data->pX, double, N); - ALLOCATE (data->pY, double, N); - ALLOCATE (data->index, int, N); - return TRUE; -} - -void PlxFitDataFree (PlxFitData *data) { - FREE (data->X); - FREE (data->Y); - FREE (data->dX); - FREE (data->dY); - FREE (data->t); - FREE (data->pX); - FREE (data->pY); - FREE (data->index); } @@ -481,4 +372,8 @@ tgt->pX[i] = src->pX[N]; tgt->pY[i] = src->pY[N]; + + // *** make this optional? + tgt->Wx[i] = src->Wx[N]; + tgt->Wy[i] = src->Wy[N]; } return TRUE; @@ -539,4 +434,7 @@ pYmin = MIN (pYmin, fitdata->pY[Nsubset]); pYmax = MAX (pYmax, fitdata->pY[Nsubset]); + + fitdata->Wx[Nsubset] = 1.0; + fitdata->Wy[Nsubset] = 1.0; fitdata->index[Nsubset] = i; Nsubset++; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm.c (revision 39266) @@ -1,21 +1,3 @@ # include "astro.h" -# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ - -# define ESCAPE(MSG,...) { \ - gprint (GP_ERR, MSG, __VA_ARGS__); \ - return FALSE; } - -typedef struct { - double Ro, dRo; - double Do, dDo; - - double uR, duR; - double uD, duD; - - double chisq; - int Nfit; -} PMFit; - -int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE); int fitpm (int argc, char **argv) { @@ -122,5 +104,5 @@ } - PMFit fit; + PlxFit fit; if (!FitPMonly (&fit, X, dX, Y, dY, t, n, VERBOSE)) { return FALSE; @@ -158,5 +140,5 @@ /* do we want an init function which does the alloc and a clear function to free? */ -int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { +int FitPMonly (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { int i; @@ -241,4 +223,7 @@ fit[0].duD = sqrt(A[3][3]); + fit[0].p = 0.0; + fit[0].dp = NAN; + // add up the chi square for the fit chisq = 0.0; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 39266) @@ -1,28 +1,3 @@ # include "astro.h" -# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ - -# define ESCAPE(MSG,...) { \ - gprint (GP_ERR, MSG, __VA_ARGS__); \ - return FALSE; } - -typedef struct { - double Ro, dRo; - double Do, dDo; - - double uR, duR; - double uD, duD; - - double chisq; - int Nfit; -} PMFit_IRLS; - -int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE); -int IRLS_converged (PMFit_IRLS *fit); -int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts, - double **A, double **B, int VERBOSE); -double weight_cauchy (double x); -double dpsi_cauchy (double x); -double MAD(double *in, int N); - int fitpm_irls (int argc, char **argv) { @@ -129,5 +104,5 @@ } - PMFit_IRLS fit; + PlxFit fit; if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) { return FALSE; @@ -165,5 +140,5 @@ /* do we want an init function which does the alloc and a clear function to free? */ -int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { +int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { int i,j; @@ -222,5 +197,5 @@ // Solve OLS equation - if (!weighted_LS(T,X,Wx,Y,Wy,Npts, + if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts, A,B,VERBOSE)) { // Handle fail case @@ -267,5 +242,5 @@ // Solve - if (!weighted_LS(T,X,Wx,Y,Wy,Npts, + if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts, A,B,VERBOSE)) { // Handle fail case @@ -284,5 +259,5 @@ u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i])); } - sigma_hat = MAD(u,Npts) / 0.6745; + sigma_hat = MedianAbsDeviation(u,Npts) / 0.6745; // Check convergence @@ -310,4 +285,6 @@ double Sum_Wx, Sum_Wy; + Sum_Wx = 0.0; + Sum_Wy = 0.0; ax = 0.0; ay = 0.0; bx = 0.0; by = 0.0; @@ -387,56 +364,5 @@ } - -double weight_cauchy (double x) { - double r = x / 2.385; - return (1.0 / (1.0 + SQ(r))); -} - -// dpsi = (d/dx) (x * weight(x)) -double dpsi_cauchy (double x) { - double r2 = SQ(x / 2.385); - return ((1.0 - r2) / (SQ(1 + r2))); -} - - -// median absolute deviation -// MAD = median(abs(x - median(x))) -double MAD(double *in, int N) { - double *x; - double median = 0.0; - int i; - - ALLOCATE(x,double,N); - for (i = 0; i < N; i++) { - x[i] = in[i]; - } - - dsort(x,N); - - if (N % 2) { - median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); - } else { - median = x[(int)(0.5*N)]; - } - - for (i = 0; i < N; i++ ) { - x[i] = fabs(x[i] - median); - } - - dsort(x,N); - - if (N % 2) { - median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); - } else { - median = x[(int)(0.5*N)]; - } - - return(median); -} - - - -int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts, - double **A, double **B, int VERBOSE) { +int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE) { int i,j; @@ -497,2 +423,50 @@ return TRUE; } + +double weight_cauchy (double x) { + double r = x / 2.385; + return (1.0 / (1.0 + SQ(r))); +} + +// dpsi = (d/dx) (x * weight(x)) +double dpsi_cauchy (double x) { + double r2 = SQ(x / 2.385); + return ((1.0 - r2) / (SQ(1 + r2))); +} + + +// median absolute deviation +// MAD = median(abs(x - median(x))) +double MedianAbsDeviation(double *in, int N) { + double *x; + double median = 0.0; + int i; + + ALLOCATE(x,double,N); + for (i = 0; i < N; i++) { + x[i] = in[i]; + } + + dsort(x,N); + + if (N % 2) { + median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); + } else { + median = x[(int)(0.5*N)]; + } + + for (i = 0; i < N; i++ ) { + x[i] = fabs(x[i] - median); + } + + dsort(x,N); + + if (N % 2) { + median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); + } else { + median = x[(int)(0.5*N)]; + } + + return(median); +} + Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fiximage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fiximage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fiximage.c (revision 39266) @@ -26,5 +26,5 @@ gfits_free_matrix (&mask[0].matrix); gfits_free_header (&mask[0].header); - CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (mask[0].file, "(empty)"); memset (mask[0].matrix.buffer, 0, Nx*Ny*sizeof(float)); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/getcoords.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/getcoords.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/getcoords.c (revision 39266) @@ -73,9 +73,9 @@ gfits_free_matrix (&outC1[0].matrix); gfits_free_header (&outC1[0].header); - CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0); + if (!CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0)) return FALSE; gfits_free_matrix (&outC2[0].matrix); gfits_free_header (&outC2[0].header); - CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0); + if (!CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0)) return FALSE; float *valC1 = (float *) outC1[0].matrix.buffer; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/imfit.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/imfit.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/imfit.c (revision 39266) @@ -141,5 +141,5 @@ strcpy (out[0].file, "(empty)"); - CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0); + if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE; /* four panels: 1) raw image. 2) fit 3) raw - fit 4) ?? */ Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c (revision 39266) @@ -9,4 +9,5 @@ int cneedles PROTO((int, char **)); int cplot PROTO((int, char **)); +int crotation PROTO((int, char **)); int csystem PROTO((int, char **)); int ctimes PROTO((int, char **)); @@ -78,4 +79,5 @@ {1, "cneedles", cneedles, "plot vectors in sky coordinates"}, {1, "cplot", cplot, "plot vectors in sky coordinates"}, + {1, "crotation", crotation, "rotate in 3D"}, {1, "csystem", csystem, "convert between coordinate systems"}, {1, "ctimes", ctimes, "convert between time formats"}, Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/shimage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/shimage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/shimage.c (revision 39266) @@ -30,6 +30,6 @@ int Ny = src[0].header.Naxis[1]; - gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0); - gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); if (!CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE; + gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); if (!CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE; // l=0 allocates space for a single value Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/vshimage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/vshimage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/vshimage.c (revision 39266) @@ -32,8 +32,8 @@ int Ny = src[0].header.Naxis[1]; - gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0); - gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0); - gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0); - gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0); + gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); if (!CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE; + gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); if (!CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0)) return FALSE; + gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); if (!CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE; + gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); if (!CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0)) return FALSE; // l=0 allocates space for a single value Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/densify.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/densify.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/densify.c (revision 39266) @@ -96,5 +96,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/fft2d.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/fft2d.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/fft2d.c (revision 39266) @@ -58,6 +58,6 @@ /* fix up output headers (real) & allocate data buffer */ - CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0); - CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0)) return FALSE; + if (!CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0)) return FALSE; gfits_copy_header (&Ire[0].header, &Ore[0].header); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gaussdeviate.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gaussdeviate.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gaussdeviate.c (revision 39266) @@ -17,7 +17,7 @@ ResetVector (vec, OPIHI_FLT, Npts); - gauss_init (2048); + ohana_gaussdev_init (); for (i = 0; i < Npts; i++) { - vec[0].elements.Flt[i] = rnd_gauss (mean, sigma); + vec[0].elements.Flt[i] = ohana_gaussdev_rnd (mean, sigma); } return (TRUE); @@ -29,8 +29,10 @@ } -double int_gauss (int i); - int gaussintegral (int argc, char **argv) { + gprint (GP_ERR, "fix gaussintegral\n"); + return (FALSE); + +/* int i, Npts; Vector *vec; @@ -46,7 +48,7 @@ ResetVector (vec, OPIHI_FLT, Npts); - gauss_init (Npts); + ohana_gaussdev_init (); for (i = 0; i < Npts; i++) { - vec[0].elements.Flt[i] = int_gauss (i); + vec[0].elements.Flt[i] = gaussian_int (i); } return (TRUE); @@ -55,4 +57,5 @@ gprint (GP_ERR, "USAGE: gaussintegral Npts mean sigma\n"); return (FALSE); +*/ } Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gridify.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gridify.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/gridify.c (revision 39266) @@ -85,5 +85,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); } @@ -111,5 +111,5 @@ gfits_free_matrix (&ct[0].matrix); gfits_free_header (&ct[0].header); - CreateBuffer (ct, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (ct, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (ct[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mcreate.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mcreate.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mcreate.c (revision 39266) @@ -27,7 +27,7 @@ if (Nz) { - CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0); + if (!CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0)) return FALSE; } else { - CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0); + if (!CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0)) return FALSE; } return (TRUE); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/medimage_commands.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/medimage_commands.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/medimage_commands.c (revision 39266) @@ -82,5 +82,5 @@ gfits_free_matrix (&output->matrix); gfits_free_header (&output->header); - CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE; float *outvalue = (float *) output->matrix.buffer; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mslice.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mslice.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/mslice.c (revision 39266) @@ -35,5 +35,5 @@ gfits_free_matrix (&out[0].matrix); gfits_free_header (&out[0].header); - CreateBuffer (out, Nx, Ny, -32, 1.0, 0.0); + if (!CreateBuffer (out, Nx, Ny, -32, 1.0, 0.0)) return FALSE; float *inF = (float *) in[0].matrix.buffer + plane*Nx*Ny; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/read_vectors.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/read_vectors.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/read_vectors.c (revision 39266) @@ -348,4 +348,5 @@ } +# undef ESCAPE # define ESCAPE(...) { \ gprint (GP_ERR, __VA_ARGS__); \ Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/reindex.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/reindex.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/reindex.c (revision 39266) @@ -1,5 +1,3 @@ # include "data.h" - -# define ESCAPE(MSG,...){ gprint (GP_ERR, MSG, __VA_ARGS__); goto error; } int reindex (int argc, char **argv) { Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/shift.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/shift.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/shift.c (revision 39266) @@ -48,5 +48,5 @@ gfits_free_matrix (&out[0].matrix); gfits_free_header (&out[0].header); - CreateBuffer (out, nx, ny, -32, 0.0, 1.0); + if (!CreateBuffer (out, nx, ny, -32, 0.0, 1.0)) return FALSE; DXin = (dx < 0) ? -dx : 0; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/squash3d.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/squash3d.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/squash3d.c (revision 39266) @@ -46,5 +46,5 @@ int oNx = iNz; int oNy = iNy; - CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0); + if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE; float *oBuf = (float *) tgt[0].matrix.buffer; @@ -70,5 +70,5 @@ int oNx = iNx; int oNy = iNz; - CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0); + if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE; float *oBuf = (float *) tgt[0].matrix.buffer; @@ -94,5 +94,5 @@ int oNx = iNx; int oNy = iNy; - CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0); + if (!CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0)) return FALSE; float *oBuf = (float *) tgt[0].matrix.buffer; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/tdhistogram.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/tdhistogram.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/tdhistogram.c (revision 39266) @@ -140,5 +140,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0); + if (!CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/vgrid.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/vgrid.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.data/vgrid.c (revision 39266) @@ -43,5 +43,5 @@ gfits_free_matrix (&bf[0].matrix); gfits_free_header (&bf[0].header); - CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (bf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dimm/camera_cmds.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dimm/camera_cmds.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dimm/camera_cmds.c (revision 39266) @@ -107,5 +107,5 @@ gfits_free_matrix (&buf[0].matrix); gfits_free_header (&buf[0].header); - CreateBuffer (buf, dx, dy, -32, 0.0, 1.0); + if (!CreateBuffer (buf, dx, dy, -32, 0.0, 1.0)) return FALSE; strcpy (buf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordimage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordimage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordimage.c (revision 39266) @@ -79,10 +79,10 @@ gfits_free_matrix (&bufX[0].matrix); gfits_free_header (&bufX[0].header); - CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + if (!CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here strcpy (bufX[0].file, "(empty)"); gfits_free_matrix (&bufY[0].matrix); gfits_free_header (&bufY[0].header); - CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + if (!CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here strcpy (bufY[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordmosaic.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordmosaic.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/coordmosaic.c (revision 39266) @@ -88,10 +88,10 @@ gfits_free_matrix (&bufX[0].matrix); gfits_free_header (&bufX[0].header); - CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + if (!CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here strcpy (bufX[0].file, "(empty)"); gfits_free_matrix (&bufY[0].matrix); gfits_free_header (&bufY[0].header); - CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + if (!CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0)) return FALSE; // initialized to 0.0 here strcpy (bufY[0].file, "(empty)"); @@ -224,4 +224,5 @@ } +# undef ESCAPE # define ESCAPE(MSG) { \ if (src) fclose (src); \ Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/fitcolors.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/fitcolors.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/fitcolors.c (revision 39266) @@ -145,5 +145,5 @@ gfits_free_matrix (&buf[0].matrix); gfits_free_header (&buf[0].header); - CreateBuffer (buf, NP, NP, -32, 0.0, 1.0); + if (!CreateBuffer (buf, NP, NP, -32, 0.0, 1.0)) return FALSE; strcpy (buf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/imphot.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/imphot.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/imphot.c (revision 39266) @@ -36,5 +36,5 @@ if (GreyScale) { if ((buf = SelectBuffer (bufname, ANYBUFFER, TRUE)) == NULL) return (FALSE); - CreateBuffer (buf, 100, 200, -32, 0.0, 1.0); + if (!CreateBuffer (buf, 100, 200, -32, 0.0, 1.0)) return FALSE; } Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/objectcoverage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/objectcoverage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/objectcoverage.c (revision 39266) @@ -144,5 +144,5 @@ gfits_free_matrix (&buf[0].matrix); gfits_free_header (&buf[0].header); - CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (buf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/skycoverage.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/skycoverage.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/dvo/skycoverage.c (revision 39266) @@ -201,5 +201,5 @@ gfits_free_matrix (&buf[0].matrix); gfits_free_header (&buf[0].header); - CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0); + if (!CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0)) return FALSE; strcpy (buf[0].file, "(empty)"); Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/astro.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/astro.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/astro.h (revision 39266) @@ -12,3 +12,59 @@ void FreeAstro (void); +typedef struct { + double *X; + double *Y; + double *t; + double *pX; + double *pY; + double *dX; + double *dY; + double *Wx; + double *Wy; + int *index; + int Npts; +} PlxFitData; + +typedef struct { + double Ro, dRo; + double Do, dDo; + + double uR, duR; + double uD, duD; + + double p, dp; + + double chisq; + int Nfit; + int getChisq; +} PlxFit; + +int VectorRobustStats (Vector *vector, double *median, double *sigma); +double VectorFractionInterpolate (double *values, float fraction, int Npts); + +int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal); +int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean); +int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE); + +int PlxFitDataAlloc (PlxFitData *data, int N); +void PlxFitDataFree (PlxFitData *data); +int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt); + +int FitPMonly (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE); +int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time); + +/***** */ + +int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE); +int 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); + +int 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); +int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE); + +double weight_cauchy (double x); +double dpsi_cauchy (double x); +double MedianAbsDeviation(double *in, int N); + # endif Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/shell.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/shell.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/include/shell.h (revision 39266) @@ -4,4 +4,6 @@ # ifndef SHELL_H # define SHELL_H + +# define ESCAPE(MSG,...) { gprint (GP_ERR, MSG, __VA_ARGS__); return FALSE; } # define ISVAR(a) (isalnum (a) || (a == ':') || (a == '_')) Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/Makefile (revision 39266) @@ -30,5 +30,4 @@ $(SDIR)/starfuncs.$(ARCH).o \ $(SDIR)/hermitian.$(ARCH).o \ -$(SDIR)/gaussian.$(ARCH).o \ $(SDIR)/graphtools.$(ARCH).o \ $(SDIR)/queues.$(ARCH).o \ @@ -36,4 +35,7 @@ $(SDIR)/open_kapa.$(ARCH).o \ $(SDIR)/style_args.$(ARCH).o + +# moved to libohana +# $(SDIR)/gaussian.$(ARCH).o \ # fix malloc Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/spline.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/spline.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.data/spline.c (revision 39266) @@ -32,4 +32,8 @@ float dx, a, b, value; + // saturate correction at high and low ends + if (X < x[0]) return y[0]; + if (X > x[N-1]) return y[N-1]; + /* find correct element in array (x must be sorted) */ lo = 0; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/BufferOps.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/BufferOps.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/BufferOps.c (revision 39266) @@ -92,9 +92,12 @@ gfits_free_matrix (&buf[0].matrix); gfits_free_header (&buf[0].header); - CreateBuffer (buf, Nx, Ny, bitpix, bzero, bscale); + if (!CreateBuffer (buf, Nx, Ny, bitpix, bzero, bscale)) return FALSE; return TRUE; } int CreateBuffer (Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale) { + + if (Nx < 0) { gprint (GP_ERR, "invalid negative matrix size Nx : %d\n", Nx); return FALSE; } + if (Ny < 0) { gprint (GP_ERR, "invalid negative matrix size Ny : %d\n", Ny); return FALSE; } /* store the default output values */ @@ -123,4 +126,8 @@ // buffer[x + Nx*y + Nx*Ny*z + ...] int CreateBuffer3D (Buffer *buf, int Nx, int Ny, int Nz, int bitpix, float bzero, float bscale) { + + if (Nx < 0) { gprint (GP_ERR, "invalid negative matrix size Nx : %d\n", Nx); return FALSE; } + if (Ny < 0) { gprint (GP_ERR, "invalid negative matrix size Ny : %d\n", Ny); return FALSE; } + if (Nz < 0) { gprint (GP_ERR, "invalid negative matrix size Nz : %d\n", Nz); return FALSE; } /* store the default output values */ Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/ListOps.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/ListOps.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/lib.shell/ListOps.c (revision 39266) @@ -212,4 +212,5 @@ if (temp != NULL) { if (!strcmp (temp, "-x")) goto escape; + if (!strcmp (temp, "-glob")) goto escape; if (!strcmp (temp, "-split")) goto escape; if (!strcmp (temp, "-splitbychar")) goto escape; Index: branches/eam_branches/ipp-20151113/Ohana/src/opihi/mana/simsignal.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/opihi/mana/simsignal.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/opihi/mana/simsignal.c (revision 39266) @@ -6,5 +6,4 @@ float *buf; double cvalue, dvalue, sigma, SN, period; - double rnd_gauss (); if (argc != 5) { @@ -28,5 +27,5 @@ ivalue = scale * 0.5 * (sin (i*2*M_PI/period) + 1) + 0.5; /* - dvalue = rnd_gauss (cvalue, sigma); + dvalue = ohana_gaussdev_rnd (cvalue, sigma); cvalue = (dvalue + range) / (2.0*range); dvalue = MAX (0, MIN (0.99999, cvalue)); @@ -44,5 +43,5 @@ 8 bit = 2^8 - gauss_init (2*scale); + ohana_gaussdev_init (); sigma = 2.0 / SN; range = 1 + 5*sigma; Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/Makefile (revision 39266) @@ -10,4 +10,6 @@ MAN = $(HOME)/doc INC = $(HOME)/include +TESTDIR = $(HOME)/test +TESTBIN = $(HOME)/test include ../../Makefile.Common @@ -22,6 +24,4 @@ install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client -# $(DESTBIN)/testparallax - RELASTRO = \ $(SRC)/ConfigInit.$(ARCH).o \ @@ -29,8 +29,11 @@ $(SRC)/FitMosaic.$(ARCH).o \ $(SRC)/FitSimple.$(ARCH).o \ -$(SRC)/FitAstromOps.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPM_IRLS.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPMandPar_IRLS.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -103,8 +106,11 @@ $(SRC)/ConfigInit.$(ARCH).o \ $(SRC)/FitSimple.$(ARCH).o \ -$(SRC)/FitAstromOps.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPM_IRLS.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPMandPar_IRLS.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -159,2 +165,19 @@ $(TESTPAR): $(INC)/relastro.h $(BIN)/testparallax.$(ARCH): $(TESTPAR) + +fitpm: $(BIN)/fitpm.$(ARCH) + +FITPM = \ +$(SRC)/fitpm.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPM_IRLS.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPMandPar_IRLS.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o + +$(FITPM): $(INC)/relastro.h +$(BIN)/fitpm.$(ARCH): $(FITPM) Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/include/relastro.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/include/relastro.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/include/relastro.h (revision 39266) @@ -134,7 +134,14 @@ double **A; double **B; + double **Cov; + double *Beta; + double *Beta_prev; int Nterms; + int getChisq; + int getError; } FitAstromData; +// XXX do we need doubles for all of these? I actually only have of order 100 of these +// allocated at a time, so size is not an issue. typedef struct { double X, dX; @@ -143,4 +150,7 @@ double D, dD; double T, dT; + double Wx, Wy; + double rx, ry; + double u; double pR; double pD; @@ -148,4 +158,5 @@ double C_red; int measure; + int mask; } FitAstromPoint; @@ -164,4 +175,5 @@ FitAstromPoint *points; FitAstromPoint *sample; + FitAstromPoint *nomask; int Npoints; int NpointsAlloc; @@ -723,4 +735,6 @@ int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode); int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints); +int BootstrapSaveUnmasked (FitAstromPoint *nomask, FitAstromPoint *points, int Npoints); + int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog); int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange); @@ -739,2 +753,13 @@ int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average); void AstromErrorSetLoop (int Nloop, int isImageMode); + +int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE); +int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE); + +double MedianAbsDeviation(FitAstromPoint *points, int Npoints); + +int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE); +int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE); + +double weight_cauchy (double x); +double dpsi_cauchy (double x); Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BootstrapOps.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BootstrapOps.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BootstrapOps.c (revision 39266) @@ -0,0 +1,116 @@ +# include "relastro.h" +# define PAR_TOOFEW 5 + +int BootstrapSaveUnmasked (FitAstromPoint *nomask, FitAstromPoint *points, int Npoints) { + int i, N; + + // I need to draw Npoints random entries from 'points' with replacement: + N = 0; + for (i = 0; i < Npoints; i++) { + if (points[i].mask) continue; + nomask[N] = points[i]; + N ++; + } + return N; +} + +int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) { + int i; + + // I need to draw Npoints random entries from 'points' with replacement: + for (i = 0; i < Npoints; i++) { + int N = Npoints * drand48(); + sample[i] = points[N]; + } + return TRUE; +} + +// calculate mean and sigma points for the 5 fit parameter +int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) { + + // generate a histogram for the selected element + double *values = NULL; + ALLOCATE (values, double, Nfit); + + int i; + + for (i = 0; i < Nfit; i++) { + switch (mode) { + case FIT_RESULT_RA: + values[i] = fit[i].Ro; + break; + case FIT_RESULT_DEC: + values[i] = fit[i].Do; + break; + case FIT_RESULT_uR: + values[i] = fit[i].uR; + break; + case FIT_RESULT_uD: + values[i] = fit[i].uD; + break; + case FIT_RESULT_PLX: + values[i] = fit[i].p; + break; + default: + myAbort ("invalid option"); + } + } + + dsort (values, Nfit); + +# if (0) + double median; + if (Nfit % 2) { + int Ncenter = Nfit / 2; + median = values[Ncenter]; + } else { + int Ncenter = Nfit / 2 - 1; + median = 0.5*(values[Ncenter] + values[Ncenter + 1]); + } +# endif + + double Slo = VectorFractionInterpolate (values, 0.158655, Nfit); + double Shi = VectorFractionInterpolate (values, 0.841345, Nfit); + double sigma = (Shi - Slo) / 2.0; + + switch (mode) { + case FIT_RESULT_RA: + // result->Ro = median; + result->dRo = sigma; + break; + case FIT_RESULT_DEC: + // result->Do = median; + result->dDo = sigma; + break; + case FIT_RESULT_uR: + // result->uR = median; + result->duR = sigma; + break; + case FIT_RESULT_uD: + // result->uD = median; + result->duD = sigma; + break; + case FIT_RESULT_PLX: + // result->p = median; + result->dp = sigma; + break; + default: + myAbort ("invalid option"); + } + + return TRUE; +} + +double VectorFractionInterpolate (double *values, float fraction, int Npts) { + + float F = fraction * Npts; + int N = fraction * Npts; + + if (N < 0 ) return NAN; + if (N >= Npts - 2) return NAN; + + // interpolate between N,N+1 + + double S = (F - N) * (values[N+1] - values[N]) + values[N]; + return S; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BrightCatalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/BrightCatalog.c (revision 39266) @@ -67,6 +67,4 @@ GET_COLUMN(R, "RA", double); GET_COLUMN(D, "DEC", double); - GET_COLUMN(RoffGAL, "ROFF_GAL", float); - GET_COLUMN(DoffGAL, "DOFF_GAL", float); GET_COLUMN(M, "MAG_SYS", float); GET_COLUMN(Mcal, "MAG_CAL", float); @@ -94,6 +92,4 @@ measure[i].R = R[i]; measure[i].D = D[i]; - measure[i].RoffGAL = RoffGAL[i]; - measure[i].DoffGAL = DoffGAL[i]; measure[i].M = M[i]; measure[i].Mcal = Mcal[i]; @@ -121,6 +117,4 @@ free (R ); free (D ); - free (RoffGAL ); - free (DoffGAL ); free (M ); free (Mcal ); @@ -362,6 +356,4 @@ gfits_define_bintable_column (&theader, "D", "RA", "ra", "degrees", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "DEC", "dec", "degrees", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "ROFF_GAL", "ra offset", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "DOFF_GAL", "dec offset", "arcsec", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); @@ -382,13 +374,9 @@ // generate the output array that carries the data - ohana_memcheck (1); gfits_create_table (&theader, &ftable); - ohana_memcheck (1); // create intermediate storage arrays double *R ; ALLOCATE (R , double, catalog->Nmeasure); double *D ; ALLOCATE (D , double, catalog->Nmeasure); - float *RoffGAL ; ALLOCATE (RoffGAL , float, catalog->Nmeasure); - float *DoffGAL ; ALLOCATE (DoffGAL , float, catalog->Nmeasure); float *M ; ALLOCATE (M , float, catalog->Nmeasure); float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); @@ -407,5 +395,4 @@ int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); - ohana_memcheck (1); // assign the storage arrays @@ -414,6 +401,4 @@ R[i] = measure[i].R ; D[i] = measure[i].D ; - RoffGAL[i] = measure[i].RoffGAL ; - DoffGAL[i] = measure[i].DoffGAL ; M[i] = measure[i].M ; Mcal[i] = measure[i].Mcal ; @@ -433,13 +418,8 @@ photcode[i] = measure[i].photcode ; } - ohana_memcheck (1); // add the columns to the output array gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Nmeasure); - ohana_memcheck (1); gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Nmeasure); - ohana_memcheck (1); - gfits_set_bintable_column (&theader, &ftable, "ROFF_GAL", RoffGAL, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "DOFF_GAL", DoffGAL, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); @@ -461,6 +441,4 @@ free (R ); free (D ); - free (RoffGAL ); - free (DoffGAL ); free (M ); free (Mcal ); Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitAstromOps.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitAstromOps.c (revision 39266) @@ -2,4 +2,64 @@ // I am modifying FitPM with an eye to (a) threaded operations and (b) bootstrap resampling tests. + +int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) { + + int k; + + int Npoints = fitStats->Npoints; + FitAstromPoint *points = fitStats->points; + + // find Tmin & Tmax from the list of accepted measurements + double Tmin = points[0].T; + double Tmax = points[0].T; + double pRmin = +2.0; + double pRmax = -2.0; + double pDmin = +2.0; + double pDmax = -2.0; + + *Tmean = 0.0; + + double Tsum = 0.0; + double Wsum = 0.0; + for (k = 0; k < Npoints; k++) { + Tmin = MIN(Tmin, points[k].T); + Tmax = MAX(Tmax, points[k].T); + + float wx = 1.0 / SQ(points[k].dX); + + Tsum += points[k].T * wx; + Wsum += wx; + + // at this point, T is in years since J2000 + ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T); + pRmin = MIN (pRmin, points[k].pR); + pRmax = MAX (pRmax, points[k].pR); + pDmin = MIN (pDmin, points[k].pD); + pDmax = MAX (pDmax, points[k].pD); + } + *Trange = Tmax - Tmin; + + // mean epoch + *Tmean = Tsum / Wsum; + + // for HIGH_SPEED, just use the center of the range + if (RELASTRO_OP == OP_HIGH_SPEED) { + *Tmean = 0.5*(Tmax - Tmin); + } + + *parRange = hypot (pRmax - pRmin, pDmax - pDmin); + + /* we need to do the fit in a locally linear space; choose a ref coordinate */ + fitStats->coords.crval1 = points[0].R; + fitStats->coords.crval2 = points[0].D; + + // project all of the R,D coordinates to a plane centered on this coordinate. set + // the times to be relative to Tmean + for (k = 0; k < Npoints; k++) { + RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords); + points[k].T -= *Tmean; + } + return TRUE; +} FitStats *FitStatsInit (int Nmax, int Nboot) { @@ -27,4 +87,5 @@ fitStats->points = NULL; fitStats->sample = NULL; + fitStats->nomask = NULL; fitStats->Npoints = 0; fitStats->NpointsAlloc = Nmax; @@ -32,4 +93,5 @@ ALLOCATE (fitStats->points, FitAstromPoint, Nmax); ALLOCATE (fitStats->sample, FitAstromPoint, Nmax); + ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax); } @@ -79,4 +141,5 @@ FREE (fitStats->points); FREE (fitStats->sample); + FREE (fitStats->nomask); FitAstromDataFree (fitStats->fitdataPos); @@ -93,7 +156,14 @@ /* do I need to do this as 2 2x2 matrix equations? */ - fit->A = array_init (Nterms, Nterms); - fit->B = array_init (Nterms, 1); + fit->B = array_init (Nterms, 1); + fit->A = array_init (Nterms, Nterms); + fit->Cov = array_init (Nterms, Nterms); + + ALLOCATE (fit->Beta, double, Nterms); + ALLOCATE (fit->Beta_prev, double, Nterms); fit->Nterms = Nterms; + + fit->getChisq = TRUE; + fit->getError = TRUE; return fit; @@ -106,4 +176,9 @@ array_free (fit->A, fit->Nterms); array_free (fit->B, fit->Nterms); + array_free (fit->Cov, fit->Nterms); + + free (fit->Beta); + free (fit->Beta_prev); + free (fit); return; @@ -126,4 +201,13 @@ object->C_red = 0.0; object->measure= -1; + + object->Wx = 1.0; + object->Wy = 1.0; + + object->rx = 0.0; + object->ry = 0.0; + object->u = 0.0; + + object->mask = 0; // keep point if mask == 0 return; } @@ -147,2 +231,48 @@ return; } + +double weight_cauchy (double x) { + double r = x / 2.385; + return (1.0 / (1.0 + SQ(r))); +} + +// dpsi = (d/dx) (x * weight(x)) +double dpsi_cauchy (double x) { + double r2 = SQ(x / 2.385); + return ((1.0 - r2) / (SQ(1 + r2))); +} + + +// median absolute deviation +// MAD = median(abs(x - median(x))) +double MedianAbsDeviation(FitAstromPoint *points, int Npoints) { + + double *x; + double median = 0.0; + int i; + + ALLOCATE(x, double, Npoints); + for (i = 0; i < Npoints; i++) { + x[i] = points[i].u; + } + dsort(x, Npoints); + + if (Npoints % 2) { + median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]); + } else { + median = x[(int)(0.5*Npoints)]; + } + + for (i = 0; i < Npoints; i++ ) { + x[i] = fabs(x[i] - median); + } + dsort(x, Npoints); + + if (Npoints % 2) { + median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]); + } else { + median = x[(int)(0.5*Npoints)]; + } + + return median; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPM_IRLS.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPM_IRLS.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPM_IRLS.c (revision 39266) @@ -0,0 +1,248 @@ +# include "relastro.h" + +// These should probably be tunable: +# define MAX_ITERATIONS 10 +# define FIT_TOLERANCE 1e-4 +# define FLT_TOLERANCE 1e-6 +# define WEIGHT_THRESHOLD 0.3 + +int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { + + int i,j; + + int Ndof = 2 * Npoints - data->Nterms; + + // Convert the measurement errors into initial weights. + for (i = 0; i < Npoints; i++) { + points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX); + points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY); + } + + // Solve OLS equation + if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) { + myAbort ("handle failures, please!"); + return(FALSE); + } + + // Calculate r vector of residuals and least squares sigma + double sigma_ols = 0.0; + for (i = 0; i < Npoints; i++) { + points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); + points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); + sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); + } + sigma_ols = sqrt(sigma_ols / (float)Ndof); + + // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) + for (i = 0; i < data->Nterms; i++) { + for (j = 0; j < data->Nterms; j++) { + data->Cov[i][j] = data->A[i][j]; + } + data->Beta[i] = data->B[i][0]; + } + + // Iteratively reweight and solve + double sigma_hat = 0.0; // save for the error model + int converged = FALSE; + int iterations = 0; + + // modify the weight based on the distance from the previous fit. try up to MAX_ITERATIONS. + // at the end "fit", has the last fit parameters + for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) { + // Save Beta. + for (i = 0; i < data->Nterms; i ++) { + data->Beta_prev[i] = data->Beta[i]; + } + + // Assign weights based on the deviation + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + } + + // Solve with the new weights + if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) { + myAbort ("handle failures, please!"); + return(FALSE); + } + + // store the new Beta. + for (i = 0; i < data->Nterms; i++) { + data->Beta[i] = data->B[i][0]; + } + + // calculate the residuals: + for (i = 0; i < Npoints; i++) { + points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); + points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); + points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); + } + sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + + // Check convergence + converged = TRUE; + for (i = 0; i < data->Nterms; i++) { + // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good + if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && + (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) { + converged = FALSE; + } + } + } + if (!converged) { + myAbort ("raise a warning on non-convergence"); + } + + // calculate the weight thresholds to mask the bad points: + double Sum_Wx = 0.0; + double Sum_Wy = 0.0; + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + + Sum_Wx += points[i].Wx; + Sum_Wy += points[i].Wy; + } + double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); + double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); + + // set a mask (which can be used by the bootstrap resampling analysis) + for (i = 0; i < Npoints; i++) { + // keep if either is above threshold? + // drop if either is below threshold? + // points are marked as keep by default + if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { + points[i].mask = 1; // keep point if mask == 0 + } + } + + // this section calculates the formal error on the weighted fit using the covariance values + // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37 + if (data->getError) { + double ax = 0.0, ay = 0.0; + double bx = 0.0, by = 0.0; + + for (i = 0; i < Npoints; i++) { + ax += dpsi_cauchy(points[i].rx / points[i].dX); + ay += dpsi_cauchy(points[i].ry / points[i].dY); + + bx += SQ(points[i].Wx); + by += SQ(points[i].Wy); + } + ax /= 1.0 * Npoints; // mean(psi_dot(r)) + ay /= 1.0 * Npoints; + bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p)) + by /= 1.0 * (Npoints - data->Nterms); + + double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax; + double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay; + + double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax; + double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay; + + // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1) + double sigma_final_x = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms))); + double sigma_final_y = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms))); + + fit[0].dRo = sqrt(data->Cov[0][0]); + fit[0].duR = sqrt(data->Cov[1][1]); + fit[0].dDo = sqrt(data->Cov[2][2]); + fit[0].duD = sqrt(data->Cov[3][3]); + + fit[0].dRo *= sigma_final_x; + fit[0].duR *= sigma_final_x; + fit[0].dDo *= sigma_final_y; + fit[0].duD *= sigma_final_y; + } + + // (optionally) add up the chi square for the fit, only counting the unmasked points + double chisq = 0.0; + fit[0].Nfit = 0; + for (i = 0; i < Npoints; i++) { + if (points[i].mask) continue; + fit[0].Nfit ++; + + if (data->getChisq) { + double Xf = fit[0].Ro + fit[0].uR*points[i].T; + double Yf = fit[0].Do + fit[0].uD*points[i].T; + double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); + double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); + chisq += SQ(points[i].X - Xf) * wx; + chisq += SQ(points[i].Y - Yf) * wy; + } + } + + // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) + fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); + + return (TRUE); +} + +int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { + + int i; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + + for (i = 0; i < Npoints; i++) { + wx = points[i].Wx; + wy = points[i].Wy; + + Wx += wx; + Wy += wy; + + double TWx = points[i].T*wx; + double TWy = points[i].T*wy; + + Tx += TWx; + Ty += TWy; + + Tx2 += points[i].T*TWx; + Ty2 += points[i].T*TWy; + + Xs += points[i].X*wx; + Ys += points[i].Y*wy; + + XT += points[i].X*TWx; + YT += points[i].Y*TWy; + } + + // X^T W X + data->A[0][0] = Wx; + data->A[0][1] = Tx; + + data->A[1][0] = Tx; + data->A[1][1] = Tx2; + data->A[2][2] = Wy; + data->A[2][3] = Ty; + data->A[3][2] = Ty; + data->A[3][3] = Ty2; + + // X^T W Y + data->B[0][0] = Xs; + data->B[1][0] = XT; + data->B[2][0] = Ys; + data->B[3][0] = YT; + + if (!dgaussjordan (data->A, data->B, 4, 1)) { +# if (DEBUG) + if (VERBOSE) fprintf (stderr, "error in fit\n"); + int j; + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + fprintf (stderr, "%e ", data->A[i][j]); + } + fprintf (stderr, " : %e\n", data->B[i][0]); + } +# endif + return FALSE; + } + + fit->Ro = data->B[0][0]; + fit->uR = data->B[1][0]; + fit->Do = data->B[2][0]; + fit->uD = data->B[3][0]; + fit->p = 0.0; + + return TRUE; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar.c (revision 39266) @@ -15,4 +15,7 @@ for (i = 0; i < Npoints; i++) { + if (points[i].mask) continue; // respect the mask if set + fit->Nfit ++; + /* handle case where dX or dY = 0.0 */ wx = 1.0 / SQ(points[i].dX); @@ -89,6 +92,4 @@ fit->dp = sqrt(data->A[4][4]); - fit->Nfit = Npoints; - return (TRUE); } Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar_IRLS.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar_IRLS.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/FitPMandPar_IRLS.c (revision 39266) @@ -0,0 +1,276 @@ +# include "relastro.h" + +// These should probably be tunable: +# define MAX_ITERATIONS 10 +# define FIT_TOLERANCE 1e-4 +# define FLT_TOLERANCE 1e-6 +# define WEIGHT_THRESHOLD 0.3 + +int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { + + int i,j; + + int Ndof = 2 * Npoints - data->Nterms; + + // Convert the measurement errors into initial weights. + for (i = 0; i < Npoints; i++) { + points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX); + points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY); + } + + // Solve OLS equation + if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) { + myAbort ("handle failures, please!"); + return(FALSE); + } + + // Calculate r vector of residuals and least squares sigma + double sigma_ols = 0.0; + for (i = 0; i < Npoints; i++) { + points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p); + points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p); + sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); + } + sigma_ols = sqrt(sigma_ols / (float)Ndof); + + // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) + for (i = 0; i < data->Nterms; i++) { + for (j = 0; j < data->Nterms; j++) { + data->Cov[i][j] = data->A[i][j]; + } + data->Beta[i] = data->B[i][0]; + } + + // Iteratively reweight and solve + double sigma_hat = 0.0; // save for the error model + int converged = FALSE; + int iterations = 0; + + // modify the weight based on the distance from the previous fit. try up to MAX_ITERATIONS. + // at the end "fit", has the last fit parameters + for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) { + // Save Beta. + for (i = 0; i < data->Nterms; i ++) { + data->Beta_prev[i] = data->Beta[i]; + } + + // Assign weights based on the deviation + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + } + + // Solve with the new weights + if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) { + myAbort ("handle failures, please!"); + return(FALSE); + } + + // store the new Beta. + for (i = 0; i < data->Nterms; i++) { + data->Beta[i] = data->B[i][0]; + } + + // calculate the residuals: + for (i = 0; i < Npoints; i++) { + points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p); + points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p); + points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); + } + sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + + // Check convergence + converged = TRUE; + for (i = 0; i < data->Nterms; i++) { + // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good + if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && + (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) { + converged = FALSE; + } + } + } + if (!converged) { + fprintf (stderr, "raise a warning on non-convergence\n"); + } + + // calculate the weight thresholds to mask the bad points: + double Sum_Wx = 0.0; + double Sum_Wy = 0.0; + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + + Sum_Wx += points[i].Wx; + Sum_Wy += points[i].Wy; + } + double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); + double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); + + // set a mask (which can be used by the bootstrap resampling analysis) + for (i = 0; i < Npoints; i++) { + // keep if either is above threshold? + // drop if either is below threshold? + // points are marked as keep by default + if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { + points[i].mask = 1; // keep point if mask == 0 + } + } + + // this section calculates the formal error on the weighted fit using the covariance values + // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37 + if (data->getError) { + double ax = 0.0, ay = 0.0; + double bx = 0.0, by = 0.0; + + for (i = 0; i < Npoints; i++) { + ax += dpsi_cauchy(points[i].rx / points[i].dX); + ay += dpsi_cauchy(points[i].ry / points[i].dY); + + bx += SQ(points[i].Wx); + by += SQ(points[i].Wy); + } + ax /= 1.0 * Npoints; // mean(psi_dot(r)) + ay /= 1.0 * Npoints; + bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p)) + by /= 1.0 * (Npoints - data->Nterms); + + double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax; + double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay; + + double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax; + double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay; + + // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1) + double sigma_final_x = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms))); + double sigma_final_y = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms))); + + fit[0].dRo = sqrt(data->Cov[0][0]); + fit[0].duR = sqrt(data->Cov[1][1]); + fit[0].dDo = sqrt(data->Cov[2][2]); + fit[0].duD = sqrt(data->Cov[3][3]); + fit[0].dp = sqrt(data->Cov[4][4]); + + fit[0].dRo *= sigma_final_x; + fit[0].duR *= sigma_final_x; + fit[0].dDo *= sigma_final_y; + fit[0].duD *= sigma_final_y; + fit[0].dp *= sqrt(sigma_final_x * sigma_final_y); + } + + // (optionally) add up the chi square for the fit, only counting the unmasked points + double chisq = 0.0; + fit[0].Nfit = 0; + for (i = 0; i < Npoints; i++) { + if (points[i].mask) continue; + fit[0].Nfit ++; + + if (data->getChisq) { + double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR; + double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD; + double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); + double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); + chisq += SQ(points[i].X - Xf) * wx; + chisq += SQ(points[i].Y - Yf) * wy; + } + } + + // the reduced chisq is divided by Ndof = (2*Nfit - Nterms) + fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); + + return (TRUE); +} + +int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { + + int i; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2; + + PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0; + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + + for (i = 0; i < Npoints; i++) { + + // if (VERBOSE == 2) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], WX[i], Y[i], WY[i], T[i], pR[i], pD[i]); + + wx = points[i].Wx; + wy = points[i].Wy; + + Wx += wx; + Wy += wy; + + Tx += points[i].T*wx; + Ty += points[i].T*wy; + + Tx2 += SQ(points[i].T)*wx; + Ty2 += SQ(points[i].T)*wy; + + PR += points[i].pR*wx; + PD += points[i].pD*wy; + + PRT += points[i].pR*points[i].T*wx; + PDT += points[i].pD*points[i].T*wy; + + PRX += points[i].pR*points[i].X*wx; + PDY += points[i].pD*points[i].Y*wy; + + PR2 += SQ(points[i].pR)*wx; + PD2 += SQ(points[i].pD)*wy; + + Xs += points[i].X*wx; + Ys += points[i].Y*wy; + + XT += points[i].X*points[i].T*wx; + YT += points[i].Y*points[i].T*wy; + } + + data->A[0][0] = Wx; + data->A[0][1] = Tx; + data->A[0][4] = PR; + + data->A[1][0] = Tx; + data->A[1][1] = Tx2; + data->A[1][4] = PRT; + + data->A[2][2] = Wy; + data->A[2][3] = Ty; + data->A[2][4] = PD; + + data->A[3][2] = Ty; + data->A[3][3] = Ty2; + data->A[3][4] = PDT; + + data->A[4][0] = PR; + data->A[4][1] = PRT; + data->A[4][2] = PD; + data->A[4][3] = PDT; + data->A[4][4] = PR2 + PD2; + + data->B[0][0] = Xs; + data->B[1][0] = XT; + data->B[2][0] = Ys; + data->B[3][0] = YT; + data->B[4][0] = PRX + PDY; + + if (!dgaussjordan (data->A, data->B, 5, 1)) { +# if (DEBUG) + if (VERBOSE) fprintf (stderr, "error in fit\n"); + int j; + for (i = 0; i < 5; i++) { + for (j = 0; j < 5; j++) { + fprintf (stderr, "%e ", data->A[i][j]); + } + fprintf (stderr, " : %e\n", data->B[i][0]); + } +# endif + return FALSE; + } + + fit->Ro = data->B[0][0]; + fit->uR = data->B[1][0]; + fit->Do = data->B[2][0]; + fit->uD = data->B[3][0]; + fit->p = data->B[4][0]; + + return TRUE; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/ImageOps.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/ImageOps.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/ImageOps.c (revision 39266) @@ -769,12 +769,4 @@ ref[i].R += dTime * catalog[c].average[n].uR / 3600.0 / cos(ref[i].D*RAD_DEG); ref[i].D += dTime * catalog[c].average[n].uD / 3600.0; - } - - // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position - if (0) { - myAssert (!isnan(measure[0].RoffGAL), "oops"); - myAssert (!isnan(measure[0].DoffGAL), "oops"); - ref[i].R += measure[0].RoffGAL / 3600.0; - ref[i].D += measure[0].DoffGAL / 3600.0; } Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/UpdateObjects.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/UpdateObjects.c (revision 39266) @@ -140,4 +140,63 @@ // use a fit with fewer parameters. + // fit the parallax + proper-motion model + // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail + // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only + if (mode == FIT_PM_AND_PAR) { + if (Trange < PM_DT_MIN) { + mode = FIT_AVERAGE; + goto justPosition; + } + if (parRange < PAR_FACTOR_MIN) { + mode = FIT_PM_ONLY; + goto skipParallax; + } + if (fitStats->Npoints <= PAR_TOOFEW) { + mode = FIT_PM_ONLY; + goto skipParallax; + } + + // we are going to use the IRLS analysis to calculate the mean solution and the masking + // then run N_BOOTSTRAP_SAMPLES to measure the errors + fitStats->fitdataPar->getError = (N_BOOTSTRAP_SAMPLES < 3); + FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, VERBOSE); + + if (N_BOOTSTRAP_SAMPLES >= 3) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + if (Nnomask < 5) { + myAbort ("handle this case, please"); + } + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask); + fitStats->Nfit ++; + } + // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + } + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); + if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + + average[0].flags |= ID_STAR_FIT_PAR; + fitStats->Npar ++; + + // XXX a hard-wired hack... + // unless there is a clear problems (below) with the parallax fit, we will use it + if ((fabs(fitPar.uR) > 4.0) || (fabs(fitPar.uD) > 4.0)) { + mode = FIT_PM_ONLY; + } else { + average[0].flags |= ID_STAR_USE_PAR; + } + } + +skipParallax: + // *** first fit for the proper motion (skip fit if Trange or Npts is too small) *** if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { @@ -151,20 +210,29 @@ } - if (fitStats->NfitAlloc == 1) { - // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: - FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints); + if (average[0].flags & ID_STAR_USE_PAR) { + if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { + myAbort ("oops"); + } } else { - fitStats->Nfit = 0; - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); - if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue; - fitStats->Nfit ++; - } - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); - } - FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY); + fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3); + FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE); + + if (N_BOOTSTRAP_SAMPLES >= 3) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + if (Nnomask < 5) { + myAbort ("handle this case, please"); + } + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue; + fitStats->Nfit ++; + } + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + } + } // project Ro, Do back to RA,DEC @@ -172,59 +240,19 @@ if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - fitPM.p = fitPM.dp = 0.0; average[0].flags |= ID_STAR_FIT_PM; fitStats->Npm ++; // XXX a hard-wired hack... - if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) { + // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit + if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) { mode = FIT_AVERAGE; average[0].flags |= ID_STAR_BAD_PM; + } else { + if (!(average[0].flags & ID_STAR_USE_PAR)) { + average[0].flags |= ID_STAR_USE_PM; + } } } - // fit the parallax + proper-motion model - // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail - // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only - if (mode == FIT_PM_AND_PAR) { - if (parRange < PAR_FACTOR_MIN) { - mode = FIT_PM_ONLY; - goto justPosition; - } - if (fitStats->Npoints <= PAR_TOOFEW) { - mode = FIT_PM_ONLY; - goto justPosition; - } - - if (fitStats->NfitAlloc == 1) { - // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: - FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints); - } else { - fitStats->Nfit = 0; - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); - FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints); - fitStats->Nfit ++; - } - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); - } - FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR); - - // project Ro, Do back to RA,DEC - XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); - if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - - average[0].flags |= ID_STAR_FIT_PAR; - fitStats->Npar ++; - - // XXX a hard-wired hack... - if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) { - mode = FIT_PM_ONLY; - } - } - justPosition: { @@ -232,4 +260,6 @@ // if we only have one point, this is silly... + // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway? + if (fitStats->NfitAlloc == 1) { FitAstromResultSetPM (&fitPos, 1, average); @@ -618,64 +648,4 @@ } -int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) { - - int k; - - int Npoints = fitStats->Npoints; - FitAstromPoint *points = fitStats->points; - - // find Tmin & Tmax from the list of accepted measurements - double Tmin = points[0].T; - double Tmax = points[0].T; - double pRmin = +2.0; - double pRmax = -2.0; - double pDmin = +2.0; - double pDmax = -2.0; - - *Tmean = 0.0; - - double Tsum = 0.0; - double Wsum = 0.0; - for (k = 0; k < Npoints; k++) { - Tmin = MIN(Tmin, points[k].T); - Tmax = MAX(Tmax, points[k].T); - - float wx = 1.0 / SQ(points[k].dX); - - Tsum += points[k].T * wx; - Wsum += wx; - - // at this point, T is in years since J2000 - ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T); - pRmin = MIN (pRmin, points[k].pR); - pRmax = MAX (pRmax, points[k].pR); - pDmin = MIN (pDmin, points[k].pD); - pDmax = MAX (pDmax, points[k].pD); - } - *Trange = Tmax - Tmin; - - // mean epoch - *Tmean = Tsum / Wsum; - - // for HIGH_SPEED, just use the center of the range - if (RELASTRO_OP == OP_HIGH_SPEED) { - *Tmean = 0.5*(Tmax - Tmin); - } - - *parRange = hypot (pRmax - pRmin, pDmax - pDmin); - - /* we need to do the fit in a locally linear space; choose a ref coordinate */ - fitStats->coords.crval1 = points[0].R; - fitStats->coords.crval2 = points[0].D; - - // project all of the R,D coordinates to a plane centered on this coordinate. set - // the times to be relative to Tmean - for (k = 0; k < Npoints; k++) { - RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords); - points[k].T -= *Tmean; - } - return TRUE; -} - int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) { @@ -689,103 +659,4 @@ } return Nmax; -} - -int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) { - int i; - - // I need to draw Npoints random entries from 'points' with replacement: - for (i = 0; i < Npoints; i++) { - int N = Npoints * drand48(); - sample[i] = points[N]; - } - return TRUE; -} - -// calculate mean and sigma points for the 5 fit parameter -int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) { - - // generate a histogram for the selected element - double *values = NULL; - ALLOCATE (values, double, Nfit); - - int i; - - for (i = 0; i < Nfit; i++) { - switch (mode) { - case FIT_RESULT_RA: - values[i] = fit[i].Ro; - break; - case FIT_RESULT_DEC: - values[i] = fit[i].Do; - break; - case FIT_RESULT_uR: - values[i] = fit[i].uR; - break; - case FIT_RESULT_uD: - values[i] = fit[i].uD; - break; - case FIT_RESULT_PLX: - values[i] = fit[i].p; - break; - default: - myAbort ("invalid option"); - } - } - - dsort (values, Nfit); - - double median; - if (Nfit % 2) { - int Ncenter = Nfit / 2; - median = values[Ncenter]; - } else { - int Ncenter = Nfit / 2 - 1; - median = 0.5*(values[Ncenter] + values[Ncenter + 1]); - } - - double Slo = VectorFractionInterpolate (values, 0.158655, Nfit); - double Shi = VectorFractionInterpolate (values, 0.841345, Nfit); - double sigma = (Shi - Slo) / 2.0; - - switch (mode) { - case FIT_RESULT_RA: - result->Ro = median; - result->dRo = sigma; - break; - case FIT_RESULT_DEC: - result->Do = median; - result->dDo = sigma; - break; - case FIT_RESULT_uR: - result->uR = median; - result->duR = sigma; - break; - case FIT_RESULT_uD: - result->uD = median; - result->duD = sigma; - break; - case FIT_RESULT_PLX: - result->p = median; - result->dp = sigma; - break; - default: - myAbort ("invalid option"); - } - - return TRUE; -} - -double VectorFractionInterpolate (double *values, float fraction, int Npts) { - - float F = fraction * Npts; - int N = fraction * Npts; - - if (N < 0 ) return NAN; - if (N >= Npts - 2) return NAN; - - // interpolate between N,N+1 - - double S = (F - N) * (values[N+1] - values[N]) + values[N]; - return S; } Index: branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/fitpm.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/fitpm.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/relastro/src/fitpm.c (revision 39266) @@ -0,0 +1,209 @@ +# include "relastro.h" + +int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); + +# define MJD_MIN_MJD 55197 /* 2010/01/01,00:00:00 */ +# define MJD_MAX_MJD 57023 /* 2015/01/01,00:00:00 */ +# define MJD_MIN_YRS 10.00 /* 2010/01/01,00:00:00 */ +# define MJD_MAX_YRS 14.9993 /* 2015/01/01,00:00:00 */ +# define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */ + +# define POS_ERROR 0.010 /* arcsec */ +# define OUT_ERROR 0.500 /* arcsec */ + +# define N_STARS 3000 +# define N_POINTS 100 +# define N_OUTLIERS 100 +# define N_BOOTSTRAP 100 + +# define dcos(THETA) cos(RAD_DEG*THETA) +# define dsin(THETA) sin(RAD_DEG*THETA) + +enum { + FIT_PM_NONE, + FIT_PM_IRLS, + FIT_PM_NOCLIP, + FIT_PLX_IRLS, + FIT_PLX_BOOT, + FIT_PLX_NOCLIP, +}; + +int main (int argc, char **argv) { + + if (argc != 2) { + fprintf (stderr, "USAGE: %s (mode)\n", argv[0]); + exit (2); + } + + int mode = FIT_NONE; + if (!strcasecmp(argv[1], "pm")) { + mode = FIT_PM_NOCLIP; + } + if (!strcasecmp(argv[1], "pm-irls")) { + mode = FIT_PM_IRLS; + } + if (!strcasecmp(argv[1], "plx")) { + mode = FIT_PLX_NOCLIP; + } + if (!strcasecmp(argv[1], "plx-irls")) { + mode = FIT_PLX_IRLS; + } + if (!strcasecmp(argv[1], "plx-boot")) { + mode = FIT_PLX_BOOT; + } + + // plan_tests (14); + + // diag ("relastro fitpm tests"); + + // init the random seed + { + struct timeval now; + gettimeofday (&now, NULL); + long A = now.tv_sec + now.tv_usec * 1000000; + srand48(A); + } + + ohana_gaussdev_init (); + + FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); + fitStats->fitdataPM->getChisq = TRUE; + fitStats->fitdataPM->getError = TRUE; + fitStats->fitdataPar->getChisq = TRUE; + fitStats->fitdataPar->getError = TRUE; + // NOTE: surprisingly, the chisq and error calculations only add ~2-3% to the exec time + + int i, k; + + int Nstars = N_STARS; + for (i = 0; i < Nstars; i++) { + + double Ro = 360.0*(drand48() - 0.5); + double Phi = 2.0*(drand48() - 0.5); + double Do = DEG_RAD * asin(Phi); + + double uR = 0.5*(drand48() - 0.5); + double uD = 0.5*(drand48() - 0.5); + + double plx; + switch (mode) { + case FIT_PM_IRLS: + case FIT_PM_NOCLIP: + plx = 0.0*(drand48() + 0.0); + break; + case FIT_PLX_IRLS: + case FIT_PLX_BOOT: + case FIT_PLX_NOCLIP: + plx = 0.5*(drand48() + 0.0); + break; + default: + myAbort("oops"); + } + + + mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS); + + double Tmean, Trange, parRange; + FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange); + + FitAstromResult fitPM; + FitAstromResultInit (&fitPM); + + switch (mode) { + case FIT_PM_NOCLIP: + FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints); + break; + case FIT_PM_IRLS: + FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0); + break; + case FIT_PLX_NOCLIP: + FitPMandPar (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints); + break; + case FIT_PLX_IRLS: + FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0); + break; + + case FIT_PLX_BOOT: + FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0); + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask); + fitStats->Nfit ++; + } + // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + break; + + default: + myAbort("oops"); + } + + XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords); + + fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %f\n", + Ro, Do, uR, uD, plx, + fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean, fitPM.dRo, fitPM.dDo, fitPM.duR, fitPM.duD, fitPM.dp, fitPM.Nfit, fitPM.chisq); + + // ok (TRUE, "fitted star"); + } + // return exit_status(); + ohana_gaussdev_free(); + + exit (0); +} + +int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad) { + + int i; + + int Ntotal = Npoints + Nbad; + + FitAstromPoint *points = fitStats->points; + + for (i = 0; i < Npoints; i++) { + FitAstromPointInit (&points[i]); + + // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) + points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + + double pR, pD; + ParFactor (&pR, &pD, Ro, Do, points[i].T); + + double dR = ohana_gaussdev_rnd(0.0, POS_ERROR); + double dD = ohana_gaussdev_rnd(0.0, POS_ERROR); + + points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do); + points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600; + + points[i].dX = POS_ERROR; + points[i].dY = POS_ERROR; + } + + for (i = Npoints; i < Ntotal; i++) { + FitAstromPointInit (&points[i]); + + // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) + points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + + double pR, pD; + ParFactor (&pR, &pD, Ro, Do, points[i].T); + + double dR = OUT_ERROR*(drand48() - 0.5); + double dD = OUT_ERROR*(drand48() - 0.5); + + points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do); + points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600; + + points[i].dX = POS_ERROR; + points[i].dY = POS_ERROR; + } + fitStats->Npoints = Ntotal; + + return TRUE; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/BrightCatalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/BrightCatalog.c (revision 39266) @@ -67,4 +67,5 @@ GET_COLUMN(M, "MAG_SYS", float); GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(Mflat, "MAG_FLAT", float); GET_COLUMN(dM, "MAG_ERR", float); GET_COLUMN(airmass, "AIRMASS", float); @@ -89,4 +90,5 @@ measure[i].M = M[i]; measure[i].Mcal = Mcal[i]; + measure[i].Mflat = Mflat[i]; measure[i].dM = dM[i]; measure[i].airmass = airmass[i]; @@ -109,4 +111,5 @@ free (M ); free (Mcal ); + free (Mflat ); free (dM ); free (airmass ); @@ -269,4 +272,5 @@ gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_FLAT", "magnitude (flat)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_ERR", "magnitude (err)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "AIRMASS", "airmass", NULL, 1.0, 0.0); @@ -290,4 +294,5 @@ float *M ; ALLOCATE (M , float, catalog->Nmeasure); float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *Mflat ; ALLOCATE (Mflat , float, catalog->Nmeasure); float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); @@ -310,4 +315,5 @@ M[i] = measure[i].M ; Mcal[i] = measure[i].Mcal ; + Mflat[i] = measure[i].Mflat ; dM[i] = measure[i].dM ; airmass[i] = measure[i].airmass ; @@ -329,4 +335,5 @@ gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_FLAT", Mflat, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); @@ -346,4 +353,5 @@ free (M ); free (Mcal ); + free (Mflat ); free (dM ); free (airmass ); Index: branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/setMrelCatalog.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/relphot/src/setMrelCatalog.c (revision 39266) @@ -1,3 +1,11 @@ # include "relphot.h" + +/* +# define TEST_OBJ_ID 0x3ae2 +# define TEST_CAT_ID 0x16c5f +*/ +# define TEST_OBJ_ID 0 +# define TEST_CAT_ID 0 + int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags); @@ -113,5 +121,5 @@ // option for a test print - if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) { + if ((averageT[0].objID == TEST_OBJ_ID) && (averageT[0].catID == TEST_CAT_ID)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); @@ -353,5 +361,5 @@ if (isSetMrelFinal) { - if (FALSE && (average->objID == 0x1) && (average->catID)) { + if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) { fprintf (stderr, "test obj\n"); } @@ -723,5 +731,5 @@ // option for a test print - if (FALSE && (average[0].objID == 0x0000049c) && (average[0].catID == 0x00016c5e)) { + if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/Makefile =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/Makefile (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/Makefile (revision 39266) @@ -1,3 +1,5 @@ -default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client +# default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client +default: uniphot setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client + help: @echo "make options: uniphot setphot setfwhm default help install default (uniphot setphot setfwhm)" @@ -34,5 +36,6 @@ ckids_client: $(BIN)/ckids_client.$(ARCH) -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 +# 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 +install: $(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 UNIPHOT = \ Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/include/setphot.h =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/include/setphot.h (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/include/setphot.h (revision 39266) @@ -17,4 +17,25 @@ unsigned int flags; } ImageSubset; + +// we have one correction (an image) for each filter and chip +typedef struct { + int Nx; // number of chips in x + int Ny; // number of chips in y + int Nfilter; // number of filters + int Nseason; // number of correction peridos + + int Nchips; // chip offset (Nx*Ny) + int Nflats; // season offset (Nx*Ny*Nfilters) + int Nvalues; // Nx*Ny*Nfilters*Nseason + + int dX; // superpixel size + int dY; // superpixel size + + int NxCCD; // number of pixels + int NyCCD; // number of pixels + + Matrix **matrix; // allocate an array of pointers + // index = ix + iy*Nx + filter*Nchips + dir*Ngroup +} CamCorrection; /* global variables set in parameter file */ Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_correction.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_correction.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_correction.c (revision 39266) @@ -29,6 +29,6 @@ int Ndir, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD; - if (!gfits_scan (&header, "NDIR", "%d", 1, &Ndir)) return FALSE; - if (!gfits_scan (&header, "NFILTER", "%d", 1, &Nfilter)) return FALSE; + if (!gfits_scan (&header, "NDIR", "%d", 1, &Ndir)) return FALSE; + if (!gfits_scan (&header, "NFILTER", "%d", 1, &Nfilter)) return FALSE; if (!gfits_scan (&header, "NX", "%d", 1, &Nx)) return FALSE; if (!gfits_scan (&header, "NY", "%d", 1, &Ny)) return FALSE; @@ -61,5 +61,4 @@ } - int found = FALSE; while (TRUE) { Header theader; @@ -90,8 +89,4 @@ // assert that cam->matrix[index] is NULL? cam->matrix[index] = matrix; - - if (!found) { - - } } return TRUE; Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_zpt_correction.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_zpt_correction.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/cam_zpt_correction.c (revision 39266) @@ -0,0 +1,272 @@ +# include "setphot.h" + +// camera systematic (photoflat) correction functions: +// load the correction set from a file + +// note that we are allocating pointers for XY00 - XY77, but only the octal elements are set (and not the 00,07,70,77 ones) + +static CamCorrection *cam = NULL; + +CamCorrection *get_cam_correction_ptr () { + return cam; +} + +void free_cam_correction (CamCorrection *cam) { + + if (!cam) return; + + for (i = 0; i < cam->Nvalues; i++) { + FREE (cam->matrix[i]); + } + + free (cam->matrix); + free (cam); +} + +CamCorrection *alloc_cam_correction (int Nx, int Ny, int Nfilter, int Nseason) { + + CamCorrection *cam; + + ALLOCATE (cam, CamCorrection, 1); + + ALLOCATE (cam->matrix, Matrix *, Nseason*Nfilter*Nx*Ny); + + cam->Nx = Nx; // for gpc1, should be 8 + cam->Ny = Ny; // for gpc1, should be 8 + cam->Nfilter = Nfilter; // for gpc1, should be 5 + cam->Nseason = Nseason; // for gpc1, should be 5 (1 for loading highres correction) + + cam->Nchips = Nx * Ny; + cam->Nflats = Nx * Ny * Nfilter; + cam->Nvalues = Nx * Ny * Nfilter * Nseason; + + ALLOCATE (cam->matrix, Matrix *, cam->Nvalues); + for (i = 0; i < cam->Nvalues; i++) { + cam->matrix[i] = NULL; + } + + return cam; +} + +int load_cam_correction (char *filename) { + + int i, ix, iy, dir, filter; + Header header; + + /* open file for input */ + FILE *f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + gprint (GP_ERR, "can't open file for read : %s\n", filename); + return FALSE; + } + + // read the PHU header + if (!gfits_load_header (f, &header)) return FALSE; + + int Nbytes = gfits_data_size (&header); + fseeko (f, Nbytes, SEEK_CUR); + + int Ndir, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD; + if (!gfits_scan (&header, "NSEASON", "%d", 1, &Nseason)) return FALSE; + if (!gfits_scan (&header, "NFILTER", "%d", 1, &Nfilter)) return FALSE; + if (!gfits_scan (&header, "NX", "%d", 1, &Nx)) return FALSE; + if (!gfits_scan (&header, "NY", "%d", 1, &Ny)) return FALSE; + if (!gfits_scan (&header, "DX", "%d", 1, &dX)) return FALSE; + if (!gfits_scan (&header, "DY", "%d", 1, &dY)) return FALSE; + + if (!gfits_scan (&header, "NX_CHIP", "%d", 1, &NxCCD)) return FALSE; + if (!gfits_scan (&header, "NY_CHIP", "%d", 1, &NyCCD)) return FALSE; + + CamCorrection *cam = alloc_cam_correction (Nx, Ny, Nfilter, Nseason); + + cam->dX = dX; // superpixel sampling, x-dir + cam->dY = dY; // superpixel sampling, y-dir + + cam->NxCCD = NxCCD; // pixels per chip + cam->NyCCD = NyCCD; // pixels per chip + + while (TRUE) { + Header theader; + + // load data for this header : if not found, assume we hit the end of the file + if (!gfits_load_header (f, &theader)) break; + + if (!gfits_scan (&theader, "FILTER", "%d", 1, &filter)) return FALSE; + if (!gfits_scan (&theader, "SEASON", "%d", 1, &season)) return FALSE; + + // XXX NOTE: astroflat.20150209.fits had ix and iy backwards in header + // double-check that the new flats are OK + + if (!gfits_scan (&theader, "X_CHIP", "%d", 1, &ix)) return FALSE; + if (!gfits_scan (&theader, "Y_CHIP", "%d", 1, &iy)) return FALSE; + + Matrix *matrix = NULL; + ALLOCATE (matrix, Matrix, 1); + if (!gfits_load_matrix (f, matrix, &theader)) break; + + int index = ix + iy*cam->Nx + filter*cam->Nchips + season*cam->Nflats; + myAssert (index >= 0, "index too small"); + myAssert (index < cam->Nvalues, "index too big"); + + myAssert (!cam->matrix[index], "entry already assigned?"); + + // assert that cam->matrix[index] is NULL? + cam->matrix[index] = matrix; + } + return TRUE; +} + +// input is: +// chipID == octal chip ID (eg 40 for XY40) +// filter == (01234 = grizy) +int get_cam_correction (int chipID, int filter, float Xccd, float Yccd, double *dX, double *dY) { + + *dX = 0.0; + *dY = 0.0; + + int index, jx, jy, NxModel; + Matrix *matrix; + float *buffer, value; + + // split out the chipID number (eg 43 for XY43) into X and Y elements + int ix = (int) (chipID / 10); + int iy = (int) (chipID % 10); + + // dX (dir == 0) + index = ix + iy*cam->Nx + filter*cam->Nchips; + + matrix = cam->matrix[index]; + NxModel = cam->matrix[index]->Naxis[0]; + + buffer = (float *) matrix->buffer; + + jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0); + jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0); + + value = buffer[jx + jy*NxModel]; + + *dX = isnan(value) ? 0.0 : value; + + // dY (dir == 1) uses the next block + index += cam->Ngroup; + + matrix = cam->matrix[index]; + NxModel = cam->matrix[index]->Naxis[0]; + + buffer = (float *) matrix->buffer; + + jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0); + jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0); + + value = buffer[jx + jy*NxModel]; + + *dY = isnan(value) ? 0.0 : value; + + return TRUE; +} + +float get_cam_correction_by_flatid (int flat_id, float Xccd, float Yccd) { + + float dM = 0.0; + + // validate the flat_id (not out of range?) + int seq = camcorr->IDtoSeq[flat_id]; + + Matrix *matrix = cam->matrix[seq]; + NxModel = cam->matrix[index]->Naxis[0]; + + buffer = (float *) matrix->buffer; + + jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0); + jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0); + + value = buffer[jx + jy*NxModel]; + + dM = isnan(value) ? 0.0 : value; + + return dM; +} + +int merge_flatcorr_and_camcorr (FlatCorrectionTable *flatcorr, CamCorrection *camcorr) { + + // this function is specifically tuned for merging the ubercal flats with the high-res flats + + // we have an ubercal FlatCorrectionTable with Nseasons, Nfilter, Nx, Ny entries + // we have a highres flat with Nfilter, Nx, Ny images + + myAssert (flatcorr->Nfilter <= camcorr->Nfilter, "unclear how to use more flatcorr filters"); + myAssert (camcorr->Nseason == 1, "no code to merge seasons flats to seasonal flats"); + + // we need to create a new CamCorrection with Nx,Ny,Nfilter from camcorr and Nseason from flatcorr: + CamCorrection *newcorr = alloc_cam_correction (camcorr->Nx, camcorr->Ny, camcorr->Nfilter, flatcorr->Nseason); + newcorr->dX = camcorr->dX; // superpixel sampling, x-dir + newcorr->dY = camcorr->dY; // superpixel sampling, y-dir + newcorr->NxCCD = camcorr->NxCCD; // pixels per chip + newcorr->NyCCD = camcorr->NyCCD; // pixels per chip + + for (Ns = 0; Ns < flatcorr->Nseason; Ns++) { + newcorr->tstart[Ns] = flatcorr->tstart[Ns]; + newcorr->tstop[Ns] = flatcorr->tstop[Ns]; + } + + // create the new flat images (1 per season) + for (Ns = 0; Ns < newcorr->Nseason; Ns++) { + for (Nf = 0; Nf < newcorr->Nfilter; Nf++) { + for (ix = 0; ix < newcorr->Nx; ix++) { + for (iy = 0; iy < newcorr->Ny; iy++) { + int index = ix + iy*newcorr->Nx + Nf*newcorr->Nchips + Ns*newcorr->Nflats; + ALLOCATE (newcorr->matrix[index], Matrix, 1); + + // camcorr only had one season + int idxOld = ix + iy*camcorr->Nx + Nf*camcorr->Nchips; + + // allocates the buffer + gfits_copy_matrix (camcorr->matrix[idxOld], newcorr->matrix[index]); + } + } + } + } + + // FlatCorrectionTable has: + // table->image[seq] with metadata about the specific image + // table->offset[seq][ix][iy] with dM values for image[seq] + + // loop over the flatcorr images and match from image[seq] to newcorr: + + for (i = 0; i < flatcorr->Nimage; i++) { + + int ID = flatcorr->image[i].ID; + myAssert (flatcorrTable->IDtoSeq[ID] == i, "oops"); + + int season = -1; + for (j = 0; j < newcorr->Nseason; j++) { + if (flatcorr->image[i].tstart != newcorr->tstart[j]) continue; + season = j; + } + myAssert (season >= 0, "oops"); + + int photcode = flatcorr->image[i].photcode; + + iy = flatcorr->image[i].photcode % 10; + ix = (int)(flatcorr->image[i].photcode / 10) % 10; + filter = (int)(flatcorr->image[i].photcode / 100) % 10; + + int seq = ix + iy * newcorr->Nx + filter * newcorr->Nchips + season * camcorr->Nflat; + + // we now have newcorr->matrix[seq] + newcorr->matrix[seq]; + + // now loop over the newcorr pixels to get the flatcorr value: + for (jx = 0; jx < newcorr->NxCCD; jx++) { + for (jy = 0; jy < newcorr->NyCCD; jy++) { + + // convert jx,jy in newcorr pixels to Jx,Jy in flatcorr pixels + newcorr->matrix[seq]->buffer[jx + jy*NxCCD] += flatcorr->offset[i][Jx][Jy]; + + } + } + } + + return TRUE; +} + Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot.c (revision 39266) @@ -33,4 +33,17 @@ int N; + + CAM_FILE = NULL; + if ((N = get_argument (argc, argv, "-CAM"))) { + remove_argument (N, &argc, argv); + char *tmpfile = strcreate (argv[N]); + CAM_FILE = abspath (tmpfile, DVO_MAX_PATH); + remove_argument (N, &argc, argv); + } + CAM_RESET = FALSE; + if ((N = get_argument (argc, argv, "-CAM-reset"))) { + remove_argument (N, &argc, argv); + CAM_RESET = TRUE; + } VERBOSE = FALSE; Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot_client.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot_client.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/initialize_setphot_client.c (revision 39266) @@ -26,4 +26,17 @@ int N; + + CAM_FILE = NULL; + if ((N = get_argument (argc, argv, "-CAM"))) { + remove_argument (N, &argc, argv); + char *tmpfile = strcreate (argv[N]); + CAM_FILE = abspath (tmpfile, DVO_MAX_PATH); + remove_argument (N, &argc, argv); + } + CAM_RESET = FALSE; + if ((N = get_argument (argc, argv, "-CAM-reset"))) { + remove_argument (N, &argc, argv); + CAM_RESET = TRUE; + } VERBOSE = FALSE; Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/match_camcorr_to_images.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/match_camcorr_to_images.c (revision 39266) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/match_camcorr_to_images.c (revision 39266) @@ -0,0 +1,53 @@ +# include "setphot.h" + +// the date/time of the image is used to find the 'season' +// the photcode is used to find the actual flat correction + +int match_camcorr_to_images (Image *image, off_t Nimage) { + + int i, j; + + CamCorrection *camcorr = get_cam_correction_ptr (); + + // we have an array of CamCorrection->matrix values, where each matrix is a flat-field image + // the sequence is seq = ix + Nx_chip*iy + filter*Nchips + season*(Nchips*Nfilter) + + // we derive ix,iy,filter from photcode : 10143 -> ix = 4, iy = 3, filter = 1 + // season comes from the date/time + + int minCode = 10000; + int maxCode = 10576; + + for (i = 0; i < Nimage; i++) { + if (!image[i].photcode) continue; // skip PHU images + + if (image[i].photcode < minCode) continue; // skip non-gpc1-chip images + if (image[i].photcode > maxCode) continue; // skip non-gpc1-chip images + + int found = FALSE; + for (j = 0; !found && (j < camcorr->Nseason); j++) { + if (image[i].tzero < camcorr->tstart[j]) continue; + if (image[i].tzero > camcorr->tstop[j]) continue; + + int iy = image[i].photcode % 10; + int ix = (int)(image[i].photcode / 10) % 10; + int filter = (int)(image[i].photcode / 100) % 10; + + int seq = ix + iy * camcorr->Nx + filter * camcorr->Nchips + j * camcorr->Ngroup; + + image[i].photom_map_id = camcorr->imageID[seq].ID; + found = TRUE; + } + if (!found) { + Nmissed ++; + fprintf (stderr, "image does not match flatcorr Table ranges: %s\n", image[i].name); + } + + if (Nmissed > 100) { + fprintf (stderr, "something is wrong: too many images do not match\n"); + exit (2); + } + + } + return TRUE; +} Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/setphot.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/setphot.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/setphot.c (revision 39266) @@ -8,5 +8,6 @@ ZptTable *zpts; Image *image; - FlatCorrectionTable flatcorrTable; + + // FlatCorrectionTable flatcorrTable; /* get configuration info, args, lockfile */ @@ -33,4 +34,22 @@ } + // if supplied, read in the camera-static flat-field correction + if (CAM_FILE) { + if (!load_cam_correction (CAM_FILE)) { + fprintf (stderr, "failed to load camera-static flat-field correction\n"); + exit (2); + } + merge_flatcorr_with_cam (flatcorrTable); + } else { + convert_flatcorr_to_cam (flatcorrTable); + } + + { + char newflatfile[DVO_MAX_PATH]; + int size = snprintf (newflatfile, DVO_MAX_PATH, "%s/flatfield.fits", CATDIR); + assert (size < DVO_MAX_PATH); + save_cam_correction (newflatfile); + } + if (!zpts) Shutdown ("failed to load zero points, or empty table"); @@ -42,9 +61,11 @@ if (UBERCAL) { - match_flatcorr_to_images (image, Nimage, &flatcorrTable); + // we are going to deprecate the flatcorr imaage lookup + // match_flatcorr_to_images (image, Nimage, &flatcorrTable); + match_camcorr_to_images (image, Nimage); } if (!IMAGES_ONLY) { - status = update_dvo_setphot (image, Nimage, &flatcorrTable); + status = update_dvo_setphot (image, Nimage); } Index: branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/update_catalog_setphot.c =================================================================== --- branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/update_catalog_setphot.c (revision 39257) +++ branches/eam_branches/ipp-20151113/Ohana/src/uniphot/src/update_catalog_setphot.c (revision 39266) @@ -1,5 +1,5 @@ # include "setphot.h" -void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) { +void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage) { off_t i, j, found; @@ -7,4 +7,6 @@ // if we are resetting, reset all flags DVOMeasureFlags photomFlags = ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_PHOTOM_UBERCAL; + + CamCorrection *camcorr = get_cam_correction_ptr (); found = 0; @@ -21,13 +23,14 @@ float Mcal = image[id].Mcal; float dMcal = image[id].dMcal; - float Mcal_offset = 0.0; + float Mflat = 0.0; // if we know about a flat-field correction, then we need to apply the sub-chip correction int flat_id = image[id].photom_map_id; if (flat_id > 0) { - Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd); + Mflat = CamCorrectionOffset (camcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd); } - catalog[0].measure[m].Mcal = Mcal - Mcal_offset; + catalog[0].measure[m].Mcal = Mcal; + catalog[0].measure[m].Mflat = Mflat; catalog[0].measure[m].dMcal = dMcal; myAssert(isfinite(catalog[0].measure[m].Mcal), "oops: ubercal made a nan");