Changeset 35966 for branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c

Timestamp:

Aug 19, 2013, 6:40:03 AM (13 years ago)

Author:

eugene

Message:

replace ad-hoc code for central pixel with pmModelCP call

File:

• branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c (modified) (3 diffs)

branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c

@@ -48 +48 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_DEV.h"
@@ -89 +90 @@
 static bool limitsApply = true;         // Apply limits?
 
-// # include "pmModel_SERSIC.CP.h"
-
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -129 +128 @@
     // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel)
     psF32 radius = hypot(X, Y);
-    if (radius < 1.0) {
-      // subdivide the central 2,3,4 pixels by Nx,Ny 
-      float Npix = 0.0;
-      float Fpix = 0.0;
-      float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS];
-      float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS];
-      # define DELTA 0.02
-      for (float ix = 0.1; ix <= 0.9; ix += DELTA) {
-        for (float iy = 0.1; iy <= 0.9; iy += DELTA) {
-          psF32 X  = Xpix + ix;
-          psF32 Y  = Ypix + iy;
-          psF32 px = X / PAR[PM_PAR_SXX];
-          psF32 py = Y / PAR[PM_PAR_SYY];
-          psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-          
-          // sqrt(z) is r
-          float q = kappa*pow(z,ALPHA);
-          psF32 f0 = exp(-q);
-          
-          psF32 f1 = PAR[PM_PAR_I0]*f0;
-          psF32 fx = PAR[PM_PAR_SKY] + f1;
-          Fpix += fx;
-          Npix += 1.0;
-        }
-      }
-      f = Fpix / Npix;
+    if (radius <= 1.5) {
+        // Nsub ~ 10*index^2 + 1
+        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV"));
+        int Nsub = 2 * ((int)(25 / axes.minor)) + 1;
+        Nsub = PS_MIN (Nsub, 121);
+        Nsub = PS_MAX (Nsub, 11);
+        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
+        f1 = PAR[PM_PAR_I0]*f0;
+        f = PAR[PM_PAR_SKY] + f1;
     }