Index: trunk/Ohana/src/addstar/src/loadgalphot_readstars.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 39121) +++ trunk/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 39122) @@ -237,6 +237,7 @@ // I need to find the minimum position (interpolated) in this 2D space + // I want to use this flag to set a bit in secfilt. is that possible? or maybe galphot.dummy becomes flags if (!FitChisqMinimum (fit, &stars[i].galphot, &GAL_CHISQ[Nsample*i], &GAL_FLUX[Nsample*i], &GAL_FLUX_ERR[Nsample*i], Nsample, FR_MAJOR_MIN[i], FR_MAJOR_MAX[i], FR_MAJOR_DEL[i], FR_MINOR_MIN[i], FR_MINOR_MAX[i], FR_MINOR_DEL[i])) { - Nbad ++; + stars[i].flag = TRUE; } @@ -298,8 +299,6 @@ // here are the steps: // find the minimum chisq grid point (set X,Y values) - // choose the nearest 9 points - // fit 2D parabola (6 free parameters) - // use this to outlier-clip remaining points - // fit 2D parabola to the rest + // rescale chisq values so minimum is == Npix + // fit 2D parabola memset (mask, 0, Npts*sizeof(char)); @@ -324,27 +323,26 @@ } if (!Nvalid) { - fprintf (stderr, "no good points\n"); return FALSE; // we cannot do anything if there is no finite chisq } - if (Nvalid < 6) { - galphot->majorAxis = Xpt[iMin]; - galphot->minorAxis = Ypt[iMin]; - galphot->majorAxisErr = MajorDel; - galphot->minorAxisErr = MinorDel; - fprintf (stderr, "too few good points: %d\n", Nvalid); - return FALSE; - } - - // re-fit all except the most extreme set - if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) { - // raise a flag? - galphot->majorAxis = Xpt[iMin]; - galphot->minorAxis = Ypt[iMin]; - galphot->majorAxisErr = MajorDel; - galphot->minorAxisErr = MinorDel; - fprintf (stderr, "failed fit\n"); - return FALSE; - } + for (i = 0; i < Npts; i++) { + chisq[i] *= (galphot->Npix / chisqMin); + } + + if (Nvalid < 6) goto bad_fit; + +# if (0) + { + FILE *f = fopen ("galphot.dump.dat", "w"); + for (i = 0; i < Npts; i++) { + fprintf (f, "%d %d %f %f %f %f\n", i, mask[i], Xpt[i], Ypt[i], chisq[i], flux[i]); + } + fclose (f); + } +# endif + + // fit with 3 iterations, 5 sigma clipping (set above after fit2d_init) + // on failure (for any reason), the stars[i].flag is set to TRUE (XXX bad choice) + if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) goto bad_fit; // get X,Y for the chisq min from the fit: @@ -353,38 +351,97 @@ float Ymin = (fit->c10*fit->c11 - 2.0*fit->c20*fit->c01) * Det; + // allow the fitting minimum to be just outside the grid, but no further + int inRange = TRUE; + inRange = inRange && (Xmin > MajorMin - MajorDel); + inRange = inRange && (Xmin < MajorMax + MajorDel); + inRange = inRange && (Ymin > MinorMin - MinorDel); + inRange = inRange && (Ymin < MinorMax + MinorDel); + + if (!inRange) goto bad_fit; + // chisqMin @ (Xmin,Ymin) - chisqMin = - fit->c00 + - fit->c10*Xmin + fit->c20*Xmin*Xmin + - fit->c01*Ymin + fit->c02*Ymin*Ymin + - fit->c11*Xmin*Ymin; - - float chisqOff = chisqMin + 1.0; - - float A, B, C, Xoff, Yoff; - - // Xoff @ chisqMin + 1.0 + float A, B, C, Q; + + // Xoff @ chisqMin + 1.0, delta Y = 0.0 A = fit->c20; - B = fit->c10 + fit->c11*Ymin; - C = fit->c00 + fit->c01*Ymin + fit->c02*Ymin*Ymin - chisqOff; - Xoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A); - - // Xoff @ chisqMin + 1.0 + B = fit->c10 + 2*fit->c20*Xmin + fit->c11*Ymin; + C = -1; + Q = B*B - 4*A*C; + if (Q < 0.0) goto bad_err; + + float dXmin = (-B + sqrt(Q)) / (2.0 * A); + + // Yoff @ chisqMin + 1.0, delta X = 0.0 A = fit->c02; - B = fit->c01 + fit->c11*Xmin; - C = fit->c00 + fit->c10*Xmin + fit->c20*Xmin*Xmin - chisqOff; - Yoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A); + B = fit->c01 + 2*fit->c02*Ymin + fit->c11*Xmin; + C = -1; + Q = B*B - 4*A*C; + if (Q < 0.0) goto bad_err; + + float dYmin = (-B + sqrt(Q)) / (2.0 * A); galphot->majorAxis = Xmin; galphot->minorAxis = Ymin; - galphot->majorAxisErr = fabs(Xoff - Xmin); - galphot->minorAxisErr = fabs(Yoff - Ymin); + galphot->majorAxisErr = dXmin; + galphot->minorAxisErr = dYmin; + + // convert (Xmin, Ymin) to (iXmin, iYmin) + float iXminF = (Xmin - MajorMin) / MajorDel; + float iYminF = (Ymin - MinorMin) / MinorDel; + + // interpolate in X and in Y + int iXmin = floor(iXminF); + iXmin = MAX(0,MIN(Nx - 2, iXmin)); // force range of iXmin to be 0,Nx-1 + float fXmin = iXminF - iXmin; + + int iYmin = floor(iYminF); + iYmin = MAX(0,MIN(Ny - 2, iYmin)); // force range of iYmin to be 0,Ny-1 + float fYmin = iYminF - iYmin; + + float V00 = flux[iXmin+0 + (iYmin+0)*Nx]; + float V01 = flux[iXmin+0 + (iYmin+1)*Nx]; + float V10 = flux[iXmin+1 + (iYmin+0)*Nx]; + float V11 = flux[iXmin+1 + (iYmin+1)*Nx]; + + float Vx0 = V10*fXmin + V00*(1.0 - fXmin); + float Vx1 = V11*fXmin + V01*(1.0 - fXmin); + + float value = fYmin * Vx1 + (1.0 - fYmin) * Vx0; + + float dV00 = fluxErr[iXmin+0 + (iYmin+0)*Nx]; + float dV01 = fluxErr[iXmin+0 + (iYmin+1)*Nx]; + float dV10 = fluxErr[iXmin+1 + (iYmin+0)*Nx]; + float dV11 = fluxErr[iXmin+1 + (iYmin+1)*Nx]; + + float dVx0 = dV10*fXmin + dV00*(1.0 - fXmin); + float dVx1 = dV11*fXmin + dV01*(1.0 - fXmin); + + float dValue = fYmin * dVx1 + (1.0 - fYmin) * dVx0; // use bilinear interpolation to choose Flux @ (Xmin,Ymin)? + // float magRaw = -2.5*log10(flux[iMin]) + ZeroPt; // correct for exptime? + // float magErrRaw = sqrt(fluxErr[iMin]) / flux[iMin]; + + galphot->mag = -2.5*log10(value) + ZeroPt; // correct for exptime? + galphot->magErr = sqrt(dValue) / value; + galphot->chisq = chisqMin; + + return TRUE; + +bad_fit: + galphot->majorAxis = Xpt[iMin]; + galphot->minorAxis = Ypt[iMin]; + galphot->majorAxisErr = MajorDel; + galphot->minorAxisErr = MinorDel; + galphot->mag = -2.5*log10(flux[iMin]) + ZeroPt; // correct for exptime? galphot->magErr = sqrt(fluxErr[iMin]) / flux[iMin]; galphot->chisq = chisqMin; - - return TRUE; + return FALSE; + +bad_err: + galphot->majorAxisErr = MajorDel; + galphot->minorAxisErr = MinorDel; + return FALSE; }