Index: /trunk/psphot/src/models/pmModel_STRAIL.c =================================================================== --- /trunk/psphot/src/models/pmModel_STRAIL.c (revision 10414) +++ /trunk/psphot/src/models/pmModel_STRAIL.c (revision 10414) @@ -0,0 +1,568 @@ + +/****************************************************************************** + 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] = length; + params->data.F32[8] = theta; +*****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_STRAIL +# define PM_MODEL_FLUX pmModelFlux_STRAIL +# define PM_MODEL_GUESS pmModelGuess_STRAIL +# define PM_MODEL_LIMITS pmModelLimits_STRAIL +# define PM_MODEL_RADIUS pmModelRadius_STRAIL +# define PM_MODEL_FROM_PSF pmModelFromPSF_STRAIL +# define PM_MODEL_FIT_STATUS pmModelFitStatus_STRAIL + +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 trailLength = PAR[7]; + psF32 theta = PAR[8]; + + psF32 x0 = PAR[2]; //streak center + psF32 y0 = PAR[3]; //streak center + + //S values (1/sigma for x and y case, sigma for xy case) + psF32 sx = PAR[4]; + psF32 sy = PAR[5]; + psF32 sxy = PAR[6]; + + psF32 sinT=sin(theta); + psF32 cosT=cos(theta); + psF32 sin2T=sin(2.0*theta); + psF32 cos2T=cos(2.0*theta); + + // printf("Trying object at %4.1f,%4.1f with length %3.1f and angle %1.3f\r", x0, y0, length, theta); + + //current location relative to trail center + psF32 X = x->data.F32[0] - x0; + psF32 Y = x->data.F32[1] - y0; + + //x' and y' location (trail-orienter coords) + psF32 xs = X*cosT + Y*sinT; + psF32 ys = -1.0*X*sinT + Y*cosT; + + //initialize variables to be changed below + psF32 x1 = 0; + psF32 y1 = 0; + psF32 px = 0; + psF32 py = 0; + psF32 z = 0; + psF32 zx = 0; + psF32 t = 0; + psF32 tx = 0; + psF32 r = 0; + psF32 rx = 0; + psF32 f = 0; + + psF32 sxrot = 0; + psF32 syrot = 0; + psF32 sxyrot = 0; + psF32 dsxrot = 0; + psF32 dsyrot = 0; + psF32 dsxyrot = 0; + + // psF32 Rx = 0; + // psF32 Ry = 0; + // psF32 Rxy = 0; + + + //calculate new S values (1/sigma) for rotated frame + psF32 sxrotsq = PS_SQR(cosT*sx) + PS_SQR(sinT*sy) + cosT*sinT*sxy; + psF32 syrotsq = PS_SQR(cosT*sy) + PS_SQR(sinT*sx) - cosT*sinT*sxy; + + // psF32 testtwo=10.1; + // psF32 testone=fabsf(testtwo); + // fprintf (stderr, "Test: %f is the absolute value of %f?\n",testone,testtwo); + if (sxrotsq<0) { + sxrot = sqrt(-(sxrotsq)); + syrot = sqrt(syrotsq); + fprintf (stderr, "error in sxrotsq: Neg, sxrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",sxrotsq,sx,sy,sxy,theta); + } else if (syrotsq<0) { + sxrot = sqrt(sxrotsq); + syrot = sqrt(-(syrotsq)); + fprintf (stderr, "error in syrotsq: Neg, syrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",syrotsq,sx,sy,sxy,theta); + } else if (sxrotsq==0){ + sxrot = 0.01; + syrot = sqrt(syrotsq); + fprintf (stderr, "error in sxrotsq: Zero, sxrotsq=%f \n",sxrotsq); + } else if (syrotsq==0) { + syrot = 0.01; + sxrot = sqrt(sxrotsq); + fprintf (stderr, "error in syrotsq: Zero, syrotsq=%f \n",syrotsq); + // return(0); + }else { + sxrot = sqrt(sxrotsq); + syrot = sqrt(syrotsq); + } + + sxyrot = sxy*cos2T + (PS_SQR(sy) - PS_SQR(sx))*sin2T; + + if (isnan(sxrot)) { + fprintf (stderr, "error in sxrot syrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",syrot,sx,sy,sxy,cosT,sinT); + } else if (isnan(syrot)) { + fprintf (stderr, "error in syrot sxrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",sxrot,sx,sy,sxy,cosT,sinT); + } + + //calculate length of pipe (not of trail motion) + psF32 length = trailLength - 2.0*2.0/sxrot; + + + if (xs < length/(-2.0)) { + x1 = (xs+length/2.0)*cosT - ys*sinT; //Endpoint1 + y1 = (xs+length/2.0)*sinT + ys*cosT; //Endpoint1 + //1.6 factor comes from by-eye testing of fits, sqrt from the later squaring of it + //1.6 ~= phi (golden mean)...coincidence? + px = sxrot*x1/sqrt(1.6); + py = syrot*y1; + + //first find out what the falloff in the x direction is... + zx = 0.5*PS_SQR(px); + tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; + rx = 1.0 / (tx + zx*zx*zx*zx/24.0); + + //...and now in the y direction + z = 0.5*PS_SQR(py)+sxyrot*x1*y1; + t = 1 + z + z*z/2.0; + r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + f = PAR[1]*rx*r + PAR[0]; + + } else if (xs > length/2.0){ + x1 = (xs-length/2.0)*cosT - ys*sinT; //Endpoint2 + y1 = (xs-length/2.0)*sinT + ys*cosT; //Endpoint2 + px = sxrot*x1/sqrt(1.6); + py = syrot*y1; + zx = 0.5*PS_SQR(px); + tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; + rx = 1.0 / (tx + zx*zx*zx*zx/24.0); + + z = 0.5*PS_SQR(py)+sxyrot*x1*y1; + t = 1 + z + z*z/2.0; + r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + f = PAR[1]*rx*r + PAR[0]; + + if (isnan(r)) { + fprintf (stderr, "error in +r t=%f z=%f\n",t,z); + } + + } else { + x1 = -ys*sinT; + y1 = ys*cosT; + px = sx*x1; + py = sy*y1; + z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + sxy*x1*y1; + t = 1 + z + z*z/2.0; + r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + rx = 1.0; //so that dF/dF0 can be generalized + f = PAR[1]*r + PAR[0]; + + if (isnan(r)) { + fprintf (stderr, "error in midr t=%f z=%f\n",t,z); + } + } + + + //ok...so this is df/dPAR[X] + if (deriv != NULL) { + psF32 q = 1; + //stable + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r * rx; + + if (isnan(deriv->data.F32[0])) { + fprintf (stderr, "error in deriv0\n"); + } else if (isnan(deriv->data.F32[1])) { + fprintf (stderr, "error in deriv1 r=%f rx=%f\n",r,rx); + } + + dsxrot = 2.0*PS_SQR(sy)*sinT*cosT - 2.0*PS_SQR(sx)*sinT*cosT + sxy*(PS_SQR(cosT)-PS_SQR(sinT)); + dsyrot = 2.0*PS_SQR(sx)*sinT*cosT - 2.0*PS_SQR(sy)*sinT*cosT - sxy*(PS_SQR(cosT)-PS_SQR(sinT)); + dsxyrot = 2.0*cos2T*(PS_SQR(sy)-PS_SQR(sx)) - 2.0*sxy*sin2T; + + if (isnan(dsxrot)) { + fprintf (stderr, "error in dsxrot\n"); + } else if (isnan(dsyrot)) { + fprintf (stderr, "error in dsyrot\n"); + } else if (isnan(dsxyrot)) { + fprintf (stderr, "error in dsxyrot\n"); + } + + //initialize variables definied in the if statements + // psF32 XXX=0; + // psF32 YYY=0; + + + // variable over piecewise func + // change the endcaps to be 4th order gaussians with sxrot_fit=1.6*sxrot + // y' direction can ge adequately modelled by a 3rd order gaussian with syrot + if (xs < length/(-2.0)) { + + q=PAR[1]*r*rx; + deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); + deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); + deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT+2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T+y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); + deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)+2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T+y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); + deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1+2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1+2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T+y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)+x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); + deriv->data.F32[7] = -q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); + deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)+y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)+x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); + + /* if (isnan(deriv->data.F32[2])) { + fprintf (stderr, "error in deriv2\n"); + } else if (isnan(deriv->data.F32[3])) { + fprintf (stderr, "error in deriv3\n"); + } else if (isnan(deriv->data.F32[4])) { + fprintf (stderr, "error in deriv4\n"); + } else if (isnan(deriv->data.F32[5])) { + fprintf (stderr, "error in deriv5\n"); + } else if (isnan(deriv->data.F32[6])) { + fprintf (stderr, "error in deriv6\n"); + } else if (isnan(deriv->data.F32[7])) { + fprintf (stderr, "error in deriv7\n"); + } else if (isnan(deriv->data.F32[8])) { + fprintf (stderr, "error in deriv8\n"); + } + */ + + + + } else if (xs > length/2.0){ + + q=PAR[1]*r*rx; + deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); + deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); + deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT-2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T-y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); + deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)-2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T-y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); + deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1-2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1-2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T-y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)-x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); + deriv->data.F32[7] = q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); + deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)-y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)-x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); + + + /* if (isnan(deriv->data.F32[2])) { + fprintf (stderr, "error in deriv2\n"); + } else if (isnan(deriv->data.F32[3])) { + fprintf (stderr, "error in deriv3\n"); + } else if (isnan(deriv->data.F32[4])) { + fprintf (stderr, "error in deriv4\n"); + } else if (isnan(deriv->data.F32[5])) { + fprintf (stderr, "error in deriv5\n"); + } else if (isnan(deriv->data.F32[6])) { + fprintf (stderr, "error in deriv6\n"); + } else if (isnan(deriv->data.F32[7])) { + fprintf (stderr, "error in deriv7\n"); + } else if (isnan(deriv->data.F32[8])) { + fprintf (stderr, "error in deriv8\n"); + } + */ + + } else { + // this does not change from before, as the y' falloff can be modelled by the standard 3rd order gaussian + // note difference from a pure gaussian: q = PAR[1]*r + q = PAR[1]*r*r*t; + deriv->data.F32[2] = q*(PS_SQR(sx*sinT)*x1 - PS_SQR(sy)*y1*cosT*sinT - sxy*x1*sinT*cosT + sxy*y1*PS_SQR(sinT)); + deriv->data.F32[3] = q*(-1*PS_SQR(sx)*x1*sinT*cosT + PS_SQR(sy*cosT)*y1 + sxy*x1*PS_SQR(cosT) - sxy*y1*sinT*cosT); + deriv->data.F32[4] = -q*sx*PS_SQR(x1); + deriv->data.F32[5] = -q*sy*PS_SQR(y1); + deriv->data.F32[6] = -q*x1*y1; + deriv->data.F32[7] = 0; + deriv->data.F32[8] = -q*( PS_SQR(sx)*x1*(xs*sinT - ys*cosT) + PS_SQR(sy)*y1*(xs*cosT - ys*sinT) + sxy*x1*(xs*cosT - ys*sinT) + sxy*y1*(xs*sinT - ys*cosT) ); + + if (isnan(deriv->data.F32[2])) { + fprintf (stderr, "error in deriv2\n"); + } else if (isnan(deriv->data.F32[3])) { + fprintf (stderr, "error in deriv3\n"); + } else if (isnan(deriv->data.F32[4])) { + fprintf (stderr, "error in deriv4\n"); + } else if (isnan(deriv->data.F32[5])) { + fprintf (stderr, "error in deriv5\n"); + } else if (isnan(deriv->data.F32[6])) { + fprintf (stderr, "error in deriv6\n"); + } else if (isnan(deriv->data.F32[7])) { + fprintf (stderr, "error in deriv7\n"); + } else if (isnan(deriv->data.F32[8])) { + fprintf (stderr, "error in deriv8\n"); + } + + + } + } + return(f); +} + +//fixed +// XXX this needs to apply the axis ratio limits to prevent avoid solutions +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) { + + float beta_lim = 0; + float params_min = 0; + float params_max = 0; + float f1, f2, q1; + float q2 = 0; + + // we need to calculate the limits for SXY specially + if (nParam == PM_PAR_SXY) { + f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + assert (q1 > 0); + q2 = 0.5*sqrt (q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: + switch (nParam) { + case PM_PAR_SKY: beta_lim = 1000; break; + case PM_PAR_I0: beta_lim = 10000; break; // too small? + case PM_PAR_XPOS: beta_lim = 50; break; + case PM_PAR_YPOS: beta_lim = 50; break; + case PM_PAR_SXX: beta_lim = 0.5; break; + case PM_PAR_SYY: beta_lim = 0.5; break; + case PM_PAR_SXY: beta_lim = 1.0; break; // set this to q2? + case 7: beta_lim = 10.0; break; + case 8: beta_lim = M_PI/6.0; break; + + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam); + } + if (fabs(beta[nParam]) > fabs(beta_lim)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); + return false; + } + return true; + case PS_MINIMIZE_PARAM_MIN: + switch (nParam) { + case PM_PAR_SKY: params_min = -1000; break; + case PM_PAR_I0: params_min = 0; break; + case PM_PAR_XPOS: params_min = -100; break; + case PM_PAR_YPOS: params_min = -100; break; + case PM_PAR_SXX: params_min = 0.5; break; + case PM_PAR_SYY: params_min = 0.5; break; + case PM_PAR_SXY: params_min = -5.0; break; // set this to -q2? + case 7: params_min = 0; break; + case 8: params_min = -1*M_PI; break; + + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam); + } + if (params[nParam] < params_min) { + params[nParam] = params_min; + return false; + } + return true; + case PS_MINIMIZE_PARAM_MAX: + switch (nParam) { + case PM_PAR_SKY: params_max = 1e5; break; + case PM_PAR_I0: params_max = 1e8; break; + case PM_PAR_XPOS: params_max = 1e4; break; + case PM_PAR_YPOS: params_max = 1e4; break; + case PM_PAR_SXX: params_max = 100; break; + case PM_PAR_SYY: params_max = 100; break; + case PM_PAR_SXY: params_max = +q2; break; + case 7: params_max = 150; break; + case 8: params_max = M_PI; break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam); + } + if (params[nParam] > params_max) { + params[nParam] = params_max; + return false; + } + return true; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits"); + } + psAbort ("psModules.pmModel_GAUSS", "should not reach here"); + return false; +} + +//fixed +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]); + psF32 Rx=2./PS_SQR(PAR[4]); + psF32 Ry=2./PS_SQR(PAR[5]); + psF32 Rxy=PAR[6]; + + + psF32 theta = PAR[8]; + psF32 sinT=sin(theta); + psF32 cosT=cos(theta); + + psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y + + psF64 A4 = Syrot*PAR[7]; + + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3) + A4; + // Area is equivalent to 2 pi sigma^2 + rectangle + + // 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 + z*z/2 + z*z*z/6); + norm += f; + } + norm *= 0.01; + + psF64 Flux = params->data.F32[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// also prevent 0 returns, and just send a v small number +// return the radius which yields the requested flux + +//fixed, but need to change how it is called to accomodate 2 radii +psF64 PM_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); + + psF32 *PAR = params->data.F32; + psF32 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + + psF32 theta = PAR[8]; + psF32 sinT=sin(theta); + psF32 cosT=cos(theta); + psF32 Rx=2./PS_SQR(PAR[4]); + psF32 Ry=2./PS_SQR(PAR[5]); + psF32 Rxy=PAR[6]; + psF32 length=PAR[7]; + + psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y + + psF64 radius = 0; + if (flux > 0){ + psF64 radius0 = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + psF64 radius1 = Syrot * sqrt (2.0 * log(params->data.F32[1] / flux)); + + if (radius0 > radius1) { + radius=radius0+length/2.0; + } else { + radius=radius1+length/2.0; + } + } else { + radius = 1000; + } + + if (radius < 0.01){ + radius = 0.01; + } + + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +//fixed I think...no good way of guessing as far as I can tell +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { + + pmMoments *Smoments = source->moments; + 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; + //sometimes these moment inputs are zero...why? + + // solve the math to go from Moments To Shape + // limit the axis ratio to 20 (a guess) + axes = psEllipseMomentsToAxes(moments, 20.0); + shape = psEllipseAxesToShape(axes); + + params[0] = Smoments->Sky; + params[1] = Smoments->Peak - Smoments->Sky; + params[2] = Smoments->x; + params[3] = Smoments->y; + params[7] = 2 * axes.major; + params[8] = axes.theta; + + if (moments.x2 == 0 || moments.y2 == 0){ + params[4] = 2.0; + params[5] = 2.0; + params[6] = 0.0; + } else { + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + } + + // printf("Who is NaN? momx: %4.3f momy: %4.3f momxy: %4.3f\n", moments.x2,moments.y2, moments.xy); + + return(true); +} + +//fixed +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]; + out[7] = in[7]; + out[8] = in[8]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params_NEW->data[i-4]; + out[i] = psPolynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +//done I think +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[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + // status &= ((dPAR[7]/PAR[7]) < 0.5); + // status &= ((dPAR[8]/PAR[8]) < 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/psphot/src/psphotModelGroupInit.c =================================================================== --- /trunk/psphot/src/psphotModelGroupInit.c (revision 10413) +++ /trunk/psphot/src/psphotModelGroupInit.c (revision 10414) @@ -5,7 +5,9 @@ # include "models/pmModel_TEST1.c" +# include "models/pmModel_STRAIL.c" static pmModelGroup userModels[] = { {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1, pmModelFlux_TEST1, pmModelRadius_TEST1, pmModelLimits_TEST1, pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelFitStatus_TEST1}, + {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL, pmModelFlux_STRAIL, pmModelRadius_STRAIL, pmModelLimits_STRAIL, pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelFitStatus_STRAIL}, };