Index: /trunk/psModules/src/objects/Makefile.am
===================================================================
--- /trunk/psModules/src/objects/Makefile.am	(revision 40490)
+++ /trunk/psModules/src/objects/Makefile.am	(revision 40491)
@@ -92,5 +92,6 @@
 	models/pmModel_EXP.c \
 	models/pmModel_DEV.c \
-	models/pmModel_TRAIL.c
+	models/pmModel_TRAIL.c \
+	models/pmModel_HSC_V1.c
 
 pkginclude_HEADERS = \
@@ -142,5 +143,6 @@
 	models/pmModel_EXP.h \
 	models/pmModel_DEV.h \
-	models/pmModel_TRAIL.h
+	models/pmModel_TRAIL.h \
+	models/pmModel_HSC_V1.h
 
 CLEANFILES = *~
Index: /trunk/psModules/src/objects/models/fwhm.sh
===================================================================
--- /trunk/psModules/src/objects/models/fwhm.sh	(revision 40490)
+++ /trunk/psModules/src/objects/models/fwhm.sh	(revision 40491)
@@ -161,8 +161,8 @@
 end
 
-# qgauss is like ps1_v1 with z^2.25
-macro find.fwhm.qgauss
+# hsc
+macro find.fwhm.hscv1
   if ($0 != 2)
-    echo "USAGE: find.qgauss (K)"
+    echo "USAGE: find.fwhm (K)"
     break
   end
@@ -180,15 +180,16 @@
   create z 0 50 $dz
 
-  # set f = 1.0 / (1.0 + $K*z + z^2.25)
-  # f = $K*z + z^2.25 - 1.0, find where f == 0.0
-  
-  set f = $K*z + z^2.25 - 1.0
-  set dfdz = $K + 2.25*z^1.25
+  # set f = 1.0 / (1.0 + $K*z + z^1.667)
+  # f = $K*z + z^1.667 - 1.0, find where f == 0.0
+  
+  set f = $K*z + z^1.8 - 1.0
+  set dfdz = $K + 1.8*z^0.8
 
   #lim -n 0 z f; clear; box; plot z f
   #lim -n 1 z dfdz; clear; box; plot z dfdz  
 
-  $nZ0 = 0
-  $nZ1 = 5 / $dz
+  # these just need to bound the solution
+  $nZ0 = 0.5 / $dz
+  $nZ1 = 5.0 / $dz
 
   $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
@@ -206,4 +207,66 @@
     $dZ = $Fg / $dFdz_g
 
+    # echo -no-return $Zg $Fg $dZ $dFdz_g
+
+    $Zg -= $dZ
+
+    $Zg = max ($Zg , 0)
+    $Zg = min ($Zg , z[-2])
+    $nZg = int($Zg / $dz)
+    
+    # echo ": $Zg $nZg"
+  end 
+  # echo $Zg $Fg $dZ $dFdz_g
+  $Zhm = $Zg
+  $FWHM = 2*sqrt(2*$Zg)
+
+  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
+end
+
+# qgauss is like ps1_v1 with z^2.25
+macro find.fwhm.qgauss
+  if ($0 != 2)
+    echo "USAGE: find.qgauss (K)"
+    break
+  end
+  
+  $K = $1
+
+  if ($K == 0.0)
+    $Zhm = 1.0
+    $FWHM = 2*sqrt(2)
+    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
+    return
+  end
+
+  $dz = 0.01
+  create z 0 50 $dz
+
+  # set f = 1.0 / (1.0 + $K*z + z^2.25)
+  # f = $K*z + z^2.25 - 1.0, find where f == 0.0
+  
+  set f = $K*z + z^2.25 - 1.0
+  set dfdz = $K + 2.25*z^1.25
+
+  #lim -n 0 z f; clear; box; plot z f
+  #lim -n 1 z dfdz; clear; box; plot z dfdz  
+
+  $nZ0 = 0
+  $nZ1 = 5 / $dz
+
+  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
+  $nZg = int($Zg / $dz)
+  $dZ = 1.0
+
+  for i 0 10
+    # interpolate between $nZg and $nZg + 1
+    $Fg = f[$nZg] + ($Zg - $nZg*$dz)*(f[$nZg+1] - f[$nZg])/$dz
+    $dFdz_g = dfdz[$nZg] + ($Zg - $nZg*$dz)*(dfdz[$nZg+1] - dfdz[$nZg])/$dz
+
+    # $Fg = f[$nZg]
+    # $dFdz_g = dfdz[$nZg]
+
+    $dZ = $Fg / $dFdz_g
+
     # echo $Zg $Fg $dZ $dFdz_g
 
