Changeset 2854

Timestamp:

Dec 29, 2004, 10:59:52 PM (22 years ago)

Author:

eugene

Message:

minor bug fixes

Location:

trunk/Ohana/src/opihi

Files:

• cmd.astro/imfit.c (modified) (7 diffs)
• cmd.data/ungridify.c (modified) (2 diffs)

trunk/Ohana/src/opihi/cmd.astro/imfit.c

@@ -2 +2 @@
 
 float fgaussTD (float, float, float *, int, float *);
+float pgaussTD (float, float, float *, int, float *);
+float vgaussTD (float, float, float *, int, float *);
+float qgaussTD (float, float, float *, int, float *);
 
 int imfit (int argc, char **argv) {
@@ -11 +14 @@
   char line[64];
   Buffer *buf;
+  float (*fitfunc)(float, float, float *, int, float *);
 
   Save = FALSE;
@@ -22 +26 @@
     remove_argument (N, &argc, argv);
     VERBOSE = TRUE;
+  }
+
+  Npar = 7;
+  fitfunc = fgaussTD;
+  if ((N = get_argument (argc, argv, "-pgauss"))) {
+    remove_argument (N, &argc, argv);
+    fitfunc = pgaussTD;
+  }
+  if ((N = get_argument (argc, argv, "-qgauss"))) {
+    remove_argument (N, &argc, argv);
+    fitfunc = qgaussTD;
+  }
+  if ((N = get_argument (argc, argv, "-vgauss"))) {
+    remove_argument (N, &argc, argv);
+    fitfunc = vgaussTD;
+    Npar = 9;
   }
 
@@ -42 +62 @@
   if (sy + ny >= Ny) goto range;
 
-  Npar = 7;
   ALLOCATE (par, float, 2*Npar);
   bzero (par, Npar*sizeof(float));
@@ -73 +92 @@
   par[5] = atof (get_variable ("Zpk"));
   par[6] = atof (get_variable ("Sg"));
+  if (Npar == 9) {
+    par[7] = par[8] = 1;
+  }
 
   /* run fit routine */
-  ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fgaussTD, VERBOSE);
+  ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
   dchisq = ochisq;
 
-  for (i = 0; (i < 20) && ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar))); i++) {
-    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fgaussTD, VERBOSE);
+  for (i = 0; (i < 20) && ((dchisq <= 0.0) || (dchisq > 0.001*(Npts - Npar))); i++) {
+    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     dchisq = ochisq - chisq;
     ochisq = chisq;
@@ -110 +132 @@
     for (j = 0; j < ny; j++) {
       for (i = 0; i < nx; i++) {
-        vf = fgaussTD ((float)(i+sx), (float)(j+sy), par, Npar, dpar);
+        vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, dpar);
         vr = Vi[(i+sx)+(j+sy)*Nx];
         Vo[(i   )+(j   )*2*nx] = vr;
@@ -172 +194 @@
   return (f);
 }
+
+/***** fix the derivatives.  q should not be in function */
+/* pars: 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) */
+  q = par[5]*r;
+  k = q*r*(1 + z + 0.5*z*z);
+  /* note difference from gaussian: q = par[5]*r */
+  f = q + 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;
+
+  return (f);
+}
+
+/* pars: 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;
+
+  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 + par[7]*0.5*z*z*(1 + par[8]*z/3)); /* ~ exp (-Z) */
+  q = par[5]*r*r*(1 + z + 0.5*z*z);
+  /* note difference from gaussian: q = par[5]*r */
+  f = q + 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] = +par[5]*r*z*z/2;
+  dpar[8] = +dpar[7]*z/3;
+
+  return (f);
+}
+
+/* pars: x, y, sx, sy, sxy, I, sky */
+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 + z + par[7]*pow(z,par[8])); /* ~ exp (-Z) */
+  q = par[5]*r*r*par[7]*par[8]*pow(z,par[8]-1.0);
+  /* note difference from gaussian: q = par[5]*r */
+  f = q + 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;
+
+  return (f);
+}
+
 /* alternate for exp(-z) */
 /* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */

trunk/Ohana/src/opihi/cmd.data/ungridify.c

@@ -10 +10 @@
   Vector *vx, *vy, *vz;
 
-  if (argc != 11) {
+  if (argc != 9) {
     fprintf (stderr, "USAGE: ungridify buffer Xmin Xmax Ymin Ymax x y z\n");
     return (FALSE);
   }
   
-  if ((bf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((bf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
   Xmin = atof (argv[2]);
   Xmax = atof (argv[3]);
@@ -26 +26 @@
   NY = bf[0].matrix.Naxis[1];
 
-  if ((vx = SelectVector (argv[6], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((vy = SelectVector (argv[7], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((vz = SelectVector (argv[8], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vx = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vy = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vz = SelectVector (argv[8], ANYVECTOR, TRUE)) == NULL) return (FALSE);
   REALLOCATE (vx[0].elements, float, Nx*Ny);
   REALLOCATE (vy[0].elements, float, Nx*Ny);