Index: /trunk/psModules/src/objects/pmModel.h =================================================================== --- /trunk/psModules/src/objects/pmModel.h (revision 8881) +++ /trunk/psModules/src/objects/pmModel.h (revision 8882) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-17 18:01:05 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -49,4 +49,18 @@ } pmModel; + +/** Symbolic names for the elements of [d]params + * Note: these are #defines not enums as a given element of [d]params + * may/will correspond to different parameters in different contexts + */ +#define PM_PAR_SKY 0 // Sky +#define PM_PAR_FLUX 1 // Flux +#define PM_PAR_XPOS 2 // X centre of object +#define PM_PAR_YPOS 3 // Y centre of object +#define PM_PAR_SXX 4 // shape X^2 moment +#define PM_PAR_SYY 5 // shape Y^2 moment +#define PM_PAR_SXY 6 // shape XY moment +#define PM_PAR_7 7 // ??? Unknown parameter +#define PM_PAR_8 8 // ??? Unknown parameter /** pmModelAlloc() Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 8881) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 8882) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -191,5 +191,5 @@ psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); - // measure the chi-square trend as a function of flux (PAR[1]) + // measure the chi-square trend as a function of flux (PAR[PM_PAR_FLUX]) // this should be linear for Poisson errors and quadratic for constant sky errors psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); @@ -206,5 +206,5 @@ mask->data.U8[i] = 1; } else { - flux->data.F64[i] = model->params->data.F32[1]; + flux->data.F64[i] = model->params->data.F32[PM_PAR_FLUX]; chisq->data.F64[i] = model[0].chisq/model[0].nDOF; mask->data.U8[i] = 0; Index: /trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 8881) +++ /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 8882) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -116,14 +116,14 @@ type = pmSourceDophotType (source); - lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); - - shape.sx = PAR[4]; - shape.sy = PAR[5]; - shape.sxy = PAR[6]; + lsky = (PAR[PM_PAR_SKY] < 1.0) ? 0.0 : log10(PAR[PM_PAR_SKY]); + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = psEllipseShapeToAxes (shape); psLineInit (line); - psLineAdd (line, "%6.1f ", PAR[2]); - psLineAdd (line, "%6.1f ", PAR[3]); + psLineAdd (line, "%6.1f ", PAR[PM_PAR_XPOS]); + psLineAdd (line, "%6.1f ", PAR[PM_PAR_YPOS]); psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->psfMag + ZERO_POINT)); psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*source->errMag))); @@ -253,7 +253,7 @@ dPAR = model->dparams->data.F32; - PAR[0] = pow (atof (array->data[5]), 10.0); - PAR[2] = atof (array->data[0]); - PAR[3] = atof (array->data[1]); + PAR[PM_PAR_SKY] = pow (atof (array->data[5]), 10.0); + PAR[PM_PAR_XPOS] = atof (array->data[0]); + PAR[PM_PAR_YPOS] = atof (array->data[1]); source->psfMag = atof (array->data[2]); source->extMag = atof (array->data[6]); Index: /trunk/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 8881) +++ /trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 8882) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -62,5 +62,5 @@ dPAR = model->dparams->data.F32; - dmag = dPAR[1] / PAR[1]; + dmag = dPAR[PM_PAR_FLUX] / PAR[PM_PAR_FLUX]; type = pmSourceDophotType (source); if ((source->apMag < 99.0) && (source->psfMag < 99.0)) { @@ -70,16 +70,16 @@ } - shape.sx = PAR[4]; - shape.sy = PAR[5]; - shape.sxy = PAR[6]; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = psEllipseShapeToAxes (shape); psLineInit (line); psLineAdd (line, "%3d", type); - psLineAdd (line, "%8.2f", PAR[2]); - psLineAdd (line, "%8.2f", PAR[3]); + psLineAdd (line, "%8.2f", PAR[PM_PAR_XPOS]); + psLineAdd (line, "%8.2f", PAR[PM_PAR_YPOS]); psLineAdd (line, "%8.3f", source->psfMag); psLineAdd (line, "%6.3f", dmag); - psLineAdd (line, "%9.2f", PAR[0]); + psLineAdd (line, "%9.2f", PAR[PM_PAR_SKY]); psLineAdd (line, "%9.3f", axes.major); psLineAdd (line, "%9.3f", axes.minor); Index: /trunk/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 8881) +++ /trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 8882) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -84,8 +84,9 @@ // dPos is positional error, dMag is mag error - dPos = hypot (dPAR[2], dPAR[3]); + dPos = hypot (dPAR[PM_PAR_XPOS], dPAR[PM_PAR_YPOS]); fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", - PAR[2], PAR[3], PAR[0], source->psfMag, source->errMag, dPos); + PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_SKY], + source->psfMag, source->errMag, dPos); for (j = 4; j < model->params->n; j++) { @@ -142,8 +143,9 @@ // dPos is shape error // XXX these are hardwired for SGAUSS - dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); + dPos = hypot ((dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]), (dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY])); fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", - PAR[2], PAR[3], PAR[0], source->extMag, source->errMag, dPos); + PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_SKY], + source->extMag, source->errMag, dPos); for (j = 4; j < model->params->n; j++) { Index: /trunk/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_SX.c (revision 8881) +++ /trunk/psModules/src/objects/pmSourceIO_SX.c (revision 8882) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-09-22 12:24:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -60,16 +60,16 @@ // pmSourceSextractType (source, &type, &flags); - shape.sx = PAR[4]; - shape.sy = PAR[5]; - shape.sxy = PAR[6]; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = psEllipseShapeToAxes (shape); psLineInit (line); psLineAdd (line, "%5.2f", 0.0); // should be type - psLineAdd (line, "%11.3f", PAR[2]); - psLineAdd (line, "%11.3f", PAR[3]); + psLineAdd (line, "%11.3f", PAR[PM_PAR_XPOS]); + psLineAdd (line, "%11.3f", PAR[PM_PAR_YPOS]); psLineAdd (line, "%9.4f", source->psfMag); psLineAdd (line, "%9.4f", source->errMag); - psLineAdd (line, "%13.4f", PAR[0]); + psLineAdd (line, "%13.4f", PAR[PM_PAR_SKY]); psLineAdd (line, "%9.2f", axes.major); psLineAdd (line, "%9.2f", axes.minor); Index: /trunk/psphot/src/models/pmModel_GAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_GAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_GAUSS.c (revision 8882) @@ -139,12 +139,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) return true; Index: /trunk/psphot/src/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 8882) @@ -16,20 +16,20 @@ 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; psF32 t = 1 + z + z*z/2.0; psF32 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ - psF32 f = PAR[1]*r + PAR[0]; + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { - // note difference from a pure gaussian: q = PAR[1]*r - psF32 q = PAR[1]*r*r*t; + // note difference from a pure gaussian: q = PAR[PM_PAR_FLUX]*r + psF32 q = PAR[PM_PAR_FLUX]*r*r*t; deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -2.0*q*px*X; deriv->data.F32[5] = -2.0*q*py*Y; @@ -78,7 +78,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 @@ -92,5 +92,5 @@ norm *= 0.01; - psF64 Flux = params->data.F32[1] * Area * norm; + psF64 Flux = params->data.F32[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -154,12 +154,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) return true; Index: /trunk/psphot/src/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 8882) @@ -26,22 +26,22 @@ 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; - psF32 r = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); - psF32 f = PAR[1]*r + PAR[0]; + psF32 r = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7] + QG_S1*pow(z, dQG_S1)); + psF32 t = PAR[PM_PAR_FLUX]*r*r; + psF32 q = t*(PAR[PM_PAR_7] + QG_S1*pow(z, dQG_S1)); deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -2.0*q*px*X; deriv->data.F32[5] = -2.0*q*py*Y; @@ -111,7 +111,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 @@ -120,10 +120,10 @@ norm = 0.0; for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); + f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); norm += f; } norm *= 0.01; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -138,15 +138,15 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) return (1.0); - if (flux >= PAR[1]) return (1.0); + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); + 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 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); psF64 dz = 1.0 / (2.0 * sigma*sigma); - psF64 limit = flux / PAR[1]; + psF64 limit = flux / PAR[PM_PAR_FLUX]; // we can do this much better with intelligent choices here for (z = 0.0; z < 20.0; z += dz) { - f = 1.0 / (1 + PAR[7]*z + pow(z, QG_S1)); + f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); if (f < limit) break; } @@ -184,12 +184,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) return true; Index: /trunk/psphot/src/models/pmModel_RGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_RGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_RGAUSS.c (revision 8882) @@ -20,25 +20,25 @@ 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; // psF32 FACT = 1 + 5*exp(-5*PAR[7]); - psF32 p = pow(z, PAR[7] - 1.0); + psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); psF32 r = 1.0 / (1 + z + z*p); - psF32 f = PAR[1]*r + PAR[0]; + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(1 + PAR[7]*p); + psF32 t = PAR[PM_PAR_FLUX]*r*r; + psF32 q = t*(1 + PAR[PM_PAR_7]*p); deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -q*px*X*2; deriv->data.F32[5] = -q*py*Y*2; @@ -60,7 +60,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 @@ -69,5 +69,5 @@ norm = 0.0; for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + z + pow(z, PAR[7])); + f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); norm += f; } @@ -87,15 +87,15 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) return (1.0); - if (flux >= PAR[1]) return (1.0); + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); + 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 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); psF64 dz = 1.0 / (2.0 * sigma*sigma); - psF64 limit = flux / PAR[1]; + psF64 limit = flux / PAR[PM_PAR_FLUX]; // we can do this much better with intelligent choices here for (z = 0.0; z < 20.0; z += dz) { - p = pow(z, PAR[7]); + p = pow(z, PAR[PM_PAR_7]); f = 1.0 / (1 + z + p); if (f < limit) break; Index: /trunk/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 8882) @@ -25,27 +25,27 @@ 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 = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8); - // note that if z -> 0, dPAR[7] -> -inf - // also z^(PAR[7]-1) -> Inf - - psF32 pr = z*PAR[8]; + 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 = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y), 1e-8); + // note that if z -> 0, dPAR[PM_PAR_7] -> -inf + // also z^(PAR[PM_PAR_7]-1) -> Inf + + psF32 pr = z*PAR[PM_PAR_8]; psF32 pr3 = pr*pr*pr; - psF32 p = pow(z, PAR[7] - 1.0); + psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); psF32 r = 1.0 / (1 + z*p + pr*pr3); - psF32 f = PAR[1]*r + PAR[0]; + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3); + psF32 t = PAR[PM_PAR_FLUX]*r*r; + psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR[PM_PAR_8]*pr3); deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -2.0*q*px*X; deriv->data.F32[5] = -2.0*q*py*Y; @@ -230,7 +230,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 @@ -239,11 +239,11 @@ norm = 0.0; for (z = 0.005; z < 50; z += 0.01) { - psF32 pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + psF32 pr = PAR[PM_PAR_8]*z; + f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); norm += f; } norm *= 0.01; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -262,21 +262,21 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) return (1.0); - if (flux >= PAR[1]) return (1.0); + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); + if (flux >= PAR[PM_PAR_FLUX]) return (1.0); // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; + shape.sx = 1.0 / PAR[PM_PAR_SXX]; + shape.sy = 1.0 / PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = EllipseShapeToAxes (shape); psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; + psF64 limit = flux / PAR[PM_PAR_FLUX]; // XXX : we can do this faster with an intelligent starting choice for (r = 0.0; r < 20.0; r += dr) { z = SQ(r); - pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + pr = PAR[PM_PAR_8]*z; + f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); if (f < limit) break; } @@ -315,26 +315,26 @@ psF32 *dPAR = model->dparams->data.F32; - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; + shape.sx = 1.0 / PAR[PM_PAR_SXX]; + shape.sy = 1.0 / PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = EllipseShapeToAxes (shape); dP = 0; - dP += PS_SQR(dPAR[4] / PAR[4]); - dP += PS_SQR(dPAR[5] / PAR[5]); - dP += PS_SQR(dPAR[7] / PAR[7]); + dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); + dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); + dP += PS_SQR(dPAR[PM_PAR_7] / PAR[PM_PAR_7]); dP = sqrt (dP); status = true; status &= (dP < 0.5); - status &= (PAR[1] > 0); - status &= ((dPAR[1]/PAR[1]) < 0.5); - status &= (fabs(PAR[8]) < 0.5); - status &= (dPAR[8] < 0.1); + status &= (PAR[PM_PAR_FLUX] > 0); + status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5); + status &= (fabs(PAR[PM_PAR_8]) < 0.5); + status &= (dPAR[PM_PAR_8] < 0.1); status &= (axes.major > 1.41); status &= (axes.minor > 1.41); status &= ((axes.major / axes.minor) < 5.0); - status &= (PAR[7] > 0.5); + status &= (PAR[PM_PAR_7] > 0.5); if (status) return true; Index: /trunk/psphot/src/models/pmModel_TAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_TAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_TAUSS.c (revision 8882) @@ -1,3 +1,2 @@ - /****************************************************************************** params->data.F32[0] = So; @@ -150,12 +149,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/psphot/src/models/pmModel_TGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 8882) @@ -27,22 +27,22 @@ 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; psF32 er = exp(TRF*z); - psF32 r = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); // (1/R) - psF32 f = PAR[1]*r + PAR[0]; + psF32 r = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S)); // (1/R) + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { - psF32 t = PAR[1]*r*r; // df/dR - psF32 q = t*(TRF*er + PAR[7] + TG_S*pow(z, dTG_S)); // (df/dR)(dR/dz) + psF32 t = PAR[PM_PAR_FLUX]*r*r; // df/dR + psF32 q = t*(TRF*er + PAR[PM_PAR_7] + TG_S*pow(z, dTG_S)); // (df/dR)(dR/dz) deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -2.0*q*px*X; deriv->data.F32[5] = -2.0*q*py*Y; @@ -112,7 +112,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 @@ -122,10 +122,10 @@ for (z = 0.005; z < 50; z += 0.01) { psF32 er = exp(TRF*z); - f = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); + f = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S)); norm += f; } norm *= 0.01; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -140,16 +140,16 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) return (1.0); - if (flux >= PAR[1]) return (1.0); + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); + 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 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); psF64 dz = 1.0 / (2.0 * sigma*sigma); - psF64 limit = flux / PAR[1]; + psF64 limit = flux / PAR[PM_PAR_FLUX]; // we can do this much better with intelligent choices here for (z = 0.0; z < 20.0; z += dz) { psF32 er = exp(TRF*z); - f = 1.0 / (er + PAR[7]*z + pow(z, TG_S)); + f = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S)); if (f < limit) break; } @@ -187,12 +187,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) return true; Index: /trunk/psphot/src/models/pmModel_ZGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_ZGAUSS.c (revision 8881) +++ /trunk/psphot/src/models/pmModel_ZGAUSS.c (revision 8882) @@ -23,24 +23,24 @@ 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; psF32 pr = PAR8*z; - psF32 p = pow(z, PAR[7] - 1.0); + psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); - psF32 f = PAR[1]*r + PAR[0]; + psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR8*pr*pr*pr); + psF32 t = PAR[PM_PAR_FLUX]*r*r; + psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR8*pr*pr*pr); deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); deriv->data.F32[4] = -q*px*X; deriv->data.F32[5] = -q*py*Y; @@ -57,7 +57,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 @@ -67,10 +67,10 @@ psF32 pr = PAR8*z; for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); norm += f; } norm *= 0.01; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; return(Flux); @@ -89,15 +89,15 @@ if (flux <= 0) return (1.0); - if (PAR[1] <= 0) return (1.0); - if (flux >= PAR[1]) return (1.0); + if (PAR[PM_PAR_FLUX] <= 0) return (1.0); + if (flux >= PAR[PM_PAR_FLUX]) return (1.0); // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; + shape.sx = 1.0 / PAR[PM_PAR_SXX]; + shape.sy = 1.0 / PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; axes = EllipseShapeToAxes (shape); psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; + psF64 limit = flux / PAR[PM_PAR_FLUX]; // XXX : we can do this faster with an intelligent starting choice @@ -105,5 +105,5 @@ z = SQ(r); pr = PAR8*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); if (f < limit) break; }