Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 42078) @@ -94,4 +94,7 @@ $(SRC)/imfit-qfgauss.$(ARCH).o \ $(SRC)/imfit-qrgauss.$(ARCH).o \ +$(SRC)/imfit-rgauss.$(ARCH).o \ +$(SRC)/imfit-fgauss-pol.$(ARCH).o \ +$(SRC)/imfit-rgauss-pol.$(ARCH).o \ $(SRC)/imfit-trail.$(ARCH).o Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss-pol.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss-pol.c (revision 42078) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss-pol.c (revision 42078) @@ -0,0 +1,84 @@ +# include "imfit.h" + +/** fgaussPol : a real 2D Gaussian **/ + +opihi_flt fgaussPolTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *); +void fgaussPolCL (); + +void fgauss_pol_setup (char *name) { + + if (strcmp(name, "fgauss-pol")) return; + + fitfunc = fgaussPolTD; + imfit_cleanup = fgaussPolCL; + Npar = 7; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, opihi_flt, MAX (Npar, 1)); + bzero (par, Npar*sizeof(opihi_flt)); + ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(opihi_flt)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[5] = get_variable_default ("Zpk", 10000); + par[6] = get_variable_default ("Sg", 0.0); + + opihi_flt Sxx = get_variable_default ("SXg", 2.0) / 2.35; // convert FWHM to sigma / sqrt(2) (Sxx) + opihi_flt Syy = get_variable_default ("SYg", 2.0) / 2.35; // convert FWHM to sigma / sqrt(2) (Syy) + opihi_flt Sxy = get_variable_default ("SXYg", 0); + + // z = (x^2 + y^2) / R^2 + (x^2 - y^2) / T^2 + x y / Q : NOTE Q is not squared to allow positive and negative values + + par[2] = 2.0/sqrt((1.0/SQ(Sxx) + 1.0/SQ(Syy))); // par[2] = R + par[3] = 2.0/sqrt((1.0/SQ(Sxx) - 1.0/SQ(Syy))); // par[3] = T + par[4] = 1.0 / Sxy; // par[4] = Q + + sky = &par[6]; +} + +void fgaussPolCL () { + opihi_flt Sxx = par[2]*par[3] / (2.0*sqrt(SQ(par[3]) + SQ(par[2]))); + opihi_flt Syy = par[2]*par[3] / (2.0*sqrt(SQ(par[3]) - SQ(par[2]))); + opihi_flt Sxy = 1.0 / par[4]; + + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("SXg", Sxx * 2.35 / sqrt(2.0)); + set_variable ("SYg", Syy * 2.35 / sqrt(2.0)); + set_variable ("SXYg", Sxy); + set_variable ("Zpk", par[5]); + set_variable ("Sg", par[6]); +} + +/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */ +opihi_flt fgaussPolTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + opihi_flt X = x - par[0]; + opihi_flt Y = y - par[1]; + + opihi_flt P_R = (SQ(X) + SQ(Y)) / SQ(par[2]); + opihi_flt P_T = (SQ(X) - SQ(Y)) / SQ(par[3]); + opihi_flt P_Q = X * Y / par[4]; + + opihi_flt z = P_R + P_T + P_Q; + + opihi_flt r = exp (-z); + opihi_flt q = par[5]*r; + opihi_flt f = q + par[6]; + + if (dpar != NULL) { + dpar[0] = +q*(2*X/SQ(par[2]) + 2*X/SQ(par[3]) + Y/par[4]); + dpar[1] = +q*(2*Y/SQ(par[2]) - 2*Y/SQ(par[3]) + X/par[4]); + + dpar[2] = 2*q*P_R/par[2]; + dpar[3] = 2*q*P_T/par[3]; + dpar[4] = q*P_Q/par[4]; + + dpar[5] = +r; + dpar[6] = +1; + } + return (f); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c (revision 42078) @@ -1,3 +1,5 @@ # include "imfit.h" + +/** fgauss : a real 2D Gaussian **/ opihi_flt fgaussTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *); Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c (revision 42078) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c (revision 42078) @@ -0,0 +1,76 @@ +# include "imfit.h" + +opihi_flt qgaussTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *); +void qgaussCL (); + +void qgauss_setup (char *name) { + + if (strcmp(name, "qgauss")) return; + + fitfunc = qgaussTD; + imfit_cleanup = qgaussCL; + Npar = 8; + Nfpar = 1; + + /* allocate free and fixed parameters */ + ALLOCATE (par, opihi_flt, MAX (Npar, 1)); + bzero (par, Npar*sizeof(opihi_flt)); + ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(opihi_flt)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[2] = 2.35 / get_variable_default ("SXg", 2.0); + par[3] = 2.35 / get_variable_default ("SYg", 2.0); + par[4] = get_variable_default ("SXYg", 0); + par[5] = get_variable_default ("Zpk", 10000); + par[6] = get_variable_default ("Sg", 0.0); + par[7] = get_variable_default ("Sr", 1.0); + fpar[0] = get_variable_default ("Npow", 2.25); + + sky = &par[6]; +} + +void qgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("SXg", 2.35 / par[2]); + set_variable ("SYg", 2.35 / par[3]); + set_variable ("SXYg", par[4]); + set_variable ("Zpk", par[5]); + set_variable ("Sg", par[6]); + set_variable ("Sr", par[7]); +} + +/* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ +opihi_flt qgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + opihi_flt X, Y, px, py; + opihi_flt z, r, q, f; + + X = x - par[0]; + Y = y - par[1]; + + px = par[2]*X; + py = par[3]*Y; + + z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y; + + r = 1.0 / (1 + par[7]*z + pow(z,fpar[0])); + f = par[5]*r + par[6]; + q = par[5]*SQ(r)*(par[7] + fpar[0]*pow(z,(fpar[0]-1))); + + if (dpar != NULL) { + dpar[0] = q*(2*px*par[2] + par[4]*Y); + dpar[1] = q*(2*py*par[3] + par[4]*X); + dpar[2] = -2*q*px*X*2; + dpar[3] = -2*q*py*Y*2; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + dpar[7] = -par[5]*SQ(r)*z; + } + return (f); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 42078) @@ -25,5 +25,5 @@ par[4] = get_variable_default ("SXYg", 0); par[5] = get_variable_default ("Zpk", 10000); - par[6] = get_variable_default ("Sg", 0.0); + par[6] = get_variable_default ("Sg", 0.0); // sky par[7] = get_variable_default ("Sr", 1.0); fpar[0] = get_variable_default ("Npow", 2.25); @@ -40,8 +40,8 @@ set_variable ("Zpk", par[5]); set_variable ("Sg", par[6]); - set_variable ("Sr", par[7]); + set_variable ("Sr", par[7]); } -/* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ +/* qgauss: (1 + Sr*z + z^Npow)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ opihi_flt qgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { OHANA_UNUSED_PARAM(Npar); @@ -60,7 +60,7 @@ r = 1.0 / (1 + par[7]*z + pow(z,fpar[0])); f = par[5]*r + par[6]; - q = par[5]*SQ(r)*(par[7] + fpar[0]*pow(z,(fpar[0]-1))); if (dpar != NULL) { + q = par[5]*SQ(r)*(par[7] + fpar[0]*pow(z,(fpar[0]-1))); dpar[0] = q*(2*px*par[2] + par[4]*Y); dpar[1] = q*(2*py*par[3] + par[4]*X); @@ -75,2 +75,23 @@ } +/****************************************************************************** + * this file defines the QGAUSS source shape model. Note that these model functions are + * loaded by pmModelClass.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 as a + * PSF representations share a few parameters, for which # define names are listed in + * pmModel.h: + + 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 (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 (k) + *****************************************************************************/ Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c (revision 42078) @@ -59,5 +59,5 @@ r = 1.0 / (1 + fpar[0]*z + pow(z,par[7])); - f = par[5]*r + par[6]; + f = par[5]*r + par[6]; // Io * f(r) + Sky q = par[5]*SQ(r)*(fpar[0] + par[7]*pow(z,(par[7]-1))); Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c (revision 42078) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c (revision 42078) @@ -0,0 +1,87 @@ +# include "imfit.h" + +opihi_flt rgaussTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *); +void rgaussCL (); + +void rgauss_setup (char *name) { + + if (strcmp(name, "rgauss")) return; + + fitfunc = rgaussTD; + cleanup = rgaussCL; + Npar = 10; + Nfpar = 1; + + /* allocate free and fixed parameters */ + ALLOCATE (par, opihi_flt, MAX (Npar, 1)); + bzero (par, Npar*sizeof(opihi_flt)); + ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(opihi_flt)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0); + par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0); + par[4] = 0.0; + par[5] = get_variable_default ("Zpk", 10000); + par[6] = get_variable_default ("Sg", 0.0); + par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0); + par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0); + par[9] = get_variable_default ("SXYf", 0.0); + + fpar[0] = get_variable_default ("Npow", 2.25); + + sky = &par[6]; +} + +/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ +opihi_flt rgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + opihi_flt X, Y, px1, py1, px2, py2; + opihi_flt z1, z2, r, q1, q2, f; + + X = x - par[0]; + Y = y - par[1]; + + px1 = par[2]*X; + py1 = par[3]*Y; + px2 = par[7]*X; + py2 = par[8]*Y; + + z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y; + z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y; + + r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,fpar[0])); + f = par[5]*r + par[6]; + + q1 = par[5]*SQ(r)*(1 + z1); + q2 = par[5]*SQ(r)*fpar[0]*pow(z2,(fpar[0]-1)); + + if (dpar != NULL) { + dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y); + dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X); + dpar[2] = -2*q1*px1*X; + dpar[3] = -2*q1*py1*Y; + dpar[4] = -q1*X*Y; + dpar[5] = +r; + dpar[6] = +1; + dpar[7] = -2*q2*px2*X; + dpar[8] = -2*q2*py2*Y; + dpar[9] = -q2*X*Y; + } + return (f); +} + +int rgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("SXg", 2.35 * sqrt(2.0) / par[2]); + set_variable ("SYg", 2.35 * sqrt(2.0) / par[3]); + set_variable ("SXYg", par[4]); + set_variable ("Zpk", par[5]); + set_variable ("Sg", par[6]); + set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]); + set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]); + set_variable ("SXYf", par[9]); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss-pol.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss-pol.c (revision 42078) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss-pol.c (revision 42078) @@ -0,0 +1,114 @@ +# include "imfit.h" + +opihi_flt rgaussPolTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *); +void rgaussPolCL (); + +void rgauss_pol_setup (char *name) { + + if (strcmp(name, "rgauss-pol")) return; + + fitfunc = rgaussPolTD; + imfit_cleanup = rgaussPolCL; + Npar = 8; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, opihi_flt, MAX (Npar, 1)); + bzero (par, Npar*sizeof(opihi_flt)); + ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(opihi_flt)); + + // NOTE 1: SXg & SYg are FWHM, SXYg is the second moment calculated by starfuncs.c:get_aperture_stats + // NOTE 2: this function (unlike imfit-rgauss.c) uses polarization parameters for the elliptical profile + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[5] = get_variable_default ("Zpk", 10000); + par[6] = get_variable_default ("Sg", 0.0); // sky + par[7] = get_variable_default ("Sr", 2.0); + + opihi_flt Sxx = get_variable_default ("SXg", 2.0) * sqrt(2.0) / 2.35; // convert FWHM to sigma / sqrt(2) (Sxx) + opihi_flt Syy = get_variable_default ("SYg", 2.0) * sqrt(2.0) / 2.35; // convert FWHM to sigma / sqrt(2) (Syy) + opihi_flt Sxy = get_variable_default ("SXYg", 0); + + // z = (x^2 + y^2) / R^2 + (x^2 - y^2) / T^2 + x y / Q : NOTE Q is not squared to allow positive and negative values + + par[2] = 2.0/sqrt((1.0/SQ(Sxx) + 1.0/SQ(Syy))); // par[2] = R + par[3] = 2.0/sqrt((1.0/SQ(Sxx) - 1.0/SQ(Syy))); // par[3] = T + par[4] = 1.0 / Sxy; // par[4] = Q + + sky = &par[6]; +} + +void rgaussPolCL () { + opihi_flt Sxx = par[2]*par[3] / sqrt(SQ(par[3]) + SQ(par[2])); + opihi_flt Syy = par[2]*par[3] / sqrt(SQ(par[3]) - SQ(par[2])); + opihi_flt Sxy = 1.0 / par[4]; + + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("SXg", Sxx * 2.35 / sqrt(2.0)); + set_variable ("SYg", Syy * 2.35 / sqrt(2.0)); + set_variable ("SXYg", Sxy); + set_variable ("Zpk", par[5]); + set_variable ("Sg", par[6]); + set_variable ("Sr", par[7]); +} + +/* rgaussPol: (1 + z + z^alpha)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ +opihi_flt rgaussPolTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + opihi_flt X = x - par[0]; + opihi_flt Y = y - par[1]; + + opihi_flt P_R = (SQ(X) + SQ(Y)) / SQ(par[2]); + opihi_flt P_T = (SQ(X) - SQ(Y)) / SQ(par[3]); + opihi_flt P_Q = X * Y / par[4]; + + opihi_flt z = P_R + P_T + P_Q; + + opihi_flt p = pow(z,par[7] - 1.0); + opihi_flt r = 1.0 / (1 + z + z*p); + opihi_flt f = par[5]*r + par[6]; + + if (dpar != NULL) { + opihi_flt t = par[5]*SQ(r); + opihi_flt q = t*(1 + par[7]*p); + + dpar[0] = +q*(2*X/SQ(par[2]) + 2*X/SQ(par[3]) + Y/par[4]); + dpar[1] = +q*(2*Y/SQ(par[2]) - 2*Y/SQ(par[3]) + X/par[4]); + + dpar[2] = 2*q*P_R/par[2]; + dpar[3] = 2*q*P_T/par[3]; + dpar[4] = q*P_Q/par[4]; + + dpar[5] = +r; + dpar[6] = +1; + + // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0 + dpar[7] = (z == 0.0) ? 0.0 : -t*log(z)*p*z; + } + return (f); +} + +/****************************************************************************** + * this file defines the RGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelClass.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 as a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: + + 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) + *****************************************************************************/ + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 42078) @@ -9,7 +9,7 @@ fitfunc = rgaussTD; - cleanup = rgaussCL; - Npar = 10; - Nfpar = 1; + imfit_cleanup = rgaussCL; + Npar = 8; + Nfpar = 0; /* allocate free and fixed parameters */ @@ -21,67 +21,80 @@ par[0] = get_variable_default ("Xg", 0); par[1] = get_variable_default ("Yg", 0); - par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0); - par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0); - par[4] = 0.0; + par[2] = 2.35 / get_variable_default ("SXg", 2.0); + par[3] = 2.35 / get_variable_default ("SYg", 2.0); + par[4] = get_variable_default ("SXYg", 0); par[5] = get_variable_default ("Zpk", 10000); - par[6] = get_variable_default ("Sg", 0.0); - par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0); - par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0); - par[9] = get_variable_default ("SXYf", 0.0); + par[6] = get_variable_default ("Sg", 0.0); // sky + par[7] = get_variable_default ("Sr", 2.0); - fpar[0] = get_variable_default ("Npow", 2.25); +// fpar[0] = get_variable_default ("Npow", 2.25); sky = &par[6]; } -/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ +void rgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("SXg", 2.35 / par[2]); + set_variable ("SYg", 2.35 / par[3]); + set_variable ("SXYg", par[4]); + set_variable ("Zpk", par[5]); + set_variable ("Sg", par[6]); + set_variable ("Sr", par[7]); +} + +/* rgauss: (1 + z + z^alpha)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ opihi_flt rgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) { OHANA_UNUSED_PARAM(Npar); - opihi_flt X, Y, px1, py1, px2, py2; - opihi_flt z1, z2, r, q1, q2, f; + opihi_flt X = x - par[0]; + opihi_flt Y = y - par[1]; + + opihi_flt px = par[2]*X; + opihi_flt py = par[3]*Y; - X = x - par[0]; - Y = y - par[1]; - - px1 = par[2]*X; - py1 = par[3]*Y; - px2 = par[7]*X; - py2 = par[8]*Y; + opihi_flt z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y; - z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y; - z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y; - - r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,fpar[0])); - f = par[5]*r + par[6]; - - q1 = par[5]*SQ(r)*(1 + z1); - q2 = par[5]*SQ(r)*fpar[0]*pow(z2,(fpar[0]-1)); + opihi_flt p = pow(z,par[7] - 1.0); + opihi_flt r = 1.0 / (1 + z + z*p); + opihi_flt f = par[5]*r + par[6]; if (dpar != NULL) { - dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y); - dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X); - dpar[2] = -2*q1*px1*X; - dpar[3] = -2*q1*py1*Y; - dpar[4] = -q1*X*Y; + opihi_flt t = par[5]*SQ(r); + opihi_flt q = t*(1 + par[7]*p); + + dpar[0] = +q*(par[2]*px + par[4]*Y); + dpar[1] = +q*(par[3]*py + par[4]*X); + dpar[2] = -q*px*X; + dpar[3] = -q*py*Y; + dpar[4] = -q*X*Y; dpar[5] = +r; dpar[6] = +1; - dpar[7] = -2*q2*px2*X; - dpar[8] = -2*q2*py2*Y; - dpar[9] = -q2*X*Y; + + // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0 + dpar[7] = (z == 0.0) ? 0.0 : -t*log(z)*p*z; } return (f); } -int rgaussCL () { - set_variable ("Xg", par[0]); - set_variable ("Yg", par[1]); - set_variable ("SXg", 2.35 * sqrt(2.0) / par[2]); - set_variable ("SYg", 2.35 * sqrt(2.0) / par[3]); - set_variable ("SXYg", par[4]); - set_variable ("Zpk", par[5]); - set_variable ("Sg", par[6]); - set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]); - set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]); - set_variable ("SXYf", par[9]); -} +/****************************************************************************** + * this file defines the RGAUSS source shape model (XXX need a better name!). Note that these + * model functions are loaded by pmModelClass.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 as a PSF + * representations share a few parameters, for which # define names are listed in pmModel.h: + + 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) + *****************************************************************************/ + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c (revision 42078) @@ -4,11 +4,14 @@ void trailCL (); +// fitted parameters: # define PAR_X 0 # define PAR_Y 1 # define PAR_THETA 2 -# define PAR_SIGMA 3 -# define PAR_LENGTH 4 -# define PAR_I0 5 -# define PAR_SKY 6 +# define PAR_LENGTH 3 +# define PAR_I0 4 +# define PAR_SKY 5 + +// fixed parameters: +# define PAR_SIGMA 0 void trail_setup (char *name) { @@ -18,6 +21,6 @@ fitfunc = trailTD; imfit_cleanup = trailCL; - Npar = 7; - Nfpar = 0; + Npar = 6; + Nfpar = 1; /* allocate free and fixed parameters */ @@ -30,9 +33,10 @@ par[PAR_Y ] = get_variable_default ("Yg", 0.0); par[PAR_THETA ] = get_variable_default ("Tg", 0.0); - par[PAR_SIGMA ] = get_variable_default ("Wg", 2.0); par[PAR_LENGTH ] = get_variable_default ("Lg", 10.0); par[PAR_I0 ] = get_variable_default ("Zpk", 10000.0); par[PAR_SKY ] = get_variable_default ("Sg", 0.0); sky = &par[PAR_SKY]; + + fpar[PAR_SIGMA ] = get_variable_default ("Wg", 2.0); } @@ -40,9 +44,10 @@ set_variable ("Xg", par[PAR_X ]); set_variable ("Yg", par[PAR_Y ]); - set_variable ("Wg", par[PAR_SIGMA ]); set_variable ("Tg", par[PAR_THETA ]); set_variable ("Lg", par[PAR_LENGTH]); set_variable ("Zpk", par[PAR_I0 ]); set_variable ("Sg", par[PAR_SKY ]); + + set_variable ("Wg", fpar[PAR_SIGMA]); } @@ -54,5 +59,5 @@ opihi_flt Y = y - par[PAR_Y]; - opihi_flt S2 = 2.0 * SQ(par[PAR_SIGMA]); + opihi_flt S2 = 2.0 * SQ(fpar[PAR_SIGMA]); opihi_flt ST = sin(RAD_DEG*par[PAR_THETA]); @@ -78,11 +83,14 @@ // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + // since X = X - par[PAR_X], dFoo/dXo = -dFoo/dX float dRdX = +ST; float dRdY = -CT; - float dRdT = -Y*ST - X*CT; + float dRdT = (-Y*ST - X*CT)*RAD_DEG; + // note PAR_THETA is in degrees float dGdX = dGdR * dRdX; float dGdY = dGdR * dRdY; float dGdT = dGdR * dRdT; + // dGdL is 0.0 because dRdL is 0.0 (R is not a function of L) // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo @@ -96,5 +104,6 @@ float dZmdL = -0.5 / sqrt(S2); - float dZpdT = (-X*ST + Y*CT) / sqrt(S2); + // note PAR_THETA is in degrees + float dZpdT = (-X*ST + Y*CT) * RAD_DEG / sqrt(S2); float dZmdT = dZpdT; // dZpdT = dZmdT @@ -119,5 +128,4 @@ // dGdL is 0.0 because dRdL is 0.0 float dPdL = Gxy * (dEpdL - dEmdL); - float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT); @@ -127,5 +135,4 @@ dpar[PAR_LENGTH] = par[PAR_I0] * dPdL; dpar[PAR_THETA] = par[PAR_I0] * dPdT; - dpar[PAR_SIGMA] = 0; // we don't actually allow this to vary, so we do not need to calculate it } Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 42078) @@ -3,24 +3,29 @@ int imfit (int argc, char **argv) { - int i, j, N, Npts, Save, VERBOSE; - int sx, sy, nx, ny, Nx, Ny; - float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold; - opihi_flt *x, *y, *z, *dz; - float *V; + int N; Buffer *buf; - Save = FALSE; + char *Save = NULL; if ((N = get_argument (argc, argv, "-save"))) { remove_argument (N, &argc, argv); - Save = TRUE; - } - - // int ShapeVariation = FALSE; - // if ((N = get_argument (argc, argv, "-shapes"))) { - // remove_argument (N, &argc, argv); - // ShapeVariation = TRUE; - // } - - SatThreshold = 0xffff; + Save = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + int Insert = FALSE; + if ((N = get_argument (argc, argv, "-insert"))) { + remove_argument (N, &argc, argv); + Insert = TRUE; + if (Save) { gprint (GP_ERR, "-save and -insert are mutually exclusive\n"); free (Save); return (FALSE); } + } + + int minIter = 5; + if ((N = get_argument (argc, argv, "-min-iter"))) { + remove_argument (N, &argc, argv); + minIter = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + int SatThreshold = 0xffff; if ((N = get_argument (argc, argv, "-sat"))) { remove_argument (N, &argc, argv); @@ -30,5 +35,5 @@ /* Gain in e/DN */ - Gain = 1.0; + float Gain = 1.0; if ((N = get_argument (argc, argv, "-gain"))) { remove_argument (N, &argc, argv); @@ -38,5 +43,5 @@ /* RD noise in DN */ - RDnoise = 0.0; + float RDnoise = 0.0; if ((N = get_argument (argc, argv, "-rdnoise"))) { remove_argument (N, &argc, argv); @@ -45,5 +50,5 @@ } - VERBOSE = FALSE; + int VERBOSE = FALSE; if ((N = get_argument (argc, argv, "-v"))) { remove_argument (N, &argc, argv); @@ -51,21 +56,27 @@ } - /* set fitting function */ + /* set fitting function : defines par, Npar, fitfunc, etc globals (imfit.h) */ fgauss_setup ("fgauss"); if ((N = get_argument (argc, argv, "-func"))) { fitfunc = NULL; remove_argument (N, &argc, argv); - fgauss_setup (argv[N]); - pgauss_setup (argv[N]); - pgauss_psf_setup (argv[N]); + fgauss_setup (argv[N]); // OK + pgauss_setup (argv[N]); // OK sgauss_setup (argv[N]); - sgauss_psf_setup (argv[N]); - qgauss_setup (argv[N]); - qgauss_psf_setup (argv[N]); - qfgauss_setup (argv[N]); + qgauss_setup (argv[N]); // OK + rgauss_setup (argv[N]); + qfgauss_setup (argv[N]); // OK qrgauss_setup (argv[N]); trail_setup (argv[N]); + pgauss_psf_setup (argv[N]); + sgauss_psf_setup (argv[N]); + qgauss_psf_setup (argv[N]); + + fgauss_pol_setup (argv[N]); + rgauss_pol_setup (argv[N]); + if (fitfunc == NULL) { gprint (GP_ERR, "unknown function %s\n", argv[N]); + FREE (Save); return (FALSE); } @@ -74,39 +85,65 @@ if (argc != 6) { - gprint (GP_ERR, "USAGE: imfit sx sy nx ny\n"); + gprint (GP_ERR, "USAGE: imfit Xo Yo dX dY\n"); + gprint (GP_ERR, "options: [-save buffer] [-insert] [-sat value] [-gain value] [-rdnoise value] [-v] [-func option]\n"); + gprint (GP_ERR, " (Xo,Yo) : center\n"); + gprint (GP_ERR, " (dX,dY) : window size\n"); + FREE (Save); return (FALSE); } /* non-optional arguments */ - if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - sx = atof (argv[2]); - sy = atof (argv[3]); - nx = atof (argv[4]); - ny = atof (argv[5]); - Nx = buf[0].matrix.Naxis[0]; - Ny = buf[0].matrix.Naxis[1]; - - /* check if region is valid */ - if (sx + 0.5*nx < 0) goto range; - if (sy + 0.5*ny < 0) goto range; - if (sx + 0.5*nx >= Nx) goto range; - if (sy + 0.5*ny >= Ny) goto range; + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) { FREE (Save); return (FALSE); } + int Xo = atof (argv[2]); + int Yo = atof (argv[3]); + int dX = atof (argv[4]); + int dY = atof (argv[5]); + int Nx = buf[0].matrix.Naxis[0]; + int Ny = buf[0].matrix.Naxis[1]; + + int sx = Xo - dX/2; + int sy = Yo - dY/2; + + /* check if region is valid (center must be in range of image pixels) */ + if (Xo < 0) goto range; + if (Yo < 0) goto range; + if (Xo >= Nx) goto range; + if (Yo >= Ny) goto range; + + // image value in DN + // rdnoise in DN + // sigma_DN^2 = sigma_e^2 / gain^2 + // sigma_e^2 = Ne + // sigma_e^2 = DN * gain + // sigma_DN^2 = DN * gain / gain^2 = DN / gain + + if (Insert) { + float *Vi = (float *)buf[0].matrix.buffer; + for (int j = 0; j < dY; j++) { + for (int i = 0; i < dX; i++) { + float vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL); + Vi[(i+sx)+(j+sy)*Nx] += vf; + } + } + return TRUE; + } /* convert array z[x,y] to x[i], y[i], z[i] */ N = 0; - Npts = nx*ny; - ALLOCATE (x, opihi_flt, 2*Npts); - ALLOCATE (y, opihi_flt, 2*Npts); - ALLOCATE (z, opihi_flt, 2*Npts); - ALLOCATE (dz, opihi_flt, 2*Npts); - for (j = 0; j < ny; j++) { + int Npts = dX*dY; + ALLOCATE_PTR (x, opihi_flt, 2*Npts); + ALLOCATE_PTR (y, opihi_flt, 2*Npts); + ALLOCATE_PTR (z, opihi_flt, 2*Npts); + ALLOCATE_PTR (dz, opihi_flt, 2*Npts); + for (int j = 0; j < dY; j++) { if (j + sy < 0) continue; if (j + sy >= Ny) continue; - V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; - for (i = 0; i < nx; i++) { + float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; + for (int i = 0; i < dX; i++) { if (i + sx < 0) continue; if (i + sx >= Nx) continue; - if (*V > SatThreshold) goto next; - dz[N] = (SQ(RDnoise) + *V/Gain); + if (*V > SatThreshold) goto next; // skip pixels above threshold + if (!isfinite(*V)) goto next; // skip nan pixels + dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise if (dz[N] <= 0) goto next; dz[N] = 1.0 / dz[N]; @@ -122,17 +159,19 @@ /* run fit routine */ - ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE); - dchisq = ochisq; - chisq = ochisq; - -//for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) { - - for (i = 0; (i < 25); i++) { + float ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE); + float dchisq = ochisq; + float chisq = ochisq; + + int Niter = 0; + // for (int i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) { + + // keep iterating if chisq is increasing or + for (Niter = 0; (Niter < 25) && ((Niter < minIter) || ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar)))); Niter++) { chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE); dchisq = ochisq - chisq; + // fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq); ochisq = chisq; - fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq); } - set_int_variable ("Niter", i); + set_int_variable ("Niter", Niter); /** create output image (keep in sky) **/ @@ -141,24 +180,25 @@ float *Vi, *Vo, vr, vf; - if ((out = SelectBuffer ("out", ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((out = SelectBuffer (Save, ANYBUFFER, TRUE)) == NULL) { free (Save); return (FALSE); } free (out[0].header.buffer); free (out[0].matrix.buffer); strcpy (out[0].file, "(empty)"); - if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE; - - /* four panels: 1) raw image. 2) fit 3) raw - fit 4) ?? */ + if (!CreateBuffer (out, 2*dX, 2*dY, -32, 0.0, 1.0)) { free (Save); return FALSE; } + + /* four panels: 1) raw image. 2) fit 3) raw - fit 4) absolute deviation */ Vi = (float *)buf[0].matrix.buffer; Vo = (float *)out[0].matrix.buffer; - for (j = 0; j < ny; j++) { - for (i = 0; i < nx; i++) { + for (int j = 0; j < dY; j++) { + for (int i = 0; i < dX; i++) { vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL); vr = Vi[(i+sx)+(j+sy)*Nx]; - Vo[(i )+(j )*2*nx] = vr; - Vo[(i+nx)+(j )*2*nx] = vf; - Vo[(i )+(j+ny)*2*nx] = vr - vf + *sky; - Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky; + Vo[(i )+(j )*2*dX] = vr; + Vo[(i+dX)+(j )*2*dX] = vf; + Vo[(i )+(j+dY)*2*dX] = vr - vf + *sky; + Vo[(i+dX)+(j+dY)*2*dX] = fabs(vr-vf) + *sky; } } + free (Save); } @@ -169,5 +209,5 @@ if (VERBOSE) { - for (i = 0; i < Npar; i++) { + for (int i = 0; i < Npar; i++) { gprint (GP_ERR, "%g ", par[i]); } Index: /trunk/Ohana/src/opihi/cmd.astro/region.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/region.c (revision 42077) +++ /trunk/Ohana/src/opihi/cmd.astro/region.c (revision 42078) @@ -159,6 +159,7 @@ if (!strcasecmp (argv[CtypeArg], "GLS")) { strcpy (graphmode.coords.ctype, "DEC--GLS"); goto got_ctype; } if (!strcasecmp (argv[CtypeArg], "PAR")) { strcpy (graphmode.coords.ctype, "DEC--PAR"); goto got_ctype; } + if (!strcasecmp (argv[CtypeArg], "MOL")) { strcpy (graphmode.coords.ctype, "DEC--MOL"); goto got_ctype; } gprint (GP_ERR, "ERROR: invalid projection type %s\n", argv[CtypeArg]); - gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR\n"); + gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR, MOL\n"); return FALSE; } Index: /trunk/Ohana/src/opihi/cmd.astro/test/imfit.sh =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/test/imfit.sh (revision 42078) +++ /trunk/Ohana/src/opihi/cmd.astro/test/imfit.sh (revision 42078) @@ -0,0 +1,358 @@ + +macro trailfit.dither.XYLT + if ($0 != 5) + echo "USAGE: testfit.trail dX dY dL dT" + break + end + + $NTEST = 10000 + + # truth values for the fake object to be fitted + $Xg_tru = 250 + $Yg_tru = 250 + $Tg_tru = 45.0 + $Wg_tru = 2.0 + $Lg_tru = 30.0 + $Zpk_tru = 500.0 + $Sg_tru = 0.0 + + set.input.pars tru + + # generate a trail to be fitted + mcreate zbuf 500 500 + set tbuf = zbuf + 100 + + # insert an object with the parameters above + imfit -func trail tbuf $Xg $Yg 50 50 -insert + + # make a normalized noise image + mgaussdev sigma 500 500 0.0 1.0 + + # generate the observed image (truth + noise) + set obuf = tbuf + sigma*sqrt(tbuf) + + delete -q Xin Yin Lin Tin Xv Yv Tv Lv Zv Sv ChiSqV + + # run a bunch of tests, modifing the centroid guess by a random amount + for i 0 $NTEST + if ($i % 100 == 0) echo -no-return . + + # reset the guesses + set.input.pars tru + $Xg = $Xg_tru + $1*(rnd(0) - 0.5) + $Yg = $Yg_tru + $2*(rnd(0) - 0.5) + $Lg = $Lg_tru + max(3, $3*(rnd(0) - 0.5)); # too small L guess is silly + $Tg = $Tg_tru + $4*(rnd(0) - 0.5) + + concat $Xg Xin + concat $Yg Yin + concat $Lg Lin + concat $Tg Tin + + imfit -func trail obuf $Xg $Yg 50 50 + # imfit -func trail obuf $Xg $Yg 50 50 -save fit + + concat $Xg Xv + concat $Yg Yv + concat $Tg Tv + concat $Lg Lv + concat $Zpk Zv + concat $Sg Sv + concat $ChiSq ChiSqV + end +end + +macro trailfit.dither.LT + if ($0 != 3) + echo "USAGE: testfit.trail dL dT" + break + end + + $NTEST = 1000 + + # truth values for the fake object to be fitted + $Xg_tru = 250 + $Yg_tru = 250 + $Tg_tru = 45.0 + $Wg_tru = 2.0 + $Lg_tru = 30.0 + $Zpk_tru = 500.0 + $Sg_tru = 0.0 + + set.input.pars tru + + # generate a trail to be fitted + mcreate zbuf 500 500 + set tbuf = zbuf + 100 + + # insert an object with the parameters above + imfit -func trail tbuf $Xg $Yg 50 50 -insert + + # make a normalized noise image + mgaussdev sigma 500 500 0.0 1.0 + + # generate the observed image (truth + noise) + set obuf = tbuf + sigma*sqrt(tbuf) + + delete -q Xin Yin Xv Yv Tv Lv Zv Sv ChiSqV + + # run a bunch of tests, modifing the centroid guess by a random amount + for i 0 $NTEST + if ($i % 100 == 0) echo -no-return . + + # reset the guesses + set.input.pars tru + $Lg = $Lg_tru + $1*(rnd(0) - 0.5) + $Tg = $Tg_tru + $2*(rnd(0) - 0.5) + + concat $Lg Lin + concat $Tg Tin + + imfit -func trail obuf $Xg $Yg 50 50 + # imfit -func trail obuf $Xg $Yg 50 50 -save fit + + concat $Xg Xv + concat $Yg Yv + concat $Tg Tv + concat $Lg Lv + concat $Zpk Zv + concat $Sg Sv + concat $ChiSq ChiSqV + end +end + +macro trailfit.dither.pos + if ($0 != 3) + echo "USAGE: testfit.trail dX dY" + break + end + + $NTEST = 1000 + + # truth values for the fake object to be fitted + $Xg_tru = 250 + $Yg_tru = 250 + $Tg_tru = 0.0 + $Wg_tru = 2.0 + $Lg_tru = 10.0 + $Zpk_tru = 500.0 + $Sg_tru = 0.0 + + set.input.pars tru + + # generate a trail to be fitted + mcreate zbuf 500 500 + set tbuf = zbuf + 100 + + # insert an object with the parameters above + imfit -func trail tbuf $Xg $Yg 50 50 -insert + + # make a normalized noise image + mgaussdev sigma 500 500 0.0 1.0 + + # generate the observed image (truth + noise) + set obuf = tbuf + sigma*sqrt(tbuf) + + delete -q Xin Yin Xv Yv Tv Lv Zv Sv ChiSqV + + # run a bunch of tests, modifing the centroid guess by a random amount + for i 0 $NTEST + if ($i % 100 == 0) echo -no-return . + + # reset the guesses + set.input.pars tru + $Xg = $Xg_tru + $1*(rnd(0) - 0.5) + $Yg = $Yg_tru + $2*(rnd(0) - 0.5) + + concat $Xg Xin + concat $Yg Yin + + imfit -func trail obuf $Xg $Yg 50 50 + # imfit -func trail obuf $Xg $Yg 50 50 -save fit + + concat $Xg Xv + concat $Yg Yv + concat $Tg Tv + concat $Lg Lv + concat $Zpk Zv + concat $Sg Sv + concat $ChiSq ChiSqV + end +end + +macro testfit.trail + if ($0 != 3) + echo "USAGE: testfit.trail dX dY" + break + end + + $myfunc = trail + + mcreate zbuf 500 500 + set tbuf = zbuf + 100 + + # note the SXg, SYg are FWHM, not sigma values + $Xg_in = 250 + $Yg_in = 250 + $Tg_in = 0.0 + $Wg_in = 2.0 + $Lg_in = 10.0 + $Zpk_in = 500.0 + $Sg_in = 0.0 + + $Xg = $Xg_in + $Yg = $Yg_in + $Tg = $Tg_in + $Wg = $Wg_in + $Lg = $Lg_in + $Zpk = $Zpk_in + $Sg = $Sg_in + + # insert an object with the parameters above + imfit -func $myfunc tbuf $Xg $Yg 50 50 -insert + + # make a normalized noise image + mgaussdev sigma 500 500 0.0 1.0 + + # generate the observed image (truth + noise) + set obuf = tbuf + sigma*sqrt(tbuf) + + # make the guess wrong by a small amount + $Xg += $1 + $Yg += $2 + imfit -func $myfunc obuf $Xg $Yg 50 50 -save fit -v + + echo $Xg : $Xg_in {$Xg - $Xg_in } + echo $Yg : $Yg_in {$Yg - $Yg_in } + echo $Tg : $Tg_in {$Tg - $Tg_in } + echo $Wg : $Wg_in {$Wg - $Wg_in } + echo $Lg : $Lg_in {$Lg - $Lg_in } + echo $Zpk : $Zpk_in {$Zpk - $Zpk_in } + echo $Sg : $Sg_in {$Sg - $Sg_in } +end + +macro testfit.func + if ($0 != 2) + echo "USAGE: testfit.func (func)" + break + end + + $myfunc = $1 + + $NX = 101 + $NY = 101 + mcreate zbuf $NX $NY + set tbuf = zero(zbuf) + + # note the SXg, SYg are FWHM, not sigma values + $Xg_in = 50 + $Yg_in = 50 + $SXg_in = 4.0 + $SYg_in = 4.0 + $SXYg_in = 0.0 + $Zpk_in = 500.0 + $Sg_in = 0.0 + $Sr_in = 1.0 + $Npow_in = 2.25 + + $Xg = $Xg_in + $Yg = $Yg_in + $SXg = $SXg_in + $SYg = $SYg_in + $SXYg = $SXYg_in + $Zpk = $Zpk_in + $Sg = $Sg_in + $Sr = $Sr_in + $Npow = $Npow_in + + # insert an object with the parameters above + imfit -func $myfunc tbuf $Xg $Yg 50 50 -insert + + mgaussdev sigma $NX $NY 0.0 1.0 + + set noise = sigma*sqrt(tbuf) + set obuf = tbuf + noise + + imfit -func $myfunc obuf $Xg $Yg 50 50 -save fit -v + + echo $Xg : $Xg_in {$Xg - $Xg_in } + echo $Yg : $Yg_in {$Yg - $Yg_in } + echo $SXg : $SXg_in {$SXg - $SXg_in } + echo $SYg : $SYg_in {$SYg - $SYg_in } + echo $SXYg : $SXYg_in {$SXYg - $SXYg_in } + echo $Zpk : $Zpk_in {$Zpk - $Zpk_in } + echo $Sg : $Sg_in {$Sg - $Sg_in } + echo $Sr : $Sr_in {$Sr - $Sr_in } + echo $Npow : $Npow_in {$Npow - $Npow_in } +end + +macro testfit.fgauss.pol + if ($0 != 1) + echo "USAGE: testfit.fgauss.pol" + break + end + + $myfunc = fgauss + + mcreate zbuf 101 101 + set tbuf_o = zero(zbuf) + set tbuf_p = zero(zbuf) + + # note the SXg, SYg are FWHM, not sigma values + $Xg_in = 50 + $Yg_in = 50 + $SXg_in = 4.0 + $SYg_in = 4.0 + $SXYg_in = 0.0 + $Zpk_in = 500.0 + $Sg_in = 0.0 + $Sr_in = 1.0 + $Npow_in = 2.25 + + $Xg = $Xg_in + $Yg = $Yg_in + $SXg = $SXg_in + $SYg = $SYg_in + $SXYg = $SXYg_in + $Zpk = $Zpk_in + $Sg = $Sg_in + $Sr = $Sr_in + $Npow = $Npow_in + + # insert an object with the parameters above + imfit -func $myfunc tbuf_o $Xg $Yg 50 50 -insert + + # insert an object with the parameters using the -pol version + imfit -func $myfunc-pol tbuf_p $Xg $Yg 50 50 -insert + + set dt = tbuf_o - tbuf_p + stat dt + + tv -ch 1 tbuf_o 0 $Zpk_in + tv -ch 2 tbuf_p 0 $Zpk_in + tv -ch 3 dt {-0.1*$Zpk_in} {0.2*$Zpk_in} +end + +macro set.input.pars + if ($0 != 2) + echo "USAGE: set.input.pars (name)" + break + end + + # copy to the + $Xg = $Xg_$1 + $Yg = $Yg_$1 + $Tg = $Tg_$1 + $Wg = $Wg_$1 + $Lg = $Lg_$1 + $Zpk = $Zpk_$1 + $Sg = $Sg_$1 +end + + # echo $Xg : $Xg_in {$Xg - $Xg_in } + # echo $Yg : $Yg_in {$Yg - $Yg_in } + # echo $Tg : $Tg_in {$Tg - $Tg_in } + # echo $Wg : $Wg_in {$Wg - $Wg_in } + # echo $Lg : $Lg_in {$Lg - $Lg_in } + # echo $Zpk : $Zpk_in {$Zpk - $Zpk_in } + # echo $Sg : $Sg_in {$Sg - $Sg_in }