Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c	(revision 35965)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c	(revision 35966)
@@ -48,4 +48,5 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_DEV.h"
@@ -89,6 +90,4 @@
 static bool limitsApply = true;         // Apply limits?
 
-// # include "pmModel_SERSIC.CP.h"
-
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -129,30 +128,13 @@
     // 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;
     }   
 
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c	(revision 35965)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c	(revision 35966)
@@ -45,4 +45,5 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_EXP.h"
@@ -124,29 +125,8 @@
     // only worry about the central 4 pixels at most
     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];
-      for (float ix = 0.1; ix < 1.0; ix += 0.2) {
-	for (float iy = 0.1; iy < 1.0; iy += 0.2) {
-	  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*sqrt(z);
-	  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) {
+	f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
+	f1 = PAR[PM_PAR_I0]*f0;
+	f = PAR[PM_PAR_SKY] + f1;
     }