Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 3143) +++ /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 3144) @@ -35,5 +35,4 @@ $(SDIR)/gauss.$(ARCH).o \ $(SDIR)/getvel.$(ARCH).o \ -$(SDIR)/imfit.$(ARCH).o \ $(SDIR)/medianmap.$(ARCH).o \ $(SDIR)/mkgauss.$(ARCH).o \ @@ -53,7 +52,20 @@ $(SDIR)/transform.$(ARCH).o +imfitstuff = \ +$(SDIR)/imsub.$(ARCH).o \ +$(SDIR)/imfit.$(ARCH).o \ +$(SDIR)/imfit-fgauss.$(ARCH).o \ +$(SDIR)/imfit-pgauss.$(ARCH).o \ +$(SDIR)/imfit-Pgauss.$(ARCH).o \ +$(SDIR)/imfit-qgauss.$(ARCH).o \ +$(SDIR)/imfit-Qgauss.$(ARCH).o \ +$(SDIR)/imfit-sgauss.$(ARCH).o \ +$(SDIR)/imfit-Sgauss.$(ARCH).o \ +$(SDIR)/imfit-qfgauss.$(ARCH).o \ +$(SDIR)/imfit-qrgauss.$(ARCH).o + libastrocmd: $(DESTLIB)/libastrocmd.a -$(LIB)/libastrocmd.$(ARCH).a: $(astrocmd) +$(LIB)/libastrocmd.$(ARCH).a: $(astrocmd) $(imfitstuff) # utilities ################################################# Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-Pgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-Pgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-Pgauss.c (revision 3144) @@ -0,0 +1,66 @@ +# include "imfit.h" + +float PgaussTD (float, float, float *, int, float *); +void PgaussCL (); + +void Pgauss_setup (char *name) { + + if (strcmp(name, "Pgauss")) return; + + fitfunc = PgaussTD; + imfit_cleanup = PgaussCL; + Npar = 4; + Nfpar = 3; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[2] = get_variable_default ("Zpk", 10000); + par[3] = get_variable_default ("Sg", 0.0); + sky = &par[3]; + + fpar[0] = 2.35 / get_variable_default ("SXg", 2.0); + fpar[1] = 2.35 / get_variable_default ("SYg", 2.0); + fpar[2] = get_variable_default ("SXYg", 0); +} + +void PgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("Zpk", par[2]); + set_variable ("Sg", par[3]); +} + +/* pseudo 2D gaussian -- x, y, (sx), (sy), (sxy), I, sky */ +float PgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, q, f; + + /* par -> fpar: (2,0), (3,1), (4,2) */ + + X = x - par[0]; + Y = y - par[1]; + + px = fpar[0]*X; + py = fpar[1]*Y; + + z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y; + r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */ + f = par[2]*r + par[3]; + q = par[2]*r*r*(1 + z + 0.5*z*z); + /* note difference from gaussian: q = par[5]*r */ + + if (dpar != NULL) { + dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y); + dpar[1] = q*(2*py*fpar[1] + fpar[2]*X); + dpar[2] = +r; + dpar[3] = +1; + } + return (f); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-Qgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-Qgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-Qgauss.c (revision 3144) @@ -0,0 +1,68 @@ +# include "imfit.h" + +float QgaussTD (float, float, float *, int, float *); +void QgaussCL (); + +void Qgauss_setup (char *name) { + + if (strcmp(name, "Qgauss")) return; + + fitfunc = QgaussTD; + imfit_cleanup = QgaussCL; + Npar = 4; + Nfpar = 5; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[2] = get_variable_default ("Zpk", 10000); + par[3] = get_variable_default ("Sg", 0.0); + + fpar[0] = 2.35 / get_variable_default ("SXg", 15.0); + fpar[1] = 2.35 / get_variable_default ("SYg", 15.0); + fpar[2] = get_variable_default ("SXYg", 0.0); + fpar[3] = get_variable_default ("Sr", 1.0); + fpar[4] = get_variable_default ("Npow", 2.25); + + sky = &par[3]; +} + +void QgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("Zpk", par[2]); + set_variable ("Sg", par[3]); +} + +/* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ +float QgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, q, f; + + X = x - par[0]; + Y = y - par[1]; + + px = fpar[0]*X; + py = fpar[1]*Y; + + z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y; + + r = 1.0 / (1 + fpar[3]*z + pow(z,fpar[4])); + f = par[2]*r + par[3]; + q = par[2]*SQ(r)*(fpar[3] + fpar[4]*pow(z,(fpar[4]-1))); + + if (dpar != NULL) { + dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y); + dpar[1] = q*(2*py*fpar[1] + fpar[2]*X); + dpar[2] = +r; + dpar[3] = +1; + } + return (f); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-Sgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-Sgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-Sgauss.c (revision 3144) @@ -0,0 +1,74 @@ +# include "imfit.h" + +float SgaussTD (float, float, float *, int, float *); +void SgaussCL (); + +void Sgauss_setup (char *name) { + + if (strcmp(name, "Sgauss")) return; + + fitfunc = SgaussTD; + imfit_cleanup = SgaussCL; + Npar = 4; + Nfpar = 7; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + par[0] = get_variable_default ("Xg", 0); + par[1] = get_variable_default ("Yg", 0); + par[2] = get_variable_default ("Zpk", 10000); + par[3] = get_variable_default ("Sg", 0.0); + + fpar[0] = 2.35 / get_variable_default ("SXg", 15.0); + fpar[1] = 2.35 / get_variable_default ("SYg", 15.0); + fpar[2] = get_variable_default ("SXYg", 0.0); + fpar[3] = 2.35 / get_variable_default ("SXf", 15.0); + fpar[4] = 2.35 / get_variable_default ("SYf", 15.0); + fpar[5] = get_variable_default ("SXYf", 0.0); + fpar[6] = get_variable_default ("Npow", 2.25); + + sky = &par[3]; +} + +void SgaussCL () { + set_variable ("Xg", par[0]); + set_variable ("Yg", par[1]); + set_variable ("Zpk", par[2]); + set_variable ("Sg", par[3]); +} + +/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ +float SgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, py1, px2, py2; + float z1, z2, r, q1, q2, f, f1, f2; + + X = x - par[0]; + Y = y - par[1]; + + px1 = fpar[0]*X; + py1 = fpar[1]*Y; + px2 = fpar[3]*X; + py2 = fpar[4]*Y; + + z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + fpar[2]*X*Y; + z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + fpar[5]*X*Y; + + r = 1.0 / (1 + z1 + pow(z2,fpar[6])); + f = par[2]*r + par[3]; + + q1 = par[2]*SQ(r); + q2 = par[2]*SQ(r)*fpar[6]*pow(z2,(fpar[6]-1)); + + if (dpar != NULL) { + dpar[0] = q1*(2*px1*fpar[0] + fpar[2]*Y) + q2*(2*px2*fpar[3] + fpar[5]*Y); + dpar[1] = q1*(2*py1*fpar[1] + fpar[2]*X) + q2*(2*py2*fpar[4] + fpar[5]*X); + dpar[2] = +r; + dpar[3] = +1; + } + return (f); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-fgauss.c (revision 3144) @@ -0,0 +1,68 @@ +# include "imfit.h" + +float fgaussTD (float, float, float *, int, float *); +void fgaussCL (); + +void fgauss_setup (char *name) { + + if (strcmp(name, "fgauss")) return; + + fitfunc = fgaussTD; + imfit_cleanup = fgaussCL; + Npar = 7; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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); + sky = &par[6]; +} + +void fgaussCL () { + 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]); +} + +/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */ +float fgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float 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 = exp (-z); + q = par[5]*r; + f = q + par[6]; + + 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; + dpar[3] = -2*q*py*Y; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + } + return (f); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-pgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-pgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-pgauss.c (revision 3144) @@ -0,0 +1,69 @@ +# include "imfit.h" + +float pgaussTD (float, float, float *, int, float *); +void pgaussCL (); + +void pgauss_setup (char *name) { + + if (strcmp(name, "pgauss")) return; + + fitfunc = pgaussTD; + imfit_cleanup = pgaussCL; + Npar = 7; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + par[0] = get_variable_default ("Xg", 0.0); + par[1] = get_variable_default ("Yg", 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.0); + par[5] = get_variable_default ("Zpk", 10000); + par[6] = get_variable_default ("Sg", 0.0); + sky = &par[6]; +} + +void pgaussCL () { + 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]); +} + +/* pseudo 2D gaussian -- x, y, sx, sy, sxy, I, sky */ +float pgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float 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 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */ + f = par[5]*r + par[6]; + q = par[5]*r*r*(1 + z + 0.5*z*z); + /* note difference from gaussian: q = par[5]*r */ + + 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; + dpar[3] = -2*q*py*Y; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + } + return (f); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-qfgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-qfgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-qfgauss.c (revision 3144) @@ -0,0 +1,74 @@ +# include "imfit.h" + +float qfgaussTD (float, float, float *, int, float *); +void qfgaussCL (); + +void qfgauss_setup (char *name) { + + if (strcmp(name, "qfgauss")) return; + + fitfunc = qfgaussTD; + imfit_cleanup = qfgaussCL; + Npar = 7; + Nfpar = 2; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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); + + fpar[0] = get_variable_default ("Npow", 2.25); + fpar[1] = get_variable_default ("Sr", 1.0); + + sky = &par[6]; +} + +void qfgaussCL () { + 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]); +} + +/* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky */ +float qfgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float 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 + fpar[1]*z + pow(z,fpar[0])); + f = par[5]*r + par[6]; + q = par[5]*SQ(r)*(fpar[1] + 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; + } + return (f); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-qgauss.c (revision 3144) @@ -0,0 +1,75 @@ +# include "imfit.h" + +float qgaussTD (float, float, float *, int, float *); +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, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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 */ +float qgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float 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-qrgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c (revision 3144) @@ -0,0 +1,78 @@ +# include "imfit.h" + +float qrgaussTD (float, float, float *, int, float *); +void qrgaussCL (); + +void qrgauss_setup (char *name) { + + if (strcmp(name, "qrgauss")) return; + + fitfunc = qrgaussTD; + imfit_cleanup = qrgaussCL; + Npar = 8; + Nfpar = 1; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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 ("Npow", 2.25); + + fpar[0] = get_variable_default ("Sr", 1.0); + + sky = &par[6]; +} + +void qrgaussCL () { + 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 ("Npow", par[7]); +} + +/* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */ +float qrgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, q, f; + + /* if (par[7] < 1.0) par[7] = 1.0; */ + + 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 + fpar[0]*z + pow(z,par[7]) + 0.05*z*z*z); + f = par[5]*r + par[6]; + q = par[5]*SQ(r)*(fpar[0] + par[7]*pow(z,(par[7]-1)) + 0.15*z*z); + + 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] = -5*par[5]*SQ(r)*log(z)*pow(z,par[7]); + } + return (f); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (revision 3144) @@ -0,0 +1,86 @@ +# include "imfit.h" + +float rgaussTD (float, float, float *, int, float *); +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, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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 */ +float rgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, py1, px2, py2; + float 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-serbulge.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-serbulge.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-serbulge.c (revision 3144) @@ -0,0 +1,102 @@ +# include "imfit.h" + +float serbulgeTD (float, float, float *, int, float *); +void serbulgeCL (); + +void serbulge_setup (char *name) { + + if (strcmp(name, "serbulge")) return; + + fitfunc = serbulgeTD; + cleanup = serbulgeCL; + Npar = 12; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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) / 2.0; + + 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[10] = get_variable_default ("Zpk", 10000) / 2.0; + + par[11] = get_variable_default ("Sr", 1.0); + + sky = &par[6]; +} + +/* 0 1 2 3 4 5 6 7 8 9 10 11 */ +/* sersic galaxy model w/ bulge: -- x, y, (sx, sy, sxy)_1, I_1, sky, (sx, sy, sxy)_2, I_2, n */ +/* exp (-b (r/r_e)^(1/n)) + pgauss (r) */ +float serbulgeTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, px2, py1, py2; + float z1, z2, r1, r2, t, 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; + + /* bulge component */ + r1 = 1.0 / (1 + z1 + 0.5*z1*z1*(1 + z1/3)); /* ~ exp (-Z) */ + + /* disk component */ + t = pow (z2, par[11]); + r2 = exp (-t); + + f = par[5]*r1 + par[10]*r2 + par[6]; + + q1 = par[5]*r1*r1*(1 + z1 + 0.5*z1*z1); + q2 = par[10]*r2*par[11]*pow(z2, par[11]-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] = +r1; + dpar[6] = +1; + dpar[7] = -2*q2*px2*X*50; + dpar[8] = -2*q2*py2*Y*50; + dpar[9] = -q2*X*Y*50; + dpar[10] = +r2*50; + dpar[11] = -q2*log(z2)*t*50; + } + return (f); +} + +int serbulgeCL () { + 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 ("Zb", 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]); + set_variable ("Zd", par[10]); + set_variable ("Sr", par[11]); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-sersic.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-sersic.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-sersic.c (revision 3144) @@ -0,0 +1,74 @@ +# include "imfit.h" + +float sersicTD (float, float, float *, int, float *); +void sersicCL (); + +void sersic_setup (char *name) { + + if (strcmp(name, "sersic")) return; + + fitfunc = sersicTD; + cleanup = sersicCL; + Npar = 8; + Nfpar = 0; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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] = get_variable_default ("Sr", 1.0); + + sky = &par[6]; +} + +/* sersic galaxy model -- x, y, sx, sy, sxy, I, sky, n */ +/* exp (-b (r/r_e)^(1/n)) */ +float sersicTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, t, 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; + t = pow (z, par[7]); + r = exp (-t); + f = par[5]*r + par[6]; + q = par[5]*r*par[7]*pow(z, par[7]-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; + dpar[3] = -2*q*py*Y; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + dpar[7] = -q*log(z)*t; + } + return (f); +} + +int sersicCL () { + 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 ("Sr", par[7]); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-sgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-sgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-sgauss.c (revision 3144) @@ -0,0 +1,98 @@ +# include "imfit.h" +# define FFACTOR 200 +# define FSCALE 1.2 + +float sgaussTD (float, float, float *, int, float *); +void sgaussCL (); + +void sgauss_setup (char *name) { + + if (strcmp(name, "sgauss")) return; + + fitfunc = sgaussTD; + imfit_cleanup = sgaussCL; + Npar = 10; + Nfpar = 2; + + /* allocate free and fixed parameters */ + ALLOCATE (par, float, MAX (Npar, 1)); + bzero (par, Npar*sizeof(float)); + ALLOCATE (fpar, float, MAX (Nfpar, 1)); + bzero (fpar, Nfpar*sizeof(float)); + + 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] = 2.35 / get_variable_default ("SXf", 15.0); + par[8] = 2.35 / get_variable_default ("SYf", 15.0); + par[9] = get_variable_default ("SXYf", 0.0); + + fpar[0] = get_variable_default ("Npow", 2.25); + fpar[1] = get_variable_default ("Npin", 1.00); // drop this? + + sky = &par[6]; +} + +void sgaussCL () { + 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 ("SXf", 2.35 / par[7]); + set_variable ("SYf", 2.35 / par[8]); + set_variable ("SXYf", par[9]); +} + +/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ +float sgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, py1, px2, py2; + float z1, z2, r, q1, q2, f, f1, f2; + + 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 + pow(z2,fpar[0])); + f = par[5]*r + par[6]; + + q1 = par[5]*SQ(r); + 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); + + /* these fudge factors impede the growth of par[2] beyond par[7] */ + f1 = fabs(par[7]) / fabs(par[2]); + f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1; + dpar[2] = -2*q1*px1*X*f2; + + f1 = fabs(par[8]) / fabs(par[3]); + f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1; + dpar[3] = -2*q1*py1*Y*f2; + + 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); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-test.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-test.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-test.c (revision 3144) @@ -0,0 +1,314 @@ + if (ShapeVariation) { + /* find dChi/dSx and dChi/dSy given by increasing shape terms by 5% */ + float tp1, tp2, chix, chiy; + chix = chiy = 0; + if (fitfunc == sgaussTD) { + tp1 = par[2]; + tp2 = par[7]; + par[2] = par[2]*1.05; + par[7] = par[7]*1.05; + chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[2] = tp1; + par[7] = tp2; + + tp1 = par[3]; + tp2 = par[8]; + par[3] = par[3]*1.05; + par[8] = par[8]*1.05; + chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[3] = tp1; + par[8] = tp2; + } + if (fitfunc == pgaussTD) { + tp1 = par[2]; + par[2] = par[2]*1.05; + chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[2] = tp1; + + tp1 = par[3]; + par[3] = par[3]*1.05; + chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[3] = tp1; + } + if (fitfunc == SgaussTD) { + tp1 = par[0]; + tp2 = par[3]; + par[0] = par[0]*1.05; + par[3] = par[3]*1.05; + chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[0] = tp1; + par[3] = tp2; + + tp1 = par[1]; + tp2 = par[4]; + par[1] = par[1]*1.05; + par[4] = par[4]*1.05; + chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[1] = tp1; + par[4] = tp2; + } + if (fitfunc == PgaussTD) { + tp1 = par[0]; + par[0] = par[0]*1.05; + chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[0] = tp1; + + tp1 = par[1]; + par[1] = par[1]*1.05; + chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; + par[1] = tp1; + } + set_variable ("dChiX", chix/chisq); + set_variable ("dChiY", chiy/chisq); + } + +# if (0) +/* pars: x, y, sx, sy, sxy, I, sky */ +float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, Z, E, F, q, R, f, p2, p3; + + X = x - par[0]; + Y = y - par[1]; + + p2 = X / par[2]; + p3 = Y / par[3]; + + Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y); /* R */ + E = 1.0 / (1 + Z); + + q = par[5] * E; + R = q*E; + F = 0.5 / Z; + + f = q + par[6]; + + dpar[0] = F*R*(p2 + par[4]*Y); + dpar[1] = F*R*(p3 + par[4]*X); + dpar[2] = F*0.5*R*p2*p2; + dpar[3] = F*0.5*R*p3*p3; + dpar[4] = -R*X*Y*F; + + dpar[5] = E; + dpar[6] = 1; + return (f); +} + +/* pars: x, y, sx, sy, sxy, I, sky */ +float fbarTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, Z, E, F, q, R, f, p2, p3; + + X = x - par[0]; + Y = y - par[1]; + + p2 = X / par[2]; + p3 = Y / par[3]; + + Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y; /* R */ + E = 1.0 / (1 + Z*Z*Z); + + q = par[5] * E; + F = 3*Z*Z; + R = q*E*F; + + f = q + par[6]; + + dpar[0] = R*(p2 + par[4]*Y); + dpar[1] = R*(p3 + par[4]*X); + dpar[2] = 0.5*R*p2*p2; + dpar[3] = 0.5*R*p3*p3; + dpar[4] = -R*X*Y; + + dpar[5] = E; + dpar[6] = 1; + return (f); +} + +/* convert from x,y to major,minor */ +void fix_ellipsegauss_pars (float *par, int Npar) { + + float p2, p4, angle, t1, t2, tmp, area; + + /* par[0], par[1] = Xo, Yo - stay the same */ + + p2 = 1/par[2]; + p4 = 1/par[3]; + + angle = 0.5 * atan2 (-2*par[4], p4 - p2); + + tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4])); + t1 = (p2 + p4 + tmp) / 2; + t2 = t1 - tmp; + + par[2] = 2.35482*sqrt(1/t2); + par[3] = 2.35482*sqrt(1/t1); + par[4] = DEG_RAD * angle; + + area = 2*M_PI/sqrt(t1*t2); + + par[5] *= area; + +} +# endif + +/*** options for later + + Subtract = FALSE; + if ((N = get_argument (argc, argv, "-sub"))) { + remove_argument (N, &argc, argv); + Subtract = TRUE; + } + + DFact = 1; + if ((N = get_argument (argc, argv, "-D"))) { + remove_argument (N, &argc, argv); + DFact = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + fitfunc = fgaussTD; + if ((N = get_argument (argc, argv, "-gal"))) { + remove_argument (N, &argc, argv); + fitfunc = fgalaxyTD; + } + if ((N = get_argument (argc, argv, "-bar"))) { + remove_argument (N, &argc, argv); + fitfunc = fbarTD; + } + + + f1 = 1; + if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c); + + f2 = 1; + if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c); + + if (Subtract) { + tmpsky = par[6]; + par[6] = 0; + for (N = j = 0; j < ny; j++) { + V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; + for (i = 0; i < nx; i++, V++, N++) { + dx = i + sx; + dy = j + sy; + *V -= fitfunc (dx, dy, par, Npar, (float *) NULL); + } + } + par[6] = tmpsky; + } + +***/ + +# if (0) + +/* these two tests were not very succcessful. the first did not model the shape well because + it could not match the change in roundness with radius. the second did not work because the + parameters were degenerate (amplitude and slope of second component) */ + +/* test: fixed, non-integer higher-order term -- x, y, sx, sy, sxy, I, sky, f1, f2 */ +float qgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, q, f, k; + + 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; + k = pow(z,1.75*par[8]); + r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */ + q = par[5]*r*r*(1 + 1.75*par[7]*par[8]*pow(z,1.75*par[8]-1)); + /* note difference from gaussian: q = par[5]*r */ + f = par[5]*r + par[6]; + + 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; + dpar[3] = -2*q*py*Y; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + dpar[7] = -10*par[5]*r*r*k; + dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z); + + return (f); +} + +/* test: two component model: inner pseudo gaussian with outer z^1.75 x, y, sx, sy, sxy, I, sky */ +float rgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, py1, px2, py2; + float z1, z2, r1, r2, q1, q2, f; + + X = x - par[0]; + Y = y - par[1]; + + px1 = par[2]*X; + py1 = par[3]*Y; + px2 = par[8]*X; + py2 = par[9]*Y; + + z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y; + z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[10]*X*Y; + + r1 = 1.0 / (1 + z1 + 0.5*SQ(z1)*(1 + z1/3)); /* ~ exp (-Z) */ + r2 = 1.0 / (1 + pow(z2,1.75)); + + f = par[5]*r1 + par[6] + par[7]*r2; + + q1 = par[5]*SQ(r1)*(1 + z1 + 0.5*SQ(z1)); + q2 = par[7]*SQ(r2)*(1.75*pow(z2,0.75)); + + dpar[ 0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[8] + par[10]*Y); + dpar[ 1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[9] + par[10]*X); + dpar[ 2] = -2*q1*px1*X; + dpar[ 3] = -2*q1*py1*Y; + dpar[ 4] = -q1*X*Y; + dpar[ 5] = +r1; + dpar[ 6] = +1; + dpar[ 7] = +r2*2; + dpar[ 8] = -2*q2*px2*X*2; + dpar[ 9] = -2*q2*py2*Y*2; + dpar[10] = -q2*X*Y; + + return (f); +} + +# endif + + /* forcing values to have a rational range + ALLOCATE (parmin, float, Npar); + ALLOCATE (parmax, float, Npar); + bzero (parmin, Npar*sizeof(float)); + bzero (parmax, Npar*sizeof(float)); + parmin[0] = parmin[1] = 0; + parmax[0] = buf[0].matrix.Naxis[0]; + parmax[1] = buf[0].matrix.Naxis[1]; + + parmin[2] = parmin[3] = 0.01; + parmax[2] = parmax[3] = 100.0; + parmin[4] = -1000; + parmax[4] = -1000; + + parmin[5] = 1; + parmax[5] = 1e5; + + parmin[6] = 0.0; + parmax[6] = 1e5; + + if (Npar == 9) { + parmin[7] = parmin[8] = 0.01; + parmax[7] = parmax[8] = 10.0; + } + if (Npar == 10) { + parmin[7] = parmin[8] = 0.01; + parmax[7] = parmax[8] = 10.0; + parmin[9] = -1000; + parmax[9] = -1000; + } + */ + Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-tgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-tgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-tgauss.c (revision 3144) @@ -0,0 +1,81 @@ +# include "imfit.h" + +float tgaussTD (float, float, float *, int, float *); +void tgaussCL (); + +void tgauss_setup (char *name) { + + if (strcmp(name, "tgauss")) return; + + fitfunc = tgaussTD; + cleanup = tgaussCL; + Npar = 10; + Nfpar = 2; + + 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); + fpar[1] = get_variable_default ("Npin", 1.00); // drop this? + + sky = &par[6]; +} + +/* two components: (1 + z_1^M + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ +float tgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px1, py1, px2, py2; + float 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 + pow(z1,fpar[1]) + pow(z2,fpar[0])); + f = par[5]*r + par[6]; + + q1 = par[5]*SQ(r)*fpar[1]*pow(z1,(fpar[1]-1)); + 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*2; + dpar[3] = -2*q1*py1*Y*2; + 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 tgaussCL () { + 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-vgauss.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit-vgauss.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit-vgauss.c (revision 3144) @@ -0,0 +1,70 @@ +# include "imfit.h" + +float vgaussTD (float, float, float *, int, float *); +void vgaussCL (); + +void vgauss_setup (char *name) { + + if (strcmp(name, "vgauss")) return; + + fitfunc = vgaussTD; + cleanup = vgaussCL; + Npar = 9; + Nfpar = 0; + + 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] = 1; + par[8] = 1; + sky = &par[6]; +} + +/* pseudo 2D gaussian with floating 2nd and 3rd order terms -- x, y, sx, sy, sxy, I, sky, f1, f2 */ +float vgaussTD (float x, float y, float *par, int Npar, float *dpar) { + + float X, Y, px, py; + float z, r, q, f, k; + + 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; + k = 0.5*z*z*(1 + par[8]*z/3); + r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */ + f = par[5]*r + par[6]; + q = par[5]*r*r*(1 + par[7]*z*(1 + par[8]*z/2)); + /* note difference from gaussian: q = par[5]*r */ + + 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; + dpar[3] = -2*q*py*Y; + dpar[4] = -q*X*Y; + dpar[5] = +r; + dpar[6] = +1; + dpar[7] = -100*par[5]*r*r*k; + dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6; + } + return (f); +} + +int vgaussCL () { + 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", par[7]); + set_variable ("SYf", par[8]); +} Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 3143) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 3144) @@ -1,55 +1,12 @@ -# include "astro.h" -# define FFACTOR 200 -# define FSCALE 1.2 - -float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts, - float *par, int Npar, - float (funcs)(float, float, float *, int, float *)); - -float fgaussTD (float, float, float *, int, float *); -float pgaussTD (float, float, float *, int, float *); -float sgaussTD (float, float, float *, int, float *); -float vgaussTD (float, float, float *, int, float *); - -float PgaussTD (float, float, float *, int, float *); -float SgaussTD (float, float, float *, int, float *); - -float tgaussTD (float, float, float *, int, float *); -float rgaussTD (float, float, float *, int, float *); - -float sersic (float, float, float *, int, float *); -float serbulge (float, float, float *, int, float *); - -static float Npow; -static float Npin; -static int MissingVariable = FALSE; -static float *fpar; -static int Nfpar; - -int mrq2dlimits (float *pmin, float *pmax, int Npar); - -float get_variable_default (char *name, float dvalue) { - - char *value; - float fvalue; - - value = get_variable (name); - if (value == NULL) { - MissingVariable = TRUE; - return (dvalue); - } - fvalue = atof (value); - return (fvalue); -} +# include "imfit.h" int imfit (int argc, char **argv) { - int i, j, N, Npts, Npar, Save, VERBOSE, ShapeVariation; + int i, j, N, Npts, Save, VERBOSE, ShapeVariation; int sx, sy, nx, ny, Nx, Ny; - float chisq, ochisq, dchisq, Gain, RDnoise; - float *x, *y, *z, *dz, *par, *V, *parmin, *parmax; + float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold; + float *x, *y, *z, *dz, *V, *parmin, *parmax; char line[64]; Buffer *buf; - float (*fitfunc)(float, float, float *, int, float *); Save = FALSE; @@ -63,4 +20,11 @@ remove_argument (N, &argc, argv); ShapeVariation = TRUE; + } + + SatThreshold = 0xffff; + if ((N = get_argument (argc, argv, "-sat"))) { + remove_argument (N, &argc, argv); + SatThreshold = atof (argv[N]); + remove_argument (N, &argc, argv); } @@ -88,50 +52,17 @@ /* set fitting function */ - Nfpar = 0; - Npar = 7; - fitfunc = fgaussTD; + fgauss_setup ("fgauss"); if ((N = get_argument (argc, argv, "-func"))) { fitfunc = NULL; remove_argument (N, &argc, argv); - if (!strcmp(argv[N], "fgauss")) { - fitfunc = fgaussTD; - } - if (!strcmp(argv[N], "pgauss")) { - fitfunc = pgaussTD; - } - if (!strcmp(argv[N], "vgauss")) { - fitfunc = vgaussTD; - Npar = 9; - } - if (!strcmp(argv[N], "sgauss")) { - fitfunc = sgaussTD; - Npar = 10; - } - if (!strcmp(argv[N], "rgauss")) { - fitfunc = rgaussTD; - Npar = 10; - } - if (!strcmp(argv[N], "tgauss")) { - fitfunc = tgaussTD; - Npar = 10; - } - if (!strcmp(argv[N], "Pgauss")) { - fitfunc = PgaussTD; - Npar = 4; - Nfpar = 3; - } - if (!strcmp(argv[N], "Sgauss")) { - fitfunc = SgaussTD; - Npar = 4; - Nfpar = 6; - } - if (!strcmp(argv[N], "sersic")) { - fitfunc = sersic; - Npar = 8; - } - if (!strcmp(argv[N], "serbulge")) { - fitfunc = serbulge; - Npar = 12; - } + fgauss_setup (argv[N]); + pgauss_setup (argv[N]); + Pgauss_setup (argv[N]); + sgauss_setup (argv[N]); + Sgauss_setup (argv[N]); + qgauss_setup (argv[N]); + Qgauss_setup (argv[N]); + qfgauss_setup (argv[N]); + qrgauss_setup (argv[N]); if (fitfunc == NULL) { fprintf (stderr, "unknown function %s\n", argv[N]); @@ -155,15 +86,9 @@ Ny = buf[0].matrix.Naxis[1]; - if (sx < 0) goto range; - if (sy < 0) goto range; - if (sx + nx >= Nx) goto range; - if (sy + ny >= Ny) goto range; - - ALLOCATE (par, float, Npar); - bzero (par, Npar*sizeof(float)); - if (Nfpar > 0) { - ALLOCATE (fpar, float, Nfpar); - bzero (fpar, Nfpar*sizeof(float)); - } + /* 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; /* convert array z[x,y] to x[i], y[i], z[i] */ @@ -175,8 +100,13 @@ ALLOCATE (dz, float, 2*Npts); for (j = 0; j < ny; 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++, V++) { + for (i = 0; i < nx; i++) { + if (i + sx < 0) continue; + if (i + sx >= Nx) continue; + if (*V > SatThreshold) goto next; dz[N] = (SQ(RDnoise) + *V/Gain); - if (dz[N] <= 0) continue; + if (dz[N] <= 0) goto next; dz[N] = 1.0 / dz[N]; x[N] = i + sx; @@ -184,66 +114,13 @@ z[N] = *V; N++; + next: + V++; } } Npts = N; - /* make initial guesses on parameters */ - if (Npar >= 7) { /* not Pgauss, Sgauss */ - 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); - } - if (Npar == 4) { /* Pgauss, Sgauss */ - par[0] = get_variable_default ("Xg", 0); - par[1] = get_variable_default ("Yg", 0); - par[2] = get_variable_default ("Zpk", 10000); - par[3] = get_variable_default ("Sg", 0.0); - } - if (fitfunc == vgaussTD) { - par[7] = par[8] = 1; - } - if (fitfunc == sersic) { - par[7] = get_variable_default ("Sr", 1.0); - } - if (fitfunc == serbulge) { - Npow = get_variable_default ("Npow", 2.25); - Npin = get_variable_default ("Npin", 1.00); - 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[11] = get_variable_default ("Sr", 1.0); - /* arbitrary guess: equal distribution between bulge and disk */ - par[10] = par[5]/2.0; - par[5] = par[10]; - } - if (Npar == 10) { /** sgauss, tgauss, rgauss **/ - Npow = get_variable_default ("Npow", 2.25); - Npin = get_variable_default ("Npin", 1.00); - 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); - } - if (fitfunc == PgaussTD) { - fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0); - fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0); - fpar[2] = get_variable_default ("SXYg", 0.0); - } - if (fitfunc == SgaussTD) { - fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0); - fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0); - fpar[2] = get_variable_default ("SXYg", 0.0); - fpar[3] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0); - fpar[4] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0); - fpar[5] = get_variable_default ("SXYf", 0.0); - } - /* run fit routine */ ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE); dchisq = ochisq; - for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) { chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE); @@ -252,74 +129,9 @@ } set_int_variable ("Niter", i); - - - if (ShapeVariation) { - /* find dChi/dSx and dChi/dSy given by increasing shape terms by 5% */ - float tp1, tp2, chix, chiy; - chix = chiy = 0; - if (fitfunc == sgaussTD) { - tp1 = par[2]; - tp2 = par[7]; - par[2] = par[2]*1.05; - par[7] = par[7]*1.05; - chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[2] = tp1; - par[7] = tp2; - - tp1 = par[3]; - tp2 = par[8]; - par[3] = par[3]*1.05; - par[8] = par[8]*1.05; - chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[3] = tp1; - par[8] = tp2; - } - if (fitfunc == pgaussTD) { - tp1 = par[2]; - par[2] = par[2]*1.05; - chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[2] = tp1; - - tp1 = par[3]; - par[3] = par[3]*1.05; - chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[3] = tp1; - } - if (fitfunc == SgaussTD) { - tp1 = par[0]; - tp2 = par[3]; - par[0] = par[0]*1.05; - par[3] = par[3]*1.05; - chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[0] = tp1; - par[3] = tp2; - - tp1 = par[1]; - tp2 = par[4]; - par[1] = par[1]*1.05; - par[4] = par[4]*1.05; - chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[1] = tp1; - par[4] = tp2; - } - if (fitfunc == PgaussTD) { - tp1 = par[0]; - par[0] = par[0]*1.05; - chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[0] = tp1; - - tp1 = par[1]; - par[1] = par[1]*1.05; - chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq; - par[1] = tp1; - } - set_variable ("dChiX", chix/chisq); - set_variable ("dChiY", chiy/chisq); - } /** create output image (keep in sky) **/ if (Save) { Buffer *out; - float *Vi, *Vo, vr, vf, sky; + float *Vi, *Vo, vr, vf; if ((out = SelectBuffer ("out", ANYBUFFER, TRUE)) == NULL) return (FALSE); @@ -335,9 +147,4 @@ out[0].header.Naxis[1] = 2*ny; CreateBuffer (out); - if (Npar == 4) { - sky = par[3]; - } else { - sky = par[6]; - } /* four panels: 1) raw image. 2) fit 3) raw - fit 4) ?? */ @@ -350,61 +157,13 @@ 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+ny)*2*nx] = vr - vf + *sky; + Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky; } } } - /* parameter renormalize */ - if (Npar >= 7) { - 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]); - } - if (fitfunc == vgaussTD) { - set_variable ("SXf", par[7]); - set_variable ("SYf", par[8]); - } - if (fitfunc == sersic) { - set_variable ("Sr", par[7]); - } - if (fitfunc == serbulge) { - set_variable ("Zb", par[5]); - 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]); - set_variable ("Zd", par[10]); - set_variable ("Sr", par[11]); - } - if (Npar == 10) { - 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]); - } - if (fitfunc == PgaussTD) { - set_variable ("Xg", par[0]); - set_variable ("Yg", par[1]); - set_variable ("Zpk", par[2]); - set_variable ("Sg", par[3]); - set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]); - set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]); - set_variable ("SXYg", fpar[2]); - } - if (fitfunc == SgaussTD) { - set_variable ("Xg", par[0]); - set_variable ("Yg", par[1]); - set_variable ("Zpk", par[2]); - set_variable ("Sg", par[3]); - set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]); - set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]); - set_variable ("SXYg", fpar[2]); - set_variable ("SXf", 2.35 * sqrt(2.0) / fpar[3]); - set_variable ("SYf", 2.35 * sqrt(2.0) / fpar[4]); - set_variable ("SXYf", fpar[5]); - } + /* save parameters to opihi variables */ + imfit_cleanup (); + set_variable ("ChiSq", chisq/(Npts - Npar)); @@ -421,5 +180,5 @@ free (dz); free (par); - if (Nfpar > 0) free (fpar); + free (fpar); mrq2dfree (Npar); @@ -431,611 +190,2 @@ } -/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */ -float fgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float 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 = exp (-z); - q = par[5]*r; - f = q + par[6]; - - 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; - dpar[3] = -2*q*py*Y; - dpar[4] = -q*X*Y; - dpar[5] = +r; - dpar[6] = +1; - } - return (f); -} - -/* sersic galaxy model -- x, y, sx, sy, sxy, I, sky, n */ -/* exp (-b (r/r_e)^(1/n)) */ -float sersic (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float z, r, t, 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; - t = pow (z, par[7]); - r = exp (-t); - f = par[5]*r + par[6]; - q = par[5]*r*par[7]*pow(z, par[7]-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; - dpar[3] = -2*q*py*Y; - dpar[4] = -q*X*Y; - dpar[5] = +r; - dpar[6] = +1; - dpar[7] = -q*log(z)*t; - } - return (f); -} - -/* 0 1 2 3 4 5 6 7 8 9 10 11 */ -/* sersic galaxy model w/ bulge: -- x, y, (sx, sy, sxy)_1, I_1, sky, (sx, sy, sxy)_2, I_2, n */ -/* exp (-b (r/r_e)^(1/n)) + pgauss (r) */ -float serbulge (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, px2, py1, py2; - float z1, z2, r1, r2, t, 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; - - /* bulge component */ - r1 = 1.0 / (1 + z1 + 0.5*z1*z1*(1 + z1/3)); /* ~ exp (-Z) */ - - /* disk component */ - t = pow (z2, par[11]); - r2 = exp (-t); - - f = par[5]*r1 + par[10]*r2 + par[6]; - - q1 = par[5]*r1*r1*(1 + z1 + 0.5*z1*z1); - q2 = par[10]*r2*par[11]*pow(z2, par[11]-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] = +r1; - dpar[6] = +1; - dpar[7] = -2*q2*px2*X*50; - dpar[8] = -2*q2*py2*Y*50; - dpar[9] = -q2*X*Y*50; - dpar[10] = +r2*50; - dpar[11] = -q2*log(z2)*t*50; - } - return (f); -} - -/* pseudo 2D gaussian -- x, y, sx, sy, sxy, I, sky */ -float pgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float 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 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */ - f = par[5]*r + par[6]; - q = par[5]*r*r*(1 + z + 0.5*z*z); - /* note difference from gaussian: q = par[5]*r */ - - 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; - dpar[3] = -2*q*py*Y; - dpar[4] = -q*X*Y; - dpar[5] = +r; - dpar[6] = +1; - } - return (f); -} - -/* pseudo 2D gaussian -- x, y, (sx), (sy), (sxy), I, sky */ -float PgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float z, r, q, f; - - /* par -> fpar: (2,0), (3,1), (4,2) */ - - X = x - par[0]; - Y = y - par[1]; - - px = fpar[0]*X; - py = fpar[1]*Y; - - z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y; - r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */ - f = par[2]*r + par[3]; - q = par[2]*r*r*(1 + z + 0.5*z*z); - /* note difference from gaussian: q = par[5]*r */ - - if (dpar != NULL) { - dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y); - dpar[1] = q*(2*py*fpar[1] + fpar[2]*X); - dpar[2] = +r; - dpar[3] = +1; - } - return (f); -} - -/* pseudo 2D gaussian with floating 2nd and 3rd order terms -- x, y, sx, sy, sxy, I, sky, f1, f2 */ -float vgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float z, r, q, f, k; - - 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; - k = 0.5*z*z*(1 + par[8]*z/3); - r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */ - f = par[5]*r + par[6]; - q = par[5]*r*r*(1 + par[7]*z*(1 + par[8]*z/2)); - /* note difference from gaussian: q = par[5]*r */ - - 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; - dpar[3] = -2*q*py*Y; - dpar[4] = -q*X*Y; - dpar[5] = +r; - dpar[6] = +1; - dpar[7] = -100*par[5]*r*r*k; - dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6; - } - return (f); -} - -/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ -float sgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, py1, px2, py2; - float z1, z2, r, q1, q2, f, f1, f2; - - 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 + pow(z2,Npow)); - f = par[5]*r + par[6]; - - q1 = par[5]*SQ(r); - q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-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); - - /* these fudge factors impede the growth of par[2] beyond par[7] */ - f1 = fabs(par[7]) / fabs(par[2]); - f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1; - dpar[2] = -2*q1*px1*X*f2; - - f1 = fabs(par[8]) / fabs(par[3]); - f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1; - dpar[3] = -2*q1*py1*Y*f2; - - 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); -} - -/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ -float SgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, py1, px2, py2; - float z1, z2, r, q1, q2, f, f1, f2; - - X = x - par[0]; - Y = y - par[1]; - - px1 = fpar[0]*X; - py1 = fpar[1]*Y; - px2 = fpar[3]*X; - py2 = fpar[4]*Y; - - z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + fpar[2]*X*Y; - z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + fpar[5]*X*Y; - - r = 1.0 / (1 + z1 + pow(z2,Npow)); - f = par[2]*r + par[3]; - - q1 = par[2]*SQ(r); - q2 = par[2]*SQ(r)*Npow*pow(z2,(Npow-1)); - - if (dpar != NULL) { - dpar[0] = q1*(2*px1*fpar[0] + fpar[2]*Y) + q2*(2*px2*fpar[3] + fpar[5]*Y); - dpar[1] = q1*(2*py1*fpar[1] + fpar[2]*X) + q2*(2*py2*fpar[4] + fpar[5]*X); - dpar[2] = +r; - dpar[3] = +1; - } - return (f); -} - -/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ -float rgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, py1, px2, py2; - float 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,Npow)); - f = par[5]*r + par[6]; - - q1 = par[5]*SQ(r)*(1 + z1); - q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-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); -} - -/* two components: (1 + z_1^M + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */ -float tgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, py1, px2, py2; - float 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 + pow(z1,Npin) + pow(z2,Npow)); - f = par[5]*r + par[6]; - - q1 = par[5]*SQ(r)*Npin*pow(z1,(Npin-1)); - q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-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*2; - dpar[3] = -2*q1*py1*Y*2; - 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); -} - -/* alternate for exp(-z) */ -/* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */ /* ~ exp (-Z) */ -/* E = 1.0 / (1 + Z + (0.5*Z*Z + 0.1666666*Z*Z*Z)); */ /* ~ exp (-Z) */ - -# if (0) -/* pars: x, y, sx, sy, sxy, I, sky */ -float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, Z, E, F, q, R, f, p2, p3; - - X = x - par[0]; - Y = y - par[1]; - - p2 = X / par[2]; - p3 = Y / par[3]; - - Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y); /* R */ - E = 1.0 / (1 + Z); - - q = par[5] * E; - R = q*E; - F = 0.5 / Z; - - f = q + par[6]; - - dpar[0] = F*R*(p2 + par[4]*Y); - dpar[1] = F*R*(p3 + par[4]*X); - dpar[2] = F*0.5*R*p2*p2; - dpar[3] = F*0.5*R*p3*p3; - dpar[4] = -R*X*Y*F; - - dpar[5] = E; - dpar[6] = 1; - return (f); -} - -/* pars: x, y, sx, sy, sxy, I, sky */ -float fbarTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, Z, E, F, q, R, f, p2, p3; - - X = x - par[0]; - Y = y - par[1]; - - p2 = X / par[2]; - p3 = Y / par[3]; - - Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y; /* R */ - E = 1.0 / (1 + Z*Z*Z); - - q = par[5] * E; - F = 3*Z*Z; - R = q*E*F; - - f = q + par[6]; - - dpar[0] = R*(p2 + par[4]*Y); - dpar[1] = R*(p3 + par[4]*X); - dpar[2] = 0.5*R*p2*p2; - dpar[3] = 0.5*R*p3*p3; - dpar[4] = -R*X*Y; - - dpar[5] = E; - dpar[6] = 1; - return (f); -} - -/* convert from x,y to major,minor */ -void fix_ellipsegauss_pars (float *par, int Npar) { - - float p2, p4, angle, t1, t2, tmp, area; - - /* par[0], par[1] = Xo, Yo - stay the same */ - - p2 = 1/par[2]; - p4 = 1/par[3]; - - angle = 0.5 * atan2 (-2*par[4], p4 - p2); - - tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4])); - t1 = (p2 + p4 + tmp) / 2; - t2 = t1 - tmp; - - par[2] = 2.35482*sqrt(1/t2); - par[3] = 2.35482*sqrt(1/t1); - par[4] = DEG_RAD * angle; - - area = 2*M_PI/sqrt(t1*t2); - - par[5] *= area; - -} -# endif - -/*** options for later - - Subtract = FALSE; - if ((N = get_argument (argc, argv, "-sub"))) { - remove_argument (N, &argc, argv); - Subtract = TRUE; - } - - DFact = 1; - if ((N = get_argument (argc, argv, "-D"))) { - remove_argument (N, &argc, argv); - DFact = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - fitfunc = fgaussTD; - if ((N = get_argument (argc, argv, "-gal"))) { - remove_argument (N, &argc, argv); - fitfunc = fgalaxyTD; - } - if ((N = get_argument (argc, argv, "-bar"))) { - remove_argument (N, &argc, argv); - fitfunc = fbarTD; - } - - - f1 = 1; - if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c); - - f2 = 1; - if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c); - - if (Subtract) { - tmpsky = par[6]; - par[6] = 0; - for (N = j = 0; j < ny; j++) { - V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; - for (i = 0; i < nx; i++, V++, N++) { - dx = i + sx; - dy = j + sy; - *V -= fitfunc (dx, dy, par, Npar, (float *) NULL); - } - } - par[6] = tmpsky; - } - -***/ - -# if (0) - -/* these two tests were not very succcessful. the first did not model the shape well because - it could not match the change in roundness with radius. the second did not work because the - parameters were degenerate (amplitude and slope of second component) */ - -/* test: fixed, non-integer higher-order term -- x, y, sx, sy, sxy, I, sky, f1, f2 */ -float qgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px, py; - float z, r, q, f, k; - - 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; - k = pow(z,1.75*par[8]); - r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */ - q = par[5]*r*r*(1 + 1.75*par[7]*par[8]*pow(z,1.75*par[8]-1)); - /* note difference from gaussian: q = par[5]*r */ - f = par[5]*r + par[6]; - - 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; - dpar[3] = -2*q*py*Y; - dpar[4] = -q*X*Y; - dpar[5] = +r; - dpar[6] = +1; - dpar[7] = -10*par[5]*r*r*k; - dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z); - - return (f); -} - -/* test: two component model: inner pseudo gaussian with outer z^1.75 x, y, sx, sy, sxy, I, sky */ -float rgaussTD (float x, float y, float *par, int Npar, float *dpar) { - - float X, Y, px1, py1, px2, py2; - float z1, z2, r1, r2, q1, q2, f; - - X = x - par[0]; - Y = y - par[1]; - - px1 = par[2]*X; - py1 = par[3]*Y; - px2 = par[8]*X; - py2 = par[9]*Y; - - z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y; - z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[10]*X*Y; - - r1 = 1.0 / (1 + z1 + 0.5*SQ(z1)*(1 + z1/3)); /* ~ exp (-Z) */ - r2 = 1.0 / (1 + pow(z2,1.75)); - - f = par[5]*r1 + par[6] + par[7]*r2; - - q1 = par[5]*SQ(r1)*(1 + z1 + 0.5*SQ(z1)); - q2 = par[7]*SQ(r2)*(1.75*pow(z2,0.75)); - - dpar[ 0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[8] + par[10]*Y); - dpar[ 1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[9] + par[10]*X); - dpar[ 2] = -2*q1*px1*X; - dpar[ 3] = -2*q1*py1*Y; - dpar[ 4] = -q1*X*Y; - dpar[ 5] = +r1; - dpar[ 6] = +1; - dpar[ 7] = +r2*2; - dpar[ 8] = -2*q2*px2*X*2; - dpar[ 9] = -2*q2*py2*Y*2; - dpar[10] = -q2*X*Y; - - return (f); -} - -# endif - - /* forcing values to have a rational range - ALLOCATE (parmin, float, Npar); - ALLOCATE (parmax, float, Npar); - bzero (parmin, Npar*sizeof(float)); - bzero (parmax, Npar*sizeof(float)); - parmin[0] = parmin[1] = 0; - parmax[0] = buf[0].matrix.Naxis[0]; - parmax[1] = buf[0].matrix.Naxis[1]; - - parmin[2] = parmin[3] = 0.01; - parmax[2] = parmax[3] = 100.0; - parmin[4] = -1000; - parmax[4] = -1000; - - parmin[5] = 1; - parmax[5] = 1e5; - - parmin[6] = 0.0; - parmax[6] = 1e5; - - if (Npar == 9) { - parmin[7] = parmin[8] = 0.01; - parmax[7] = parmax[8] = 10.0; - } - if (Npar == 10) { - parmin[7] = parmin[8] = 0.01; - parmax[7] = parmax[8] = 10.0; - parmin[9] = -1000; - parmax[9] = -1000; - } - */ - Index: /trunk/Ohana/src/opihi/cmd.astro/imsub.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imsub.c (revision 3144) +++ /trunk/Ohana/src/opihi/cmd.astro/imsub.c (revision 3144) @@ -0,0 +1,79 @@ +# include "imfit.h" + +int imsub (int argc, char **argv) { + + int i, j, N, VERBOSE; + int sx, sy, nx, ny, Nx, Ny; + float value; + float *V; + Buffer *buf; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + + /* set fitting function */ + 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_setup (argv[N]); + sgauss_setup (argv[N]); + Sgauss_setup (argv[N]); + qgauss_setup (argv[N]); + Qgauss_setup (argv[N]); + qfgauss_setup (argv[N]); + qrgauss_setup (argv[N]); + if (fitfunc == NULL) { + fprintf (stderr, "unknown function %s\n", argv[N]); + return (FALSE); + } + remove_argument (N, &argc, argv); + } + + if (argc != 6) { + fprintf (stderr, "USAGE: imfit sx sy nx ny\n"); + 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 >= Nx) goto range; + if (sy >= Ny) goto range; + if (sx + nx < 0) goto range; + if (sy + ny < 0) goto range; + + /* subtract model fit, but not local sky */ + for (j = 0; j < ny; 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++, V++) { + if (i + sx < 0) continue; + if (i + sx >= Nx) continue; + value = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL); + *V -= value; + } + } + + free (par); + free (fpar); + return (TRUE); + +range: + fprintf (stderr, "region out of range\n"); + return (FALSE); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/init.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 3143) +++ /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 3144) @@ -15,4 +15,5 @@ int getvel PROTO((int, char **)); int imfit PROTO((int, char **)); +int imsub PROTO((int, char **)); int medianmap PROTO((int, char **)); int mkgauss PROTO((int, char **)); @@ -46,4 +47,5 @@ {"getvel", getvel, "rotcurve to velocities"}, {"imfit", imfit, "fit function"}, + {"imsub", imsub, "subtract function"}, {"medianmap", medianmap, "small median image"}, {"mkgauss", mkgauss, "generate a 2-D gaussian centered in image"}, Index: /trunk/Ohana/src/opihi/include/imfit.h =================================================================== --- /trunk/Ohana/src/opihi/include/imfit.h (revision 3144) +++ /trunk/Ohana/src/opihi/include/imfit.h (revision 3144) @@ -0,0 +1,10 @@ +# include "astro.h" + +int Npar; +int Nfpar; +float *par; +float *fpar; +float *sky; + +float (*fitfunc)(float, float, float *, int, float *); +void (*imfit_cleanup)();