IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36085 for trunk/psModules


Ignore:
Timestamp:
Aug 31, 2013, 5:55:16 AM (13 years ago)
Author:
eugene
Message:

merge changes from EAM branch (pmModel_CentralPixel functions for sersic-like model support; add interactive mode to PCM LMM fitting; harder PCM LMM damping (nu scaled by 3 not 2); add EXT_AND SKY and SHAPE fitting modes, include sky in INDEX-only fits; ifdefs for trace-only values; allow positions to go to 100,000; major re-work of central pixel & flux for sersic-like models; add pmPCMMakeModel function to generate the flux for an arbitrary model

Location:
trunk/psModules
Files:
23 edited
4 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules

  • trunk/psModules/src/detrend/pmPattern.c

    r35270 r36085  
    12221222    }
    12231223    double max_XX = 0;
     1224#if (PS_TRACE_ON)
    12241225    double solution_V = 0;
    12251226    int i_peak = -1;
     1227#endif
    12261228    for (int i = 0; i < numChips + 4; i++) { // If any cells have no value of themself, set the matrix to 1.0.
    12271229      if (XX->data.F64[i][i] == 0.0) {
    … …  
    12301232      if (XX->data.F64[i][i] > max_XX) {
    12311233        max_XX = XX->data.F64[i][i];
     1234#if (PS_TRACE_ON)
    12321235        solution_V = solution->data.F64[i];
    12331236        i_peak = i;
     1237#endif
    12341238      }
    12351239    }

  • trunk/psModules/src/objects/Makefile.am

    r34823 r36085  
    2020        pmModelClass.c \
    2121        pmModelUtils.c \
     22        pmModel_CentralPixel.c \
    2223        pmSource.c \
    2324        pmPhotObj.c \
    … …  
    9798        pmModelClass.h \
    9899        pmModelUtils.h \
     100        pmModel_CentralPixel.h \
    99101        pmSource.h \
    100102        pmPhotObj.h \
    … …  
    110112        pmSourceOutputs.h \
    111113        pmSourceIO.h \
    112         pmSourceSatstar.h \ 
     114        pmSourceSatstar.h \
    113115        pmSourcePlots.h \
    114116        pmSourceVisual.h \

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

    r35768 r36085  
    1616   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
    1717   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    18    * PM_PAR_7   7   - normalized dev parameter
    1918
    2019   note that a standard dev model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    … …  
    4948#include "pmPSFtry.h"
    5049#include "pmDetections.h"
     50#include "pmModel_CentralPixel.h"
    5151
    5252#include "pmModel_DEV.h"
    … …  
    6363# define PM_MODEL_SET_LIMITS      pmModelSetLimits_DEV
    6464
    65 // f = exp(-z^0.125)
     65// f = exp(-kappa*r^(1/index))
     66// f = exp(-kappa*z^(0.5/index))
     67// index = 4, 0.5/index = 0.125
    6668# define ALPHA 0.125
    67 // # define ALPHA 0.25
    6869
    6970// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    … …  
    7374// Lax parameter limits
    7475static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0 };
    75 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     76static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    7677
    7778// Moderate parameter limits
    … …  
    8687static float *paramsMinUse = paramsMinLax;
    8788static float *paramsMaxUse = paramsMaxLax;
    88 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     89static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 };
    8990
    9091static bool limitsApply = true;         // Apply limits?
    91 
    92 # include "pmModel_SERSIC.CP.h"
    9392
    9493psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    109108    psAssert (z >= 0, "do not allow negative z values in model");
    110109
    111     float index = 0.5 / ALPHA;
    112     float par7 = ALPHA;
    113     float bn = 1.9992*index - 0.3271;
    114     float Io = exp(bn);
    115 
    116     psF32 f2 = bn*pow(z,ALPHA);
    117     psF32 f1 = Io*exp(-f2);
    118 
    119     psF32 radius = hypot(X, Y);
    120     if (radius < 1.0) {
    121 
    122         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    123 
    124         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    125         psEllipseAxes axes;
    126         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    127 
    128         // get the central pixel flux from the lookup table
    129         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    130         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    131         float yPix = (index - centralPixelYo) / centralPixeldY;
    132         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    133 
    134         // the integral of a Sersic has an analytical form as follows:
    135         float logGamma = lgamma(2.0*index);
    136         float bnFactor = pow(bn, 2.0*index);
    137         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    138 
    139         // XXX interpolate to get the value
    140         // XXX for the moment, just integerize
    141         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    142         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    143        
    144         float px1 = 1.0 / PAR[PM_PAR_SXX];
    145         float py1 = 1.0 / PAR[PM_PAR_SYY];
    146         float z10 = PS_SQR(px1);
    147         float z01 = PS_SQR(py1);
    148 
    149         // which pixels do we need for this interpolation?
    150         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    151         if ((X >= 0) && (Y >= 0)) {
    152             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    153             float V00 = Vcenter;
    154             float V10 = Io*exp(-bn*pow(z10,par7));
    155             float V01 = Io*exp(-bn*pow(z01,par7));
    156             float V11 = Io*exp(-bn*pow(z11,par7));
    157             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    158         }
    159         if ((X < 0) && (Y >= 0)) {
    160             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    161             float V00 = Io*exp(-bn*pow(z10,par7));
    162             float V10 = Vcenter;
    163             float V01 = Io*exp(-bn*pow(z11,par7));
    164             float V11 = Io*exp(-bn*pow(z01,par7));
    165             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    166         }
    167         if ((X >= 0) && (Y < 0)) {
    168             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    169             float V00 = Io*exp(-bn*pow(z01,par7));
    170             float V10 = Io*exp(-bn*pow(z11,par7));
    171             float V01 = Vcenter;
    172             float V11 = Io*exp(-bn*pow(z10,par7));
    173             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    174         }
    175         if ((X < 0) && (Y < 0)) {
    176             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    177             float V00 = Io*exp(-bn*pow(z11,par7));
    178             float V10 = Io*exp(-bn*pow(z10,par7));
    179             float V01 = Io*exp(-bn*pow(z01,par7));
    180             float V11 = Vcenter;
    181             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    182         }
     110    // for DEV, we can hard-wire kappa(4):
     111    // float index = 4.0;
     112    float kappa = 7.670628;
     113
     114    // r = sqrt(z)
     115    float q = kappa*pow(z,ALPHA);
     116    float f0 = exp(-q);
     117
     118    assert (isfinite(q));
     119    assert (isfinite(f0));
     120
     121    // only worry about the central pixels at most
     122    float radius = hypot(X, Y);
     123    if (radius <= 1.5) {
     124        // Nsub ~ 10*index^2 + 1
     125        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV"));
     126        int Nsub = 2 * ((int)(25 / axes.minor)) + 1;
     127        Nsub = PS_MIN (Nsub, 121);
     128        Nsub = PS_MAX (Nsub, 11);
     129        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
    183130    }   
    184131
    185     psF32 z0 = PAR[PM_PAR_I0]*f1;
    186     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    187 
    188     assert (isfinite(f2));
     132    float f1 = PAR[PM_PAR_I0]*f0;
     133    float f = PAR[PM_PAR_SKY] + f1;
     134
    189135    assert (isfinite(f1));
    190     assert (isfinite(z0));
    191     assert (isfinite(f0));
     136    assert (isfinite(f));
    192137
    193138    if (deriv != NULL) {
    … …  
    195140
    196141        dPAR[PM_PAR_SKY]  = +1.0;
    197         dPAR[PM_PAR_I0]   = +2.0*f1; // XXX extra damping..
    198 
    199         // gradient is infinite for z = 0; saturate at z = 0.01
    200         psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0);
    201 
    202         assert (isfinite(z1));
    203 
    204         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    205         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    206         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
    207         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    208         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    209     }
    210     return (f0);
     142        dPAR[PM_PAR_I0]   = +f0;
     143
     144        if (z > 0.01) {
     145          float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0);
     146          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     147          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     148          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     149          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     150          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     151        } else {
     152          // gradient -> 0 for z -> 0, but has undef form
     153          float z1 = f1*kappa*ALPHA*pow(z,ALPHA);
     154          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     155          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     156          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     157          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     158          dPAR[PM_PAR_SXY]  = -1.0*z1;
     159        }
     160    }
     161    return (f);
    211162}
    212163
    … …  
    302253    }
    303254
    304     // the normalization is modified by the slope
    305     float index = 0.5 / ALPHA;
    306     float bn = 1.9992*index - 0.3271;
    307     float Io = exp(0.5*bn);
    308 
    309255    // set the model normalization
    310256    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    311257      return false;
    312258    }
    313     PAR[PM_PAR_I0] /= Io;
    314259
    315260    // set the model position
    … …  
    328273    psEllipseAxes axes;
    329274    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    330     float AspectRatio = axes.minor / axes.major;
    331 
    332     float index = 4.0;
    333     float bn = 1.9992*index - 0.3271;
    334 
    335     // the integral of a Sersic has an analytical form as follows:
    336     float logGamma = lgamma(2.0*index);
    337     float bnFactor = pow(bn, 2.0*index);
    338     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    339    
    340     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    341 
    342     return(Flux);
     275
     276    float norm = 0.00168012;
     277    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     278
     279    return(flux);
    343280}
    344281
    … …  
    359296    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    360297
    361     // f = Io exp(-z^n) -> z^n = ln(Io/f)
    362     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    363     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
     298    // static value for DEV:
     299    float kappa = 7.670628;
     300
     301    // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
     302    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     303    psF64 radius = axes.major * pow(zn, 0.5 / ALPHA);
    364304
    365305    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",

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

    r35768 r36085  
    4545#include "pmPSFtry.h"
    4646#include "pmDetections.h"
     47#include "pmModel_CentralPixel.h"
    4748
    4849#include "pmModel_EXP.h"
    … …  
    6566// Lax parameter limits
    6667static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 };
    67 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     68static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    6869
    6970// Moderate parameter limits
    … …  
    7879static float *paramsMinUse = paramsMinLax;
    7980static float *paramsMaxUse = paramsMaxLax;
    80 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5};
     81static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5};
    8182
    8283static bool limitsApply = true;         // Apply limits?
    8384
    84 # include "pmModel_SERSIC.CP.h"
     85// # include "pmModel_SERSIC.CP.h"
     86
     87// the problems I'm having with the SERSIC-like functions are:
     88// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius)
     89// 2) getting the central pixel right
     90// 3) getting the derivaties right.
    8591
    8692psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    101107    psAssert (z >= 0, "do not allow negative z values in model");
    102108
    103     float index = 1.0;
    104     float par7 = 0.5;
    105     float bn = 1.9992*index - 0.3271;
    106     float Io = exp(bn);
    107 
    108     psF32 f2 = bn*sqrt(z);
    109     psF32 f1 = Io*exp(-f2);
    110 
     109    // for EXP, we can hard-wire kappa(1):
     110    // float index = 1.0;
     111    float kappa = 1.70056;
     112
     113    // sqrt(z) is r
     114    float q = kappa*sqrt(z);
     115    psF32 f0 = exp(-q);
     116
     117    assert (isfinite(q));
     118
     119    // only worry about the central 4 pixels at most
    111120    psF32 radius = hypot(X, Y);
    112     if (radius < 1.0) {
    113 
    114         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    115 
    116         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    117         psEllipseAxes axes;
    118         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    119 
    120         // get the central pixel flux from the lookup table
    121         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    122         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    123         float yPix = (index - centralPixelYo) / centralPixeldY;
    124         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    125 
    126         // the integral of a Sersic has an analytical form as follows:
    127         float logGamma = lgamma(2.0*index);
    128         float bnFactor = pow(bn, 2.0*index);
    129         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    130 
    131         // XXX interpolate to get the value
    132         // XXX for the moment, just integerize
    133         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    134         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    135        
    136         float px1 = 1.0 / PAR[PM_PAR_SXX];
    137         float py1 = 1.0 / PAR[PM_PAR_SYY];
    138         float z10 = PS_SQR(px1);
    139         float z01 = PS_SQR(py1);
    140 
    141         // which pixels do we need for this interpolation?
    142         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    143         if ((X >= 0) && (Y >= 0)) {
    144             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    145             float V00 = Vcenter;
    146             float V10 = Io*exp(-bn*pow(z10,par7));
    147             float V01 = Io*exp(-bn*pow(z01,par7));
    148             float V11 = Io*exp(-bn*pow(z11,par7));
    149             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    150         }
    151         if ((X < 0) && (Y >= 0)) {
    152             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    153             float V00 = Io*exp(-bn*pow(z10,par7));
    154             float V10 = Vcenter;
    155             float V01 = Io*exp(-bn*pow(z11,par7));
    156             float V11 = Io*exp(-bn*pow(z01,par7));
    157             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    158         }
    159         if ((X >= 0) && (Y < 0)) {
    160             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    161             float V00 = Io*exp(-bn*pow(z01,par7));
    162             float V10 = Io*exp(-bn*pow(z11,par7));
    163             float V01 = Vcenter;
    164             float V11 = Io*exp(-bn*pow(z10,par7));
    165             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    166         }
    167         if ((X < 0) && (Y < 0)) {
    168             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    169             float V00 = Io*exp(-bn*pow(z11,par7));
    170             float V10 = Io*exp(-bn*pow(z10,par7));
    171             float V01 = Io*exp(-bn*pow(z01,par7));
    172             float V11 = Vcenter;
    173             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    174         }
    175     }   
    176 
    177     psF32 z0 = PAR[PM_PAR_I0]*f1;
    178     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    179 
    180     assert (isfinite(f2));
     121    if (radius <= 1.5) {
     122        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
     123    }
     124    assert (isfinite(f0));
     125
     126    psF32 f1 = PAR[PM_PAR_I0]*f0;
     127    psF32 f = PAR[PM_PAR_SKY] + f1;
     128
    181129    assert (isfinite(f1));
    182     assert (isfinite(z0));
    183     assert (isfinite(f0));
     130    assert (isfinite(f));
    184131
    185132    if (deriv != NULL) {
    … …  
    187134
    188135        dPAR[PM_PAR_SKY]  = +1.0;
    189         dPAR[PM_PAR_I0]   = +f1;
    190 
    191         // gradient is infinite for z = 0; saturate at z = 0.01
    192         // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
    193         psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
    194 
    195         // XXX dampen SXX and SYY as in GAUSS?
    196         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    197         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    198         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
    199         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    200         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    201     }
    202     return (f0);
     136        dPAR[PM_PAR_I0]   = +f0;
     137
     138        if (z > 0.01) {
     139          float z1 = 0.5*f1*kappa/sqrt(z);
     140          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     141          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     142          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     143          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     144          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     145        } else {
     146          // gradient -> 0 for z -> 0, but has undef form
     147          float z1 = 0.5*f1*kappa;
     148          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     149          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     150          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     151          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     152          dPAR[PM_PAR_SXY]  = -1.0*z1;
     153        }
     154    }
     155    return (f);
    203156}
    204157
    … …  
    314267    psEllipseAxes axes;
    315268    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    316     float AspectRatio = axes.minor / axes.major;
    317 
    318     float index = 1.0;
    319     float bn = 1.9992*index - 0.3271;
    320 
    321     // the integral of a Sersic has an analytical form as follows:
    322     float logGamma = lgamma(2.0*index);
    323     float bnFactor = pow(bn, 2.0*index);
    324     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    325    
    326     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    327 
    328     return(Flux);
     269
     270    // static value for EXP:
     271    float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr
     272
     273    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     274
     275    return(flux);
    329276}
    330277
    … …  
    345292    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    346293
    347     // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
    348     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    349     psF64 radius = axes.major * sqrt (2.0) * zn;
     294    // static value for EXP:
     295    float kappa = 1.70056;
     296
     297    // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
     298    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     299    psF64 radius = axes.major * zn;
    350300
    351301    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",
    … …  
    501451    return;
    502452}
     453
     454# if (0)
     455void bilin_inter_function () {
     456        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     457        psEllipseAxes axes;
     458        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     459
     460        // get the central pixel flux from the lookup table
     461        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     462        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     463        float yPix = (index - centralPixelYo) / centralPixeldY;
     464        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     465
     466        // the integral of a Sersic has an analytical form as follows:
     467        float logGamma = lgamma(2.0*index);
     468        float bnFactor = pow(bn, 2.0*index);
     469        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     470
     471        // XXX interpolate to get the value
     472        // XXX for the moment, just integerize
     473        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     474        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     475       
     476        float px1 = 1.0 / PAR[PM_PAR_SXX];
     477        float py1 = 1.0 / PAR[PM_PAR_SYY];
     478        float z10 = PS_SQR(px1);
     479        float z01 = PS_SQR(py1);
     480
     481        // which pixels do we need for this interpolation?
     482        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     483        if ((X >= 0) && (Y >= 0)) {
     484            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     485            float V00 = Vcenter;
     486            float V10 = Io*exp(-bn*pow(z10,par7));
     487            float V01 = Io*exp(-bn*pow(z01,par7));
     488            float V11 = Io*exp(-bn*pow(z11,par7));
     489            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     490        }
     491        if ((X < 0) && (Y >= 0)) {
     492            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     493            float V00 = Io*exp(-bn*pow(z10,par7));
     494            float V10 = Vcenter;
     495            float V01 = Io*exp(-bn*pow(z11,par7));
     496            float V11 = Io*exp(-bn*pow(z01,par7));
     497            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     498        }
     499        if ((X >= 0) && (Y < 0)) {
     500            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     501            float V00 = Io*exp(-bn*pow(z01,par7));
     502            float V10 = Io*exp(-bn*pow(z11,par7));
     503            float V01 = Vcenter;
     504            float V11 = Io*exp(-bn*pow(z10,par7));
     505            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     506        }
     507        if ((X < 0) && (Y < 0)) {
     508            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     509            float V00 = Io*exp(-bn*pow(z11,par7));
     510            float V10 = Io*exp(-bn*pow(z10,par7));
     511            float V01 = Io*exp(-bn*pow(z01,par7));
     512            float V11 = Vcenter;
     513            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     514        }
     515}
     516# endif

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

    r35768 r36085  
    6161// Lax parameter limits
    6262static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 };
    63 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     63static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    6464
    6565// Moderate parameter limits

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

    r35768 r36085  
    6161// Lax parameter limits
    6262static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 };
    63 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     63static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    6464
    6565// Moderate parameter limits

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

    r35768 r36085  
    7070// Lax parameter limits
    7171static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
    72 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     72static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    7373
    7474// Moderate parameter limits
    7575// Tolerate a small divot (k < 0)
    7676static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
    77 static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     77static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    7878
    7979// Strict parameter limits
    8080// k = PAR_7 < 0 is very undesirable (big divot in the middle)
    8181static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
    82 static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     82static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    8383
    8484// Parameter limits to use

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

    r35768 r36085  
    7070// Lax parameter limits
    7171static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
    72 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     72static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    7373
    7474// Moderate parameter limits
    7575// Tolerate a small divot (k < 0)
    7676static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
    77 static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     77static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    7878
    7979// Strict parameter limits
    8080// k = PAR_7 < 0 is very undesirable (big divot in the middle)
    8181static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
    82 static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     82static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
    8383
    8484// Parameter limits to use

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

    r35768 r36085  
    6666// Lax parameter limits
    6767static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
    68 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     68static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 };
    6969
    7070// Moderate parameter limits

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

    r35768 r36085  
    2020   * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
    2121     construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
    22      Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
     22     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using:
    2323     shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
    2424     psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
    … …  
    5555#include "pmPSFtry.h"
    5656#include "pmDetections.h"
     57#include "pmModel_CentralPixel.h"
    5758
    5859#include "pmModel_SERSIC.h"
    … …  
    7475
    7576// Lax parameter limits
    76 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
    77 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     77static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 };
     78static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 };
    7879
    7980// Moderate parameter limits
    … …  
    8889static float *paramsMinUse = paramsMinLax;
    8990static float *paramsMaxUse = paramsMaxLax;
    90 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     91static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0};
    9192
    9293static bool limitsApply = true;         // Apply limits?
    9394
    94 # include "pmModel_SERSIC.CP.h"
     95// # include "pmModel_SERSIC.CP.h"
    9596
    9697psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    111112    psAssert (z >= 0, "do not allow negative z values in model");
    112113
    113     float index = 0.5 / PAR[PM_PAR_7];
    114     float par7 = PAR[PM_PAR_7];
    115     float bn = 1.9992*index - 0.3271;
    116     float Io = exp(bn);
    117 
    118     psF32 f2 = bn*pow(z,par7);
    119     psF32 f1 = Io*exp(-f2);
     114    float Sindex = 0.5 / PAR[PM_PAR_7];
     115    float kappa = pmSersicKappa (Sindex);
     116
     117    float q = kappa*pow(z,PAR[PM_PAR_7]);
     118    psF32 f0 = exp(-q);
     119
     120    assert (isfinite(q));
    120121
    121122    psF32 radius = hypot(X, Y);
    122     if (radius < 1.0) {
    123 
    124         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    125 
    126         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    127         psEllipseAxes axes;
    128         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    129 
    130         // get the central pixel flux from the lookup table
    131         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    132         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    133         float yPix = (index - centralPixelYo) / centralPixeldY;
    134         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    135 
    136         // the integral of a Sersic has an analytical form as follows:
    137         float logGamma = lgamma(2.0*index);
    138         float bnFactor = pow(bn, 2.0*index);
    139         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    140 
    141         // XXX interpolate to get the value
    142         // XXX for the moment, just integerize
    143         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    144         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    145        
    146         float px1 = 1.0 / PAR[PM_PAR_SXX];
    147         float py1 = 1.0 / PAR[PM_PAR_SYY];
    148         float z10 = PS_SQR(px1);
    149         float z01 = PS_SQR(py1);
    150 
    151         // which pixels do we need for this interpolation?
    152         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    153         if ((X >= 0) && (Y >= 0)) {
    154             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    155             float V00 = Vcenter;
    156             float V10 = Io*exp(-bn*pow(z10,par7));
    157             float V01 = Io*exp(-bn*pow(z01,par7));
    158             float V11 = Io*exp(-bn*pow(z11,par7));
    159             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    160         }
    161         if ((X < 0) && (Y >= 0)) {
    162             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    163             float V00 = Io*exp(-bn*pow(z10,par7));
    164             float V10 = Vcenter;
    165             float V01 = Io*exp(-bn*pow(z11,par7));
    166             float V11 = Io*exp(-bn*pow(z01,par7));
    167             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    168         }
    169         if ((X >= 0) && (Y < 0)) {
    170             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    171             float V00 = Io*exp(-bn*pow(z01,par7));
    172             float V10 = Io*exp(-bn*pow(z11,par7));
    173             float V01 = Vcenter;
    174             float V11 = Io*exp(-bn*pow(z10,par7));
    175             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    176         }
    177         if ((X < 0) && (Y < 0)) {
    178             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    179             float V00 = Io*exp(-bn*pow(z11,par7));
    180             float V10 = Io*exp(-bn*pow(z10,par7));
    181             float V01 = Io*exp(-bn*pow(z01,par7));
    182             float V11 = Vcenter;
    183             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    184         }
    185     }   
    186 
    187     psF32 z0 = PAR[PM_PAR_I0]*f1;
    188     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    189 
    190     if (!isfinite(z0)) {
    191         fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
     123    if (radius <= 1.5) {
     124        // Nsub ~ 10*index^2 + 1
     125        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC"));
     126        int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1;
     127        Nsub = PS_MIN (Nsub, 121);
     128        Nsub = PS_MAX (Nsub, 11);
     129        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub);
     130    }
     131    if (!isfinite(f0)) {
     132        fprintf(stderr, "f0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, q);
    192133        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
    193134            PAR[5], PAR[6], PAR[7]);
    194135    }
    195 
    196     assert (isfinite(f2));
     136    assert (isfinite(f0));
     137
     138    psF32 f1 = PAR[PM_PAR_I0]*f0;
     139    psF32 f = PAR[PM_PAR_SKY] + f1;
     140
    197141    assert (isfinite(f1));
    198     assert (isfinite(z0));
    199     assert (isfinite(f0));
     142    assert (isfinite(f));
    200143
    201144    if (deriv != NULL) {
    … …  
    203146
    204147        dPAR[PM_PAR_SKY]  = +1.0;
    205         dPAR[PM_PAR_I0]   = +f1;
    206 
    207         // gradient is infinite for z = 0; saturate at z = 0.01
    208         psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
    209 
    210         dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
    211         dPAR[PM_PAR_7]   *= 3.0;
    212 
    213         assert (isfinite(z1));
     148        dPAR[PM_PAR_I0]   = +f0;
     149
     150        if (z > 0.01) {
     151          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0);
     152          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     153          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     154          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     155          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     156          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     157          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z);
     158        } else {
     159          // gradient -> 0 for z -> 0, but has undef form
     160          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]);
     161          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     162          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     163          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     164          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     165          dPAR[PM_PAR_SXY]  = -1.0*z1;
     166          // dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001);
     167          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001); // factor of 16 to reduce the gain
     168        }
    214169        assert (isfinite(dPAR[PM_PAR_7]));
    215 
    216         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    217         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    218         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    219         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    220         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    221     }
    222     return (f0);
     170    }
     171    return (f);
    223172}
    224173
    … …  
    370319    psEllipseAxes axes;
    371320    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    372     float AspectRatio = axes.minor / axes.major;
    373 
    374     float index = 0.5 / PAR[PM_PAR_7];
    375     float bn = 1.9992*index - 0.3271;
    376 
    377     // the integral of a Sersic has an analytical form as follows:
    378     float logGamma = lgamma(2.0*index);
    379     float bnFactor = pow(bn, 2.0*index);
    380     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    381    
    382     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    383 
    384     return(Flux);
     321
     322    float Sindex = 0.5 / PAR[PM_PAR_7];
     323    float norm = pmSersicNorm (Sindex);
     324
     325    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     326
     327    return(flux);
    385328}
    386329
    … …  
    401344    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    402345
     346    float Sindex = 0.5 / PAR[PM_PAR_7];
     347    float kappa = pmSersicKappa (Sindex);
     348
    403349    // f = Io exp(-z^n) -> z^n = ln(Io/f)
    404     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    405     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     350    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     351    psF64 radius = axes.major * pow(zn, Sindex);
    406352
    407353    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f",

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

    r35768 r36085  
    6161// Lax parameter limits
    6262static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2,   0.5, -3.3, -0.5 };
    63 static float paramsMaxLax[] = {  1.0e5, 1.0e+8, +1.0e4, +1.0e4, 150.0, +3.3 , 5.0 };
     63static float paramsMaxLax[] = {  1.0e5, 1.00+9, +1.0e5, +1.0e5, 150.0, +3.3 , 5.0 };
    6464
    6565// Moderate parameter limits

  • trunk/psModules/src/objects/pmFootprintCullPeaks.c

    r34800 r36085  
    2525bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp);
    2626
    27  /*
    28   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    29   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    30   * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
    31   * starting point, discard the peak.
    32   */
     27/*
     28 * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
     29 * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
     30 * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
     31 * starting point, discard the peak.
     32 */
    3333
    3434# define IN_PEAK 1
    … …  
    4848
    4949    if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do
    50         return PS_ERR_NONE;
     50        return PS_ERR_NONE;
    5151    }
    5252
    … …  
    9191
    9292        // max flux is above threshold for brightest peak
    93       pmPeak *maxPeak = NULL;
    94       for (int i = 0; i < fp->peaks->n; i++) {
    95         pmPeak *testPeak = fp->peaks->data[i];
    96         float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
     93        pmPeak *maxPeak = NULL;
     94        for (int i = 0; i < fp->peaks->n; i++) {
     95            pmPeak *testPeak = fp->peaks->data[i];
     96            float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
    9797       
    98         if (isfinite(this_peak)) {
    99           maxPeak = fp->peaks->data[i];
    100           break;
    101         }
    102       }
    103       psAssert(maxPeak,"maxPeak was not set in these peaks");
    104       //      = fp->peaks->data[0];
     98            if (isfinite(this_peak)) {
     99                maxPeak = fp->peaks->data[i];
     100                break;
     101            }
     102        }
     103        psAssert(maxPeak,"maxPeak was not set in these peaks");
     104        //      = fp->peaks->data[0];
    105105        float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;
    106106
    … …  
    130130        }
    131131#if (0)
    132         if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) {
    133             psWarning ("upper limit: %f does not include max flux: %f",
    134                     threshbounds->data.F32[threshbounds->n-1], maxFlux);
    135         }
     132        if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) {
     133            psWarning ("upper limit: %f does not include max flux: %f",
     134                       threshbounds->data.F32[threshbounds->n-1], maxFlux);
     135        }
    136136#endif
    137137        psHistogram *threshist = psHistogramAllocGeneric(threshbounds);

  • trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

    r35768 r36085  
    4242#include "pmPCMdata.h"
    4343
     44# define SAVE_IMAGES 0
     45# if (SAVE_IMAGES)
     46int psphotSaveImage (psMetadata *header, psImage *image, char *filename);
     47# endif
     48
    4449# define FACILITY "psModules.objects"
    4550
    … …  
    9196    psF32 lambda = 0.001;
    9297    psF32 dLinear = 0.0;
    93     psF32 nu = 2.0;
     98    psF32 nu = 3.0;
    9499
    95100# if (USE_FFT && PRE_CONVOLVE)
    … …  
    130135        }
    131136
     137        char key[10]; // used for interactive responses
     138        bool testValue = false;
     139
    132140        // set a new guess for Alpha, Beta, Params
    133141        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
     142            if (min->isInteractive) {
     143                fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] ");
     144                if (!fgets(key, 8, stdin)) {
     145                    psWarning("Unable to read option");
     146                }
     147                switch (key[0]) {
     148                  case 'n':
     149                  case 'N':
     150                    done = true;
     151                    break;
     152                  case 'y':
     153                  case 'Y':
     154                  case '\n':
     155                    lambda *= 10.0;
     156                    continue;
     157                  default:
     158                    lambda *= 10.0;
     159                    continue;
     160                }
     161                if (done) break;
     162            }
    134163            min->iter ++;
    135164            if (min->iter >=  min->maxIter) break;
    … …  
    138167        }
    139168
     169        if (min->isInteractive) {
     170            p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: ");
     171            fprintf (stdout, "last chisq : %f\n", min->value);
     172            bool getOptions = true;
     173            while (getOptions) {
     174                fprintf (stdout, "options: (m)odify, (g)o, (q)uit: ");
     175                if (!fgets(key, 8, stdin)) {
     176                    psWarning("Unable to read option");
     177                }
     178                switch (key[0]) {
     179                  case 'm':
     180                  case 'M':
     181                    testValue = TRUE;
     182                    fprintf (stdout, "enter (Npar) (value): ");
     183                    int Npar = 0;
     184                    float value= 0;
     185                    int Nscan = fscanf (stdin, "%d %f", &Npar, &value);
     186                    if (Nscan != 2) {
     187                      fprintf (stderr, "scan failure\n");
     188                    }
     189                    Params->data.F32[Npar] = value;
     190                    break;
     191                  case 'g':
     192                  case 'G':
     193                  case '\n':
     194                    getOptions = false;
     195                    break;
     196                  default:
     197                    done = true;
     198                    break;
     199                }
     200                fprintf (stderr, "foo\n");
     201            }
     202            if (done) break;
     203        }
     204           
    140205        // dump some useful info if trace is defined
    141206        if (psTraceGetLevel(FACILITY) >= 6) {
    … …  
    202267        // XXX : Madsen gives suggestion for better use of rho
    203268        // rho is positive if the new chisq is smaller
    204         if (rho >= -1e-6) {
     269        if (testValue || (rho >= -1e-6)) {
    205270            min->value = Chisq;
    206271            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
    … …  
    215280          case 0:
    216281            if (rho >= -1e-6) {
    217                 lambda *= 0.25;
     282                lambda *= 0.1;
    218283            } else {
    219284                lambda *= 10.0;
    … …  
    234299            if (rho > 0.0) {
    235300                lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
    236                 nu = 2.0;
     301                nu = 3.0;
    237302            } else {
    238303                lambda *= nu;
    239                 nu *= 2.0;
     304                nu *= 3.0;
    240305            }
    241306            break;
    … …  
    474539    // XXX TEST : SAVE IMAGES
    475540# if (SAVE_IMAGES)
    476     psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
    477     psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
    478     psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
    479     psphotSaveImage (NULL, source->pixels, "obj.fits");
    480     psphotSaveImage (NULL, source->maskObj, "mask.fits");
    481     psphotSaveImage (NULL, source->variance, "variance.fits");
     541    static int Npass = 0;
     542    char name[128];
     543    snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name);
     544    snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name);
     545    snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name);
     546    snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name);
     547    snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name);
     548    snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name);
     549    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
     550        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
     551        if (!dmodelConv) continue;
     552        snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name);
     553    }
     554    Npass ++;
    482555# endif
    483556

  • trunk/psModules/src/objects/pmPCMdata.c

    r35768 r36085  
    173173}
    174174
     175int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) {
     176
     177    // set parameter mask based on fitting mode
     178    int nParams = 0;
     179
     180    switch (mode) {
     181      case PM_SOURCE_FIT_NORM:
     182        // fits only source normalization (Io)
     183        nParams = 1;
     184        psVectorInit (constraint->paramMask, 1);
     185        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     186        break;
     187
     188      case PM_SOURCE_FIT_PSF:
     189        // fits only x,y,Io
     190        nParams = 3;
     191        psVectorInit (constraint->paramMask, 1);
     192        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     193        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
     194        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
     195        break;
     196
     197      case PM_SOURCE_FIT_EXT:
     198        // fits all params except sky
     199        nParams = params->n - 1;
     200        psVectorInit (constraint->paramMask, 0);
     201        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     202        break;
     203
     204      case PM_SOURCE_FIT_EXT_AND_SKY:
     205        // fits all params including sky
     206        nParams = params->n;
     207        psVectorInit (constraint->paramMask, 0);
     208        break;
     209
     210      case PM_SOURCE_FIT_SHAPE:
     211        // fits shape (Sxx, Sxy, Syy) and Io
     212        nParams = 5;
     213        psVectorInit (constraint->paramMask, 1);
     214        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 0;
     215        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     216        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXX] = 0;
     217        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXY] = 0;
     218        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SYY] = 0;
     219        break;
     220
     221      case PM_SOURCE_FIT_INDEX:
     222        // fits only Io, index (PAR7) -- only Io for models with < 8 params
     223        psVectorInit (constraint->paramMask, 1);
     224        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     225        if (params->n == 7) {
     226            nParams = 1;
     227        } else {
     228            nParams = 2;
     229            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
     230        }
     231        break;
     232
     233      case PM_SOURCE_FIT_NO_INDEX:
     234        // fits all but index (PAR7) including sky
     235        psVectorInit (constraint->paramMask, 0);
     236        if (params->n == 7) {
     237            nParams = params->n;
     238        } else {
     239            nParams = params->n - 1;
     240            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
     241        }
     242        break;
     243      default:
     244        psAbort("invalid fitting mode");
     245    }
     246    return nParams;
     247}
     248
    175249pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
    176250
    … …  
    219293    constraint->checkLimits = model->modelLimits;
    220294
    221     // set parameter mask based on fitting mode
    222     int nParams = 0;
    223     switch (fitOptions->mode) {
    224       case PM_SOURCE_FIT_NORM:
    225         // NORM-only model fits only source normalization (Io)
    226         nParams = 1;
    227         psVectorInit (constraint->paramMask, 1);
    228         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    229         break;
    230       case PM_SOURCE_FIT_PSF:
    231         // PSF model only fits x,y,Io
    232         nParams = 3;
    233         psVectorInit (constraint->paramMask, 1);
    234         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    235         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
    236         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
    237         break;
    238       case PM_SOURCE_FIT_EXT:
    239         // EXT model fits all params (except sky)
    240         nParams = params->n - 1;
    241         psVectorInit (constraint->paramMask, 0);
    242         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    243         break;
    244       case PM_SOURCE_FIT_INDEX:
    245         // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
    246         psVectorInit (constraint->paramMask, 1);
    247         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    248         if (params->n == 7) {
    249             nParams = 1;
    250         } else {
    251             nParams = 2;
    252             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
    253         }
    254         break;
    255       case PM_SOURCE_FIT_NO_INDEX:
    256         // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
    257         psVectorInit (constraint->paramMask, 0);
    258         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    259         if (params->n == 7) {
    260             nParams = params->n - 1;
    261         } else {
    262             nParams = params->n - 2;
    263             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
    264         }
    265         break;
    266       default:
    267         psAbort("invalid fitting mode");
    268     }
     295    int nParams = pmPCMsetParams (constraint, fitOptions->mode);
    269296
    270297    if (nPix <  nParams + 1) {
    … …  
    341368    }
    342369
    343     // if we changed the fit mode, we need to update nDOF
    344     int nParams = 0;
    345     // set parameter mask based on fitting mode
    346     switch (fitOptions->mode) {
    347       case PM_SOURCE_FIT_NORM:
    348         // NORM-only model fits only source normalization (Io)
    349         nParams = 1;
    350         psVectorInit (pcm->constraint->paramMask, 1);
    351         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    352         break;
    353       case PM_SOURCE_FIT_PSF:
    354         // PSF model only fits x,y,Io
    355         nParams = 3;
    356         psVectorInit (pcm->constraint->paramMask, 1);
    357         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    358         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
    359         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
    360         break;
    361       case PM_SOURCE_FIT_EXT:
    362         // EXT model fits all params (except sky)
    363         nParams = model->params->n - 1;
    364         psVectorInit (pcm->constraint->paramMask, 0);
    365         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    366         break;
    367       case PM_SOURCE_FIT_INDEX:
    368         // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
    369         psVectorInit (pcm->constraint->paramMask, 1);
    370         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    371         if (model->params->n == 7) {
    372             nParams = 1;
    373         } else {
    374             nParams = 2;
    375             pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
    376         }
    377         break;
    378       case PM_SOURCE_FIT_NO_INDEX:
    379         // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
    380         psVectorInit (pcm->constraint->paramMask, 0);
    381         pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    382         if (model->params->n == 7) {
    383             nParams = model->params->n - 1;
    384         } else {
    385             nParams = model->params->n - 2;
    386             pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
    387         }
    388         break;
    389       default:
    390         psAbort("invalid fitting mode");
    391     }
     370    int nParams = pmPCMsetParams (pcm->constraint, fitOptions->mode);
    392371
    393372    if (pcm->nPix <  nParams + 1) {
    … …  
    478457}
    479458
     459// construct a realization of the source model
     460bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) {
     461
     462    PS_ASSERT_PTR_NON_NULL(source, false);
     463
     464    // if we already have a cached image, re-use that memory
     465    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
     466    psImageInit (source->modelFlux, 0.0);
     467
     468    // modelFlux always has unity normalization (I0 = 1.0)
     469    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     470
     471    // convolve the model image with the PSF
     472    if (USE_1D_GAUSS) {
     473        // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     474        // * the model flux is not masked
     475        // * threading takes place above this level
     476       
     477        // define the Gauss parameters from the psf
     478        pmModel *modelPSF = source->modelPSF;
     479        psAssert (modelPSF, "psf model must be defined");
     480   
     481        psEllipseAxes axes;
     482        bool useReff = pmModelUseReff (modelPSF->type);
     483        psF32 *PAR = modelPSF->params->data.F32;
     484        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     485   
     486        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     487        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     488
     489        float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     490        float nsigma = 2.0;
     491
     492        psImageSmooth (source->modelFlux, sigma, nsigma);
     493    } else {
     494        // make sure we save a cached copy of the psf flux
     495        pmSourceCachePSF (source, maskVal);
     496
     497        // convert the cached cached psf model for this source to a psKernel
     498        psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
     499        if (!psf) {
     500            // NOTE: this only happens if the source is too close to an edge
     501            model->flags |= PM_MODEL_STATUS_BADARGS;
     502            return NULL;
     503        }
     504
     505        // XXX not sure if I can place the output on top of the input
     506        psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf);
     507    }
     508    return true;
     509}
     510

  • trunk/psModules/src/objects/pmPCMdata.h

    r32725 r36085  
    9898bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize);
    9999
     100bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize);
     101
    100102/// @}
    101103# endif /* PM_PCM_DATA_H */

  • trunk/psModules/src/objects/pmSourceFitModel.c

    r35768 r36085  
    6666    opt->gainFactorMode = 0;
    6767    opt->chisqConvergence = true;
     68    opt->isInteractive = false;
    6869
    6970    return opt;
    … …  
    247248    myMin->gainFactorMode = options->gainFactorMode;
    248249    myMin->chisqConvergence = options->chisqConvergence;
     250    myMin->isInteractive = options->isInteractive;
    249251
    250252    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

  • trunk/psModules/src/objects/pmSourceFitModel.h

    r35768 r36085  
    2121    PM_SOURCE_FIT_EXT_AND_SKY,
    2222    PM_SOURCE_FIT_INDEX,
     23    PM_SOURCE_FIT_SHAPE,
    2324    PM_SOURCE_FIT_NO_INDEX,
    2425    PM_SOURCE_FIT_TRAIL,
    … …  
    3738    int gainFactorMode;
    3839    bool chisqConvergence;
     40    bool isInteractive;
    3941} pmSourceFitOptions;
    4042

  • trunk/psModules/src/objects/pmSourceFitPCM.c

    r35768 r36085  
    6868    myMin->chisqConvergence = fitOptions->chisqConvergence;
    6969    myMin->gainFactorMode = fitOptions->gainFactorMode;
     70    myMin->isInteractive = fitOptions->isInteractive;
    7071
    7172    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

  • trunk/psModules/src/objects/pmSourceFitSet.c

    r35768 r36085  
    570570    myMin->gainFactorMode = options->gainFactorMode;
    571571    myMin->chisqConvergence = options->chisqConvergence;
     572    myMin->isInteractive = options->isInteractive;
    572573
    573574    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

  • trunk/psModules/src/psmodules.h

    r35768 r36085  
    129129#include <pmModelFuncs.h>
    130130#include <pmModel.h>
     131#include <pmModel_CentralPixel.h>
    131132
    132133#include <pmSourceMasks.h>

  • trunk/psModules/test/objects

    • Property svn:ignore

      •  

        old new  
        88tap_pmSourceFitModel_Delta
        99tap_pmSourcePhotometry
         10tap_pmModel_CentralPixel
         11tap_pmModel_CentralPixel_v2

  • trunk/psModules/test/objects/Makefile.am

    r29547 r36085  
    1919        tap_pmModelUtils \
    2020        tap_pmModelClass \
     21        tap_pmModel_CentralPixel \
     22        tap_pmModel_CentralPixel_v2 \
    2123        tap_pmPSF \
    2224        tap_pmTrend2D \
Note: See TracChangeset for help on using the changeset viewer.