IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3144


Ignore:
Timestamp:
Feb 5, 2005, 6:57:15 PM (21 years ago)
Author:
eugene
Message:

added imsub, new imfit functions

Location:
trunk/Ohana/src/opihi
Files:
17 added
3 edited

Legend:

Unmodified
Added
Removed

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

    r2843 r3144  
    3535$(SDIR)/gauss.$(ARCH).o            \
    3636$(SDIR)/getvel.$(ARCH).o           \
    37 $(SDIR)/imfit.$(ARCH).o            \
    3837$(SDIR)/medianmap.$(ARCH).o        \
    3938$(SDIR)/mkgauss.$(ARCH).o          \
    … …  
    5352$(SDIR)/transform.$(ARCH).o
    5453
     54imfitstuff = \
     55$(SDIR)/imsub.$(ARCH).o            \
     56$(SDIR)/imfit.$(ARCH).o            \
     57$(SDIR)/imfit-fgauss.$(ARCH).o     \
     58$(SDIR)/imfit-pgauss.$(ARCH).o     \
     59$(SDIR)/imfit-Pgauss.$(ARCH).o     \
     60$(SDIR)/imfit-qgauss.$(ARCH).o     \
     61$(SDIR)/imfit-Qgauss.$(ARCH).o     \
     62$(SDIR)/imfit-sgauss.$(ARCH).o     \
     63$(SDIR)/imfit-Sgauss.$(ARCH).o     \
     64$(SDIR)/imfit-qfgauss.$(ARCH).o    \
     65$(SDIR)/imfit-qrgauss.$(ARCH).o   
     66
    5567libastrocmd:    $(DESTLIB)/libastrocmd.a
    5668
    57 $(LIB)/libastrocmd.$(ARCH).a:           $(astrocmd)
     69$(LIB)/libastrocmd.$(ARCH).a:           $(astrocmd) $(imfitstuff)
    5870
    5971# utilities #################################################

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

    r3105 r3144  
    1 # include "astro.h"
    2 # define FFACTOR 200
    3 # define FSCALE 1.2
    4 
    5 float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts,
    6                 float *par, int Npar,
    7                 float (funcs)(float, float, float *, int, float *));
    8 
    9 float fgaussTD (float, float, float *, int, float *);
    10 float pgaussTD (float, float, float *, int, float *);
    11 float sgaussTD (float, float, float *, int, float *);
    12 float vgaussTD (float, float, float *, int, float *);
    13 
    14 float PgaussTD (float, float, float *, int, float *);
    15 float SgaussTD (float, float, float *, int, float *);
    16 
    17 float tgaussTD (float, float, float *, int, float *);
    18 float rgaussTD (float, float, float *, int, float *);
    19 
    20 float sersic (float, float, float *, int, float *);
    21 float serbulge (float, float, float *, int, float *);
    22 
    23 static float Npow;
    24 static float Npin;
    25 static int MissingVariable = FALSE;
    26 static float *fpar;
    27 static int Nfpar;
    28 
    29 int mrq2dlimits (float *pmin, float *pmax, int Npar);
    30 
    31 float get_variable_default (char *name, float dvalue) {
    32 
    33   char *value;
    34   float fvalue;
    35 
    36   value = get_variable (name);
    37   if (value == NULL) {
    38     MissingVariable = TRUE;
    39     return (dvalue);
    40   }
    41   fvalue = atof (value);
    42   return (fvalue);
    43 }
     1# include "imfit.h"
    442
    453int imfit (int argc, char **argv) {
    464
    47   int i, j, N, Npts, Npar, Save, VERBOSE, ShapeVariation;
     5  int i, j, N, Npts, Save, VERBOSE, ShapeVariation;
    486  int sx, sy, nx, ny, Nx, Ny;
    49   float chisq, ochisq, dchisq, Gain, RDnoise;
    50   float *x, *y, *z, *dz, *par, *V, *parmin, *parmax;
     7  float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold;
     8  float *x, *y, *z, *dz, *V, *parmin, *parmax;
    519  char line[64];
    5210  Buffer *buf;
    53   float (*fitfunc)(float, float, float *, int, float *);
    5411
    5512  Save = FALSE;
    … …  
    6320    remove_argument (N, &argc, argv);
    6421    ShapeVariation = TRUE;
     22  }
     23
     24  SatThreshold = 0xffff;
     25  if ((N = get_argument (argc, argv, "-sat"))) {
     26    remove_argument (N, &argc, argv);
     27    SatThreshold = atof (argv[N]);
     28    remove_argument (N, &argc, argv);
    6529  }
    6630
    … …  
    8852
    8953  /* set fitting function */
    90   Nfpar = 0;
    91   Npar = 7;
    92   fitfunc = fgaussTD;
     54  fgauss_setup ("fgauss");
    9355  if ((N = get_argument (argc, argv, "-func"))) {
    9456    fitfunc = NULL;
    9557    remove_argument (N, &argc, argv);
    96     if (!strcmp(argv[N], "fgauss")) {
    97       fitfunc = fgaussTD;
    98     }
    99     if (!strcmp(argv[N], "pgauss")) {
    100       fitfunc = pgaussTD;
    101     }
    102     if (!strcmp(argv[N], "vgauss")) {
    103       fitfunc = vgaussTD;
    104       Npar = 9;
    105     }
    106     if (!strcmp(argv[N], "sgauss")) {
    107       fitfunc = sgaussTD;
    108       Npar = 10;
    109     }
    110     if (!strcmp(argv[N], "rgauss")) {
    111       fitfunc = rgaussTD;
    112       Npar = 10;
    113     }
    114     if (!strcmp(argv[N], "tgauss")) {
    115       fitfunc = tgaussTD;
    116       Npar = 10;
    117     }
    118     if (!strcmp(argv[N], "Pgauss")) {
    119       fitfunc = PgaussTD;
    120       Npar = 4;
    121       Nfpar = 3;
    122     }
    123     if (!strcmp(argv[N], "Sgauss")) {
    124       fitfunc = SgaussTD;
    125       Npar = 4;
    126       Nfpar = 6;
    127     }
    128     if (!strcmp(argv[N], "sersic")) {
    129       fitfunc = sersic;
    130       Npar = 8;
    131     }
    132     if (!strcmp(argv[N], "serbulge")) {
    133       fitfunc = serbulge;
    134       Npar = 12;
    135     }
     58    fgauss_setup (argv[N]);
     59    pgauss_setup (argv[N]);
     60    Pgauss_setup (argv[N]);
     61    sgauss_setup (argv[N]);
     62    Sgauss_setup (argv[N]);
     63    qgauss_setup (argv[N]);
     64    Qgauss_setup (argv[N]);
     65    qfgauss_setup (argv[N]);
     66    qrgauss_setup (argv[N]);
    13667    if (fitfunc == NULL) {
    13768      fprintf (stderr, "unknown function %s\n", argv[N]);
    … …  
    15586  Ny = buf[0].matrix.Naxis[1];
    15687
    157   if (sx < 0) goto range;
    158   if (sy < 0) goto range;
    159   if (sx + nx >= Nx) goto range;
    160   if (sy + ny >= Ny) goto range;
    161 
    162   ALLOCATE (par, float, Npar);
    163   bzero (par, Npar*sizeof(float));
    164   if (Nfpar > 0) {
    165     ALLOCATE (fpar, float, Nfpar);
    166     bzero (fpar, Nfpar*sizeof(float));
    167   }   
     88  /* check if region is valid */
     89  if (sx + 0.5*nx < 0) goto range;
     90  if (sy + 0.5*ny < 0) goto range;
     91  if (sx + 0.5*nx >= Nx) goto range;
     92  if (sy + 0.5*ny >= Ny) goto range;
    16893
    16994  /* convert array z[x,y] to x[i], y[i], z[i] */
    … …  
    175100  ALLOCATE (dz, float, 2*Npts);
    176101  for (j = 0; j < ny; j++) {
     102    if (j + sy < 0) continue;
     103    if (j + sy >= Ny) continue;
    177104    V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    178     for (i = 0; i < nx; i++, V++) {
     105    for (i = 0; i < nx; i++) {
     106      if (i + sx < 0) continue;
     107      if (i + sx >= Nx) continue;
     108      if (*V > SatThreshold) goto next;
    179109      dz[N] = (SQ(RDnoise) + *V/Gain);
    180       if (dz[N] <= 0) continue;
     110      if (dz[N] <= 0) goto next;
    181111      dz[N] = 1.0 / dz[N];
    182112      x[N] = i + sx;
    … …  
    184114      z[N] = *V;
    185115      N++;
     116    next:
     117      V++;
    186118    }
    187119  }
    188120  Npts = N;
    189121
    190   /* make initial guesses on parameters */
    191   if (Npar >= 7) { /* not Pgauss, Sgauss */
    192     par[0] = get_variable_default ("Xg", 0);
    193     par[1] = get_variable_default ("Yg", 0);
    194     par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0);
    195     par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0);
    196     par[4] = 0.0;
    197     par[5] = get_variable_default ("Zpk", 10000);
    198     par[6] = get_variable_default ("Sg", 0.0);
    199   }
    200   if (Npar == 4) { /* Pgauss, Sgauss */
    201     par[0] = get_variable_default ("Xg", 0);
    202     par[1] = get_variable_default ("Yg", 0);
    203     par[2] = get_variable_default ("Zpk", 10000);
    204     par[3] = get_variable_default ("Sg", 0.0);
    205   }
    206   if (fitfunc == vgaussTD) {
    207     par[7] = par[8] = 1;
    208   }
    209   if (fitfunc == sersic) {
    210     par[7] = get_variable_default ("Sr", 1.0);
    211   }
    212   if (fitfunc == serbulge) {
    213     Npow = get_variable_default ("Npow", 2.25);
    214     Npin = get_variable_default ("Npin", 1.00);
    215     par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
    216     par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
    217     par[9] = get_variable_default ("SXYf", 0.0);
    218     par[11] = get_variable_default ("Sr", 1.0);
    219     /* arbitrary guess: equal distribution between bulge and disk */
    220     par[10] = par[5]/2.0;
    221     par[5]  = par[10];
    222   }
    223   if (Npar == 10) { /** sgauss, tgauss, rgauss **/
    224     Npow = get_variable_default ("Npow", 2.25);
    225     Npin = get_variable_default ("Npin", 1.00);
    226     par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
    227     par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
    228     par[9] = get_variable_default ("SXYf", 0.0);
    229   }
    230   if (fitfunc == PgaussTD) {
    231     fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
    232     fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
    233     fpar[2] = get_variable_default ("SXYg", 0.0);
    234   }   
    235   if (fitfunc == SgaussTD) {
    236     fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
    237     fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
    238     fpar[2] = get_variable_default ("SXYg", 0.0);
    239     fpar[3] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
    240     fpar[4] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
    241     fpar[5] = get_variable_default ("SXYf", 0.0);
    242   }   
    243 
    244122  /* run fit routine */
    245123  ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    246124  dchisq = ochisq;
    247 
    248125  for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
    249126    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    … …  
    252129  } 
    253130  set_int_variable ("Niter",  i);
    254  
    255 
    256   if (ShapeVariation) {
    257     /* find dChi/dSx and dChi/dSy given by increasing shape terms by 5% */
    258     float tp1, tp2, chix, chiy;
    259     chix = chiy = 0;
    260     if (fitfunc == sgaussTD) {
    261       tp1 = par[2];
    262       tp2 = par[7];
    263       par[2] = par[2]*1.05;
    264       par[7] = par[7]*1.05;
    265       chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    266       par[2] = tp1;
    267       par[7] = tp2;
    268 
    269       tp1 = par[3];
    270       tp2 = par[8];
    271       par[3] = par[3]*1.05;
    272       par[8] = par[8]*1.05;
    273       chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    274       par[3] = tp1;
    275       par[8] = tp2;
    276     }
    277     if (fitfunc == pgaussTD) {
    278       tp1 = par[2];
    279       par[2] = par[2]*1.05;
    280       chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    281       par[2] = tp1;
    282 
    283       tp1 = par[3];
    284       par[3] = par[3]*1.05;
    285       chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    286       par[3] = tp1;
    287     }
    288     if (fitfunc == SgaussTD) {
    289       tp1 = par[0];
    290       tp2 = par[3];
    291       par[0] = par[0]*1.05;
    292       par[3] = par[3]*1.05;
    293       chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    294       par[0] = tp1;
    295       par[3] = tp2;
    296 
    297       tp1 = par[1];
    298       tp2 = par[4];
    299       par[1] = par[1]*1.05;
    300       par[4] = par[4]*1.05;
    301       chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    302       par[1] = tp1;
    303       par[4] = tp2;
    304     }
    305     if (fitfunc == PgaussTD) {
    306       tp1 = par[0];
    307       par[0] = par[0]*1.05;
    308       chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    309       par[0] = tp1;
    310 
    311       tp1 = par[1];
    312       par[1] = par[1]*1.05;
    313       chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
    314       par[1] = tp1;
    315     }
    316     set_variable ("dChiX", chix/chisq);
    317     set_variable ("dChiY", chiy/chisq);
    318   }
    319131
    320132  /** create output image (keep in sky) **/
    321133  if (Save) {
    322134    Buffer *out;
    323     float *Vi, *Vo, vr, vf, sky;
     135    float *Vi, *Vo, vr, vf;
    324136
    325137    if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
    … …  
    335147    out[0].header.Naxis[1] = 2*ny;
    336148    CreateBuffer (out);
    337     if (Npar == 4) {
    338       sky = par[3];
    339     } else {
    340       sky = par[6];
    341     }
    342149
    343150    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
    … …  
    350157        Vo[(i   )+(j   )*2*nx] = vr;
    351158        Vo[(i+nx)+(j   )*2*nx] = vf;
    352         Vo[(i   )+(j+ny)*2*nx] = vr - vf + sky;
    353         Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + sky;
     159        Vo[(i   )+(j+ny)*2*nx] = vr - vf + *sky;
     160        Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky;
    354161      }
    355162    }
    356163  }
    357164
    358   /* parameter renormalize */
    359   if (Npar >= 7) {
    360     set_variable ("Xg",  par[0]);
    361     set_variable ("Yg",  par[1]);
    362     set_variable ("SXg", 2.35 * sqrt(2.0) / par[2]);
    363     set_variable ("SYg", 2.35 * sqrt(2.0) / par[3]);
    364     set_variable ("SXYg", par[4]);
    365     set_variable ("Zpk", par[5]);
    366     set_variable ("Sg",  par[6]);
    367   }
    368   if (fitfunc == vgaussTD) {
    369     set_variable ("SXf", par[7]);
    370     set_variable ("SYf", par[8]);
    371   }
    372   if (fitfunc == sersic) {
    373     set_variable ("Sr", par[7]);
    374   }
    375   if (fitfunc == serbulge) {
    376     set_variable ("Zb", par[5]);
    377     set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]);
    378     set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
    379     set_variable ("SXYf", par[9]);
    380     set_variable ("Zd", par[10]);
    381     set_variable ("Sr", par[11]);
    382   }
    383   if (Npar == 10) {
    384     set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]);
    385     set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
    386     set_variable ("SXYf", par[9]);
    387   }
    388   if (fitfunc == PgaussTD) {
    389     set_variable ("Xg",  par[0]);
    390     set_variable ("Yg",  par[1]);
    391     set_variable ("Zpk", par[2]);
    392     set_variable ("Sg",  par[3]);
    393     set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
    394     set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
    395     set_variable ("SXYg", fpar[2]);
    396   }
    397   if (fitfunc == SgaussTD) {
    398     set_variable ("Xg",  par[0]);
    399     set_variable ("Yg",  par[1]);
    400     set_variable ("Zpk", par[2]);
    401     set_variable ("Sg",  par[3]);
    402     set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
    403     set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
    404     set_variable ("SXYg", fpar[2]);
    405     set_variable ("SXf", 2.35 * sqrt(2.0) / fpar[3]);
    406     set_variable ("SYf", 2.35 * sqrt(2.0) / fpar[4]);
    407     set_variable ("SXYf", fpar[5]);
    408   }   
     165  /* save parameters to opihi variables */
     166  imfit_cleanup ();
     167
    409168  set_variable ("ChiSq", chisq/(Npts - Npar));
    410169
    … …  
    421180  free (dz);
    422181  free (par);
    423   if (Nfpar > 0) free (fpar);
     182  free (fpar);
    424183
    425184  mrq2dfree (Npar);
    … …  
    431190}
    432191
    433 /* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */
    434 float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    435 
    436   float X, Y, px, py;
    437   float z, r, q, f;
    438 
    439   X = x - par[0];
    440   Y = y - par[1];
    441  
    442   px = par[2]*X;
    443   py = par[3]*Y;
    444 
    445   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    446   r = exp (-z);
    447   q = par[5]*r;
    448   f = q + par[6];
    449 
    450   if (dpar != NULL) {
    451     dpar[0] = q*(2*px*par[2] + par[4]*Y);
    452     dpar[1] = q*(2*py*par[3] + par[4]*X);
    453     dpar[2] = -2*q*px*X;
    454     dpar[3] = -2*q*py*Y;
    455     dpar[4] = -q*X*Y;
    456     dpar[5] = +r;
    457     dpar[6] = +1;
    458   }
    459   return (f);
    460 }
    461 
    462 /* sersic galaxy model -- x, y, sx, sy, sxy, I, sky, n */
    463 /* exp (-b (r/r_e)^(1/n)) */
    464 float sersic (float x, float y, float *par, int Npar, float *dpar) {
    465 
    466   float X, Y, px, py;
    467   float z, r, t, q, f;
    468 
    469   X = x - par[0];
    470   Y = y - par[1];
    471  
    472   px = par[2]*X;
    473   py = par[3]*Y;
    474 
    475   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    476   t = pow (z, par[7]);
    477   r = exp (-t);
    478   f = par[5]*r + par[6];
    479   q = par[5]*r*par[7]*pow(z, par[7]-1);
    480 
    481   if (dpar != NULL) {
    482     dpar[0] = q*(2*px*par[2] + par[4]*Y);
    483     dpar[1] = q*(2*py*par[3] + par[4]*X);
    484     dpar[2] = -2*q*px*X;
    485     dpar[3] = -2*q*py*Y;
    486     dpar[4] = -q*X*Y;
    487     dpar[5] = +r;
    488     dpar[6] = +1;
    489     dpar[7] = -q*log(z)*t;
    490   }
    491   return (f);
    492 }
    493 
    494 /*                                  0  1    2   3   4      5    6     7   8   9       10  11 */
    495 /* sersic galaxy model w/ bulge: -- x, y, (sx, sy, sxy)_1, I_1, sky, (sx, sy, sxy)_2, I_2, n */
    496 /* exp (-b (r/r_e)^(1/n)) + pgauss (r) */
    497 float serbulge (float x, float y, float *par, int Npar, float *dpar) {
    498 
    499   float X, Y, px1, px2, py1, py2;
    500   float z1, z2, r1, r2, t, q1, q2, f;
    501 
    502   X = x - par[0];
    503   Y = y - par[1];
    504  
    505   px1 = par[2]*X;
    506   py1 = par[3]*Y;
    507   px2 = par[7]*X;
    508   py2 = par[8]*Y;
    509 
    510   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    511   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    512 
    513   /* bulge component */
    514   r1 = 1.0 / (1 + z1 + 0.5*z1*z1*(1 + z1/3)); /* ~ exp (-Z) */
    515 
    516   /* disk component */
    517   t = pow (z2, par[11]);
    518   r2 = exp (-t);
    519 
    520   f = par[5]*r1 + par[10]*r2 + par[6];
    521 
    522   q1 = par[5]*r1*r1*(1 + z1 + 0.5*z1*z1);
    523   q2 = par[10]*r2*par[11]*pow(z2, par[11]-1);
    524 
    525   if (dpar != NULL) {
    526     dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    527     dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    528     dpar[2] = -2*q1*px1*X;
    529     dpar[3] = -2*q1*py1*Y;
    530     dpar[4] = -q1*X*Y;
    531     dpar[5] = +r1;
    532     dpar[6] = +1;
    533     dpar[7] = -2*q2*px2*X*50;
    534     dpar[8] = -2*q2*py2*Y*50;
    535     dpar[9] = -q2*X*Y*50;
    536     dpar[10] = +r2*50;
    537     dpar[11] = -q2*log(z2)*t*50;
    538   }
    539   return (f);
    540 }
    541 
    542 /* pseudo 2D gaussian -- x, y, sx, sy, sxy, I, sky */
    543 float pgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    544 
    545   float X, Y, px, py;
    546   float z, r, q, f;
    547 
    548   X = x - par[0];
    549   Y = y - par[1];
    550  
    551   px = par[2]*X;
    552   py = par[3]*Y;
    553 
    554   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    555   r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
    556   f = par[5]*r + par[6];
    557   q = par[5]*r*r*(1 + z + 0.5*z*z);
    558   /* note difference from gaussian: q = par[5]*r */
    559 
    560   if (dpar != NULL) {
    561     dpar[0] = q*(2*px*par[2] + par[4]*Y);
    562     dpar[1] = q*(2*py*par[3] + par[4]*X);
    563     dpar[2] = -2*q*px*X;
    564     dpar[3] = -2*q*py*Y;
    565     dpar[4] = -q*X*Y;
    566     dpar[5] = +r;
    567     dpar[6] = +1;
    568   }
    569   return (f);
    570 }
    571 
    572 /* pseudo 2D gaussian -- x, y, (sx), (sy), (sxy), I, sky */
    573 float PgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    574 
    575   float X, Y, px, py;
    576   float z, r, q, f;
    577 
    578   /* par -> fpar: (2,0), (3,1), (4,2) */
    579 
    580   X = x - par[0];
    581   Y = y - par[1];
    582  
    583   px = fpar[0]*X;
    584   py = fpar[1]*Y;
    585 
    586   z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y;
    587   r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
    588   f = par[2]*r + par[3];
    589   q = par[2]*r*r*(1 + z + 0.5*z*z);
    590   /* note difference from gaussian: q = par[5]*r */
    591 
    592   if (dpar != NULL) {
    593     dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y);
    594     dpar[1] = q*(2*py*fpar[1] + fpar[2]*X);
    595     dpar[2] = +r;
    596     dpar[3] = +1;
    597   }
    598   return (f);
    599 }
    600 
    601 /* pseudo 2D gaussian with floating 2nd and 3rd order terms -- x, y, sx, sy, sxy, I, sky, f1, f2 */
    602 float vgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    603 
    604   float X, Y, px, py;
    605   float z, r, q, f, k;
    606 
    607   X = x - par[0];
    608   Y = y - par[1];
    609  
    610   px = par[2]*X;
    611   py = par[3]*Y;
    612 
    613   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    614   k = 0.5*z*z*(1 + par[8]*z/3);
    615   r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
    616   f = par[5]*r + par[6];
    617   q = par[5]*r*r*(1 + par[7]*z*(1 + par[8]*z/2));
    618   /* note difference from gaussian: q = par[5]*r */
    619 
    620   if (dpar != NULL) {
    621     dpar[0] = q*(2*px*par[2] + par[4]*Y);
    622     dpar[1] = q*(2*py*par[3] + par[4]*X);
    623     dpar[2] = -2*q*px*X;
    624     dpar[3] = -2*q*py*Y;
    625     dpar[4] = -q*X*Y;
    626     dpar[5] = +r;
    627     dpar[6] = +1;
    628     dpar[7] = -100*par[5]*r*r*k;
    629     dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
    630   }
    631   return (f);
    632 }
    633 
    634 /* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
    635 float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    636 
    637   float X, Y, px1, py1, px2, py2;
    638   float z1, z2, r, q1, q2, f, f1, f2;
    639 
    640   X = x - par[0];
    641   Y = y - par[1];
    642  
    643   px1 = par[2]*X;
    644   py1 = par[3]*Y;
    645   px2 = par[7]*X;
    646   py2 = par[8]*Y;
    647 
    648   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    649   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    650 
    651   r = 1.0 / (1 + z1 + pow(z2,Npow));
    652   f = par[5]*r + par[6];
    653 
    654   q1 = par[5]*SQ(r);
    655   q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
    656 
    657   if (dpar != NULL) {
    658     dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    659     dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    660 
    661     /* these fudge factors impede the growth of par[2] beyond par[7] */
    662     f1 = fabs(par[7]) / fabs(par[2]);
    663     f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
    664     dpar[2] = -2*q1*px1*X*f2;
    665 
    666     f1 = fabs(par[8]) / fabs(par[3]);
    667     f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
    668     dpar[3] = -2*q1*py1*Y*f2;
    669 
    670     dpar[4] = -q1*X*Y;
    671     dpar[5] = +r;
    672     dpar[6] = +1;
    673     dpar[7] = -2*q2*px2*X;
    674     dpar[8] = -2*q2*py2*Y;
    675     dpar[9] = -q2*X*Y;
    676   }
    677   return (f);
    678 }
    679 
    680 /* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
    681 float SgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    682 
    683   float X, Y, px1, py1, px2, py2;
    684   float z1, z2, r, q1, q2, f, f1, f2;
    685 
    686   X = x - par[0];
    687   Y = y - par[1];
    688  
    689   px1 = fpar[0]*X;
    690   py1 = fpar[1]*Y;
    691   px2 = fpar[3]*X;
    692   py2 = fpar[4]*Y;
    693 
    694   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + fpar[2]*X*Y;
    695   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + fpar[5]*X*Y;
    696 
    697   r = 1.0 / (1 + z1 + pow(z2,Npow));
    698   f = par[2]*r + par[3];
    699 
    700   q1 = par[2]*SQ(r);
    701   q2 = par[2]*SQ(r)*Npow*pow(z2,(Npow-1));
    702 
    703   if (dpar != NULL) {
    704     dpar[0] = q1*(2*px1*fpar[0] + fpar[2]*Y) + q2*(2*px2*fpar[3] + fpar[5]*Y);
    705     dpar[1] = q1*(2*py1*fpar[1] + fpar[2]*X) + q2*(2*py2*fpar[4] + fpar[5]*X);
    706     dpar[2] = +r;
    707     dpar[3] = +1;
    708   }
    709   return (f);
    710 }
    711 
    712 /* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
    713 float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    714 
    715   float X, Y, px1, py1, px2, py2;
    716   float z1, z2, r, q1, q2, f;
    717 
    718   X = x - par[0];
    719   Y = y - par[1];
    720  
    721   px1 = par[2]*X;
    722   py1 = par[3]*Y;
    723   px2 = par[7]*X;
    724   py2 = par[8]*Y;
    725 
    726   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    727   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    728 
    729   r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,Npow));
    730   f = par[5]*r + par[6];
    731 
    732   q1 = par[5]*SQ(r)*(1 + z1);
    733   q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
    734 
    735   if (dpar != NULL) {
    736     dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    737     dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    738     dpar[2] = -2*q1*px1*X;
    739     dpar[3] = -2*q1*py1*Y;
    740     dpar[4] = -q1*X*Y;
    741     dpar[5] = +r;
    742     dpar[6] = +1;
    743     dpar[7] = -2*q2*px2*X;
    744     dpar[8] = -2*q2*py2*Y;
    745     dpar[9] = -q2*X*Y;
    746   }
    747   return (f);
    748 }
    749 
    750 /* two components: (1 + z_1^M + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
    751 float tgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    752 
    753   float X, Y, px1, py1, px2, py2;
    754   float z1, z2, r, q1, q2, f;
    755 
    756   X = x - par[0];
    757   Y = y - par[1];
    758  
    759   px1 = par[2]*X;
    760   py1 = par[3]*Y;
    761   px2 = par[7]*X;
    762   py2 = par[8]*Y;
    763 
    764   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    765   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    766 
    767   r = 1.0 / (1 + pow(z1,Npin) + pow(z2,Npow));
    768   f = par[5]*r + par[6];
    769 
    770   q1 = par[5]*SQ(r)*Npin*pow(z1,(Npin-1));
    771   q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
    772 
    773   if (dpar != NULL) {
    774     dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    775     dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    776     dpar[2] = -2*q1*px1*X*2;
    777     dpar[3] = -2*q1*py1*Y*2;
    778     dpar[4] = -q1*X*Y;
    779     dpar[5] = +r;
    780     dpar[6] = +1;
    781     dpar[7] = -2*q2*px2*X;
    782     dpar[8] = -2*q2*py2*Y;
    783     dpar[9] = -q2*X*Y;
    784   }
    785   return (f);
    786 }
    787 
    788 /* alternate for exp(-z) */
    789 /* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
    790 /* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
    791 
    792 # if (0)
    793 /* pars: x, y, sx, sy, sxy, I, sky */
    794 float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
    795 
    796   float X, Y, Z, E, F, q, R, f, p2, p3;
    797 
    798   X = x - par[0];
    799   Y = y - par[1];
    800  
    801   p2 = X / par[2];
    802   p3 = Y / par[3];
    803 
    804   Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
    805   E = 1.0 / (1 + Z);   
    806 
    807   q = par[5] * E;
    808   R = q*E;
    809   F = 0.5 / Z;
    810  
    811   f = q + par[6];
    812 
    813   dpar[0] = F*R*(p2 + par[4]*Y);
    814   dpar[1] = F*R*(p3 + par[4]*X);
    815   dpar[2] = F*0.5*R*p2*p2;
    816   dpar[3] = F*0.5*R*p3*p3;
    817   dpar[4] = -R*X*Y*F;
    818    
    819   dpar[5] = E;
    820   dpar[6] = 1;
    821   return (f);
    822 }
    823 
    824 /* pars: x, y, sx, sy, sxy, I, sky */
    825 float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
    826 
    827   float X, Y, Z, E, F, q, R, f, p2, p3;
    828 
    829   X = x - par[0];
    830   Y = y - par[1];
    831  
    832   p2 = X / par[2];
    833   p3 = Y / par[3];
    834 
    835   Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
    836   E = 1.0 / (1 + Z*Z*Z);   
    837 
    838   q = par[5] * E;
    839   F = 3*Z*Z;
    840   R = q*E*F;
    841  
    842   f = q + par[6];
    843 
    844   dpar[0] = R*(p2 + par[4]*Y);
    845   dpar[1] = R*(p3 + par[4]*X);
    846   dpar[2] = 0.5*R*p2*p2;
    847   dpar[3] = 0.5*R*p3*p3;
    848   dpar[4] = -R*X*Y;
    849    
    850   dpar[5] = E;
    851   dpar[6] = 1;
    852   return (f);
    853 }
    854 
    855 /* convert from x,y to major,minor */
    856 void fix_ellipsegauss_pars (float *par, int Npar) {
    857 
    858   float p2, p4, angle, t1, t2, tmp, area;
    859 
    860   /* par[0], par[1] = Xo, Yo - stay the same */
    861 
    862   p2 = 1/par[2];
    863   p4 = 1/par[3];
    864 
    865   angle = 0.5 * atan2 (-2*par[4], p4 - p2);
    866 
    867   tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
    868   t1 = (p2 + p4 + tmp) / 2;
    869   t2 = t1 - tmp;
    870 
    871   par[2] = 2.35482*sqrt(1/t2);
    872   par[3] = 2.35482*sqrt(1/t1);
    873   par[4] = DEG_RAD * angle;
    874 
    875   area = 2*M_PI/sqrt(t1*t2);
    876 
    877   par[5] *= area;
    878 
    879 }
    880 # endif
    881 
    882 /***  options for later
    883 
    884   Subtract = FALSE;
    885   if ((N = get_argument (argc, argv, "-sub"))) {
    886     remove_argument (N, &argc, argv);
    887     Subtract  = TRUE;
    888   }
    889 
    890   DFact = 1;
    891   if ((N = get_argument (argc, argv, "-D"))) {
    892     remove_argument (N, &argc, argv);
    893     DFact  = atof(argv[N]);
    894     remove_argument (N, &argc, argv);
    895   }
    896 
    897   fitfunc = fgaussTD;
    898   if ((N = get_argument (argc, argv, "-gal"))) {
    899     remove_argument (N, &argc, argv);
    900     fitfunc = fgalaxyTD;
    901   }
    902   if ((N = get_argument (argc, argv, "-bar"))) {
    903     remove_argument (N, &argc, argv);
    904     fitfunc = fbarTD;
    905   }
    906 
    907 
    908   f1 = 1;
    909   if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
    910 
    911   f2 = 1;
    912   if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
    913 
    914   if (Subtract) {
    915     tmpsky = par[6];
    916     par[6] = 0;
    917     for (N = j = 0; j < ny; j++) {
    918       V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    919       for (i = 0; i < nx; i++, V++, N++) {
    920         dx = i + sx;
    921         dy = j + sy;
    922         *V -= fitfunc (dx, dy, par, Npar, (float *) NULL);
    923       }
    924     }
    925     par[6] = tmpsky;
    926   }
    927 
    928 ***/
    929 
    930 # if (0)
    931 
    932 /* these two tests were not very succcessful.  the first did not model the shape well because
    933    it could not match the change in roundness with radius.  the second did not work because the
    934    parameters were degenerate (amplitude and slope of second component) */
    935 
    936 /* test: fixed, non-integer higher-order term -- x, y, sx, sy, sxy, I, sky, f1, f2 */
    937 float qgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    938 
    939   float X, Y, px, py;
    940   float z, r, q, f, k;
    941 
    942   X = x - par[0];
    943   Y = y - par[1];
    944  
    945   px = par[2]*X;
    946   py = par[3]*Y;
    947 
    948   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
    949   k = pow(z,1.75*par[8]);
    950   r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
    951   q = par[5]*r*r*(1 + 1.75*par[7]*par[8]*pow(z,1.75*par[8]-1));
    952   /* note difference from gaussian: q = par[5]*r */
    953   f = par[5]*r + par[6];
    954 
    955   dpar[0] = q*(2*px*par[2] + par[4]*Y);
    956   dpar[1] = q*(2*py*par[3] + par[4]*X);
    957   dpar[2] = -2*q*px*X;
    958   dpar[3] = -2*q*py*Y;
    959   dpar[4] = -q*X*Y;
    960   dpar[5] = +r;
    961   dpar[6] = +1;
    962   dpar[7] = -10*par[5]*r*r*k;
    963   dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z);
    964 
    965   return (f);
    966 }
    967 
    968 /* test: two component model: inner pseudo gaussian with outer z^1.75 x, y, sx, sy, sxy, I, sky */
    969 float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    970 
    971   float X, Y, px1, py1, px2, py2;
    972   float z1, z2, r1, r2, q1, q2, f;
    973 
    974   X = x - par[0];
    975   Y = y - par[1];
    976  
    977   px1 = par[2]*X;
    978   py1 = par[3]*Y;
    979   px2 = par[8]*X;
    980   py2 = par[9]*Y;
    981 
    982   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    983   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[10]*X*Y;
    984 
    985   r1 = 1.0 / (1 + z1 + 0.5*SQ(z1)*(1 + z1/3)); /* ~ exp (-Z) */
    986   r2 = 1.0 / (1 + pow(z2,1.75));
    987 
    988   f = par[5]*r1 + par[6] + par[7]*r2;
    989 
    990   q1 = par[5]*SQ(r1)*(1 + z1 + 0.5*SQ(z1));
    991   q2 = par[7]*SQ(r2)*(1.75*pow(z2,0.75));
    992 
    993   dpar[ 0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[8] + par[10]*Y);
    994   dpar[ 1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[9] + par[10]*X);
    995   dpar[ 2] = -2*q1*px1*X;
    996   dpar[ 3] = -2*q1*py1*Y;
    997   dpar[ 4] = -q1*X*Y;
    998   dpar[ 5] = +r1;
    999   dpar[ 6] = +1;
    1000   dpar[ 7] = +r2*2;
    1001   dpar[ 8] = -2*q2*px2*X*2;
    1002   dpar[ 9] = -2*q2*py2*Y*2;
    1003   dpar[10] = -q2*X*Y;
    1004 
    1005   return (f);
    1006 }
    1007 
    1008 # endif
    1009 
    1010   /* forcing values to have a rational range
    1011   ALLOCATE (parmin, float, Npar);
    1012   ALLOCATE (parmax, float, Npar);
    1013   bzero (parmin, Npar*sizeof(float));
    1014   bzero (parmax, Npar*sizeof(float));
    1015   parmin[0] = parmin[1] = 0;
    1016   parmax[0] = buf[0].matrix.Naxis[0];
    1017   parmax[1] = buf[0].matrix.Naxis[1];
    1018 
    1019   parmin[2] = parmin[3] = 0.01;
    1020   parmax[2] = parmax[3] = 100.0;
    1021   parmin[4] = -1000;
    1022   parmax[4] = -1000;
    1023  
    1024   parmin[5] = 1;
    1025   parmax[5] = 1e5;
    1026 
    1027   parmin[6] = 0.0;
    1028   parmax[6] = 1e5;
    1029 
    1030   if (Npar == 9) {
    1031     parmin[7] = parmin[8] = 0.01;
    1032     parmax[7] = parmax[8] = 10.0;
    1033   }
    1034   if (Npar == 10) {
    1035     parmin[7] = parmin[8] = 0.01;
    1036     parmax[7] = parmax[8] = 10.0;
    1037     parmin[9] = -1000;
    1038     parmax[9] = -1000;
    1039   }
    1040   */
    1041 

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

    r2843 r3144  
    1515int getvel                  PROTO((int, char **));
    1616int imfit                   PROTO((int, char **));
     17int imsub                   PROTO((int, char **));
    1718int medianmap               PROTO((int, char **));
    1819int mkgauss                 PROTO((int, char **));
    … …  
    4647  {"getvel",      getvel,       "rotcurve to velocities"},
    4748  {"imfit",       imfit,        "fit function"},
     49  {"imsub",       imsub,        "subtract function"},
    4850  {"medianmap",   medianmap,    "small median image"},
    4951  {"mkgauss",     mkgauss,      "generate a 2-D gaussian centered in image"},
Note: See TracChangeset for help on using the changeset viewer.