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; }