IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42078


Ignore:
Timestamp:
Feb 28, 2022, 12:10:28 PM (4 years ago)
Author:
eugene
Message:

add Molweide projection; improvements to source fitting functions (imfit) for testing

Location:
trunk/Ohana/src/opihi/cmd.astro
Files:
8 edited
5 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/Makefile

    r41666 r42078  
    9494$(SRC)/imfit-qfgauss.$(ARCH).o     \
    9595$(SRC)/imfit-qrgauss.$(ARCH).o     \
     96$(SRC)/imfit-rgauss.$(ARCH).o      \
     97$(SRC)/imfit-fgauss-pol.$(ARCH).o          \
     98$(SRC)/imfit-rgauss-pol.$(ARCH).o          \
    9699$(SRC)/imfit-trail.$(ARCH).o
    97100

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

    r39457 r42078  
    11# include "imfit.h"
     2
     3/** fgauss : a real 2D Gaussian **/
    24
    35opihi_flt fgaussTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *);

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

    r39457 r42078  
    2525  par[4]  = get_variable_default ("SXYg", 0);
    2626  par[5]  = get_variable_default ("Zpk", 10000);
    27   par[6]  = get_variable_default ("Sg", 0.0);
     27  par[6]  = get_variable_default ("Sg", 0.0); // sky
    2828  par[7]  = get_variable_default ("Sr", 1.0);
    2929  fpar[0] = get_variable_default ("Npow", 2.25);
    … …  
    4040  set_variable ("Zpk",  par[5]);
    4141  set_variable ("Sg",   par[6]);
    42   set_variable ("Sr", par[7]);
     42  set_variable ("Sr",   par[7]);
    4343}
    4444
    45 /* one component, two slopes: (1 + z^M + z^N)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */
     45/* qgauss: (1 + Sr*z + z^Npow)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */
    4646opihi_flt qgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) {
    4747  OHANA_UNUSED_PARAM(Npar);
    … …  
    6060  r = 1.0 / (1 + par[7]*z + pow(z,fpar[0]));
    6161  f = par[5]*r + par[6];
    62   q = par[5]*SQ(r)*(par[7] + fpar[0]*pow(z,(fpar[0]-1)));
    6362
    6463  if (dpar != NULL) {
     64    q = par[5]*SQ(r)*(par[7] + fpar[0]*pow(z,(fpar[0]-1)));
    6565    dpar[0] = q*(2*px*par[2] + par[4]*Y);
    6666    dpar[1] = q*(2*py*par[3] + par[4]*X);
    … …  
    7575}
    7676
     77/******************************************************************************
     78 * this file defines the QGAUSS source shape model.  Note that these model functions are
     79 * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
     80 * their own.  The models use a psVector to represent the set of parameters, with the
     81 * sequence used to specify the meaning of the parameter.  The meaning of the parameters
     82 * may thus vary depending on the specifics of the model.  All models which are used as a
     83 * PSF representations share a few parameters, for which # define names are listed in
     84 * pmModel.h:
     85
     86   power-law with fitted linear term
     87   1 / (1 + kz + z^2.25)
     88
     89   * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
     90   * PM_PAR_I0 1    - central intensity
     91   * PM_PAR_XPOS 2  - X center of object
     92   * PM_PAR_YPOS 3  - Y center of object
     93   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
     94   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
     95   * PM_PAR_SXY 6   - X*Y term of elliptical contour
     96   * PM_PAR_7   7   - amplitude of the linear component (k)
     97   *****************************************************************************/

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

    r39457 r42078  
    5959
    6060  r = 1.0 / (1 + fpar[0]*z + pow(z,par[7]));
    61   f = par[5]*r + par[6];
     61  f = par[5]*r + par[6]; // Io * f(r) + Sky
    6262  q = par[5]*SQ(r)*(fpar[0] + par[7]*pow(z,(par[7]-1)));
    6363

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

    r39457 r42078  
    99
    1010  fitfunc = rgaussTD;
    11   cleanup = rgaussCL;
    12   Npar = 10;
    13   Nfpar = 1;
     11  imfit_cleanup = rgaussCL;
     12  Npar =  8;
     13  Nfpar = 0;
    1414
    1515  /* allocate free and fixed parameters */
    … …  
    2121  par[0] = get_variable_default ("Xg", 0);
    2222  par[1] = get_variable_default ("Yg", 0);
    23   par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0);
    24   par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0);
    25   par[4] = 0.0;
     23  par[2] = 2.35 / get_variable_default ("SXg", 2.0);
     24  par[3] = 2.35 / get_variable_default ("SYg", 2.0);
     25  par[4] = get_variable_default ("SXYg", 0);
    2626  par[5] = get_variable_default ("Zpk", 10000);
    27   par[6] = get_variable_default ("Sg", 0.0);
    28   par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
    29   par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
    30   par[9] = get_variable_default ("SXYf", 0.0);
     27  par[6] = get_variable_default ("Sg", 0.0); // sky
     28  par[7] = get_variable_default ("Sr", 2.0);
    3129
    32   fpar[0] = get_variable_default ("Npow", 2.25);
     30//  fpar[0] = get_variable_default ("Npow", 2.25);
    3331
    3432  sky = &par[6];
    3533}
    3634
    37 /* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
     35void rgaussCL () {
     36  set_variable ("Xg",   par[0]);
     37  set_variable ("Yg",   par[1]);
     38  set_variable ("SXg",  2.35 / par[2]);
     39  set_variable ("SYg",  2.35 / par[3]);
     40  set_variable ("SXYg", par[4]);
     41  set_variable ("Zpk",  par[5]);
     42  set_variable ("Sg",   par[6]);
     43  set_variable ("Sr",   par[7]);
     44}
     45
     46/* rgauss: (1 + z + z^alpha)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */
    3847opihi_flt rgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) {
    3948  OHANA_UNUSED_PARAM(Npar);
    4049
    41   opihi_flt X, Y, px1, py1, px2, py2;
    42   opihi_flt z1, z2, r, q1, q2, f;
     50  opihi_flt X = x - par[0];
     51  opihi_flt Y = y - par[1];
     52 
     53  opihi_flt px = par[2]*X;
     54  opihi_flt py = par[3]*Y;
    4355
    44   X = x - par[0];
    45   Y = y - par[1];
    46  
    47   px1 = par[2]*X;
    48   py1 = par[3]*Y;
    49   px2 = par[7]*X;
    50   py2 = par[8]*Y;
     56  opihi_flt z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    5157
    52   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    53   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    54 
    55   r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,fpar[0]));
    56   f = par[5]*r + par[6];
    57 
    58   q1 = par[5]*SQ(r)*(1 + z1);
    59   q2 = par[5]*SQ(r)*fpar[0]*pow(z2,(fpar[0]-1));
     58  opihi_flt p = pow(z,par[7] - 1.0);
     59  opihi_flt r = 1.0 / (1 + z + z*p);
     60  opihi_flt f = par[5]*r + par[6];
    6061
    6162  if (dpar != NULL) {
    62     dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    63     dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    64     dpar[2] = -2*q1*px1*X;
    65     dpar[3] = -2*q1*py1*Y;
    66     dpar[4] = -q1*X*Y;
     63    opihi_flt t = par[5]*SQ(r);
     64    opihi_flt q = t*(1 + par[7]*p);
     65
     66    dpar[0] = +q*(par[2]*px + par[4]*Y);
     67    dpar[1] = +q*(par[3]*py + par[4]*X);
     68    dpar[2] = -q*px*X;
     69    dpar[3] = -q*py*Y;
     70    dpar[4] = -q*X*Y;
    6771    dpar[5] = +r;
    6872    dpar[6] = +1;
    69     dpar[7] = -2*q2*px2*X;
    70     dpar[8] = -2*q2*py2*Y;
    71     dpar[9] = -q2*X*Y;
     73
     74    // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0
     75    dpar[7] = (z == 0.0) ? 0.0 : -t*log(z)*p*z;
    7276  }
    7377  return (f);
    7478}
    7579
    76 int rgaussCL () {
    77   set_variable ("Xg",   par[0]);
    78   set_variable ("Yg",   par[1]);
    79   set_variable ("SXg",  2.35 * sqrt(2.0) / par[2]);
    80   set_variable ("SYg",  2.35 * sqrt(2.0) / par[3]);
    81   set_variable ("SXYg", par[4]);
    82   set_variable ("Zpk",  par[5]);
    83   set_variable ("Sg",   par[6]);
    84   set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]);
    85   set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
    86   set_variable ("SXYf", par[9]);
    87 }
     80/******************************************************************************
     81 * this file defines the RGAUSS source shape model (XXX need a better name!).  Note that these
     82 * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
     83 * statements of their own.  The models use a psVector to represent the set of parameters, with
     84 * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
     85 * may thus vary depending on the specifics of the model.  All models which are used as a PSF
     86 * representations share a few parameters, for which # define names are listed in pmModel.h:
     87
     88   power-law with fitted slope
     89   1 / (1 + z + z^alpha)
     90
     91 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
     92 * PM_PAR_I0 1    - central intensity
     93 * PM_PAR_XPOS 2  - X center of object
     94 * PM_PAR_YPOS 3  - Y center of object
     95 * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
     96 * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     97 * PM_PAR_SXY 6   - X*Y term of elliptical contour
     98 * PM_PAR_7   7   - power-law slope (alpha)
     99 *****************************************************************************/
     100

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

    r41666 r42078  
    44void  trailCL ();
    55
     6// fitted parameters:
    67# define PAR_X      0
    78# define PAR_Y      1
    89# define PAR_THETA  2
    9 # define PAR_SIGMA  3
    10 # define PAR_LENGTH 4
    11 # define PAR_I0     5
    12 # define PAR_SKY    6
     10# define PAR_LENGTH 3
     11# define PAR_I0     4
     12# define PAR_SKY    5
     13
     14// fixed parameters:
     15# define PAR_SIGMA  0
    1316
    1417void trail_setup (char *name) {
    … …  
    1821  fitfunc = trailTD;
    1922  imfit_cleanup = trailCL;
    20   Npar = 7;
    21   Nfpar = 0;
     23  Npar  = 6;
     24  Nfpar = 1;
    2225
    2326  /* allocate free and fixed parameters */
    … …  
    3033  par[PAR_Y      ] = get_variable_default ("Yg",      0.0);
    3134  par[PAR_THETA  ] = get_variable_default ("Tg",      0.0);
    32   par[PAR_SIGMA  ] = get_variable_default ("Wg",      2.0);
    3335  par[PAR_LENGTH ] = get_variable_default ("Lg",     10.0);
    3436  par[PAR_I0     ] = get_variable_default ("Zpk", 10000.0);
    3537  par[PAR_SKY    ] = get_variable_default ("Sg",      0.0);
    3638  sky = &par[PAR_SKY];
     39
     40  fpar[PAR_SIGMA ] = get_variable_default ("Wg",      2.0);
    3741}
    3842
    … …  
    4044  set_variable ("Xg",  par[PAR_X     ]);
    4145  set_variable ("Yg",  par[PAR_Y     ]);
    42   set_variable ("Wg",  par[PAR_SIGMA ]);
    4346  set_variable ("Tg",  par[PAR_THETA ]);
    4447  set_variable ("Lg",  par[PAR_LENGTH]);
    4548  set_variable ("Zpk", par[PAR_I0    ]);
    4649  set_variable ("Sg",  par[PAR_SKY   ]);
     50
     51  set_variable ("Wg",  fpar[PAR_SIGMA]);
    4752}
    4853
    … …  
    5459  opihi_flt Y = y - par[PAR_Y];
    5560 
    56   opihi_flt S2 = 2.0 * SQ(par[PAR_SIGMA]);
     61  opihi_flt S2 = 2.0 * SQ(fpar[PAR_SIGMA]);
    5762
    5863  opihi_flt ST = sin(RAD_DEG*par[PAR_THETA]);
    … …  
    7883
    7984    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
     85    // since X = X - par[PAR_X], dFoo/dXo = -dFoo/dX
    8086    float dRdX = +ST;
    8187    float dRdY = -CT;
    82     float dRdT = -Y*ST - X*CT;
     88    float dRdT = (-Y*ST - X*CT)*RAD_DEG;
     89    // note PAR_THETA is in degrees
    8390
    8491    float dGdX = dGdR * dRdX;
    8592    float dGdY = dGdR * dRdY;
    8693    float dGdT = dGdR * dRdT;
     94    // dGdL is 0.0 because dRdL is 0.0 (R is not a function of L)
    8795
    8896    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
    … …  
    96104    float dZmdL = -0.5 / sqrt(S2);
    97105
    98     float dZpdT = (-X*ST + Y*CT) / sqrt(S2);
     106    // note PAR_THETA is in degrees
     107    float dZpdT = (-X*ST + Y*CT) * RAD_DEG / sqrt(S2);
    99108    float dZmdT = dZpdT; // dZpdT = dZmdT
    100109
    … …  
    119128    // dGdL is 0.0 because dRdL is 0.0
    120129    float dPdL = Gxy * (dEpdL - dEmdL);
    121 
    122130    float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT);
    123131
    … …  
    127135    dpar[PAR_LENGTH] = par[PAR_I0] * dPdL;
    128136    dpar[PAR_THETA]  = par[PAR_I0] * dPdT;
    129     dpar[PAR_SIGMA]  = 0;       // we don't actually allow this to vary, so we do not need to calculate it
    130137  }
    131138

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

    r41666 r42078  
    33int imfit (int argc, char **argv) {
    44
    5   int i, j, N, Npts, Save, VERBOSE;
    6   int sx, sy, nx, ny, Nx, Ny;
    7   float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold;
    8   opihi_flt *x, *y, *z, *dz;
    9   float *V;
     5  int N;
    106  Buffer *buf;
    117
    12   Save = FALSE;
     8  char *Save = NULL;
    139  if ((N = get_argument (argc, argv, "-save"))) {
    1410    remove_argument (N, &argc, argv);
    15     Save = TRUE;
    16   }
    17 
    18   // int ShapeVariation = FALSE;
    19   // if ((N = get_argument (argc, argv, "-shapes"))) {
    20   //   remove_argument (N, &argc, argv);
    21   //   ShapeVariation = TRUE;
    22   // }
    23 
    24   SatThreshold = 0xffff;
     11    Save = strcreate (argv[N]);
     12    remove_argument (N, &argc, argv);
     13  }
     14
     15  int Insert = FALSE;
     16  if ((N = get_argument (argc, argv, "-insert"))) {
     17    remove_argument (N, &argc, argv);
     18    Insert = TRUE;
     19    if (Save) { gprint (GP_ERR, "-save and -insert are mutually exclusive\n"); free (Save); return (FALSE); }
     20  }
     21
     22  int minIter = 5;
     23  if ((N = get_argument (argc, argv, "-min-iter"))) {
     24    remove_argument (N, &argc, argv);
     25    minIter = atoi (argv[N]);
     26    remove_argument (N, &argc, argv);
     27  }
     28
     29  int SatThreshold = 0xffff;
    2530  if ((N = get_argument (argc, argv, "-sat"))) {
    2631    remove_argument (N, &argc, argv);
    … …  
    3035
    3136  /* Gain in e/DN */
    32   Gain = 1.0;
     37  float Gain = 1.0;
    3338  if ((N = get_argument (argc, argv, "-gain"))) {
    3439    remove_argument (N, &argc, argv);
    … …  
    3843
    3944  /* RD noise in DN */
    40   RDnoise = 0.0;
     45  float RDnoise = 0.0;
    4146  if ((N = get_argument (argc, argv, "-rdnoise"))) {
    4247    remove_argument (N, &argc, argv);
    … …  
    4550  }
    4651
    47   VERBOSE = FALSE;
     52  int VERBOSE = FALSE;
    4853  if ((N = get_argument (argc, argv, "-v"))) {
    4954    remove_argument (N, &argc, argv);
    … …  
    5156  }
    5257
    53   /* set fitting function */
     58  /* set fitting function : defines par, Npar, fitfunc, etc globals (imfit.h) */
    5459  fgauss_setup ("fgauss");
    5560  if ((N = get_argument (argc, argv, "-func"))) {
    5661    fitfunc = NULL;
    5762    remove_argument (N, &argc, argv);
    58     fgauss_setup (argv[N]);
    59     pgauss_setup (argv[N]);
    60     pgauss_psf_setup (argv[N]);
     63    fgauss_setup (argv[N]); // OK
     64    pgauss_setup (argv[N]); // OK
    6165    sgauss_setup (argv[N]);
    62     sgauss_psf_setup (argv[N]);
    63     qgauss_setup (argv[N]);
    64     qgauss_psf_setup (argv[N]);
    65     qfgauss_setup (argv[N]);
     66    qgauss_setup (argv[N]); // OK
     67    rgauss_setup (argv[N]);
     68    qfgauss_setup (argv[N]); // OK
    6669    qrgauss_setup (argv[N]);
    6770    trail_setup (argv[N]);
     71    pgauss_psf_setup (argv[N]);
     72    sgauss_psf_setup (argv[N]);
     73    qgauss_psf_setup (argv[N]);
     74
     75    fgauss_pol_setup (argv[N]);
     76    rgauss_pol_setup (argv[N]);
     77
    6878    if (fitfunc == NULL) {
    6979      gprint (GP_ERR, "unknown function %s\n", argv[N]);
     80      FREE (Save);
    7081      return (FALSE);
    7182    }
    … …  
    7485
    7586  if (argc != 6) {
    76     gprint (GP_ERR, "USAGE: imfit <buffer> sx sy nx ny\n");
     87    gprint (GP_ERR, "USAGE: imfit <buffer> Xo Yo dX dY\n");
     88    gprint (GP_ERR, "options: [-save buffer] [-insert] [-sat value] [-gain value] [-rdnoise value] [-v] [-func option]\n");
     89    gprint (GP_ERR, "   (Xo,Yo) : center\n");
     90    gprint (GP_ERR, "   (dX,dY) : window size\n");
     91    FREE (Save);
    7792    return (FALSE);
    7893  }
    7994
    8095  /* non-optional arguments */
    81   if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    82   sx = atof (argv[2]);
    83   sy = atof (argv[3]);
    84   nx = atof (argv[4]);
    85   ny = atof (argv[5]);
    86   Nx = buf[0].matrix.Naxis[0];
    87   Ny = buf[0].matrix.Naxis[1];
    88 
    89   /* check if region is valid */
    90   if (sx + 0.5*nx < 0) goto range;
    91   if (sy + 0.5*ny < 0) goto range;
    92   if (sx + 0.5*nx >= Nx) goto range;
    93   if (sy + 0.5*ny >= Ny) goto range;
     96  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) { FREE (Save); return (FALSE); }
     97  int Xo = atof (argv[2]);
     98  int Yo = atof (argv[3]);
     99  int dX = atof (argv[4]);
     100  int dY = atof (argv[5]);
     101  int Nx = buf[0].matrix.Naxis[0];
     102  int Ny = buf[0].matrix.Naxis[1];
     103
     104  int sx = Xo - dX/2;
     105  int sy = Yo - dY/2;
     106
     107  /* check if region is valid (center must be in range of image pixels) */
     108  if (Xo < 0) goto range;
     109  if (Yo < 0) goto range;
     110  if (Xo >= Nx) goto range;
     111  if (Yo >= Ny) goto range;
     112
     113  // image value in DN
     114  // rdnoise in DN
     115  // sigma_DN^2 = sigma_e^2 / gain^2
     116  // sigma_e^2  = Ne
     117  // sigma_e^2  = DN * gain
     118  // sigma_DN^2 = DN * gain / gain^2 = DN / gain
     119
     120  if (Insert) {
     121    float *Vi = (float *)buf[0].matrix.buffer;
     122    for (int j = 0; j < dY; j++) {
     123      for (int i = 0; i < dX; i++) {
     124        float vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
     125        Vi[(i+sx)+(j+sy)*Nx] += vf;
     126      }
     127    }
     128    return TRUE;
     129  }
    94130
    95131  /* convert array z[x,y] to x[i], y[i], z[i] */
    96132  N = 0;
    97   Npts = nx*ny;
    98   ALLOCATE (x,  opihi_flt, 2*Npts);
    99   ALLOCATE (y,  opihi_flt, 2*Npts);
    100   ALLOCATE (z,  opihi_flt, 2*Npts);
    101   ALLOCATE (dz, opihi_flt, 2*Npts);
    102   for (j = 0; j < ny; j++) {
     133  int Npts = dX*dY;
     134  ALLOCATE_PTR (x,  opihi_flt, 2*Npts);
     135  ALLOCATE_PTR (y,  opihi_flt, 2*Npts);
     136  ALLOCATE_PTR (z,  opihi_flt, 2*Npts);
     137  ALLOCATE_PTR (dz, opihi_flt, 2*Npts);
     138  for (int j = 0; j < dY; j++) {
    103139    if (j + sy < 0) continue;
    104140    if (j + sy >= Ny) continue;
    105     V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    106     for (i = 0; i < nx; i++) {
     141    float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
     142    for (int i = 0; i < dX; i++) {
    107143      if (i + sx < 0) continue;
    108144      if (i + sx >= Nx) continue;
    109       if (*V > SatThreshold) goto next;
    110       dz[N] = (SQ(RDnoise) + *V/Gain);
     145      if (*V > SatThreshold) goto next; // skip pixels above threshold
     146      if (!isfinite(*V)) goto next; // skip nan pixels
     147      dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
    111148      if (dz[N] <= 0) goto next;
    112149      dz[N] = 1.0 / dz[N];
    … …  
    122159
    123160  /* run fit routine */
    124   ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    125   dchisq = ochisq;
    126   chisq  = ochisq;
    127 
    128 //for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
    129 
    130   for (i = 0; (i < 25); i++) {
     161  float ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     162  float dchisq = ochisq;
     163  float chisq  = ochisq;
     164
     165  int Niter = 0;
     166  // for (int i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     167
     168  // keep iterating if chisq is increasing or
     169  for (Niter = 0; (Niter < 25) && ((Niter < minIter) || ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar)))); Niter++) {
    131170    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    132171    dchisq = ochisq - chisq;
     172    // fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    133173    ochisq = chisq;
    134     fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    135174  } 
    136   set_int_variable ("Niter",  i);
     175  set_int_variable ("Niter",  Niter);
    137176
    138177  /** create output image (keep in sky) **/
    … …  
    141180    float *Vi, *Vo, vr, vf;
    142181
    143     if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
     182    if ((out = SelectBuffer (Save, ANYBUFFER, TRUE)) == NULL) { free (Save); return (FALSE); }
    144183    free (out[0].header.buffer);
    145184    free (out[0].matrix.buffer);
    146185
    147186    strcpy (out[0].file, "(empty)");
    148     if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE;
    149 
    150     /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
     187    if (!CreateBuffer (out, 2*dX, 2*dY, -32, 0.0, 1.0)) { free (Save); return FALSE; }
     188
     189    /* four panels: 1) raw image. 2) fit  3) raw - fit  4) absolute deviation */
    151190    Vi = (float *)buf[0].matrix.buffer;
    152191    Vo = (float *)out[0].matrix.buffer;
    153     for (j = 0; j < ny; j++) {
    154       for (i = 0; i < nx; i++) {
     192    for (int j = 0; j < dY; j++) {
     193      for (int i = 0; i < dX; i++) {
    155194        vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
    156195        vr = Vi[(i+sx)+(j+sy)*Nx];
    157         Vo[(i   )+(j   )*2*nx] = vr;
    158         Vo[(i+nx)+(j   )*2*nx] = vf;
    159         Vo[(i   )+(j+ny)*2*nx] = vr - vf + *sky;
    160         Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky;
     196        Vo[(i   )+(j   )*2*dX] = vr;
     197        Vo[(i+dX)+(j   )*2*dX] = vf;
     198        Vo[(i   )+(j+dY)*2*dX] = vr - vf + *sky;
     199        Vo[(i+dX)+(j+dY)*2*dX] = fabs(vr-vf) + *sky;
    161200      }
    162201    }
     202    free (Save);
    163203  }
    164204
    … …  
    169209
    170210  if (VERBOSE) {
    171     for (i = 0; i < Npar; i++) {
     211    for (int i = 0; i < Npar; i++) {
    172212      gprint (GP_ERR, "%g ", par[i]);
    173213    }

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

    r41515 r42078  
    159159    if (!strcasecmp (argv[CtypeArg], "GLS")) { strcpy (graphmode.coords.ctype, "DEC--GLS"); goto got_ctype; }
    160160    if (!strcasecmp (argv[CtypeArg], "PAR")) { strcpy (graphmode.coords.ctype, "DEC--PAR"); goto got_ctype; }
     161    if (!strcasecmp (argv[CtypeArg], "MOL")) { strcpy (graphmode.coords.ctype, "DEC--MOL"); goto got_ctype; }
    161162    gprint (GP_ERR, "ERROR: invalid projection type %s\n", argv[CtypeArg]);
    162     gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR\n");
     163    gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR, MOL\n");
    163164    return FALSE;
    164165  }
Note: See TracChangeset for help on using the changeset viewer.