Index: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/Makefile (revision 41964) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/Makefile (revision 41965) @@ -94,4 +94,5 @@ $(SRC)/imfit-qfgauss.$(ARCH).o \ $(SRC)/imfit-qrgauss.$(ARCH).o \ +$(SRC)/imfit-rgauss.$(ARCH).o \ $(SRC)/imfit-trail.$(ARCH).o Index: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c (revision 41965) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-q2gauss.c (revision 41965) @@ -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: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-qgauss.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 41964) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 41965) @@ -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: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c (revision 41965) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-r2gauss.c (revision 41965) @@ -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: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-rgauss.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 41964) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 41965) @@ -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: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit.c (revision 41964) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit.c (revision 41965) @@ -65,4 +65,5 @@ sgauss_setup (argv[N]); qgauss_setup (argv[N]); // OK + rgauss_setup (argv[N]); qfgauss_setup (argv[N]); // OK qrgauss_setup (argv[N]); Index: /branches/eam_branches/ipp-20211108/Ohana/src/opihi/include/imfit.h =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/opihi/include/imfit.h (revision 41964) +++ /branches/eam_branches/ipp-20211108/Ohana/src/opihi/include/imfit.h (revision 41965) @@ -14,4 +14,5 @@ void sgauss_setup (char *name); void qgauss_setup (char *name); +void rgauss_setup (char *name); void qfgauss_setup (char *name); void qrgauss_setup (char *name);