Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 3029) +++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c (revision 3030) @@ -17,4 +17,7 @@ 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; @@ -122,4 +125,12 @@ Npar = 4; Nfpar = 6; + } + if (!strcmp(argv[N], "sersic")) { + fitfunc = sersic; + Npar = 8; + } + if (!strcmp(argv[N], "serbulge")) { + fitfunc = serbulge; + Npar = 12; } if (fitfunc == NULL) { @@ -178,5 +189,5 @@ /* make initial guesses on parameters */ - if (Npar >= 7) { + if (Npar >= 7) { /* not Pgauss, Sgauss */ par[0] = get_variable_default ("Xg", 0); par[1] = get_variable_default ("Yg", 0); @@ -187,5 +198,5 @@ par[6] = get_variable_default ("Sg", 0.0); } - if (Npar == 4) { + if (Npar == 4) { /* Pgauss, Sgauss */ par[0] = get_variable_default ("Xg", 0); par[1] = get_variable_default ("Yg", 0); @@ -193,6 +204,20 @@ par[3] = get_variable_default ("Sg", 0.0); } - if (Npar == 9) { /** vgauss **/ + 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 **/ @@ -201,18 +226,18 @@ 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] = 0.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", 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", 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", 15.0); + fpar[5] = get_variable_default ("SXYf", 0.0); } @@ -339,7 +364,18 @@ set_variable ("Sg", par[6]); } - if (Npar == 9) { + 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) { @@ -371,4 +407,9 @@ set_variable ("ChiSq", chisq/(Npts - Npar)); + for (i = 0; i < Npar; i++) { + fprintf (stderr, "%g ", par[i]); + } + fprintf (stderr, "\n"); + free (x); free (y); @@ -411,4 +452,84 @@ 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); Index: /trunk/Ohana/src/opihi/mana/Makefile =================================================================== --- /trunk/Ohana/src/opihi/mana/Makefile (revision 3029) +++ /trunk/Ohana/src/opihi/mana/Makefile (revision 3030) @@ -32,4 +32,6 @@ manacmds = \ $(SDIR)/rawstars.$(ARCH).o \ +$(SDIR)/fitcontour.$(ARCH).o \ +$(SDIR)/starcontour.$(ARCH).o \ $(SDIR)/findpeaks.$(ARCH).o Index: /trunk/Ohana/src/opihi/mana/fitcontour.c =================================================================== --- /trunk/Ohana/src/opihi/mana/fitcontour.c (revision 3030) +++ /trunk/Ohana/src/opihi/mana/fitcontour.c (revision 3030) @@ -0,0 +1,102 @@ +# include "mana.h" +# define NTERM 3 + +int fitcontour (int argc, char **argv) { + + int i; + double **C, **B; + float cs1, sn1, cs, sn, x, y, r, xo, yo, t; + float dR, Rmin, Rmaj, Theta, Rx, Ry, Rxy, R1, R2, R3; + Vector *vecx, *vecy; + + /* USAGE fitcontour x y Xo Yo */ + if (argc < 5) goto usage; + + if ((vecx = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((vecy = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + xo = atof (argv[3]); + yo = atof (argv[4]); + + ALLOCATE (B, double *, NTERM); + ALLOCATE (C, double *, NTERM); + for (i = 0; i < NTERM; i++) { + ALLOCATE (C[i], double, NTERM); + bzero (C[i], NTERM*sizeof(double)); + ALLOCATE (B[i], double, 1); + bzero (B[i], sizeof(double)); + } + + /* we are fitting r = ro + rs*sin(2theta) + rc*cos(2theta) */ + /* sin(2t) = 2cos(t)sin(t) + cos(2t) = cos^2(t) - sin^2(t) + */ + for (i = 0; i < vecx[0].Nelements; i++) { + x = vecx[0].elements[i] - xo; + y = vecy[0].elements[i] - yo; + r = hypot (x, y); + + /* + t = 2*atan2(y,x); + sn = sin(t); + cs = cos(t); + */ + sn1 = y / r; + cs1 = x / r; + sn = 2*sn1*cs1; + cs = cs1*cs1 - sn1*sn1; + + C[0][0] += 1.0; + C[1][0] += sn; + C[1][1] += SQ(sn); + C[2][0] += cs; + C[2][1] += cs*sn; + C[2][2] += SQ(cs); + + B[0][0] += r; + B[1][0] += r*sn; + B[2][0] += r*cs; + } + C[0][1] = C[1][0]; + C[0][2] = C[2][0]; + C[1][2] = C[2][1]; + + gaussj (C, NTERM, B, 1); + + dR = hypot (B[1][0], B[2][0]); + Rmaj = B[0][0] + dR; + Rmin = B[0][0] - dR; + Theta = DEG_RAD*atan2 (B[1][0], B[2][0]) / 2; + + sn = B[1][0] / dR; + cs = B[2][0] / dR; + + R1 = SQ(Rmaj) + SQ(Rmin); + R2 = SQ(Rmaj) - SQ(Rmin); + R3 = Rmaj*Rmin; + + Rx = R3 / sqrt (R1 - R2*cs); + Ry = R3 / sqrt (R1 + R2*cs); + Rxy = -sn*R2 / SQ(R3); + + set_variable ("Rx", Rx); + set_variable ("Ry", Ry); + set_variable ("Rxy", Rxy); + + set_variable ("Rmin", Rmin); + set_variable ("Rmaj", Rmaj); + set_variable ("Theta", Theta); + + for (i = 0; i < NTERM; i++) { + free (B[i]); + free (C[i]); + } + free (B); + free (C); + + return (TRUE); + + usage: + fprintf (stderr, "fitcontour x y (xo) (yo)\n"); + return (FALSE); +} + Index: /trunk/Ohana/src/opihi/mana/init.c =================================================================== --- /trunk/Ohana/src/opihi/mana/init.c (revision 3029) +++ /trunk/Ohana/src/opihi/mana/init.c (revision 3030) @@ -2,8 +2,12 @@ int findpeaks PROTO((int, char **)); +int fitcontour PROTO((int, char **)); +int starcontour PROTO((int, char **)); int rawstars PROTO((int, char **)); static Command cmds[] = { {"findpeaks", findpeaks, "find image peaks"}, + {"fitcontour", fitcontour, "fit ellipse contour"}, + {"starcontour", starcontour, "object contour"}, {"rawstars", rawstars, "find raw star stats"}, }; Index: /trunk/Ohana/src/opihi/mana/starcontour.c =================================================================== --- /trunk/Ohana/src/opihi/mana/starcontour.c (revision 3030) +++ /trunk/Ohana/src/opihi/mana/starcontour.c (revision 3030) @@ -0,0 +1,129 @@ +# include "mana.h" + +int starcontour (int argc, char **argv) { + + int x, y, xs, xp, yp; + int Nx, Ny, N, Npts, min, max; + float *v; + float zt, zo, xmin, xmax; + Vector *vecx, *vecy; + Buffer *buf; + + if (argc < 5) goto usage; + + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + xp = atof (argv[2]); + yp = atof (argv[3]); + zo = atof (argv[4]); + + if ((vecx = SelectVector ("xo", ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((vecy = SelectVector ("yo", ANYVECTOR, TRUE)) == NULL) return (FALSE); + + N = 0; + Npts = 100; + REALLOCATE (vecx[0].elements, float, MAX (Npts, 1)); + REALLOCATE (vecy[0].elements, float, MAX (Npts, 1)); + + Nx = buf[0].matrix.Naxis[0]; + Ny = buf[0].matrix.Naxis[1]; + v = (float *)buf[0].matrix.buffer; + + /* find transition below (limit range?) */ + xs = xp; + for (y = yp; (y >= 0) && (v[xs + y*Nx] > zo); y--) { + /* find transition below (limit range?) */ + min = max = FALSE; + for (x = xs; (x >= 0) && !min; x--) { + min = FALSE; + zt = v[x + y*Nx]; + if (zt < zo) { + min = TRUE; + xmin = x + (zo - zt)/(v[x + 1 + y*Nx] - zt); + vecx[0].elements[N] = xmin; + vecy[0].elements[N] = y; + N ++; + if (N >= Npts) { + Npts += 100; + REALLOCATE (vecx[0].elements, float, MAX (Npts, 1)); + REALLOCATE (vecy[0].elements, float, MAX (Npts, 1)); + } + } + /* ignore edge cases? */ + } + /* find transition above (limit range?) */ + for (x = xs; (x < Nx) && !max; x++) { + max = FALSE; + zt = v[x + y*Nx]; + if (zt < zo) { + max = TRUE; + xmax = x - (zo - zt)/(v[x - 1 + y*Nx] - zt); + vecx[0].elements[N] = xmax; + vecy[0].elements[N] = y; + N ++; + if (N >= Npts) { + Npts += 100; + REALLOCATE (vecx[0].elements, float, MAX (Npts, 1)); + REALLOCATE (vecy[0].elements, float, MAX (Npts, 1)); + } + } + /* ignore edge cases? */ + } + if (min && max) { + xs = 0.5*(xmin + xmax); + } + } + + /* find transition above (limit range?) */ + xs = xp; + for (y = yp; (y < Ny) && (v[xs + y*Nx] > zo); y++) { + /* find transition below (limit range?) */ + min = max = FALSE; + for (x = xs; (x >= 0) && !min; x--) { + min = FALSE; + zt = v[x + y*Nx]; + if (zt < zo) { + min = TRUE; + xmin = x + (zo - zt)/(v[x + 1 + y*Nx] - zt); + vecx[0].elements[N] = xmin; + vecy[0].elements[N] = y; + N ++; + if (N >= Npts) { + Npts += 100; + REALLOCATE (vecx[0].elements, float, MAX (Npts, 1)); + REALLOCATE (vecy[0].elements, float, MAX (Npts, 1)); + } + } + /* ignore edge cases? */ + } + /* find transition above (limit range?) */ + for (x = xs; (x < Nx) && !max; x++) { + max = FALSE; + zt = v[x + y*Nx]; + if (zt < zo) { + max = TRUE; + xmax = x - (zo - zt)/(v[x - 1 + y*Nx] - zt); + vecx[0].elements[N] = xmax; + vecy[0].elements[N] = y; + N ++; + if (N >= Npts) { + Npts += 100; + REALLOCATE (vecx[0].elements, float, MAX (Npts, 1)); + REALLOCATE (vecy[0].elements, float, MAX (Npts, 1)); + } + } + /* ignore edge cases? */ + } + if (min && max) { + xs = 0.5*(xmin + xmax); + } + } + vecx[0].Nelements = N; + vecy[0].Nelements = N; + + return (TRUE); + + usage: + fprintf (stderr, "starcontour (buffer) (xpeak) (ypeak) (level)\n"); + return (FALSE); +} +