IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2876


Ignore:
Timestamp:
Jan 3, 2005, 4:03:21 PM (22 years ago)
Author:
eugene
Message:

minor bug fixes, imfit tests

Location:
trunk/Ohana/src/opihi
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r2854 r2876  
    55float vgaussTD (float, float, float *, int, float *);
    66float qgaussTD (float, float, float *, int, float *);
     7float rgaussTD (float, float, float *, int, float *);
     8float sgaussTD (float, float, float *, int, float *);
     9
     10static float Npow;
    711
    812int imfit (int argc, char **argv) {
    … …  
    1014  int i, j, N, Npts, Npar, Save, VERBOSE;
    1115  int sx, sy, nx, ny, Nx, Ny;
    12   float chisq, ochisq, dchisq;
     16  float chisq, ochisq, dchisq, Gain, RDnoise;
    1317  float *x, *y, *z, *dz, *par, *V;
    1418  char line[64];
    … …  
    2226  }
    2327
     28  /* Gain in e/DN */
     29  Gain = 1.0;
     30  if ((N = get_argument (argc, argv, "-gain"))) {
     31    remove_argument (N, &argc, argv);
     32    Gain = atof(argv[N]);
     33    remove_argument (N, &argc, argv);
     34  }
     35
     36  /* RD noise in DN */
     37  RDnoise = 0.0;
     38  if ((N = get_argument (argc, argv, "-rdnoise"))) {
     39    remove_argument (N, &argc, argv);
     40    RDnoise = atof(argv[N]);
     41    remove_argument (N, &argc, argv);
     42  }
     43
    2444  VERBOSE = FALSE;
    2545  if ((N = get_argument (argc, argv, "-v"))) {
    … …  
    2848  }
    2949
     50  /* set fitting function */
    3051  Npar = 7;
    3152  fitfunc = fgaussTD;
    32   if ((N = get_argument (argc, argv, "-pgauss"))) {
    33     remove_argument (N, &argc, argv);
    34     fitfunc = pgaussTD;
    35   }
    36   if ((N = get_argument (argc, argv, "-qgauss"))) {
    37     remove_argument (N, &argc, argv);
    38     fitfunc = qgaussTD;
    39   }
    40   if ((N = get_argument (argc, argv, "-vgauss"))) {
    41     remove_argument (N, &argc, argv);
    42     fitfunc = vgaussTD;
    43     Npar = 9;
     53  if ((N = get_argument (argc, argv, "-func"))) {
     54    fitfunc = NULL;
     55    remove_argument (N, &argc, argv);
     56    if (!strcmp(argv[N], "fgauss")) {
     57      fitfunc = fgaussTD;
     58    }
     59    if (!strcmp(argv[N], "pgauss")) {
     60      fitfunc = pgaussTD;
     61    }
     62    if (!strcmp(argv[N], "qgauss")) {
     63      fitfunc = qgaussTD;
     64      Npar = 9;
     65    }
     66    if (!strcmp(argv[N], "vgauss")) {
     67      fitfunc = vgaussTD;
     68      Npar = 9;
     69    }
     70    if (!strcmp(argv[N], "rgauss")) {
     71      fitfunc = rgaussTD;
     72      Npar = 11;
     73    }
     74    if (!strcmp(argv[N], "sgauss")) {
     75      fitfunc = sgaussTD;
     76      Npar = 10;
     77    }
     78    if (fitfunc == NULL) {
     79      fprintf (stderr, "unknown function %s\n", argv[N]);
     80      return (FALSE);
     81    }
     82    remove_argument (N, &argc, argv);
    4483  }
    4584
    … …  
    4988  }
    5089
     90  /* non-optional arguments */
    5191  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    5292  sx = atof (argv[2]);
    … …  
    66106
    67107  /* convert array z[x,y] to x[i], y[i], z[i] */
     108  N = 0;
    68109  Npts = nx*ny;
    69110  ALLOCATE (x,  float, 2*Npts);
    … …  
    71112  ALLOCATE (z,  float, 2*Npts);
    72113  ALLOCATE (dz, float, 2*Npts);
    73 
    74   N = 0;
    75114  for (j = 0; j < ny; j++) {
    76115    V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    77     for (i = 0; i < nx; i++, V++, N++) {
     116    for (i = 0; i < nx; i++, V++) {
     117      dz[N] = (SQ(RDnoise) + *V/Gain);
     118      if (dz[N] <= 0) continue;
     119      dz[N] = 1.0 / dz[N];
    78120      x[N] = i + sx;
    79121      y[N] = j + sy;
    80122      z[N] = *V;
    81       dz[N] = 1.0 / *V;  /* temporarily no gain (or, gain = 1) */
    82     }
    83   }
     123      N++;
     124    }
     125  }
     126  Npts = N;
    84127
    85128  /* make initial guesses on parameters */
    86   /* fitting X^2 / (2*Sx), Y^2 / (2*Sy) */
    87129  par[0] = atof (get_variable ("Xg"));
    88130  par[1] = atof (get_variable ("Yg"));
    … …  
    95137    par[7] = par[8] = 1;
    96138  }
     139  if (Npar == 11) {
     140    par[7] = 0.2*par[5];
     141    par[8] = 0.25*par[2];
     142    par[9] = 0.25*par[3];
     143    par[10] = 0.0;
     144  }
     145  if (Npar == 10) {
     146    Npow = atof (get_variable ("Npow"));
     147    par[7] = 2.35 * sqrt(2.0) / atof (get_variable ("SXf"));
     148    par[8] = 2.35 * sqrt(2.0) / atof (get_variable ("SYf"));
     149    par[9] = 0.0;
     150  }
     151
    97152
    98153  /* run fit routine */
    … …  
    100155  dchisq = ochisq;
    101156
    102   for (i = 0; (i < 20) && ((dchisq <= 0.0) || (dchisq > 0.001*(Npts - Npar))); i++) {
     157  for (i = 0; (i < 50) && ((dchisq <= 0.0) || (dchisq > 0.001*(Npts - Npar))); i++) {
    103158    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    104159    dchisq = ochisq - chisq;
    … …  
    107162  } 
    108163
    109   /** create output image **/
     164  /** create output image (keep in sky : par[6]?) **/
    110165  if (Save) {
    111166    Buffer *out;
    … …  
    137192        Vo[(i+nx)+(j   )*2*nx] = vf;
    138193        Vo[(i   )+(j+ny)*2*nx] = vr - vf;
    139         Vo[(i+nx)+(j+ny)*2*nx] = 0;
     194        Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf);
    140195      }
    141196    }
    … …  
    146201  par[2] = 2.35*sqrt(2.0) / par[2];
    147202  par[3] = 2.35*sqrt(2.0) / par[3];
     203  if (Npar == 11) {
     204    par[8] = 2.35*sqrt(2.0) / par[8];
     205    par[9] = 2.35*sqrt(2.0) / par[9];
     206  }
     207  if (Npar == 10
     208) {
     209    par[7] = 2.35*sqrt(2.0) / par[7];
     210    par[8] = 2.35*sqrt(2.0) / par[8];
     211  }
    148212  for (i = 0; i < Npar; i++) {
    149213    sprintf (line, "PAR%d", i);
    … …  
    152216  }
    153217  if (VERBOSE) fprintf (stderr, "  %f\n", chisq/Npts);
     218  set_variable ("ChiSq", (double) (chisq/Npts));
    154219
    155220  free (par);
    … …  
    210275  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    211276  r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
    212   q = par[5]*r;
    213   k = q*r*(1 + z + 0.5*z*z);
     277  f = par[5]*r + par[6];
     278  q = par[5]*r*r*(1 + z + 0.5*z*z);
    214279  /* note difference from gaussian: q = par[5]*r */
    215   f = q + par[6];
    216280
    217281  dpar[0] = q*(2*px*par[2] + par[4]*Y);
    … …  
    230294
    231295  float X, Y, px, py;
    232   float z, r, q, f;
     296  float z, r, q, f, k;
    233297
    234298  X = x - par[0];
    … …  
    239303
    240304  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    241   r = 1.0 / (1 + z + par[7]*0.5*z*z*(1 + par[8]*z/3)); /* ~ exp (-Z) */
    242   q = par[5]*r*r*(1 + z + 0.5*z*z);
     305  k = 0.5*z*z*(1 + par[8]*z/3);
     306  r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
     307  f = par[5]*r + par[6];
     308  q = par[5]*r*r*(1 + par[7]*z*(1 + par[8]*z/2));
    243309  /* note difference from gaussian: q = par[5]*r */
    244   f = q + par[6];
    245310
    246311  dpar[0] = q*(2*px*par[2] + par[4]*Y);
    … …  
    251316  dpar[5] = +r;
    252317  dpar[6] = +1;
    253   dpar[7] = +par[5]*r*z*z/2;
    254   dpar[8] = +dpar[7]*z/3;
     318  dpar[7] = -100*par[5]*r*r*k;
     319  dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
    255320
    256321  return (f);
    … …  
    261326
    262327  float X, Y, px, py;
    263   float z, r, q, f;
     328  float z, r, q, f, k;
    264329
    265330  X = x - par[0];
    … …  
    270335
    271336  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    272   r = 1.0 / (1 + z + par[7]*pow(z,par[8])); /* ~ exp (-Z) */
    273   q = par[5]*r*r*par[7]*par[8]*pow(z,par[8]-1.0);
     337  k = pow(z,1.75*par[8]);
     338  r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
     339  q = par[5]*r*r*(1 + 1.75*par[7]*par[8]*pow(z,1.75*par[8]-1));
    274340  /* note difference from gaussian: q = par[5]*r */
    275   f = q + par[6];
     341  f = par[5]*r + par[6];
    276342
    277343  dpar[0] = q*(2*px*par[2] + par[4]*Y);
    … …  
    282348  dpar[5] = +r;
    283349  dpar[6] = +1;
     350  dpar[7] = -10*par[5]*r*r*k;
     351  dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z);
     352
     353  return (f);
     354}
     355
     356/* pars: x, y, sx, sy, sxy, I, sky */
     357float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     358
     359  float X, Y, px1, py1, px2, py2;
     360  float z1, z2, r1, r2, q1, q2, f;
     361
     362  X = x - par[0];
     363  Y = y - par[1];
     364 
     365  px1 = par[2]*X;
     366  py1 = par[3]*Y;
     367  px2 = par[8]*X;
     368  py2 = par[9]*Y;
     369
     370  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
     371  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[10]*X*Y;
     372
     373  r1 = 1.0 / (1 + z1 + 0.5*SQ(z1)*(1 + z1/3)); /* ~ exp (-Z) */
     374  r2 = 1.0 / (1 + pow(z2,1.75));
     375
     376  f = par[5]*r1 + par[6] + par[7]*r2;
     377
     378  q1 = par[5]*SQ(r1)*(1 + z1 + 0.5*SQ(z1));
     379  q2 = par[7]*SQ(r2)*(1.75*pow(z2,0.75));
     380
     381  dpar[ 0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[8] + par[10]*Y);
     382  dpar[ 1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[9] + par[10]*X);
     383  dpar[ 2] = -2*q1*px1*X;
     384  dpar[ 3] = -2*q1*py1*Y;
     385  dpar[ 4] = -q1*X*Y;
     386  dpar[ 5] = +r1;
     387  dpar[ 6] = +1;
     388  dpar[ 7] = +r2*2;
     389  dpar[ 8] = -2*q2*px2*X*2;
     390  dpar[ 9] = -2*q2*py2*Y*2;
     391  dpar[10] = -q2*X*Y;
     392
     393  return (f);
     394}
     395
     396/* pars: x, y, sx, sy, sxy, I, sky */
     397float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     398
     399  float X, Y, px1, py1, px2, py2;
     400  float z1, z2, r, q1, q2, f;
     401
     402  X = x - par[0];
     403  Y = y - par[1];
     404 
     405  px1 = par[2]*X;
     406  py1 = par[3]*Y;
     407  px2 = par[7]*X;
     408  py2 = par[8]*Y;
     409
     410  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
     411  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
     412
     413  r = 1.0 / (1 + z1 + pow(z2,Npow));
     414  f = par[5]*r + par[6];
     415
     416  q1 = par[5]*SQ(r);
     417  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
     418
     419  dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
     420  dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
     421  dpar[2] = -2*q1*px1*X;
     422  dpar[3] = -2*q1*py1*Y;
     423  dpar[4] = -q1*X*Y;
     424  dpar[5] = +r;
     425  dpar[6] = +1;
     426  dpar[7] = -2*q2*px2*X;
     427  dpar[8] = -2*q2*py2*Y;
     428  dpar[9] = -q2*X*Y;
    284429
    285430  return (f);

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

    r2854 r2876  
    55  int i, j, n;
    66  int Nx, Ny, NX, NY;
    7   float Xmin, Xmax, Ymin, Ymax;
     7  int Xmin, Xmax, Ymin, Ymax;
    88  float *v, *x, *y, *z;
    99  Buffer *bf;
    … …  
    4949    }
    5050  }
     51  if (n != Nx*Ny) {
     52    fprintf (stderr, "error in ungrid: %d vs %d (%d x %d)\n", n, Nx*Ny, Nx, Ny);
     53  }
    5154  vx[0].Nelements = vy[0].Nelements = vz[0].Nelements = n;
    5255  return (TRUE);
Note: See TracChangeset for help on using the changeset viewer.