@@ -223,5 +286,5 @@
   if ($0 != 4)
     echo "USAGE: fwhm.trend (model) (struct) (output)"
-    echo "  model: pgauss, rgauss, ps1v1, qgauss"
+    echo "  model: pgauss, rgauss, ps1v1, qgauss, hscv1"
     break
   end
@@ -229,5 +292,5 @@
   delete fwhm_v k_v
 
-  $SAVE = 0
+  $SAVE = 1
   if ($SAVE) exec rm -f $3
   if ($SAVE) output $3
@@ -240,5 +303,5 @@
   if ($SAVE) echo "# define MIN_FWHM_BIN $minK"
   if ($SAVE) echo "# define N_FWHM_BIN {int(($maxK - $minK) / $delK) + 1}"
-  if ($SAVE) echo "static float $2[] = \{"
+#  if ($SAVE) echo "static float $2[] = \{"
 
   for k $minK $maxK $delK -incl
@@ -247,8 +310,8 @@
     concat $FWHM fwhm_v
   end
-  if ($SAVE) echo "\}@"
+#  if ($SAVE) echo "\}@"
   if ($SAVE) output stdout
 
-  lim k_v fwhm_v; clear; box; plot k_v fwhm_v -pt 10 -sz 1.0
-end
-
+#  lim k_v fwhm_v; clear; box; plot k_v fwhm_v -pt 10 -sz 1.0
+end
+
Index: /trunk/psModules/src/objects/models/hscv1.fwhm.h
===================================================================
--- /trunk/psModules/src/objects/models/hscv1.fwhm.h	(revision 40491)
+++ /trunk/psModules/src/objects/models/hscv1.fwhm.h	(revision 40491)
@@ -0,0 +1,111 @@
+# define FWHM_BIN 0.2
+# define MIN_FWHM_BIN -1
+# define N_FWHM_BIN 106
+static float HSC_V1_Scale[] = {
+3.75, // -1.0, 1.76
+3.55, // -0.8, 1.57
+3.35, // -0.6, 1.40
+3.17, // -0.4, 1.25
+2.99, // -0.2, 1.12
+2.83, //  0.0, 1.00
+2.68, //  0.2, 0.90
+2.54, //  0.4, 0.81
+2.41, //  0.6, 0.73
+2.30, //  0.8, 0.66
+2.19, //  1.0, 0.60
+2.10, //  1.2, 0.55
+2.01, //  1.4, 0.51
+1.93, //  1.6, 0.47
+1.86, //  1.8, 0.43
+1.79, //  2.0, 0.40
+1.74, //  2.2, 0.38
+1.68, //  2.4, 0.35
+1.63, //  2.6, 0.33
+1.58, //  2.8, 0.31
+1.54, //  3.0, 0.30
+1.50, //  3.2, 0.28
+1.46, //  3.4, 0.27
+1.43, //  3.6, 0.25
+1.39, //  3.8, 0.24
+1.36, //  4.0, 0.23
+1.33, //  4.2, 0.22
+1.31, //  4.4, 0.21
+1.28, //  4.6, 0.20
+1.26, //  4.8, 0.20
+1.23, //  5.0, 0.19
+1.21, //  5.2, 0.18
+1.19, //  5.4, 0.18
+1.17, //  5.6, 0.17
+1.15, //  5.8, 0.17
+1.13, //  6.0, 0.16
+1.12, //  6.2, 0.16
+1.10, //  6.4, 0.15
+1.08, //  6.6, 0.15
+1.07, //  6.8, 0.14
+1.05, //  7.0, 0.14
+1.04, //  7.2, 0.14
+1.03, //  7.4, 0.13
+1.01, //  7.6, 0.13
+1.00, //  7.8, 0.13
+0.99, //  8.0, 0.12
+0.98, //  8.2, 0.12
+0.97, //  8.4, 0.12
+0.95, //  8.6, 0.11
+0.94, //  8.8, 0.11
+0.93, //  9.0, 0.11
+0.92, //  9.2, 0.11
+0.91, //  9.4, 0.10
+0.91, //  9.6, 0.10
+0.90, //  9.8, 0.10
+0.89, // 10.0, 0.10
+0.88, // 10.2, 0.10
+0.87, // 10.4, 0.09
+0.86, // 10.6, 0.09
+0.85, // 10.8, 0.09
+0.85, // 11.0, 0.09
+0.84, // 11.2, 0.09
+0.83, // 11.4, 0.09
+0.83, // 11.6, 0.09
+0.82, // 11.8, 0.08
+0.81, // 12.0, 0.08
+0.81, // 12.2, 0.08
+0.80, // 12.4, 0.08
+0.79, // 12.6, 0.08
+0.79, // 12.8, 0.08
+0.78, // 13.0, 0.08
+0.77, // 13.2, 0.08
+0.77, // 13.4, 0.07
+0.76, // 13.6, 0.07
+0.76, // 13.8, 0.07
+0.75, // 14.0, 0.07
+0.75, // 14.2, 0.07
+0.74, // 14.4, 0.07
+0.74, // 14.6, 0.07
+0.73, // 14.8, 0.07
+0.73, // 15.0, 0.07
+0.72, // 15.2, 0.07
+0.72, // 15.4, 0.06
+0.71, // 15.6, 0.06
+0.71, // 15.8, 0.06
+0.70, // 16.0, 0.06
+0.70, // 16.2, 0.06
+0.70, // 16.4, 0.06
+0.69, // 16.6, 0.06
+0.69, // 16.8, 0.06
+0.68, // 17.0, 0.06
+0.68, // 17.2, 0.06
+0.68, // 17.4, 0.06
+0.67, // 17.6, 0.06
+0.67, // 17.8, 0.06
+0.66, // 18.0, 0.06
+0.66, // 18.2, 0.05
+0.66, // 18.4, 0.05
+0.65, // 18.6, 0.05
+0.65, // 18.8, 0.05
+0.65, // 19.0, 0.05
+0.64, // 19.2, 0.05
+0.64, // 19.4, 0.05
+0.64, // 19.6, 0.05
+0.63, // 19.8, 0.05
+0.63, // 20.0, 0.05
+};
Index: /trunk/psModules/src/objects/models/pmModel_HSC_V1.c
===================================================================
--- /trunk/psModules/src/objects/models/pmModel_HSC_V1.c	(revision 40491)
+++ /trunk/psModules/src/objects/models/pmModel_HSC_V1.c	(revision 40491)
@@ -0,0 +1,505 @@
+/******************************************************************************
+ * this file defines the HSC_V1 source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
+
+   power-law with fitted linear term
+   1 / (1 + kz + z^(3.33/2))  (z = r^2, so r^3.33)
+
+   * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
+   * PM_PAR_I0 1    - central intensity
+   * PM_PAR_XPOS 2  - X center of object
+   * PM_PAR_YPOS 3  - Y center of object
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
+   * PM_PAR_SXY 6   - X*Y term of elliptical contour
+   * PM_PAR_7   7   - amplitude of the linear component (k)
+   *****************************************************************************/
+
+#include <stdio.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
+#include "pmModelClass.h"
+#include "pmModel.h"
+#include "pmModelUtils.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
+#include "pmSource.h"
+#include "pmSourceFitModel.h"
+#include "pmPSF.h"
+#include "pmPSFtry.h"
+#include "pmDetections.h"
+
+#include "pmModel_HSC_V1.h"
+
+# define PM_MODEL_NPARAM          8
+# define PM_MODEL_FUNC            pmModelFunc_HSC_V1
+# define PM_MODEL_FLUX            pmModelFlux_HSC_V1
+# define PM_MODEL_GUESS           pmModelGuess_HSC_V1
+# define PM_MODEL_LIMITS          pmModelLimits_HSC_V1
+# define PM_MODEL_RADIUS          pmModelRadius_HSC_V1
+# define PM_MODEL_SET_FWHM        pmModelSetFWHM_HSC_V1
+# define PM_MODEL_FROM_PSF        pmModelFromPSF_HSC_V1
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_HSC_V1
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_HSC_V1
+# define PM_MODEL_SET_LIMITS      pmModelSetLimits_HSC_V1
+
+# define ALPHA   1.8
+# define ALPHA_M 0.8
+
+// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
+// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
+// values need to be pixel coords
+
+// Lax parameter limits
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
+
+// Moderate parameter limits
+// Tolerate a small divot (k < 0)
+static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
+
+// Strict parameter limits
+// k = PAR_7 < 0 is very undesirable (big divot in the middle)
+static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
+
+// Parameter limits to use
+static float *paramsMinUse = paramsMinLax;
+static float *paramsMaxUse = paramsMaxLax;
+static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
+
+static bool limitsApply = true;         // Apply limits?
+
+psF32 PM_MODEL_FUNC (psVector *deriv,
+                     const psVector *params,
+                     const psVector *pixcoord)
+{
+    psF32 *PAR = params->data.F32;
+
+    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
+    psF32 Y  = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS];
+    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;
+
+    // XXX if the elliptical contour is defined in valid way, this step should not be required.
+    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
+    // for negative values of z
+    // XXX use an assert here to force the elliptical parameters to be correctly determined
+    // if (z < 0) z = 0;
+    assert (z >= 0);
+
+    psF32 zp = pow(z,ALPHA_M);
+    psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
+
+    psF32 r1 = PAR[PM_PAR_I0]*r;
+    psF32 f  = r1 + PAR[PM_PAR_SKY];
+
+    if (deriv != NULL) {
+        psF32 *dPAR = deriv->data.F32;
+
+        // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
+        psF32 t = r1*r;
+        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
+
+        dPAR[PM_PAR_SKY]  = +1.0;
+        dPAR[PM_PAR_I0]   = +r;
+        dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
+        dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
+        // the extra factor of 2 below is needed to avoid excessive swings
+        dPAR[PM_PAR_SXX]  = +4.0*q*px*px/PAR[PM_PAR_SXX];
+        dPAR[PM_PAR_SYY]  = +4.0*q*py*py/PAR[PM_PAR_SYY];
+        dPAR[PM_PAR_SXY]  = -q*X*Y;
+        dPAR[PM_PAR_7]    = -t*z;
+    }
+    return(f);
+}
+
+// define the parameter limits
+// AR_MAX is the maximum allowed axis ratio
+// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2)
+# define AR_MAX 20.0
+# define AR_RATIO 0.99
+bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
+{
+    if (!limitsApply) {
+        return true;
+    }
+    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
+
+    // we need to calculate the limits for SXY specially
+    float q2 = NAN;
+    if (nParam == PM_PAR_SXY) {
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
+        q1 = (q1 < 0.0) ? 0.0 : q1;
+        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
+        // angle and let f2,f1 fight it out
+        q2 = 0.5*sqrtf(q1);
+    }
+
+    switch (mode) {
+      case PS_MINIMIZE_BETA_LIMIT: {
+          psAssert(beta, "Require beta to limit beta");
+          float limit = betaUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (fabs(beta[nParam]) > fabs(limit)) {
+              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
+              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
+                      nParam, beta[nParam], limit);
+              return false;
+          }
+          return true;
+      }
+      case PS_MINIMIZE_PARAM_MIN: {
+          psAssert(params, "Require parameters to limit parameters");
+          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
+          float limit = paramsMinUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (params[nParam] < limit) {
+              params[nParam] = limit;
+              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
+                      nParam, params[nParam], limit);
+              return false;
+          }
+          return true;
+      }
+      case PS_MINIMIZE_PARAM_MAX: {
+          psAssert(params, "Require parameters to limit parameters");
+          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
+          float limit = paramsMaxUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (params[nParam] > limit) {
+              params[nParam] = limit;
+              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
+                      nParam, params[nParam], limit);
+              return false;
+          }
+          return true;
+      }
+      default:
+        psAbort("invalid choice for limits");
+    }
+    psAbort("should not reach here");
+    return false;
+}
+
+
+// make an initial guess for parameters
+// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
+bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
+{
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
+    PAR[PM_PAR_SKY]  = 0.0;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // extra parameter
+    PAR[PM_PAR_7]    = 0.5;
+
+    return(true);
+}
+
+// integrate the model to get the full flux
+psF64 PM_MODEL_FLUX (const psVector *params)
+{
+    float z, norm;
+
+    psF32 *PAR = params->data.F32;
+
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
+    norm = 0.0;
+
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
+    float f1, f2;
+    for (float r = DR; r < 150; r += DR) {
+	z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        r += DR;
+	z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        norm += f0 + 4*f1 + f2;
+        f0 = f2;
+    }
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
+
+    return(Flux);
+}
+
+// define this function so it never returns Inf or NaN
+// return the radius which yields the requested flux
+psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
+{
+    psF64 z;
+
+    psF32 *PAR = params->data.F32;
+
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+    psF64 sigma = axes.major;
+
+    // we can solve directly if PAR_7 goes to 0.0
+    if (PAR[PM_PAR_7] == 0.0) { 
+        psF32 z =  powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
+        return ( sigma * sqrt (2.0 * z) );
+    }
+    psF64 limit = flux / PAR[PM_PAR_I0];
+
+    // use the fact that f is monotonically decreasing
+    z = 0;
+
+    // choose a z value guaranteed to be beyond our limit
+    float z0 = 0.0;
+    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+	      PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
+    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
+
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+	float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
+	z = PS_MAX(z, 0.0);
+    }
+    psF64 radius = sigma * sqrt (2.0 * z);
+
+    if (isnan(radius))
+        psAbort("error in code: radius is NaN");
+
+    return (radius);
+}
+
+// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value:
+# include "hscv1.fwhm.h"
+
+psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
+
+    psF32 *PAR = params->data.F32;
+
+    float core = PAR[PM_PAR_7];
+
+    if (!isfinite(core)) return (2.0*M_SQRT2*sigma);
+
+    int binCore = MAX(0, MIN (N_FWHM_BIN - 2, (int)((core - MIN_FWHM_BIN)/FWHM_BIN)));
+    float coreBin = binCore * FWHM_BIN + MIN_FWHM_BIN;
+
+    float scale = (core - coreBin) * (HSC_V1_Scale[binCore + 1] - HSC_V1_Scale[binCore]) / FWHM_BIN + HSC_V1_Scale[binCore];
+
+    // fprintf (stderr, "core: %f, binCore: %d, scale: %f\n", core, binCore, scale);
+
+    return (scale * sigma);
+}
+
+bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf)
+{
+    psF32 *out = modelPSF->params->data.F32;
+    psF32 *in  = modelFLT->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
+            out[i] = in[i];
+        } else {
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+        }
+    }
+
+    // the 2D PSF model fits polarization terms (E0,E1,E2)
+    // convert to shape terms (SXX,SYY,SXY)
+    bool useReff = modelPSF->class->useReff;
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
+        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
+        return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+        bool status = true;
+        status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL);
+        status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL);
+        if (!status) {
+            psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)",
+                     in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            modelPSF->flags |= PM_MODEL_STATUS_LIMITS;
+        }
+    }
+
+    return true;
+}
+
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+        if (i == PM_PAR_SKY) continue;
+        if (i == PM_PAR_I0) continue;
+        if (i == PM_PAR_XPOS) continue;
+        if (i == PM_PAR_YPOS) continue;
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
+    }
+
+    // the 2D PSF model fits polarization terms (E0,E1,E2)
+    // convert to shape terms (SXX,SYY,SXY)
+    // XXX user-defined value for limit?
+    bool useReff = model->class->useReff;
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
+        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+        return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+        bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+        if (!status) {
+            psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+            model->flags |= PM_MODEL_STATUS_LIMITS;
+        }
+    }
+    return(true);
+}
+
+bool PM_MODEL_FIT_STATUS (pmModel *model)
+{
+    bool  status;
+
+    psF32 *PAR  = model->params->data.F32;
+    psF32 *dPAR = model->dparams->data.F32;
+
+    status = true;
+    status &= (PAR[PM_PAR_I0] > 0);
+    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
+
+    return status;
+}
+
+
+void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
+{
+    switch (type) {
+      case PM_MODEL_LIMITS_NONE:
+        paramsMinUse = NULL;
+        paramsMaxUse = NULL;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_IGNORE:
+        paramsMinUse = NULL;
+        paramsMaxUse = NULL;
+        limitsApply = false;
+        break;
+      case PM_MODEL_LIMITS_LAX:
+        paramsMinUse = paramsMinLax;
+        paramsMaxUse = paramsMaxLax;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_MODERATE:
+        paramsMinUse = paramsMinModerate;
+        paramsMaxUse = paramsMaxModerate;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_STRICT:
+        paramsMinUse = paramsMinStrict;
+        paramsMaxUse = paramsMaxStrict;
+        limitsApply = true;
+        break;
+      default:
+        psAbort("Unrecognised model limits type: %x", type);
+    }
+    return;
+}
Index: /trunk/psModules/src/objects/models/pmModel_HSC_V1.h
===================================================================
--- /trunk/psModules/src/objects/models/pmModel_HSC_V1.h	(revision 40491)
+++ /trunk/psModules/src/objects/models/pmModel_HSC_V1.h	(revision 40491)
@@ -0,0 +1,16 @@
+#ifndef PM_MODEL_HSC_V1_H
+
+#include "pmModel.h"
+
+psF32 pmModelFunc_HSC_V1(psVector *deriv, const psVector *params, const psVector *pixcoord);
+bool pmModelLimits_HSC_V1(psMinConstraintMode mode, int nParam, float *params, float *beta);
+bool pmModelGuess_HSC_V1(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
+psF64 pmModelFlux_HSC_V1(const psVector *params);
+psF64 pmModelRadius_HSC_V1(const psVector *params, psF64 flux);
+psF64 pmModelSetFWHM_HSC_V1(const psVector *params, psF64 flux);
+bool pmModelFromPSF_HSC_V1(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf);
+bool  pmModelParamsFromPSF_HSC_V1(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io);
+bool pmModelFitStatus_HSC_V1(pmModel *model);
+void pmModelSetLimits_HSC_V1(pmModelLimitsType type);
+
+#endif
Index: /trunk/psModules/src/objects/pmModelClass.c
===================================================================
--- /trunk/psModules/src/objects/pmModelClass.c	(revision 40490)
+++ /trunk/psModules/src/objects/pmModelClass.c	(revision 40491)
@@ -47,4 +47,5 @@
 # include "models/pmModel_QGAUSS.h"
 # include "models/pmModel_PS1_V1.h"
+# include "models/pmModel_HSC_V1.h"
 # include "models/pmModel_RGAUSS.h"
 # include "models/pmModel_SERSIC.h"
@@ -63,4 +64,5 @@
     {"PS_MODEL_DEV",          7, 1, (pmModelFunc)pmModelFunc_DEV,     (pmModelFlux)pmModelFlux_DEV,     (pmModelRadius)pmModelRadius_DEV,     (pmModelSetFWHM)pmModelSetFWHM_DEV,     (pmModelLimits)pmModelLimits_DEV,     (pmModelGuessFunc)pmModelGuess_DEV,    (pmModelFromPSFFunc)pmModelFromPSF_DEV,    (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV,    (pmModelFitStatusFunc)pmModelFitStatus_DEV,    (pmModelSetLimitsFunc)pmModelSetLimits_DEV    },
     {"PS_MODEL_TRAIL",        7, 0, (pmModelFunc)pmModelFunc_TRAIL,   (pmModelFlux)pmModelFlux_TRAIL,   (pmModelRadius)pmModelRadius_TRAIL,   (pmModelSetFWHM)pmModelSetFWHM_TRAIL,   (pmModelLimits)pmModelLimits_TRAIL,   (pmModelGuessFunc)pmModelGuess_TRAIL,  (pmModelFromPSFFunc)pmModelFromPSF_TRAIL,  (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL,  (pmModelFitStatusFunc)pmModelFitStatus_TRAIL,  (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL  },
+    {"PS_MODEL_HSC_V1",       8, 0, (pmModelFunc)pmModelFunc_HSC_V1,  (pmModelFlux)pmModelFlux_HSC_V1,  (pmModelRadius)pmModelRadius_HSC_V1,  (pmModelSetFWHM)pmModelSetFWHM_HSC_V1,  (pmModelLimits)pmModelLimits_HSC_V1,  (pmModelGuessFunc)pmModelGuess_HSC_V1, (pmModelFromPSFFunc)pmModelFromPSF_HSC_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_HSC_V1, (pmModelFitStatusFunc)pmModelFitStatus_HSC_V1, (pmModelSetLimitsFunc)pmModelSetLimits_HSC_V1 },
 };
 
Index: /trunk/psModules/src/objects/pmSourceOutputs.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceOutputs.c	(revision 40490)
+++ /trunk/psModules/src/objects/pmSourceOutputs.c	(revision 40491)
@@ -135,4 +135,7 @@
 		outputs->psfCore = PAR[PM_PAR_7];
 	    }
+	    if (model->type == pmModelClassGetType ("PS_MODEL_HSC_V1")) {
+		outputs->psfCore = PAR[PM_PAR_7];
+	    }
 
 	    outputs->psfMajorFWHM = model->class->modelSetFWHM(model->params, axes.major);