Index: /trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 9775) @@ -1,4 +1,4 @@ /****************************************************************************** - * this file defines the PGAUSS source shape model. Note that these model functions are loaded + * this file defines the GAUSS source shape model. Note that these model functions are loaded * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own. The * models use a psVector to represent the set of parameters, with the sequence used to specify @@ -6,4 +6,7 @@ * specifics of the model. All models which are used a PSF representations share a few * parameters, for which # define names are listed in pmModel.h: + + pure Gaussian: + exp(-z) * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future @@ -114,5 +117,5 @@ shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); - shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; // Area is equivalent to 2 pi sigma^2 @@ -131,19 +134,19 @@ psEllipseShape shape; + psF32 *PAR = params->data.F32; + if (flux <= 0) return (1.0); - if (params->data.F32[PM_PAR_I0] <= 0) + if (PAR[PM_PAR_I0] <= 0) return (1.0); - if (flux >= params->data.F32[PM_PAR_I0]) + if (flux >= PAR[PM_PAR_I0]) return (1.0); - - psF32 *PAR = params->data.F32; shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); - shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; psEllipseAxes axes = psEllipseShapeToAxes (shape); - psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); return (radius); } Index: /trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 9775) @@ -6,4 +6,7 @@ * specifics of the model. All models which are used a PSF representations share a few * parameters, for which # define names are listed in pmModel.h: + + Gaussian taylor expansion + 1 / (1 + z + z^2/2 + z^3/6) * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future @@ -115,5 +118,5 @@ shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); - shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; // Area is equivalent to 2 pi sigma^2 @@ -146,25 +149,27 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + if (flux <= 0) return (1.0); - if (params->data.F32[PM_PAR_I0] <= 0) + if (PAR[PM_PAR_I0] <= 0) return (1.0); - if (flux >= params->data.F32[PM_PAR_I0]) + if (flux >= PAR[PM_PAR_I0]) return (1.0); - - psF32 *PAR = params->data.F32; shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); - shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; // this estimates the radius assuming f(z) is roughly exp(-z) psEllipseAxes axes = psEllipseShapeToAxes (shape); - psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); if (isnan(radius)) - psAbort ("psphot.model", "error in code: never return invalid radius"); + psAbort ("psphot.model", "error in code: radius is NaN"); if (radius < 0) - psAbort ("psphot.model", "error in code: never return invalid radius"); + psAbort ("psphot.model", "error in code: radius is negative"); return (radius); @@ -215,7 +220,5 @@ status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); - if (status) - return true; - return false; + return status; } Index: /trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 9775) @@ -1,32 +1,46 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - one component, two slopes: - 1 / (1 + z^M + z^N) + * this file defines the QGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: - params->data.F32[PM_PAR_SKY] = So; - params->data.F32[PM_PAR_I0] = 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] = -*****************************************************************************/ - -psF32 pmModelFunc_QGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) + power-law with fitted linear term + 1 / (1 + kz + z^2.25) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - amplitude of the linear component (k) + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_QGAUSS +# define PM_MODEL_FLUX pmModelFlux_QGAUSS +# define PM_MODEL_GUESS pmModelGuess_QGAUSS +# define PM_MODEL_LIMITS pmModelLimits_QGAUSS +# define PM_MODEL_RADIUS pmModelRadius_QGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_QGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) { psF32 *PAR = params->data.F32; - 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 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; + + // XXX if the elliptical contour is defined in valid way, this step should not be required. + // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values + // for negative values of z if (z < 0) z = 0; @@ -45,17 +59,17 @@ psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp); - dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +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; + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; + dPAR[PM_PAR_7] = -t*z; } return(f); } -bool pmModelLimits_QGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) { @@ -77,6 +91,6 @@ 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_SXX] = 0.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; params_min[0][0].data.F32[PM_PAR_7] = 0.1; @@ -86,6 +100,6 @@ 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_SXX] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; params_max[0][0].data.F32[PM_PAR_7] = 10.0; @@ -95,35 +109,37 @@ // make an initial guess for parameters -bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source) +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { pmMoments *moments = source->moments; pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - - params[PM_PAR_SKY] = moments->Sky; - params[PM_PAR_I0] = 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; + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = peak->x; + PAR[PM_PAR_YPOS] = peak->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = 0.0; + PAR[PM_PAR_7] = 1.0; return(true); } -psF64 pmModelFlux_QGAUSS(const psVector *params) -{ - float z; - float norm; +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + psEllipseShape shape; psF32 *PAR = params->data.F32; - 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); + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; // the area needs to be multiplied by the integral of f(z) @@ -150,9 +166,11 @@ // define this function so it never returns Inf or NaN // return the radius which yields the requested flux -psF64 pmModelRadius_QGAUSS (const psVector *params, psF64 flux) +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { psF64 z, f; + int Nstep = 0; + psEllipseShape shape; + psF32 *PAR = params->data.F32; - int Nstep = 0; if (flux <= 0) @@ -163,26 +181,15 @@ return (1.0); - // if Sx == Sy, sigma = Sx == Sy - psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 sigma = axes.major; + psF64 limit = flux / PAR[PM_PAR_I0]; - # if 0 - /* test example will just use both, printing both */ - psF64 dz = 1.0 / (2.0 * sigma*sigma); - - // we can do this much better with intelligent choices here - for (z = 0.0; z < 30.0; z += dz) { - Nstep ++; - 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; - } - // fprintf (stderr, "rad 1: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep); - - # else - - /* use the fact that f is monotonically decreasing */ - z = 0; + // use the fact that f is monotonically decreasing + z = 0; Nstep = 0; @@ -193,4 +200,5 @@ z0 = 0.0; + // perform a type of bisection to find the value 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)); @@ -198,5 +206,4 @@ z = 0.5*(z0 + z1); 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) { z0 = z; @@ -208,15 +215,13 @@ Nstep ++; } - // fprintf (stderr, "rad 2: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep); - # endif - psF64 radius = sigma * sqrt (2.0 * z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } + + if (isnan(radius)) + psAbort ("psphot.model", "error in code: radius is NaN"); + return (radius); } -bool pmModelFromPSF_QGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { @@ -224,20 +229,24 @@ psF32 *in = modelFLT->params->data.F32; - out[PM_PAR_SKY] = in[PM_PAR_SKY]; - out[PM_PAR_I0] = in[PM_PAR_I0]; - out[PM_PAR_XPOS] = in[PM_PAR_XPOS]; - out[PM_PAR_YPOS] = in[PM_PAR_YPOS]; - - assert (PM_PAR_YPOS + 1 == 4); // so starting at 4 is correct - 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[PM_PAR_XPOS], out[PM_PAR_YPOS]); - out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]); - } + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + return(true); } -bool pmModelFitStatus_QGAUSS (pmModel *model) +bool PM_MODEL_FIT_STATUS (pmModel *model) { @@ -258,6 +267,12 @@ status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); - if (!status) - return false; - return true; -} + return status; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 9775) @@ -1,81 +1,164 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - one component, two slopes: - 1 / (1 + z + z^Npow) + * this file defines the RGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = 1 / SigmaX; - params->data.F32[5] = 1 / SigmaY; - params->data.F32[6] = Sxy; - params->data.F32[7] = Npow -*****************************************************************************/ - -psF64 psModelFunc_RGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) + power-law with fitted slope + 1 / (1 + z + z^alpha) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - power-law slope (alpha) + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_RGAUSS +# define PM_MODEL_FLUX pmModelFlux_RGAUSS +# define PM_MODEL_GUESS pmModelGuess_RGAUSS +# define PM_MODEL_LIMITS pmModelLimits_RGAUSS +# define PM_MODEL_RADIUS pmModelRadius_RGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_RGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS + +psF64 PS_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) { 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 FACT = 1 + 5*exp(-5*PAR[7]); - - psF32 p = pow(z, PAR[7] - 1.0); + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + X*Y*PAR[PM_PAR_SXY]; + + 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_I0]*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); - - 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[4] = -q*px*X*2; - deriv->data.F32[5] = -q*py*Y*2; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -5.0*t*log(z)*p*z; - - // deriv->data.F32[4] = -1.8*q*px*X*2; - // deriv->data.F32[5] = -1.8*q*py*Y*2; - // deriv->data.F32[6] = -1.5*q*X*Y; - // deriv->data.F32[7] = -5.0*t*log(z)*p*z; + psF32 *dPAR = deriv->data.F32; + + // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r + psF32 t = PAR[PM_PAR_I0]*r*r; + psF32 q = t*(1 + PAR[PM_PAR_7]*p); + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; + dPAR[PM_PAR_7] = -5.0*t*log(z)*p*z; } return(f); } -psF64 psModelFlux_RGAUSS(const psVector *params) -{ - float f, norm, z; +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) +{ + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 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_I0] = 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.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; + params_min[0][0].data.F32[PM_PAR_7] = 1.25; + + params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; + params_max[0][0].data.F32[PM_PAR_I0] = 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] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; + params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; + params_max[0][0].data.F32[PM_PAR_7] = 4.0; + + return (TRUE); +} + +bool PS_MODEL_GUESS (psModel *model, psSource *source) +{ + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *PAR = model->params->data.F32; + + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments tmpMoments; + + // XXX fix this stuff : should be using correct ellipse relationships... + tmpMoments.x2 = PS_SQR(source->moments->Sx); + tmpMoments.y2 = PS_SQR(source->moments->Sy); + tmpMoments.xy = source->moments->Sxy; + + axes = psEllipseMomentsToAxes(tmpMoments); + shape = psEllipseAxesToShape(axes); + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = peak->x; + PAR[PM_PAR_YPOS] = peak->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = shape.sxy; + PAR[PM_PAR_7] = 2.0; + return(true); +} + +psF64 PS_MODEL_FLUX (const psVector *params) +{ + float norm, z; + psEllipseShape shape; psF32 *PAR = params->data.F32; - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; // the area needs to be multiplied by the integral of f(z) norm = 0.0; - for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + z + pow(z, PAR[7])); - norm += f; - } - norm *= 0.01; - - psF64 Flux = params->data.F32[1] * Area * norm; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); + z += DZ; + f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; return(Flux); @@ -84,63 +167,62 @@ // define this function so it never returns Inf or NaN // return the radius which yields the requested flux -psF64 psModelRadius_RGAUSS (const psVector *params, psF64 flux) +psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) { psF64 z, f, p; + psEllipseShape shape; + psF32 *PAR = params->data.F32; if (flux <= 0) return (1.0); - if (PAR[1] <= 0) + if (PAR[PM_PAR_I0] <= 0) return (1.0); - if (flux >= PAR[1]) + if (flux >= PAR[PM_PAR_I0]) return (1.0); - // if Sx == Sy, sigma = Sx == Sy - psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 sigma = axes.major; + psF64 dz = 1.0 / (2.0 * sigma*sigma); - psF64 limit = flux / PAR[1]; - - // we can do this much better with intelligent choices here - for (z = 0.0; z < 20.0; z += dz) { - p = pow(z, PAR[7]); - f = 1.0 / (1 + z + p); - if (f < limit) - break; + psF64 limit = flux / PAR[PM_PAR_I0]; + + // use the fact that f is monotonically decreasing + z = 0; + Nstep = 0; + + // 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[PM_PAR_7]; + z1 = PS_MAX (z0, z1); + z0 = 0.0; + + // perform a type of bisection to find the value + float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7])); + float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7])); + while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) { + z = 0.5*(z0 + z1); + f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); + if (f > limit) { + z0 = z; + f0 = f; + } else { + z1 = z; + f1 = f; + } + Nstep ++; } psF64 radius = sigma * sqrt (2.0 * z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } + + if (isnan(radius)) + psAbort ("psphot.model", "error in code: radius is NaN"); + return (radius); } -bool psModelGuess_RGAUSS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(source->moments->Sx); - moments.y2 = PS_SQR(source->moments->Sy); - moments.xy = source->moments->Sxy; - - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = 1.0 / shape.sx; - params->data.F32[5] = 1.0 / shape.sy; - params->data.F32[6] = shape.sxy; - params->data.F32[7] = 2.0; - return(true); -} - -bool psModelFromPSF_RGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) +bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { @@ -148,13 +230,50 @@ psF32 *in = modelFLT->params->data.F32; - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 8; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + return(true); } + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + 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[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psModules/src/objects/models/pmModel_SGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_SGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_SGAUSS.c (revision 9775) @@ -1,29 +1,46 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - - /****************************************************************************** - one component, two slopes: - 1 / (1 + z^Npow + St*z^Ntide) + * this file defines the SGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = 1 / SigmaX; - params->data.F32[5] = 1 / SigmaY; - params->data.F32[6] = Sxy; - params->data.F32[7] = Npow - params->data.F32[8] = St -*****************************************************************************/ - -# define SQ(A)((A)*(A)) + power-law with fitted slope and outer tidal radius + 1 / (1 + z^N + kz^4) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - slope of power-law component (N) + * PM_PAR_8 8 - amplitude of the tidal cutoff (k) + *****************************************************************************/ + +/*** + XXXX the model in this file needs to be tested more carefully. + the code for guessing the power-law slope based on the radial profile + is either too slow or does not work well. + fix up the code to follow conventions in the other model function files. +***/ + +XXX broken code + +# define PM_MODEL_FUNC pmModelFunc_SGAUSS +# define PM_MODEL_FLUX pmModelFlux_SGAUSS +# define PM_MODEL_GUESS pmModelGuess_SGAUSS +# define PM_MODEL_LIMITS pmModelLimits_SGAUSS +# define PM_MODEL_RADIUS pmModelRadius_SGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_SGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_SGAUSS + psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image); -psF64 p_psImageGetElementF64(psImage *a, int i, int j); - -psF32 pmModelFunc_SGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *x) { psF32 *PAR = params->data.F32; @@ -61,5 +78,5 @@ } -bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) { @@ -99,4 +116,179 @@ return (TRUE); +} + +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +{ + + pmMoments *sMoments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments moments; + + moments.x2 = PS_SQR(sMoments->Sx); + moments.y2 = PS_SQR(sMoments->Sy); + moments.xy = sMoments->Sxy; + + // solve the math to go from Moments To Shape + axes = psEllipseMomentsToAxes(moments); + shape = psEllipseAxesToShape(axes); + + params[0] = sMoments->Sky; + params[1] = sMoments->Peak - sMoments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + + # if (0) + + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + axes.major *= 2.0; + axes.minor *= 2.0; + f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + + if (f1 > f2) { + params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); + } else { + params[7] = 0.5; + } + # endif + + params[7] = 1.8; + params[8] = 0.1; + + + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + 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]) + PS_SQR(PS_SQR(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF64 r, z = 0.0, pr, f; + psF32 *PAR = params->data.F32; + + psEllipseAxes axes; + psEllipseShape shape; + + if (flux <= 0) + return (1.0); + if (PAR[1] <= 0) + return (1.0); + if (flux >= PAR[1]) + 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]; + + axes = psEllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = PS_SQR(r); + pr = PAR[8]*z; + f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr))); + if (f < limit) + break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 9; 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]); + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + psEllipseAxes axes; + psEllipseShape shape; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = psEllipseShapeToAxes (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 = 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 &= (axes.major > 1.41); + status &= (axes.minor > 1.41); + status &= ((axes.major / axes.minor) < 5.0); + status &= (PAR[7] > 0.5); + + if (status) + return true; + return false; } @@ -142,5 +334,6 @@ } -bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source) +// XXX EAM : test version using flux contours to guess slope +bool PM_MODEL_GUESS_HARD (pmModel *model, pmSource *source) { @@ -148,4 +341,5 @@ pmPeak *peak = source->peak; psF32 *params = model->params->data.F32; + float f1, f2; psEllipseAxes axes; @@ -169,7 +363,5 @@ params[6] = shape.sxy; - # if (0) - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); axes.major *= 2.0; axes.minor *= 2.0; @@ -181,181 +373,14 @@ params[7] = 0.5; } - # endif - - params[7] = 1.8; params[8] = 0.1; - return(true); } -// XXX EAM : test version using flux contours to guess slope -bool pmModelGuess_SGAUSS_HARD (pmModel *model, pmSource *source) -{ - - pmMoments *sMoments = source->moments; - pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - float f1, f2; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(sMoments->Sx); - moments.y2 = PS_SQR(sMoments->Sy); - moments.xy = sMoments->Sxy; - - // solve the math to go from Moments To Shape - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params[0] = sMoments->Sky; - params[1] = sMoments->Peak - sMoments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = 1.0 / shape.sx; - params[5] = 1.0 / shape.sy; - params[6] = shape.sxy; - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - axes.major *= 2.0; - axes.minor *= 2.0; - f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - - if (f1 > f2) { - params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); - } else { - params[7] = 0.5; - } - params[8] = 0.1; - - return(true); -} - -psF64 pmModelFlux_SGAUSS(const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - 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))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 pmModelRadius_SGAUSS (const psVector *params, psF64 flux) -{ - psF64 r, z = 0.0, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - 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]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // 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))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool pmModelFromPSF_SGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 9; 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]); - } - return(true); -} - -bool pmModelFitStatus_SGAUSS (pmModel *model) -{ - - psF32 dP; - bool status; - psEllipseAxes axes; - psEllipseShape shape; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (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 = 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 &= (axes.major > 1.41); - status &= (axes.minor > 1.41); - status &= ((axes.major / axes.minor) < 5.0); - status &= (PAR[7] > 0.5); - - if (status) - return true; - return false; -} +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psModules/src/objects/models/pmModel_TGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_TGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_TGAUSS.c (revision 9775) @@ -1,23 +1,41 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - one component, two slopes: - 1 / (1 + z^M + z^N) + * this file defines the TGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: - 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] = -*****************************************************************************/ - -psF64 psModelFunc_TGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) + power-law with fixed slope and fitted amplitude + 1 / (1 + z + kz^2.2) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - amplitude of high-order component (k) + *****************************************************************************/ + +/*** + XXXX the model in this file needs to be tested more carefully. + fix up the code to follow conventions in the other model function files. +***/ + +XXX broken code + +# define PM_MODEL_FUNC pmModelFunc_TGAUSS +# define PM_MODEL_FLUX pmModelFlux_TGAUSS +# define PM_MODEL_GUESS pmModelGuess_TGAUSS +# define PM_MODEL_LIMITS pmModelLimits_TGAUSS +# define PM_MODEL_RADIUS pmModelRadius_TGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_TGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_TGAUSS + +psF64 PS_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *x) { psF32 *PAR = params->data.F32; @@ -50,37 +68,42 @@ } -psF64 psModelFlux_TGAUSS(const psVector *params) -{ - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - - psF64 Flux = params->data.F32[1] * Area; - - return(Flux); -} - -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 psModelRadius_TGAUSS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1]) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool psModelGuess_TGAUSS (psModel *model, psSource *source) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) +{ + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 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_I0] = 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.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + 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_I0] = 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] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.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); +} + +bool PS_MODEL_GUESS (psModel *model, psSource *source) { @@ -99,5 +122,37 @@ } -bool psModelFromPSF_TGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) +psF64 PS_MODEL_FLUX (const psVector *params) +{ + psF64 A1 = 1 / PS_SQR(params->data.F32[4]); + psF64 A2 = 1 / PS_SQR(params->data.F32[5]); + psF64 A3 = params->data.F32[6]; + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + if (flux <= 0) + return (1.0); + if (params->data.F32[1] <= 0) + return (1.0); + if (flux >= params->data.F32[1]) + return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { @@ -116,2 +171,32 @@ return(true); } + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + 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[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psModules/src/objects/models/pmModel_WAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_WAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_WAUSS.c (revision 9775) @@ -1,19 +1,42 @@ -#ifdef HAVE_CONFIG_H -#include -#endif +/****************************************************************************** + * this file defines the WAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: + + power-law with fitted linear term + 1 / (1 + Az + Bz^2 + z^3/6) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - amplitude of the linear component (A) + * PM_PAR_8 8 - amplitude of the quadratic component (B) + *****************************************************************************/ -/****************************************************************************** - 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; -*****************************************************************************/ +/*** + XXXX the model in this file needs to be tested more carefully. + fix up the code to follow conventions in the other model function files. +***/ -psF64 psModelFunc_WAUSS(psVector *deriv, - const psVector *params, - const psVector *x) +XXX broken code + +# define PM_MODEL_FUNC pmModelFunc_WAUSS +# define PM_MODEL_FLUX pmModelFlux_WAUSS +# define PM_MODEL_GUESS pmModelGuess_WAUSS +# define PM_MODEL_LIMITS pmModelLimits_WAUSS +# define PM_MODEL_RADIUS pmModelRadius_WAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_WAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS + +psF64 PS_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *x) { psF32 X = x->data.F32[0] - params->data.F32[2]; @@ -43,34 +66,42 @@ } -// this is probably wrong since it uses the gauss integral 2 pi sigma^2 -psF64 psModelFlux_WAUSS(const psVector *params) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) { - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - psF64 Flux = params->data.F32[1] * Area; + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); - return(Flux); + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 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_I0] = 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.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + 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_I0] = 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] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.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); } -// return the radius which yields the requested flux -psF64 psModelRadius_WAUSS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1] <= 0) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - return (radius); -} - -bool psModelGuess_WAUSS (psModel *model, psSource *source) +bool PS_MODEL_GUESS (psModel *model, psSource *source) { @@ -90,5 +121,34 @@ } -bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) +// this is probably wrong since it uses the gauss integral 2 pi sigma^2 +psF64 PS_MODEL_FLUX (const psVector *params) +{ + psF64 A1 = 1 / PS_SQR(params->data.F32[4]); + psF64 A2 = 1 / PS_SQR(params->data.F32[5]); + psF64 A3 = params->data.F32[6]; + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + if (flux <= 0) + return (1.0); + if (params->data.F32[1] <= 0) + return (1.0); + if (flux >= params->data.F32[1] <= 0) + return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { @@ -107,2 +167,32 @@ return(true); } + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + 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[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psModules/src/objects/models/pmModel_ZGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 9774) +++ /trunk/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 9775) @@ -1,26 +1,43 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - one component, two slopes: - 1 / (1 + z^Npow + PAR8*z^4) + * this file defines the ZGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * statements of their own. The models use a psVector to represent the set of parameters, with + * the sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = 1 / SigmaX; - params->data.F32[5] = 1 / SigmaY; - params->data.F32[6] = Sxy; - params->data.F32[7] = Npow -*****************************************************************************/ - -# define SQ(A)((A)*(A)) + power-law with fitted slope and tidal cutoff + 1 / (1 + z^N + (Az)^4) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - slope of power-law component (N) + *****************************************************************************/ + +/*** + XXXX the model in this file needs to be tested more carefully. + fix up the code to follow conventions in the other model function files. +***/ + +XXX broken code + +# define PM_MODEL_FUNC pmModelFunc_ZGAUSS +# define PM_MODEL_FLUX pmModelFlux_ZGAUSS +# define PM_MODEL_GUESS pmModelGuess_ZGAUSS +# define PM_MODEL_LIMITS pmModelLimits_ZGAUSS +# define PM_MODEL_RADIUS pmModelRadius_ZGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_ZGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_ZGAUSS + # define PAR8 0.1 -psF64 psModelFunc_ZGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) +psF64 PS_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *x) { psF32 *PAR = params->data.F32; @@ -54,73 +71,42 @@ } -psF64 psModelFlux_ZGAUSS(const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - psF32 pr = PAR8*z; - for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 psModelRadius_ZGAUSS (const psVector *params, psF64 flux) -{ - psF64 r, z, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - 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]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // XXX : we can do this faster with an intelligent starting choice - for (r = 0.0; r < 20.0; r += dr) { - z = SQ(r); - pr = PAR8*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool psModelGuess_ZGAUSS (psModel *model, psSource *source) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) +{ + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 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_I0] = 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.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + 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_I0] = 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] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.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); +} + +bool PS_MODEL_GUESS (psModel *model, psSource *source) { @@ -149,5 +135,73 @@ } -bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) +psF64 PS_MODEL_FLUX (const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + psF32 pr = PAR8*z; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF64 r, z, pr, f; + psF32 *PAR = params->data.F32; + + psEllipseAxes axes; + psEllipseShape shape; + + if (flux <= 0) + return (1.0); + if (PAR[1] <= 0) + return (1.0); + if (flux >= PAR[1]) + 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]; + + axes = psEllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = SQ(r); + pr = PAR8*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + if (f < limit) + break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { @@ -166,2 +220,32 @@ return(true); } + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + 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[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS