Index: trunk/psModules/test/objects/Makefile.am
===================================================================
--- trunk/psModules/test/objects/Makefile.am	(revision 31451)
+++ trunk/psModules/test/objects/Makefile.am	(revision 36085)
@@ -19,4 +19,6 @@
 	tap_pmModelUtils \
 	tap_pmModelClass \
+	tap_pmModel_CentralPixel \
+	tap_pmModel_CentralPixel_v2 \
 	tap_pmPSF \
 	tap_pmTrend2D \
Index: trunk/psModules/test/objects/tap_pmModel_CentralPixel.c
===================================================================
--- trunk/psModules/test/objects/tap_pmModel_CentralPixel.c	(revision 36085)
+++ trunk/psModules/test/objects/tap_pmModel_CentralPixel.c	(revision 36085)
@@ -0,0 +1,202 @@
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "tap.h"
+#include "pstap.h"
+
+int main (int argc, char **argv)
+{
+    if (argc != 3) {
+	fprintf (stderr, "USAGE: tap_pmModel_CentralPixel (model) (set)\n");
+	exit (2);
+    }
+
+    int set = atoi(argv[2]);
+
+    psMemId id = psMemGetId();
+
+    plan_tests(240);
+
+    pmModelCP *t0 = pmModelCP_Alloc();
+    ok (t0, "allocated pmModelCP");
+    psFree (t0);
+
+    pmModelCPset *t1 = pmModelCPset_Alloc();
+    ok (t1, "allocated pmModelCPset");
+    psFree (t1);
+    
+    pmModelCPset *cpset = pmModelCP_Load (argv[1]);
+    ok (cpset, "loaded pmModelCPset from file");
+
+    ok (cpset->RmajorNitem ==  4, "correct number of Rmajor values");
+    ok (cpset->AratioNitem ==  7, "correct number of Aratio values");
+    ok (cpset->SindexNitem == 10, "correct number of Sindex values");
+    ok (cpset->images->n == 280, "correct number of CP images");
+
+    pmModelCP *cp = NULL;
+    
+    if (1) {
+	cp = pmModelCP_GetImage (cpset, 0.0, 1.0, 1.0);
+	ok (cp, "returned a cp image");
+	if (cp) {
+	    ok_float_tol (cp->Rmajor, 0.0, 0.001, "got image with correct Rmajor");
+	    ok_float_tol (cp->Aratio, 1.0, 0.001, "got image with correct Aratio");
+	    ok_float_tol (cp->Sindex, 1.0, 0.001, "got image with correct Sindex");
+	}
+
+	cp = pmModelCP_GetImage (cpset, 1.0, 1.0, 1.0);
+	ok (cp, "returned a cp image");
+	if (cp) {
+	    ok_float_tol (cp->Rmajor, 1.0, 0.001, "got image with correct Rmajor");
+	    ok_float_tol (cp->Aratio, 1.0, 0.001, "got image with correct Aratio");
+	    ok_float_tol (cp->Sindex, 1.0, 0.001, "got image with correct Sindex");
+	}
+    
+	cp = pmModelCP_GetImage (cpset, 0.0, 0.4, 1.0);
+	ok (cp, "returned a cp image");
+	if (cp) {
+	    ok_float_tol (cp->Rmajor, 0.0, 0.001, "got image with correct Rmajor");
+	    ok_float_tol (cp->Aratio, 0.4, 0.001, "got image with correct Aratio");
+	    ok_float_tol (cp->Sindex, 1.0, 0.001, "got image with correct Sindex");
+	}
+    
+	cp = pmModelCP_GetImage (cpset, 0.0, 0.4, 3.5);
+	ok (cp, "returned a cp image");
+	if (cp) {
+	    ok_float_tol (cp->Rmajor, 0.0, 0.001, "got image with correct Rmajor");
+	    ok_float_tol (cp->Aratio, 0.4, 0.001, "got image with correct Aratio");
+	    ok_float_tol (cp->Sindex, 3.5, 0.001, "got image with correct Sindex");
+	}
+    }
+    
+    float valuePixel = NAN;
+    float valueModel = NAN;
+
+    switch (set) {
+      case 0:
+	cp = pmModelCP_GetImage (cpset, 1.0, 1.0, 1.0);
+	valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.0, 0.0);
+	valueModel = pmModelCP_FullSersic (0.0, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	fprintf (stdout, "%f / %f = %f\n", valuePixel, valueModel, valuePixel / valueModel);
+	break;
+
+      case 1:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.0, 0.0);
+	    valueModel = pmModelCP_FullSersic (0.0, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 2:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.0, 30.0);
+	    valueModel = pmModelCP_FullSersic (0.0, 0.0, 30.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 3:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.5, 0.0);
+	    valueModel = pmModelCP_FullSersic (0.0, 0.5, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 4:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.5, 0.5, 0.0);
+	    valueModel = pmModelCP_FullSersic (0.5, 0.5, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 5:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, -0.5, 0.5, 30.0);
+	    valueModel = pmModelCP_FullSersic (-0.5, 0.5, 30.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 6:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.0, -0.5, 0.0);
+	    valueModel = pmModelCP_FullSersic (0.0, -0.5, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 7:
+	cp = cpset->images->data[224];
+	float delta = 1.0 / 11.0;
+	float offset = -10*delta;
+	for (float dx = offset; dx <= 1.0; dx += delta) {
+	    valuePixel = pmModelCP_GetFlux (cp, dx, 0.0, 0.0);
+	    valueModel = pmModelCP_FullSersic (dx, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 8:
+	cp = cpset->images->data[224];
+	valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.0, 45.0);
+	valueModel = pmModelCP_FullSersic (0.0, 0.0, 45.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	break;
+
+      case 9:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, -0.5, -0.5, 0.0);
+	    valueModel = pmModelCP_FullSersic (-0.5, -0.5, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+      case 10:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.5, 0.0, 0.0);
+	    valueModel = pmModelCP_FullSersic (0.5, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+
+      case 11:
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, -0.5, 0.0, 0.0);
+	    valueModel = pmModelCP_FullSersic (-0.5, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	    fprintf (stdout, "%f / %f = %f | %f  %f  %f\n", valuePixel, valueModel, valuePixel / valueModel, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	break;
+      case 12:
+	psTimerStart ("test");
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valuePixel = pmModelCP_GetFlux (cp, 0.0, 0.0, 0.0);
+	}
+	fprintf (stderr, "CP code: %f\n", psTimerMark ("test"));
+	psTimerStart ("test");
+	for (int i = 0; i < cpset->images->n; i++) {
+	    cp = cpset->images->data[i];
+	    valueModel = pmModelCP_FullSersic (0.0, 0.0, 0.0, pow(10.0, cp->Rmajor), cp->Aratio, cp->Sindex);
+	}
+	fprintf (stderr, "Full code: %f\n", psTimerMark ("test"));
+	break;
+    }
+
+    psFree (cpset);
+    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
+
+    return exit_status();
+}
Index: trunk/psModules/test/objects/tap_pmModel_CentralPixel_v2.c
===================================================================
--- trunk/psModules/test/objects/tap_pmModel_CentralPixel_v2.c	(revision 36085)
+++ trunk/psModules/test/objects/tap_pmModel_CentralPixel_v2.c	(revision 36085)
@@ -0,0 +1,52 @@
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "tap.h"
+#include "pstap.h"
+
+int main (int argc, char **argv)
+{
+    if (argc != 8) {
+	fprintf (stderr, "USAGE: tap_pmModel_CentralPixel (model) (dx) (dy) (theta) (Reff) (Arat) (Sidx)\n");
+	exit (2);
+    }
+
+    float dx = atof(argv[2]);
+    float dy = atof(argv[3]);
+    float theta = atof(argv[4]);
+
+    float Reff = atof(argv[5]);
+    float Arat = atof(argv[6]);
+    float Sidx = atof(argv[7]);
+
+    psMemId id = psMemGetId();
+
+    plan_tests(6);
+
+    pmModelCPset *cpset = pmModelCP_Load (argv[1]);
+    ok (cpset, "loaded pmModelCPset from file");
+
+    ok (cpset->RmajorNitem ==  4, "correct number of Rmajor values");
+    ok (cpset->AratioNitem ==  7, "correct number of Aratio values");
+    ok (cpset->SindexNitem == 10, "correct number of Sindex values");
+    ok (cpset->images->n == 280, "correct number of CP images");
+
+    pmModelCP *cp = NULL;
+    
+    float valuePixel = NAN;
+    float valueModel = NAN;
+
+    cp = pmModelCP_GetImage (cpset, log10(Reff), Arat, Sidx);
+    valuePixel = pmModelCP_GetFlux (cp, dx, dy, theta);
+    valueModel = pmModelCP_FullSersic (dx, dy, theta, Reff, Arat, Sidx);
+    fprintf (stdout, "%f / %f = %f\n", valuePixel, valueModel, valuePixel / valueModel);
+
+    fprintf (stdout, "%f, %f, %f\n", cp->Rmajor, cp->Aratio, cp->Sindex);
+
+    psFree (cpset);
+    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
+
+    return exit_status();
+}