Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 36085) @@ -16,5 +16,4 @@ * 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 - normalized dev parameter note that a standard dev model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, @@ -49,4 +48,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_DEV.h" @@ -63,7 +63,8 @@ # define PM_MODEL_SET_LIMITS pmModelSetLimits_DEV -// f = exp(-z^0.125) +// f = exp(-kappa*r^(1/index)) +// f = exp(-kappa*z^(0.5/index)) +// index = 4, 0.5/index = 0.125 # define ALPHA 0.125 -// # define ALPHA 0.25 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) @@ -73,5 +74,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits @@ -86,9 +87,7 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 }; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 }; static bool limitsApply = true; // Apply limits? - -# include "pmModel_SERSIC.CP.h" psF32 PM_MODEL_FUNC (psVector *deriv, @@ -109,85 +108,31 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 0.5 / ALPHA; - float par7 = ALPHA; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*pow(z,ALPHA); - psF32 f1 = Io*exp(-f2); - - psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } + // for DEV, we can hard-wire kappa(4): + // float index = 4.0; + float kappa = 7.670628; + + // r = sqrt(z) + float q = kappa*pow(z,ALPHA); + float f0 = exp(-q); + + assert (isfinite(q)); + assert (isfinite(f0)); + + // only worry about the central pixels at most + float radius = hypot(X, Y); + if (radius <= 1.5) { + // Nsub ~ 10*index^2 + 1 + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV")); + int Nsub = 2 * ((int)(25 / axes.minor)) + 1; + Nsub = PS_MIN (Nsub, 121); + Nsub = PS_MAX (Nsub, 11); + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub); } - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - assert (isfinite(f2)); + float f1 = PAR[PM_PAR_I0]*f0; + float f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -195,18 +140,24 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +2.0*f1; // XXX extra damping.. - - // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0); - - assert (isfinite(z1)); - - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = f1*kappa*ALPHA*pow(z,ALPHA); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + } + } + return (f); } @@ -302,14 +253,8 @@ } - // the normalization is modified by the slope - float index = 0.5 / ALPHA; - float bn = 1.9992*index - 0.3271; - float Io = exp(0.5*bn); - // set the model normalization if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { return false; } - PAR[PM_PAR_I0] /= Io; // set the model position @@ -328,17 +273,9 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 4.0; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + float norm = 0.00168012; + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -359,7 +296,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - // f = Io exp(-z^n) -> z^n = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); + // static value for DEV: + float kappa = 7.670628; + + // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * pow(zn, 0.5 / ALPHA); psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 36085) @@ -45,4 +45,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_EXP.h" @@ -65,5 +66,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits @@ -78,9 +79,14 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5}; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5}; static bool limitsApply = true; // Apply limits? -# include "pmModel_SERSIC.CP.h" +// # include "pmModel_SERSIC.CP.h" + +// the problems I'm having with the SERSIC-like functions are: +// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius) +// 2) getting the central pixel right +// 3) getting the derivaties right. psF32 PM_MODEL_FUNC (psVector *deriv, @@ -101,85 +107,26 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 1.0; - float par7 = 0.5; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*sqrt(z); - psF32 f1 = Io*exp(-f2); - + // for EXP, we can hard-wire kappa(1): + // float index = 1.0; + float kappa = 1.70056; + + // sqrt(z) is r + float q = kappa*sqrt(z); + psF32 f0 = exp(-q); + + assert (isfinite(q)); + + // only worry about the central 4 pixels at most psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } - } - - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - assert (isfinite(f2)); + if (radius <= 1.5) { + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51); + } + assert (isfinite(f0)); + + psF32 f1 = PAR[PM_PAR_I0]*f0; + psF32 f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -187,18 +134,24 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +f1; - - // gradient is infinite for z = 0; saturate at z = 0.01 - // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i] - psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z); - - // XXX dampen SXX and SYY as in GAUSS? - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = 0.5*f1*kappa/sqrt(z); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = 0.5*f1*kappa; + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + } + } + return (f); } @@ -314,17 +267,11 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 1.0; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + // static value for EXP: + float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr + + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -345,7 +292,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * zn; + // static value for EXP: + float kappa = 1.70056; + + // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * zn; psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", @@ -501,2 +451,66 @@ return; } + +# if (0) +void bilin_inter_function () { + // first, use Rmajor and index to find the central pixel flux (fraction of total flux) + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + + // get the central pixel flux from the lookup table + float xPix = (axes.major - centralPixelXo) / centralPixeldX; + xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); + float yPix = (index - centralPixelYo) / centralPixeldY; + yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); + + // the integral of a Sersic has an analytical form as follows: + float logGamma = lgamma(2.0*index); + float bnFactor = pow(bn, 2.0*index); + float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; + + // XXX interpolate to get the value + // XXX for the moment, just integerize + // XXX I need to multiply by the integrated flux to get the flux in the central pixel + float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; + + float px1 = 1.0 / PAR[PM_PAR_SXX]; + float py1 = 1.0 / PAR[PM_PAR_SYY]; + float z10 = PS_SQR(px1); + float z01 = PS_SQR(py1); + + // which pixels do we need for this interpolation? + // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) + if ((X >= 0) && (Y >= 0)) { + float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive + float V00 = Vcenter; + float V10 = Io*exp(-bn*pow(z10,par7)); + float V01 = Io*exp(-bn*pow(z01,par7)); + float V11 = Io*exp(-bn*pow(z11,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, X, Y); + } + if ((X < 0) && (Y >= 0)) { + float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative + float V00 = Io*exp(-bn*pow(z10,par7)); + float V10 = Vcenter; + float V01 = Io*exp(-bn*pow(z11,par7)); + float V11 = Io*exp(-bn*pow(z01,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); + } + if ((X >= 0) && (Y < 0)) { + float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative + float V00 = Io*exp(-bn*pow(z01,par7)); + float V10 = Io*exp(-bn*pow(z11,par7)); + float V01 = Vcenter; + float V11 = Io*exp(-bn*pow(z10,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); + } + if ((X < 0) && (Y < 0)) { + float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive + float V00 = Io*exp(-bn*pow(z11,par7)); + float V10 = Io*exp(-bn*pow(z10,par7)); + float V01 = Io*exp(-bn*pow(z01,par7)); + float V11 = Vcenter; + f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); + } +} +# endif Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 36085) @@ -70,15 +70,15 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Moderate parameter limits // Tolerate a small divot (k < 0) static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 }; -static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Strict parameter limits // k = PAR_7 < 0 is very undesirable (big divot in the middle) static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; -static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Parameter limits to use Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 36085) @@ -70,15 +70,15 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Moderate parameter limits // Tolerate a small divot (k < 0) static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 }; -static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Strict parameter limits // k = PAR_7 < 0 is very undesirable (big divot in the middle) static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; -static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Parameter limits to use Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 36085) @@ -66,5 +66,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 36085) @@ -20,5 +20,5 @@ * note that a Sersic model is usually defined in terms of R_e, the half-light radius. This construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian. - Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting: + Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using: shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2) @@ -55,4 +55,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_SERSIC.h" @@ -74,6 +75,6 @@ // Lax parameter limits -static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 }; // Moderate parameter limits @@ -88,9 +89,9 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 }; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0}; static bool limitsApply = true; // Apply limits? -# include "pmModel_SERSIC.CP.h" +// # include "pmModel_SERSIC.CP.h" psF32 PM_MODEL_FUNC (psVector *deriv, @@ -111,91 +112,33 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 0.5 / PAR[PM_PAR_7]; - float par7 = PAR[PM_PAR_7]; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*pow(z,par7); - psF32 f1 = Io*exp(-f2); + float Sindex = 0.5 / PAR[PM_PAR_7]; + float kappa = pmSersicKappa (Sindex); + + float q = kappa*pow(z,PAR[PM_PAR_7]); + psF32 f0 = exp(-q); + + assert (isfinite(q)); psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } - } - - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - if (!isfinite(z0)) { - fprintf(stderr, "z0 is not finite for %f %f %f %f %f. Parameters: \n", X, Y, radius, z, f1); + if (radius <= 1.5) { + // Nsub ~ 10*index^2 + 1 + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC")); + int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1; + Nsub = PS_MIN (Nsub, 121); + Nsub = PS_MAX (Nsub, 11); + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub); + } + if (!isfinite(f0)) { + fprintf(stderr, "f0 is not finite for %f %f %f %f %f. Parameters: \n", X, Y, radius, z, q); fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[7]); } - - assert (isfinite(f2)); + assert (isfinite(f0)); + + psF32 f1 = PAR[PM_PAR_I0]*f0; + psF32 f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -203,22 +146,28 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +f1; - - // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0); - - dPAR[PM_PAR_7] = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z); - dPAR[PM_PAR_7] *= 3.0; - - assert (isfinite(z1)); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + dPAR[PM_PAR_7] = -1.0*f1*q*log(z); + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + // dPAR[PM_PAR_7] = -1.0*f1*q*log(z + 0.0001); + dPAR[PM_PAR_7] = -1.0*f1*q*log(z + 0.0001); // factor of 16 to reduce the gain + } assert (isfinite(dPAR[PM_PAR_7])); - - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag? - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + } + return (f); } @@ -370,17 +319,11 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 0.5 / PAR[PM_PAR_7]; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + float Sindex = 0.5 / PAR[PM_PAR_7]; + float norm = pmSersicNorm (Sindex); + + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -401,7 +344,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + float Sindex = 0.5 / PAR[PM_PAR_7]; + float kappa = pmSersicKappa (Sindex); + // f = Io exp(-z^n) -> z^n = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]); + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * pow(zn, Sindex); psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", Index: trunk/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 35768) +++ trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2, 0.5, -3.3, -0.5 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e+8, +1.0e4, +1.0e4, 150.0, +3.3 , 5.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.00+9, +1.0e5, +1.0e5, 150.0, +3.3 , 5.0 }; // Moderate parameter limits