Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 9513) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 9526) @@ -4,13 +4,13 @@ 1 / (1 + z^M + z^N) - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2.0) / SigmaX; - params->data.F32[5] = sqrt(2.0) / SigmaY; - params->data.F32[6] = Sxy; - params->data.F32[7] = - params->data.F32[8] = + params->data.F32[PM_PAR_SKY] = So; + params->data.F32[PM_PAR_FLUX] = Zo; + params->data.F32[PM_PAR_XPOS] = Xo; + params->data.F32[PM_PAR_YPOS] = Yo; + params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX; + params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY; + params->data.F32[PM_PAR_SXY] = Sxy; + params->data.F32[PM_PAR_7] = + params->data.F32[PM_PAR_8] = *****************************************************************************/ @@ -21,33 +21,33 @@ psF32 *PAR = params->data.F32; - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + psF32 X = x->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = x->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = PAR[PM_PAR_SXX]*X; + psF32 py = PAR[PM_PAR_SYY]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; if (z < 0) z = 0; psF32 zp = pow(z,1.25); - psF32 r = 1.0 / (1 + PAR[7]*z + z*zp); - - psF32 r1 = PAR[1]*r; - psF32 f = r1 + PAR[0]; + psF32 r = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp); + + psF32 r1 = PAR[PM_PAR_FLUX]*r; + psF32 f = r1 + PAR[PM_PAR_SKY]; if (deriv != NULL) { psF32 *dPAR = deriv->data.F32; - // note difference from a pure gaussian: q = params->data.F32[1]*r + // note difference from a pure gaussian: q = params->data.F32[PM_PAR_FLUX]*r psF32 t = r1*r; - psF32 q = t*(PAR[7] + 2.25*zp); - - dPAR[0] = +1.0; - dPAR[1] = +r; - dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - dPAR[4] = -2.0*q*px*X; - dPAR[5] = -2.0*q*py*Y; - dPAR[6] = -q*X*Y; - dPAR[7] = -t*z; + psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp); + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_FLUX] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + dPAR[PM_PAR_YPOS] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); + dPAR[PM_PAR_SXX] = -2.0*q*px*X; + dPAR[PM_PAR_SYY] = -2.0*q*py*Y; + dPAR[PM_PAR_SXY] = -q*X*Y; + dPAR[PM_PAR_7] = -t*z; } return(f); @@ -64,30 +64,30 @@ (*params_max)->n = (*params_max)->nalloc; - beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 3e6; - beta_lim[0][0].data.F32[2] = 5; - beta_lim[0][0].data.F32[3] = 5; - beta_lim[0][0].data.F32[4] = 0.5; - beta_lim[0][0].data.F32[5] = 0.5; - beta_lim[0][0].data.F32[6] = 0.5; - beta_lim[0][0].data.F32[7] = 0.5; - - params_min[0][0].data.F32[0] = -1000; - params_min[0][0].data.F32[1] = 0; - params_min[0][0].data.F32[2] = -100; - params_min[0][0].data.F32[3] = -100; - params_min[0][0].data.F32[4] = 0.01; - params_min[0][0].data.F32[5] = 0.01; - params_min[0][0].data.F32[6] = -5.0; - params_min[0][0].data.F32[7] = 0.1; - - params_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e8; - params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[4] = 2.0; - params_max[0][0].data.F32[5] = 2.0; - params_max[0][0].data.F32[6] = +5.0; - params_max[0][0].data.F32[7] = 10.0; + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_FLUX] = 3e6; + beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_7] = 0.5; + + params_min[0][0].data.F32[PM_PAR_SKY] = -1000; + params_min[0][0].data.F32[PM_PAR_FLUX] = 0; + params_min[0][0].data.F32[PM_PAR_XPOS] = -100; + params_min[0][0].data.F32[PM_PAR_YPOS] = -100; + params_min[0][0].data.F32[PM_PAR_SXX] = 0.01; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.01; + params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; + params_min[0][0].data.F32[PM_PAR_7] = 0.1; + + params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; + params_max[0][0].data.F32[PM_PAR_FLUX] = 1e8; + params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_SXX] = 2.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 2.0; + params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; + params_max[0][0].data.F32[PM_PAR_7] = 10.0; return (TRUE); @@ -102,12 +102,12 @@ psF32 *params = model->params->data.F32; - params[0] = moments->Sky; - params[1] = moments->Peak - moments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx); - params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy); - params[6] = 0.0; - params[7] = 1.0; + params[PM_PAR_SKY] = moments->Sky; + params[PM_PAR_FLUX] = moments->Peak - moments->Sky; + params[PM_PAR_XPOS] = peak->x; + params[PM_PAR_YPOS] = peak->y; + params[PM_PAR_SXX] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx); + params[PM_PAR_SYY] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy); + params[PM_PAR_SXY] = 0.0; + params[PM_PAR_7] = 1.0; return(true); @@ -121,7 +121,7 @@ psF32 *PAR = params->data.F32; - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); + psF64 A1 = PS_SQR(PAR[PM_PAR_SXX]); + psF64 A2 = PS_SQR(PAR[PM_PAR_SYY]); + psF64 A3 = PS_SQR(PAR[PM_PAR_SXY]); psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); // Area is equivalent to 2 pi sigma^2 @@ -135,7 +135,7 @@ float f1, f2; for (z = DZ; z < 50; z += DZ) { - f1 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); z += DZ; - f2 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); norm += f0 + 4*f1 + f2; f0 = f2; @@ -143,5 +143,5 @@ norm *= DZ / 3.0; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -158,12 +158,12 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); - if (flux >= PAR[1]) + if (flux >= PAR[PM_PAR_FLUX]) return (1.0); // if Sx == Sy, sigma = Sx == Sy - psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); - psF64 limit = flux / PAR[1]; + psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); + psF64 limit = flux / PAR[PM_PAR_FLUX]; # if 0 @@ -174,6 +174,6 @@ for (z = 0.0; z < 30.0; z += dz) { Nstep ++; - f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); - // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z)); + f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + // test: f = 1.0 / (1 + PAR[PM_PAR_7]*z + PS_SQR(z)); if (f < limit) break; @@ -189,13 +189,13 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / 2.25)); - float z1 = (1.0 / limit) / PAR[7]; + float z1 = (1.0 / limit) / PAR[PM_PAR_7]; z1 = PS_MAX (z0, z1); z0 = 0.0; - float f0 = 1.0 / (1 + PAR[7]*z0 + pow(z0, 2.25)); - float f1 = 1.0 / (1 + PAR[7]*z1 + pow(z1, 2.25)); + float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25)); + float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25)); while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) { z = 0.5*(z0 + z1); - f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); // fprintf (stderr, "%f %f %f : %f %f %f\n", f0, f, f1, z0, z, z1); if (f > limit) { @@ -224,14 +224,15 @@ psF32 *in = modelFLT->params->data.F32; - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - + out[PM_PAR_SKY] = in[PM_PAR_SKY]; + out[PM_PAR_FLUX] = in[PM_PAR_FLUX]; + out[PM_PAR_XPOS] = in[PM_PAR_XPOS]; + out[PM_PAR_YPOS] = in[PM_PAR_YPOS]; + + assert (PM_PAR_YPOS == 4); // so that we copy the rest of the params for (int i = 4; i < 8; i++) { psPolynomial2D *poly = psf->params->data[i-4]; // XXX: Verify this (from EAM change) - //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]); - out[i] = psPolynomial2DEval(poly, out[2], out[3]); + //out[i] = Polynomial2DEval_EAM(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]); + out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]); } return(true); @@ -248,12 +249,12 @@ dP = 0; - dP += PS_SQR(dPAR[4] / PAR[4]); - dP += PS_SQR(dPAR[5] / PAR[5]); + dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); + dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); dP = sqrt (dP); status = true; status &= (dP < 0.5); - status &= (PAR[1] > 0); - status &= ((dPAR[1]/PAR[1]) < 0.5); + status &= (PAR[PM_PAR_FLUX] > 0); + status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5); if (!status) Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 9513) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 9526) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-13 02:23:55 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -136,12 +136,12 @@ // NORM-only model fits only source normalization (Io) nParams = 1; - paramMask->data.U8[1] = 0; + paramMask->data.U8[PM_PAR_FLUX] = 0; break; case PM_SOURCE_FIT_PSF: // PSF model only fits x,y,Io nParams = 3; - paramMask->data.U8[1] = 0; - paramMask->data.U8[2] = 0; - paramMask->data.U8[3] = 0; + paramMask->data.U8[PM_PAR_FLUX] = 0; + paramMask->data.U8[PM_PAR_XPOS] = 0; + paramMask->data.U8[PM_PAR_YPOS] = 0; break; case PM_SOURCE_FIT_EXT: @@ -149,5 +149,5 @@ nParams = params->n - 1; psVectorInit (paramMask, 0); - paramMask->data.U8[0] = 1; + paramMask->data.U8[PM_PAR_SKY] = 1; break; case PM_SOURCE_FIT_PSF_AND_SKY: @@ -226,8 +226,8 @@ // models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + onPic &= (params->data.F32[PM_PAR_XPOS] >= source->pixels->col0); + onPic &= (params->data.F32[PM_PAR_XPOS] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0); + onPic &= (params->data.F32[PM_PAR_YPOS] < source->pixels->row0 + source->pixels->numRows); if (!onPic) { model->status = PM_MODEL_OFFIMAGE; Index: trunk/psModules/src/objects/pmSourceIO_CMF.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c (revision 9513) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c (revision 9526) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-22 12:24:17 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-13 02:23:55 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -63,5 +63,17 @@ type = pmSourceDophotType (source); - lsky = (PAR[PM_PAR_SKY] < 1.0) ? 0.0 : log10(PAR[PM_PAR_SKY]); + #if 1 + + float sky = PAR[PM_PAR_SKY]; + #else // we eventually need to interpolate here to get the sky background. XXX + + file = psMetadataLookupPtr(&status, header, "PSPHOT.BACKGND"); + if (file != NULL) { + background = pmFPAviewThisReadout (view, file->fpa); + } + float sky = source->moments->Sky + 0; // etc. + #endif + + lsky = (sky < 1.0) ? 0.0 : log10(sky); shape.sx = PAR[PM_PAR_SXX];