IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 20, 2010, 1:14:11 PM (16 years ago)
Author:
eugene
Message:

Many changes:

  • psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)

  • ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:
trunk/psModules
Files:
7 edited
4 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/objects/models/pmModel_GAUSS.c

    r26916 r29004  
    2121#include <stdio.h>
    2222#include <pslib.h>
    23 
     23#include "pmHDU.h"
     24#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
     29#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
     31#include "pmFootprint.h"
     32#include "pmPeaks.h"
    2433#include "pmMoments.h"
    25 #include "pmPeaks.h"
     34#include "pmModelFuncs.h"
     35#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    2641#include "pmSource.h"
    27 #include "pmModel.h"
     42#include "pmSourceFitModel.h"
     43#include "pmPSF.h"
     44#include "pmPSFtry.h"
     45#include "pmDetections.h"
     46
    2847#include "pmModel_GAUSS.h"
    2948
     49# define PM_MODEL_NPARAM          7
    3050# define PM_MODEL_FUNC            pmModelFunc_GAUSS
    3151# define PM_MODEL_FLUX            pmModelFlux_GAUSS
    … …  
    83103        dPAR[PM_PAR_XPOS] = q*(2*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    84104        dPAR[PM_PAR_YPOS] = q*(2*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     105
    85106        // the extra factor of 2 below is needed to avoid excessive swings
    86107        dPAR[PM_PAR_SXX]  = +4.0*q*px*px/PAR[PM_PAR_SXX];
    … …  
    102123        return true;
    103124    }
    104     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     125    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    105126
    106127    // we need to calculate the limits for SXY specially
    … …  
    347368// this test is invalid if the parameters are derived
    348369// from the PSF model
     370// XXX how is this used?  it prevents forced photometry from ever being 'successful'
    349371bool PM_MODEL_FIT_STATUS (pmModel *model)
    350372{
    … …  
    394416    return;
    395417}
    396 
    397 # undef PM_MODEL_FUNC
    398 # undef PM_MODEL_FLUX
    399 # undef PM_MODEL_GUESS
    400 # undef PM_MODEL_LIMITS
    401 # undef PM_MODEL_RADIUS
    402 # undef PM_MODEL_FROM_PSF
    403 # undef PM_MODEL_PARAMS_FROM_PSF
    404 # undef PM_MODEL_FIT_STATUS
    405 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_PGAUSS.c

    r27565 r29004  
    2121#include <stdio.h>
    2222#include <pslib.h>
    23 
     23#include "pmHDU.h"
     24#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
     29#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
     31#include "pmFootprint.h"
     32#include "pmPeaks.h"
    2433#include "pmMoments.h"
    25 #include "pmPeaks.h"
     34#include "pmModelFuncs.h"
     35#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    2641#include "pmSource.h"
    27 #include "pmModel.h"
     42#include "pmSourceFitModel.h"
     43#include "pmPSF.h"
     44#include "pmPSFtry.h"
     45#include "pmDetections.h"
     46
    2847#include "pmModel_PGAUSS.h"
    2948
     49# define PM_MODEL_NPARAM          7
    3050# define PM_MODEL_FUNC            pmModelFunc_PGAUSS
    3151# define PM_MODEL_FLUX            pmModelFlux_PGAUSS
    … …  
    103123        return true;
    104124    }
    105     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     125    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    106126
    107127    // we need to calculate the limits for SXY specially
    … …  
    448468    return;
    449469}
    450 
    451 # undef PM_MODEL_FUNC
    452 # undef PM_MODEL_FLUX
    453 # undef PM_MODEL_GUESS
    454 # undef PM_MODEL_LIMITS
    455 # undef PM_MODEL_RADIUS
    456 # undef PM_MODEL_FROM_PSF
    457 # undef PM_MODEL_PARAMS_FROM_PSF
    458 # undef PM_MODEL_FIT_STATUS
    459 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_PS1_V1.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_PS1_V1.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_PS1_V1
    3252# define PM_MODEL_FLUX            pmModelFlux_PS1_V1
    … …  
    122142        return true;
    123143    }
    124     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     144    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    125145
    126146    // we need to calculate the limits for SXY specially
    … …  
    222242    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
    223243    PAR[PM_PAR_SXY]  = shape.sxy;
    224     PAR[PM_PAR_7]    = 1.0;
     244    PAR[PM_PAR_7]    = 0.5;
    225245
    226246    return(true);
    … …  
    468488    return;
    469489}
    470 
    471 # undef PM_MODEL_FUNC
    472 # undef PM_MODEL_FLUX
    473 # undef PM_MODEL_GUESS
    474 # undef PM_MODEL_LIMITS
    475 # undef PM_MODEL_RADIUS
    476 # undef PM_MODEL_FROM_PSF
    477 # undef PM_MODEL_PARAMS_FROM_PSF
    478 # undef PM_MODEL_FIT_STATUS
    479 # undef PM_MODEL_SET_LIMITS
    480 # undef ALPHA
    481 # undef ALPHA_M

  • trunk/psModules/src/objects/models/pmModel_QGAUSS.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_QGAUSS.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_QGAUSS
    3252# define PM_MODEL_FLUX            pmModelFlux_QGAUSS
    … …  
    123143        return true;
    124144    }
    125     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     145    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    126146
    127147    // we need to calculate the limits for SXY specially
    … …  
    469489    return;
    470490}
    471 
    472 # undef PM_MODEL_FUNC
    473 # undef PM_MODEL_FLUX
    474 # undef PM_MODEL_GUESS
    475 # undef PM_MODEL_LIMITS
    476 # undef PM_MODEL_RADIUS
    477 # undef PM_MODEL_FROM_PSF
    478 # undef PM_MODEL_PARAMS_FROM_PSF
    479 # undef PM_MODEL_FIT_STATUS
    480 # undef PM_MODEL_SET_LIMITS
    481 # undef ALPHA
    482 # undef ALPHA_M

  • trunk/psModules/src/objects/models/pmModel_RGAUSS.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_RGAUSS.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_RGAUSS
    3252# define PM_MODEL_FLUX            pmModelFlux_RGAUSS
    … …  
    112132        return true;
    113133    }
    114     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     134    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    115135
    116136    // we need to calculate the limits for SXY specially
    … …  
    209229    PAR[PM_PAR_XPOS] = peak->xf;
    210230    PAR[PM_PAR_YPOS] = peak->yf;
    211     PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
    212     PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
     231    PAR[PM_PAR_SXX]  = PS_MAX(0.5, shape.sx);
     232    PAR[PM_PAR_SYY]  = PS_MAX(0.5, shape.sy);
    213233    PAR[PM_PAR_SXY]  = shape.sxy;
    214     PAR[PM_PAR_7]    = 2.25;
     234    PAR[PM_PAR_7]    = 1.5;
    215235
    216236    return(true);
    … …  
    463483    return;
    464484}
    465 
    466 # undef PM_MODEL_FUNC
    467 # undef PM_MODEL_FLUX
    468 # undef PM_MODEL_GUESS
    469 # undef PM_MODEL_LIMITS
    470 # undef PM_MODEL_RADIUS
    471 # undef PM_MODEL_FROM_PSF
    472 # undef PM_MODEL_PARAMS_FROM_PSF
    473 # undef PM_MODEL_FIT_STATUS
    474 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_SERSIC.c

    r26916 r29004  
    1818   * PM_PAR_7   7   - normalized sersic parameter
    1919
    20    note that a standard sersic model uses exp(-K*(z^(1/n) - 1).  the additional elements (K,
     20   note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    2121   the -1 offset) are absorbed in this model by the normalization, the exponent, and the
    2222   radial scale.  We fit the elements in this form, then re-normalize them on output.
    … …  
    2525#include <stdio.h>
    2626#include <pslib.h>
    27 
     27#include "pmHDU.h"
     28#include "pmFPA.h"
     29
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
     33#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
     35#include "pmFootprint.h"
     36#include "pmPeaks.h"
    2837#include "pmMoments.h"
    29 #include "pmPeaks.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
    3045#include "pmSource.h"
    31 #include "pmModel.h"
     46#include "pmSourceFitModel.h"
     47#include "pmPSF.h"
     48#include "pmPSFtry.h"
     49#include "pmDetections.h"
     50
    3251#include "pmModel_SERSIC.h"
    3352
     53# define PM_MODEL_NPARAM          8
    3454# define PM_MODEL_FUNC            pmModelFunc_SERSIC
    3555# define PM_MODEL_FLUX            pmModelFlux_SERSIC
    … …  
    4767
    4868// Lax parameter limits
    49 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
     69static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
    5070static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
    5171
    … …  
    84104    assert (z >= 0);
    85105
    86     psF32 f2 = pow(z,PAR[PM_PAR_7]);
    87     psF32 f1 = exp(-f2);
     106    float index = 0.5 / PAR[PM_PAR_7];
     107    float bn = 1.9992*index - 0.3271;
     108    float Io = exp(bn);
     109
     110    psF32 f2 = bn*pow(z,PAR[PM_PAR_7]);
     111    psF32 f1 = Io*exp(-f2);
    88112    psF32 z0 = PAR[PM_PAR_I0]*f1;
    89113    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    98122
    99123        // gradient is infinite for z = 0; saturate at z = 0.01
    100         psF32 z1 = (z < 0.01) ? z0*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
     124        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
    101125
    102126        dPAR[PM_PAR_SKY]  = +1.0;
    103127        dPAR[PM_PAR_I0]   = +f1;
    104         dPAR[PM_PAR_7]    = (z == 0.0) ? 0.0 : -z0*f2*log(z);
     128        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
     129        dPAR[PM_PAR_7]   *= 3.0;
    105130
    106131        assert (isfinite(z1));
    … …  
    109134        dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    110135        dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    111         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     136        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    112137        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    113         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    114138        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    115139    }
    … …  
    127151        return true;
    128152    }
    129     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     153    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    130154
    131155    // we need to calculate the limits for SXY specially
    … …  
    201225    psF32     *PAR  = model->params->data.F32;
    202226
     227    // the other parameters depend on the guess for PAR_7
     228    if (!isfinite(PAR[PM_PAR_7])) {
     229        PAR[PM_PAR_7]    = 0.25;
     230    }   
     231    float index = 0.5 / PAR[PM_PAR_7];
     232
     233    // the scale-length is a function of the moments and the index:
     234    // Rmajor_guess = Rmajor_moments * Scale * Zero
     235    float Scale = 0.70 + 0.053 * PAR[PM_PAR_7];
     236    float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
     237
    203238    psEllipseMoments emoments;
    204239    emoments.x2 = moments->Mxx;
    … …  
    213248    if (!isfinite(axes.theta)) return false;
    214249
     250    // set a lower limit to avoid absurd solutions..
     251    float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero);
     252    float Rminor = Rmajor * (axes.minor / axes.major);
     253
     254    fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor);
     255
     256    axes.major = Rmajor;
     257    axes.minor = Rminor;
    215258    psEllipseShape shape = psEllipseAxesToShape (axes);
    216259
    … …  
    219262    if (!isfinite(shape.sxy)) return false;
    220263
     264    float bn = 1.9992*index - 0.3271;
     265    // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
     266    float Io = exp(0.5*bn);
     267
     268    // XXX do we need this factor of sqrt(2)?
     269    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
     270    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
     271
     272    float Sxx = PS_MAX(0.5, shape.sx);
     273    float Syy = PS_MAX(0.5, shape.sy);
     274
    221275    PAR[PM_PAR_SKY]  = 0.0;
    222     PAR[PM_PAR_I0]   = peak->flux;
     276    PAR[PM_PAR_I0]   = peak->flux / Io;
    223277    PAR[PM_PAR_XPOS] = peak->xf;
    224278    PAR[PM_PAR_YPOS] = peak->yf;
    225     PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
    226     PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
     279    PAR[PM_PAR_SXX]  = Sxx;
     280    PAR[PM_PAR_SYY]  = Syy;
    227281    PAR[PM_PAR_SXY]  = shape.sxy;
    228     PAR[PM_PAR_7]    = 0.5;
    229282
    230283    return(true);
    … …  
    254307    float f1, f2;
    255308    for (z = DZ; z < 50; z += DZ) {
    256         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     309        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     310        f1 = exp(-pow(z,PAR[PM_PAR_7]));
    257311        z += DZ;
    258         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     312        f2 = exp(-pow(z,PAR[PM_PAR_7]));
    259313        norm += f0 + 4*f1 + f2;
    260314        f0 = f2;
    … …  
    287341
    288342    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    289     psF64 sigma = axes.major;
    290 
    291     psF64 limit = flux / PAR[PM_PAR_I0];
    292 
    293     psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7]));
    294     psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
    295 
    296     psF64 radius = sigma * sqrt (2.0 * z);
    297     psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma);
    298 
    299     if (isnan(radius))
    300         psAbort("error in code: radius is NaN");
    301 
     343
     344    // f = Io exp(-z^n) -> z^n = ln(Io/f)
     345    psF64 zn = log(PAR[PM_PAR_I0] / flux);
     346    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     347
     348    fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
     349
     350    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
    302351    return (radius);
    303352}
    … …  
    448497    return;
    449498}
    450 
    451 # undef PM_MODEL_FUNC
    452 # undef PM_MODEL_FLUX
    453 # undef PM_MODEL_GUESS
    454 # undef PM_MODEL_LIMITS
    455 # undef PM_MODEL_RADIUS
    456 # undef PM_MODEL_FROM_PSF
    457 # undef PM_MODEL_PARAMS_FROM_PSF
    458 # undef PM_MODEL_FIT_STATUS
    459 # undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.