Fixed detection threshold missing in XML metadata output.

Fixed TFORM error message in fitshead.c (thanks to S.Guieu).
Fixed model-fitting error variance/covariance issue (bug in matrix inversion).
Added ELLIP1ERRMODEL_IMAGE, ELLIP2ERRMODEL_IMAGE, ELLIPCORRMODEL_IMAGE,
POLAR1ERRMODEL_IMAGE, POLAR2ERRMODEL_IMAGE, POLARCORRMODEL_IMAGE (more testing is needed and WORLD versions are not yet functional).
Narrowed the tolerance constraint for stopping windowed centering iterations.
Fixed incorrect reading of the DATE-OBS FITS keyword when EPOCH and MJD-OBS are missing.
Pushed version number to 2.11.11.

configure

 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for sextractor 2.11.8.
+# Generated by GNU Autoconf 2.63 for sextractor 2.11.11.
 #
 # Report bugs to <bertin@iap.fr>.
 #
@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
 # Identity of this package.
 PACKAGE_NAME='sextractor'
 PACKAGE_TARNAME='sextractor'
-PACKAGE_VERSION='2.11.8'
-PACKAGE_STRING='sextractor 2.11.8'
+PACKAGE_VERSION='2.11.11'
+PACKAGE_STRING='sextractor 2.11.11'
 PACKAGE_BUGREPORT='bertin@iap.fr'
 ac_unique_file="src/makeit.c"
@@ -1505,7 +1505,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures sextractor 2.11.8 to adapt to many kinds of systems.
+\`configure' configures sextractor 2.11.11 to adapt to many kinds of systems.
 Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1575,7 +1575,7 @@ fi
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of sextractor 2.11.8:";;
+     short | recursive ) echo "Configuration of sextractor 2.11.11:";;
    esac
   cat <<\_ACEOF
@@ -1706,7 +1706,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-sextractor configure 2.11.8
+sextractor configure 2.11.11
 generated by GNU Autoconf 2.63
 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1720,7 +1720,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
-It was created by sextractor $as_me 2.11.8, which was
+It was created by sextractor $as_me 2.11.11, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   $ $0 $@
@@ -2423,7 +2423,7 @@ fi
 # Define the identity of the package.
  PACKAGE='sextractor'
- VERSION='2.11.8'
+ VERSION='2.11.11'
 cat >>confdefs.h <<_ACEOF
@@ -28535,7 +28535,7 @@ exec 6>&1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by sextractor $as_me 2.11.8, which was
+This file was extended by sextractor $as_me 2.11.11, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   CONFIG_FILES    = $CONFIG_FILES
@@ -28598,7 +28598,7 @@ Report bugs to <bug-autoconf@gnu.org>."
 _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_version="\\
-sextractor config.status 2.11.8
+sextractor config.status 2.11.11
 configured by $0, generated by GNU Autoconf 2.63,
   with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"

configure.ac

@@ -6,7 +6,7 @@ define([AC_CACHE_LOAD],)
 define([AC_CACHE_SAVE],)
 
 # This is your standard Bertin source code...
-AC_INIT(sextractor, 2.11.8, [bertin@iap.fr])
+AC_INIT(sextractor, 2.11.11, [bertin@iap.fr])
 AC_CONFIG_SRCDIR(src/makeit.c)
 AC_CONFIG_AUX_DIR(autoconf)
 AM_CONFIG_HEADER(config.h)

man/sex.1

-.TH SEXTRACTOR "1" "July 2010" "SWarp 2.11.8" "User Commands"
+.TH SEXTRACTOR "1" "August 2010" "SWarp 2.11.11" "User Commands"
 .SH NAME
 sex \- extract a source catalog from an astronomical FITS image
 .SH SYNOPSIS

src/astrom.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	Astrometrical computations.
 *
-*	Last modify:	15/12/2009
+*	Last modify:	03/08/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -1093,12 +1093,21 @@ void	astrom_profshapeparam(picstruct *field, objstruct *obj)
       obj2->prof_cxyw = (float)(-2*xym/temp);
       }
   
-    if (FLAG(obj2.prof_e1w))
+    if (FLAG(obj2.prof_pol1w))
       {
       if (xm2+ym2 > 1.0/BIG)
         {
         obj2->prof_pol1w = (xm2 - ym2) / (xm2+ym2);
         obj2->prof_pol2w = 2.0*xym / (xm2 + ym2);
+        }
+      else
+        obj2->prof_pol1w = obj2->prof_pol2w = 0.0;
+      }
+
+    if (FLAG(obj2.prof_e1w))
+      {
+      if (xm2+ym2 > 1.0/BIG)
+        {
         temp = xm2*ym2-xym*xym;
         if (temp>=0.0)
           temp = xm2+ym2+2.0*sqrt(temp);
@@ -1108,8 +1117,7 @@ void	astrom_profshapeparam(picstruct *field, objstruct *obj)
         obj2->prof_e2w = 2.0*xym / temp;
         }
       else
-        obj2->prof_pol1w = obj2->prof_pol2w
-	= obj2->prof_e1w = obj2->prof_e2w = 0.0;
+        obj2->prof_e1w = obj2->prof_e2w = 0.0;
       }
     }
 

src/catout.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	functions for output of catalog data.
 *
-*	Last modify:	12/07/2010
+*	Last modify:	03/08/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -233,8 +233,8 @@ void	updateparamflags()
   FLAG(obj2.prof_bar_aspect) |= FLAG(obj2.prof_bar_aspecterr);
   FLAG(obj2.prof_bar_theta) |= FLAG(obj2.prof_bar_thetaerr);
   FLAG(obj2.prof_arms_mag) |= FLAG(obj2.prof_arms_magerr);
-  FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w)
-			| FLAG(obj2.prof_pol1w) | FLAG(obj2.prof_pol2w);
+  FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w);
+  FLAG(obj2.prof_pol1w) |= FLAG(obj2.prof_pol2w);
   FLAG(obj2.prof_aw) |= FLAG(obj2.prof_bw);
   FLAG(obj2.prof_cxxw) |= FLAG(obj2.prof_cyyw) | FLAG(obj2.prof_cxyw);
   FLAG(obj2.prof_thetas) |= FLAG(obj2.prof_theta1950)
@@ -245,7 +245,7 @@ void	updateparamflags()
 			| FLAG(obj2.prof_mxyw)
 			| FLAG(obj2.prof_thetaw) | FLAG(obj2.prof_aw)
 			| FLAG(obj2.prof_cxxw)
-			| FLAG(obj2.prof_e1w);
+			| FLAG(obj2.prof_e1w) | FLAG(obj2.prof_pol1w);
 
   FLAG(obj2.dtheta1950) |= FLAG(obj2.prof_theta1950)
 			| FLAG(obj2.prof_spheroid_theta1950)
@@ -313,10 +313,6 @@ void	updateparamflags()
 			| FLAG(obj2.prof_spheroid_aspecterrw)
 			| FLAG(obj2.prof_spheroid_thetaerrw);
 
- FLAG(obj2.prof_mx2w) |= FLAG(obj2.prof_my2w) | FLAG(obj2.prof_mxyw)
-			| FLAG(obj2.prof_e1w) |FLAG(obj2.prof_e2w)
-			| FLAG(obj2.prof_pol1w) |FLAG(obj2.prof_pol2w);
-
   FLAG(obj2.prof_spheroid_fluxmean) |= FLAG(obj2.prof_spheroid_mumean);
   FLAG(obj2.prof_spheroid_fluxeff) |= FLAG(obj2.prof_spheroid_mueff);
   FLAG(obj2.prof_spheroid_peak) |= FLAG(obj2.prof_spheroid_mumax)
@@ -328,16 +324,21 @@ void	updateparamflags()
 			| FLAG(obj2.prof_bar_lengthw)
 			| FLAG(obj2.prof_arms_scalew)
 			| FLAG(obj2.prof_mx2w);
+
+  FLAG(obj2.proferr_e1) |= FLAG(obj2.proferr_e2) | FLAG(obj2.profcorr_e12);
 			
   FLAG(obj2.prof_e1) |= FLAG(obj2.prof_e2)
-			| FLAG(obj2.prof_pol1) | FLAG(obj2.prof_pol2)
-			| FLAG(obj2.prof_e1w);
+			| FLAG(obj2.proferr_e1) | FLAG(obj2.prof_e1w);
+  FLAG(obj2.proferr_pol1) |= FLAG(obj2.proferr_pol2)
+			| FLAG(obj2.profcorr_pol12);
+  FLAG(obj2.prof_pol1) |= FLAG(obj2.prof_pol2)
+			| FLAG(obj2.proferr_pol1) | FLAG(obj2.prof_pol1w);
   FLAG(obj2.prof_a) |= FLAG(obj2.prof_b) | FLAG(obj2.prof_theta)
 			| FLAG(obj2.prof_aw);
   FLAG(obj2.prof_cxx) |= FLAG(obj2.prof_cyy)
 			| FLAG(obj2.prof_cxy) | FLAG(obj2.prof_cxxw);
   FLAG(obj2.prof_mx2) |= FLAG(obj2.prof_my2) | FLAG(obj2.prof_mxy)
-			| FLAG(obj2.prof_e1)
+			| FLAG(obj2.prof_e1) | FLAG(obj2.prof_pol1)
 			| FLAG(obj2.prof_a) | FLAG(obj2.prof_cxx)
 			| FLAG(obj2.prof_mx2w);
   FLAG(obj2.fluxmean_prof) |= FLAG(obj2.mumean_prof);

src/fits/fitshead.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	general functions for handling FITS file headers.
 *
-*	Last modify:	28/10/2009
+*	Last modify:	20/07/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -187,7 +187,7 @@ OUTPUT	RETURN_OK if a binary table was found and mapped, RETURN_ERROR
 	otherwise.
 NOTES	-.
 AUTHOR	E. Bertin (IAP & Leiden observatory)
-VERSION	28/10/2009
+VERSION	20/07/2010
  ***/
 int	readbintabparam_head(tabstruct *tab)
 
@@ -271,7 +271,7 @@ int	readbintabparam_head(tabstruct *tab)
         key->htype = H_STRING;
         break;
       default:
-        error(EXIT_FAILURE, "*Internal Error*: Unkwown T_TYPE for ", str);
+        error(EXIT_FAILURE, "*Error*: Unknown TFORM in ", cat->filename);
       }
 
 /*--handle the special case of multimensional arrays*/

src/fitswcs.c

@@ -9,7 +9,7 @@
 *
 *	Contents:       Read and write WCS header info.
 *
-*	Last modify:	20/02/2009
+*	Last modify:	30/07/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -315,7 +315,7 @@ INPUT	tab structure.
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	20/02/2009
+VERSION	30/07/2010
  ***/
 wcsstruct	*read_wcs(tabstruct *tab)
 
@@ -433,7 +433,7 @@ wcsstruct	*read_wcs(tabstruct *tab)
     else
       {
 /*---- Search for an observation date expressed in "civilian" format */
-      FITSREADS(buf, "DATE-OBS ", str, "");
+      FITSREADS(buf, "DATE-OBS", str, "");
       if (*str)
         {
 /*------ Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */

src/levmar/misc_core.c

@@ -543,6 +543,16 @@ LM_REAL *a, *x, *work, max, sum, tmp;
 }
 #endif /* HAVE_LAPACK */
 
+/* precision-specific definitions */
+/* Added by EB */
+#ifdef HAVE_CONFIG_H
+#include        "config.h"
+#endif
+#include        ATLAS_LAPACK_H
+#define ATLAS_POTRF	LM_CAT_(clapack_, LM_ADD_PREFIX(potrf))
+#define ATLAS_POTRI	LM_CAT_(clapack_, LM_ADD_PREFIX(potri))
+/* End added by EB */
+
 /*
  * This function computes in C the covariance matrix corresponding to a least
  * squares fit. JtJ is the approximate Hessian at the solution (i.e. J^T*J, where
@@ -564,13 +574,13 @@ LM_REAL *a, *x, *work, max, sum, tmp;
  */
 int LEVMAR_COVAR(LM_REAL *JtJ, LM_REAL *C, LM_REAL sumsq, int m, int n)
 {
-register int i;
+register int i,j;
 int rnk;
 LM_REAL fact;
 
 #ifdef HAVE_LAPACK
-   rnk=LEVMAR_PSEUDOINVERSE(JtJ, C, m);
-   if(!rnk) return 0;
+  rnk=LEVMAR_PSEUDOINVERSE(JtJ, C, m);
+  if(!rnk) return 0;
 #else
 #ifdef _MSC_VER
 #pragma message("LAPACK not available, LU will be used for matrix inversion when computing the covariance; this might be unstable at times")
@@ -581,23 +591,27 @@ LM_REAL fact;
 #endif // _MSC_VER
 
 //   rnk=LEVMAR_LUINVERSE(JtJ, C, m);
-   rnk = SVDINV(JtJ, C, m);
-   if(!rnk) return 0;
 
-//   rnk=m; /* assume full rank */
-#endif /* HAVE_LAPACK */
+  rnk = SVDINV(JtJ, C, m);
 
-//   fact=sumsq/(LM_REAL)(n-rnk);
-//   for(i=0; i<m*m; ++i)
-//     C[i]*=fact;
+  if (!rnk)
+    return 0;
 
-    fact=(LM_REAL)n/(LM_REAL)(n-rnk);
-   for(i=0; i<m*m; ++i)
-     C[i]*=fact;
+//  rnk=m; /* assume full rank */
+#endif /* HAVE_LAPACK */
+  
+  fact=(LM_REAL)n/(LM_REAL)(n-rnk);
+  for(i=0; i<m*m; ++i)
+    C[i]*=fact;
 
-   for(i=0; i<m; ++i)
-     if (C[i*(m+1)]<0.0)
-       C[i*(m+1)] = 0.0;
+  for(i=0; i<m; ++i)
+    if (C[i*(m+1)]<0.0 || fabs(C[i*(m+1)])>1e30)
+      {
+      for(j=0; j<m; ++j)
+        C[i*m+j] = 0.0;
+      for(j=0; j<m; ++j)
+        C[j*m+i] = 0.0;
+      }
 
    return rnk;
 }
@@ -832,6 +846,8 @@ VERSION	30/05/2008
 
 static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
   {
+#define MIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) < (maxarg2) ?\
+        (maxarg1) : (maxarg2))
 #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
         (maxarg1) : (maxarg2))
 #define PYTHAG(a,b)     ((at=fabs(a)) > (bt=fabs(b)) ? \
@@ -840,10 +856,9 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
 #define TOL             1.0e-6
 
-   int                  flag,i,its,j,jj,k,l,mmi,nm, nml, rank;
+   int                  flag,i,its,j,jj,k,l,nm, rank;
    LM_REAL              *vmat,*wmat,
-			*w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
-                        *rv1,*tmp,
+			*w,*rv1,*tmp,
 			c,f,h,s,x,y,z,
                         anorm, g, scale,
                         at,bt,ct,maxarg1,maxarg2,
@@ -855,249 +870,232 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
   vmat=(LM_REAL *)malloc(m*m*sizeof(LM_REAL));
   wmat=(LM_REAL *)malloc(m*sizeof(LM_REAL));
 
-  l = nm = nml = 0;	/* To avoid gcc -Wall warnings */
-  for (i=0;i<m;i++)
+  l = nm = 0;	/* To avoid gcc -Wall warnings */
+
+  for (i=0; i<m; i++)
     {
     l = i+1;
-    nml = m-l;
     rv1[i] = scale*g;
     g = s = scale = 0.0;
-    mmi = m - i;
-    ap = ap0 = a+i*(m+1);
-    for (k=mmi;k--;)
-      scale += fabs(*(ap++));
+    for (k=i; k<m; k++)
+      scale += fabs(a[k*m+i]);
     if (scale)
       {
-      for (ap=ap0,k=mmi; k--; ap++)
+      for (k=i; k<m; k++)
         {
-        *ap /= scale;
-        s += *ap**ap;
+        a[k*m+i] /= scale;
+        s += a[k*m+i]*a[k*m+i];
         }
-      f = *ap0;
+      f = a[i*m+i];
       g = -SIGN(sqrt(s),f);
-      h = f*g-s;
-      *ap0 = f-g;
-      ap10 = a+l*m+i;
-      for (j=nml;j--; ap10+=m)
+      h = f*g - s;
+      a[i*m+i] = f-g;
+      for (j=l; j<m; j++)
         {
         s = 0.0;
-        for (ap=ap0,ap1=ap10,k=mmi; k--;)
-          s += *(ap1++)**(ap++);
+        for (k=i; k<m; k++)
+          s += a[k*m+i] * a[k*m+j];
         f = s/h;
-        for (ap=ap0,ap1=ap10,k=mmi; k--;)
-          *(ap1++) += f**(ap++);
+        for (k=i; k<m; k++)
+          a[k*m+j] += f * a[k*m+i];
         }
-      for (ap=ap0,k=mmi; k--;)
-        *(ap++) *= scale;
+      for (k=i; k<m; k++)
+        a[k*m+i] *= scale;
       }
-    wmat[i] = scale*g;
+    wmat[i]= scale * g;
     g = s = scale = 0.0;
-    if (i+1 != m)
+    if (i != m-1)
       {
-      ap = ap0 = a+i+m*l;
-      for (k=nml;k--; ap+=m)
-        scale += fabs(*ap);
+      for (k=l; k<m; k++)
+        scale += fabs(a[i*m+k]);
       if (scale)
         {
-        for (ap=ap0,k=nml;k--; ap+=m)
+        for (k=l; k<m; k++)
           {
-          *ap /= scale;
-          s += *ap**ap;
+          a[i*m+k] /= scale;
+          s += a[i*m+k] * a[i*m+k];
           }
-        f=*ap0;
+        f = a[i*m+l];
         g = -SIGN(sqrt(s),f);
-        h=f*g-s;
-        *ap0=f-g;
-        rv1p = rv1+l;
-        for (ap=ap0,k=nml;k--; ap+=m)
-          *(rv1p++) = *ap/h;
-        ap10 = a+l+m*l;
-        for (j=m-l; j--; ap10++)
+        h = f*g - s;
+        a[i*m+l] = f - g;
+        for (k=l; k<m; k++)
+          rv1[k] = a[i*m+k] / h;
+        for (j=l; j<m; j++)
           {
-          for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
-            s += *ap1**ap;
-          rv1p = rv1+l;
-          for (ap1=ap10,k=nml;k--; ap1+=m)
-            *ap1 += s**(rv1p++);
+          s = 0.0;
+          for (k=l; k<m; k++)
+            s += a[j*m+k] * a[i*m+k];
+          for (k=l; k<m; k++)
+            a[j*m+k] += s * rv1[k];
           }
-        for (ap=ap0,k=nml;k--; ap+=m)
-          *ap *= scale;
+        for (k=l; k<m; k++)
+          a[i*m+k] *= scale;
         }
       }
     anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
     }
 
-  for (i=m-1;i>=0;i--)
+  for (i=m;i--;)
     {
     if (i < m-1)
       {
       if (g)
         {
-        ap0 = a+l*m+i;
-        vp0 = vmat+i*m+l;
-        vp10 = vmat+l*m+l;
-        g *= *ap0;
-        for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
-          *(vp++) = *ap/g;
-        for (j=nml; j--; vp10+=m)
+        for (j=l; j<m; j++)
+          vmat[j*m+i]=(a[i*m+j]/a[i*m+l]) / g;
+        for (j=l; j<m; j++)
           {
-          for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
-            s += *ap**(vp1++);
-          for (vp=vp0,vp1=vp10,k=nml; k--;)
-            *(vp1++) += s**(vp++);
+           s = 0.0;
+          for (k=l; k<m; k++)
+            s += a[i*m+k] * vmat[k*m+j];
+          for (k=l; k<m; k++)
+            vmat[k*m+j] += s * vmat[k*m+i];
           }
         }
-      vp = vmat+l*m+i;
-      vp1 = vmat+i*m+l;
-      for (j=nml; j--; vp+=m)
-        *vp = *(vp1++) = 0.0;
+      for (j=l; j<m; j++)
+        vmat[i*m+j] = vmat[j*m+i] =0.0;
       }
-    vmat[i*m+i]=1.0;
-    g=rv1[i];
-    l=i;
-    nml = m-l;
+    vmat[i*m+i] = 1.0;
+    g = rv1[i];
+    l = i;
     }
 
-  for (i=m; --i>=0;)
+  for (i=m; i--;)
     {
-    l=i+1;
-    nml = m-l;
-    mmi=m-i;
-    g=wmat[i];
-    ap0 = a+i*m+i;
-    ap10 = ap0 + m;
-    for (ap=ap10,j=nml;j--;ap+=m)
-      *ap=0.0;
+    l = i+1;
+    g = wmat[i];
+    for (j=l; j<m;j++)
+      a[i*m+j] = 0.0;
     if (g)
       {
-      g=1.0/g;
-      for (j=nml;j--; ap10+=m)
+      g = 1.0/g;
+      for (j=l; j<m; j++)
         {
-        for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
-              s += *(++ap)**(++ap1);
-        f = (s/(*ap0))*g;
-        for (ap=ap0,ap1=ap10,k=mmi;k--;)
-          *(ap1++) += f**(ap++);
+        s = 0.0;
+        for (k=l; k<m; k++)
+          s += a[k*m+i] * a[k*m+j];
+        f = (s / a[i*m+i]) * g;
+        for (k=i; k<m; k++)
+          a[k*m+j] += f * a[k*m+i];
         }
-      for (ap=ap0,j=mmi;j--;)
-        *(ap++) *= g;
+      for (j=i; j<m; j++)
+        a[j*m+i] *= g;
       }
     else
-      for (ap=ap0,j=mmi;j--;)
-        *(ap++)=0.0;
-    ++(*ap0);
+      for (j=i; j<m; j++)
+        a[j*m+i] = 0.0;
+      ++a[i*m+i];
     }
 
-  for (k=m; --k>=0;)
+  for (k=m; k--;)
+    {
+    for (its=0; its<30; its++)
       {
-      for (its=0;its<100;its++)
+      flag = 1;
+      for (l=k+1; l--; )
         {
-        flag=1;
-        for (l=k;l>=0;l--)
+        nm=l-1;
+        if ((double)(fabs(rv1[l])+anorm) == anorm)
           {
-          nm=l-1;
-          if (!l || fabs(rv1[l]) <= anorm*TOL)
-            {
-            flag=0;
-            break;
-            }
-          if (fabs(wmat[nm]) <= anorm*TOL)
-            break;
+          flag=0;
+          break;
           }
-        if (flag)
+        if ((double)(fabs(wmat[nm])+anorm) == anorm)
+          break;
+        }
+      if (flag)
+        {
+        c = 0.0;
+        s = 1.0;
+        for (i=l; i<=k; i++)
           {
-          c=0.0;
-          s=1.0;
-          ap0 = a+nm*m;
-          ap10 = a+l*m;
-          for (i=l; i<=k; i++,ap10+=m)
+          f = s * rv1[i];
+          rv1[i] = c * rv1[i];
+          if ((double)(fabs(f) + anorm) == anorm)
+            break;
+          g = wmat[i];
+          h = PYTHAG(f,g);
+          wmat[i] = h;
+          h = 1.0 / h;
+          c = g * h;
+          s = -f * h;
+          for (j=0; j<m; j++)
             {
-            f=s*rv1[i];
-            if (fabs(f) <= anorm*TOL)
-              break;
-            g=wmat[i];
-            h=PYTHAG(f,g);
-            wmat[i]=h;
-            h=1.0/h;
-            c=g*h;
-            s=(-f*h);
-            for (ap=ap0,ap1=ap10,j=m; j--;)
-              {
-              z = *ap1;
-              y = *ap;
-              *(ap++) = y*c+z*s;
-              *(ap1++) = z*c-y*s;
-              }
+            y = a[j*m+nm];
+            z = a[j*m+i];
+            a[j*m+nm] = y*c + z*s;
+            a[j*m+i] = z*c - y*s;
             }
           }
-        z=wmat[k];
-        if (l == k)
+        }
+      z = wmat[k];
+      if (l==k)
+        {
+        if (z < 0.0)
           {
-          if (z < 0.0)
-            {
-            wmat[k] = -z;
-            vp = vmat+k*m;
-            for (j=m; j--; vp++)
-              *vp = (-*vp);
-            }
-          break;
+          wmat[k] = -z;
+          for (j=0; j<m; j++)
+            vmat[j*m+k] = -vmat[j*m+k];
+          }
+        break;
+        }
+      if (its == 29)
+        return 0;
+      x = wmat[l];
+      nm = k-1;
+      y = wmat[nm];
+      g = rv1[nm];
+      h = rv1[k];
+      f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
+      g = PYTHAG(f, 1.0);
+      f = ((x-z)*(x+z) + h * ((y / (f+SIGN(g,f))) - h)) / x;
+      c = s = 1.0;
+      for (j=l; j<=nm; j++)
+        {
+        i = j+1;
+        g = rv1[i];
+        y = wmat[i];
+        h = s * g;
+        g = c * g;
+        z = PYTHAG(f, h);
+        rv1[j] = z;
+        c = f / z;
+        s = h / z;
+        f = x*c + g*s;
+        g = g*c - x*s;
+        h = y*s;
+        y *= c;
+        for (jj=0; jj<m; jj++)
+          {
+          x = vmat[jj*m+j];
+          z = vmat[jj*m+i];
+          vmat[jj*m+j] = x*c + z*s;
+          vmat[jj*m+i] = z*c - x*s;
           }
-        x=wmat[l];
-        nm=k-1;
-        y=wmat[nm];
-        g=rv1[nm];
-        h=rv1[k];
-        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
-        g=PYTHAG(f,1.0);
-        f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
-        c=s=1.0;
-        ap10 = a+l*m;
-        vp10 = vmat+l*m;
-        for (j=l;j<=nm;j++,ap10+=m,vp10+=m)
+        z = PYTHAG(f, h);
+        wmat[j] = z;
+        if (z)
           {
-          i=j+1;
-          g=rv1[i];
-          y=wmat[i];
-          h=s*g;
-          g=c*g;
-          z=PYTHAG(f,h);
-          rv1[j]=z;
-          c=f/z;
-          s=h/z;
-          f=x*c+g*s;
-          g=g*c-x*s;
-          h=y*s;
-          y=y*c;
-          for (vp=(vp1=vp10)+m,jj=m; jj--;)
-            {
-            z = *vp;
-            x = *vp1;
-            *(vp1++) = x*c+z*s;
-            *(vp++) = z*c-x*s;
-            }
-          z=PYTHAG(f,h);
-          wmat[j]=z;
-          if (z)
-            {
-            z=1.0/z;
-            c=f*z;
-            s=h*z;
-            }
-          f=c*g+s*y;
-          x=c*y-s*g;
-          for (ap=(ap1=ap10)+m,jj=m; jj--;)
-            {
-            z = *ap;
-            y = *ap1;
-            *(ap1++) = y*c+z*s;
-            *(ap++) = z*c-y*s;
-            }
+          z = 1.0 / z;
+          c = f*z;
+          s = h*z;
+          }
+        f = c*g + s*y;
+        x = c*y - s*g;
+        for (jj=0; jj<m; jj++)
+          {
+          y = a[jj*m+j];
+          z = a[jj*m+i];
+          a[jj*m+j] = y*c + z*s;
+          a[jj*m+i] = z*c - y*s;
           }
-        rv1[l]=0.0;
-        rv1[k]=f;
-        wmat[k]=x;
         }
+      rv1[l] = 0.0;
+      rv1[k] = f;
+      wmat[k] = x;
       }
+    }
 
   wmax=0.0;
   w = wmat;
@@ -1121,7 +1119,7 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
       {
       s = 0.0;
       for (k=0; k<m; k++)
-        s += vmat[j+k*m]*a[i+k*m]*wmat[k];
+        s += vmat[k+j*m]*a[k+i*m]*wmat[k];
       b[j+i*m] = s;
       }
 

src/paramprofit.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	Model-fitting parameter list for catalog data.
 *
-*	Last modify:	25/05/2010
+*	Last modify:	03/08/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -197,6 +197,26 @@
 	&outobj2.prof_pol2, H_FLOAT, T_FLOAT, "%10.6f", "",
 	"src.ellipticity;stat.fit;instr.det", ""},
 
+  {"ELLIP1ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
+	&outobj2.proferr_e1, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.error;src.ellipticity;stat.fit;instr.det", ""},
+  {"ELLIP2ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
+	&outobj2.proferr_e2, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.error;src.ellipticity;stat.fit;instr.det", ""},
+  {"ELLIPCORRMODEL_IMAGE", "Corr.coeff between ellip.components from model-fitting",
+	&outobj2.profcorr_e12, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
+
+  {"POLAR1ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
+	&outobj2.proferr_pol1, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.error;src.ellipticity;stat.fit;instr.det", ""},
+  {"POLAR2ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
+	&outobj2.proferr_pol2, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.error;src.ellipticity;stat.fit;instr.det", ""},
+  {"POLARCORRMODEL_IMAGE", "Corr.coeff between polar. components from fitting",
+	&outobj2.profcorr_pol12, H_FLOAT, T_FLOAT, "%10.6f", "",
+	"stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
+
   {"X2MODEL_WORLD", "Variance along X-WORLD (alpha) from model-fitting",
 	&outobj2.prof_mx2w, H_EXPO, T_DOUBLE, "%18.10e", "deg**2",
 	"src.impactParam;stat.fit", "deg2"},
@@ -212,10 +232,10 @@
   {"ELLIP2MODEL_WORLD", "Ellipticity component from model-fitting",
 	&outobj2.prof_e2w, H_FLOAT, T_FLOAT, "%10.6f", "",
 	"src.ellipticity;stat.fit", ""},
-  {"POLAR1MODEL_WORLD", "Ellipticity component (quadratic) from model-fitting",
+  {"POLAR1MODEL_WORLD", "Polarisation component from model-fitting",
 	&outobj2.prof_pol1w, H_FLOAT, T_FLOAT, "%10.6f", "",
 	"src.ellipticity;stat.fit", ""},
-  {"POLAR2MODEL_WORLD", "Ellipticity component (quadratic) from model-fitting",
+  {"POLAR2MODEL_WORLD", "Polarisation component from model-fitting",
 	&outobj2.prof_pol2w, H_FLOAT, T_FLOAT, "%10.6f", "",
 	"src.ellipticity;stat.fit", ""},
 

src/profit.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	Fit an arbitrary profile combination to a detection.
 *
-*	Last modify:	07/07/2010
+*	Last modify:	03/08/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -2027,13 +2027,18 @@ INPUT	Profile-fitting structure,
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	07/07/2010
+VERSION	03/08/2010
  ***/
 void	 profit_moments(profitstruct *profit, obj2struct *obj2)
   {
    profstruct	*prof;
-   double	m0, mx2,my2,mxy, den, temp,temp2,pmx2,theta;
-   int		p;
+   double	*jac,*jact, *pjac,*pjact,
+		*dmx2,*dmx2t,*dmy2,*dmy2t,*dmxy,*dmxyt, *de1,*de1t,*de2,*de2t,
+		m0,invm0, mx2,my2,mxy, invden, temp,temp2,invstemp2,temp3,
+		pmx2,theta, flux, var1,var2,cov, dval;
+   float	 *covart;
+   int		findex[PROF_NPROF],
+		i,j,p, nparam;
 
 /*  hw = (float)(profit->modnaxisn[0]/2);*/
 /*  hh = (float)(profit->modnaxisn[1]/2);*/
@@ -2068,19 +2073,57 @@ void	 profit_moments(profitstruct *profit, obj2struct *obj2)
 /*  obj2->prof_my2 = my2 = my2/sum - my*my;*/
 /*  obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
 
+  nparam = profit->nparam;
+  if (FLAG(obj2.proferr_e1) || FLAG(obj2.proferr_pol1))
+    {
+/*-- Set up Jacobian matrices */
+    QCALLOC(jac, double, nparam*3);
+    QMALLOC(pjac, double, nparam*3);
+    dmx2 = jac;
+    dmy2 = jac+nparam;
+    dmxy = jac+2*nparam;
+    }
+  else
+    jac = pjac = NULL;
+
   m0 = mx2 = my2 = mxy = 0.0;
   for (p=0; p<profit->nprof; p++)
     {
     prof = profit->prof[p];
-    prof_moments(profit, prof);
-    m0 += *prof->flux;
-    mx2 += prof->mx2**prof->flux;
-    my2 += prof->my2**prof->flux;
-    mxy += prof->mxy**prof->flux;
+    findex[p] = prof_moments(profit, prof, pjac);
+    flux = *prof->flux;
+    m0 += flux;
+    mx2 += prof->mx2*flux;
+    my2 += prof->my2*flux;
+    mxy += prof->mxy*flux;
+    if (jac)
+      {
+      jact = jac;
+      pjact = pjac;
+      for (j=nparam*3; j--;)
+        *(jact++) += flux * *(pjact++);
+      dmx2[findex[p]] += prof->mx2;
+      dmy2[findex[p]] += prof->my2;
+      dmxy[findex[p]] += prof->mxy;
+      }
+    }
+  invm0 = 1.0 / m0;
+  obj2->prof_mx2 = (mx2 *= invm0);
+  obj2->prof_my2 = (my2 *= invm0);
+  obj2->prof_mxy = (mxy *= invm0);
+/* Complete the flux derivative of moments */
+  if (jac)
+    {
+    for (p=0; p<profit->nprof; p++)
+      {
+      dmx2[findex[p]] -= mx2;
+      dmy2[findex[p]] -= my2;
+      dmxy[findex[p]] -= mxy;
+      }
+    jact = jac;
+    for (j=nparam*3; j--;)
+      *(jact++) *= invm0;
     }
-  obj2->prof_mx2 = mx2 / m0;
-  obj2->prof_my2 = my2 / m0;
-  obj2->prof_mxy = mxy / m0;
 
 /* Handle fully correlated profiles (which cause a singularity...) */
   if ((temp2=mx2*my2-mxy*mxy)<0.00694)
@@ -2090,21 +2133,125 @@ void	 profit_moments(profitstruct *profit, obj2struct *obj2)
     temp2 = mx2*my2-mxy*mxy;
     }
 
-  if (FLAG(obj2.prof_e1))
+  if (FLAG(obj2.prof_pol1))
     {
     if (mx2+my2 > 1.0/BIG)
       {
       obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
       obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
+      if (jac)
+        {
+/*------ Compute Jacobian of polarisation and ellipticity */
+        if (FLAG(obj2.proferr_pol1))
+          {
+          invden = 2.0/((mx2+my2)*(mx2+my2));
+/*-------- Ellipticity */
+          dmx2t = dmx2;
+          dmy2t = dmy2;
+          dmxyt = dmxy;
+          de1t = de1 = pjac;			/* We re-use pjac */
+          de2t = de2 = pjac + nparam;		/* We re-use pjac */
+          for (j=nparam; j--;)
+            {
+            *(de1t++) = invden*(*dmx2t*my2-*dmy2t*mx2);
+            *(de2t++) = invden*(*(dmxyt++)*(mx2+my2)
+			- mxy*(*(dmx2t++)-*(dmy2t++)));
+            }
+/*-------- Use the Jacobians to compute errors */
+          covart = profit->covar;
+          var1 = 0.0;
+          for (j=0; j<nparam; j++)
+            {
+            dval = 0.0;
+            de1t = de1;
+            for (i=nparam; i--;)
+              dval += *(covart++)**(de1t++);
+            var1 += dval*de1[j];
+            }
+          obj2->proferr_pol1 = (float)sqrt(var1<0.0? 0.0: var1);
+          covart = profit->covar;
+          cov = var2 = 0.0;
+          for (j=0; j<nparam; j++)
+            {
+            dval = 0.0;
+            de2t = de2;
+            for (i=nparam; i--;)
+              dval += *(covart++)**(de2t++);
+            cov += dval*de1[j];
+            var2 += dval*de2[j];
+            }
+          obj2->proferr_pol2 = (float)sqrt(var2<0.0? 0.0: var2);
+          obj2->profcorr_pol12 = (dval=var1*var2) > 0.0? (float)(cov/sqrt(dval))
+							: 0.0;
+          }
+        }
+      }
+    else
+      obj2->prof_pol1 = obj2->prof_pol2
+	= obj2->proferr_pol1 = obj2->proferr_pol2 = obj2->profcorr_pol12 = 0.0;
+    }
+
+  if (FLAG(obj2.prof_e1))
+    {
+    if (mx2+my2 > 1.0/BIG)
+      {
       if (temp2>=0.0)
-        den = mx2+my2+2.0*sqrt(temp2);
+        invden = 1.0/(mx2+my2+2.0*sqrt(temp2));
       else
-        den = mx2+my2;
-      obj2->prof_e1 = (mx2 - my2) / den;
-      obj2->prof_e2 = 2.0*mxy / den;
+        invden = 1.0/(mx2+my2);
+      obj2->prof_e1 = (float)(invden * (mx2 - my2));
+      obj2->prof_e2 = (float)(2.0 * invden * mxy);
+      if (jac)
+        {
+/*------ Compute Jacobian of polarisation and ellipticity */
+        invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0;
+        if (FLAG(obj2.proferr_e1))
+          {
+/*-------- Ellipticity */
+          dmx2t = dmx2;
+          dmy2t = dmy2;
+          dmxyt = dmxy;
+          de1t = de1 = pjac;			/* We re-use pjac */
+          de2t = de2 = pjac + nparam;		/* We re-use pjac */
+          for (j=nparam; j--;)
+            {
+            temp3 = invden*(*dmx2t+*dmy2t
+			+ invstemp2 * (*dmx2t*my2+mx2**dmy2t-2.0*mxy**dmxyt));
+            *(de1t++) = invden*(*(dmx2t++) - *(dmy2t++) - (mx2-my2)*temp3);
+            *(de2t++) = 2.0*invden*(*(dmxyt++) - mxy*temp3);
+            }
+/*-------- Use the Jacobians to compute errors */
+          covart = profit->covar;
+          var1 = 0.0;
+          for (j=0; j<nparam; j++)
+            {
+            dval = 0.0;
+            de1t = de1;
+            for (i=nparam; i--;)
+              dval += *(covart++)**(de1t++);
+            var1 += dval*de1[j];
+            }
+          obj2->proferr_e1 = (float)sqrt(var1<0.0? 0.0: var1);
+          covart = profit->covar;
+          cov = var2 = 0.0;
+          for (j=0; j<nparam; j++)
+            {
+            dval = 0.0;
+            de2t = de2;
+            for (i=nparam; i--;)
+              dval += *(covart++)**(de2t++);
+            cov += dval*de1[j];
+            var2 += dval*de2[j];
+            }
+          obj2->proferr_e2 = (float)sqrt(var2<0.0? 0.0: var2);
+          obj2->profcorr_e12 = (dval=var1*var2) > 0.0? (float)(cov/sqrt(dval))
+							: 0.0;
+          }
+        }
       }
     else
-      obj2->prof_pol1 = obj2->prof_pol2 = obj2->prof_e1 = obj2->prof_e2 = 0.0;
+      obj2->prof_e1 = obj2->prof_e2
+	= obj2->proferr_e1 = obj2->proferr_e2 = obj2->profcorr_e12 = 0.0;
     }
 
   if (FLAG(obj2.prof_cxx))
@@ -2128,6 +2275,10 @@ void	 profit_moments(profitstruct *profit, obj2struct *obj2)
     obj2->prof_theta = theta*180.0/PI;
     }
 
+/* Free memory used by Jacobians */
+  free(jac);
+  free(pjac);
+
   return;
   }
 
@@ -2695,13 +2846,12 @@ INPUT	Pointer to the profit structure,
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	08/07/2010
+VERSION	27/07/2010
  ***/
 void	profit_covarunboundtobound(profitstruct *profit,
 				double *dcovar, float *covar)
   {
-   double	*dxdy,
-		dxmin,dxmax;
+   double	*dxdy;
    float	*x,*xmin,*xmax;
    int		*fflag,
 		f,f1,f2, nfree, p,p1,p2, nparam;
@@ -2718,12 +2868,7 @@ void	profit_covarunboundtobound(profitstruct *profit,
     if (fflag[p])
       {
       if (xmin[p]!=xmax[p])
-        {
-        dxmin = x[p] - xmin[p];
-        dxmax= xmax[p] - x[p];
-        dxdy[f++] = (fabs(dxmin) < 1.0/BIG && fabs(dxmax) < 1.0/BIG) ?
-		0.0 : dxmin*dxmax/(dxmin+dxmax);
-        }
+        dxdy[f++] = (x[p] - xmin[p]) * (xmax[p] - x[p]) / (xmax[p] - xmin[p]);
       else
         dxdy[f++] = xmax[p];
       }
@@ -3546,61 +3691,160 @@ width = 3.0;
 
 
 /****** prof_moments **********************************************************
-PROTO	void	prof_moments(profitstruct *profit, profstruct *prof)
-PURPOSE	Computes (analytically or numerically) the 2nd moments of a profile
-	model.
+PROTO	int	prof_moments(profitstruct *profit, profstruct *prof)
+PURPOSE	Computes (analytically or numerically) the 2nd moments of a profile.
 INPUT	Profile-fitting structure,
-	profile structure.
-OUTPUT	-.
-NOTES	-.
+	profile structure,
+	optional pointer to 3xnparam Jacobian matrix.
+OUTPUT	Index to the profile flux for further processing.
+NOTES	.
 AUTHOR	E. Bertin (IAP)
-VERSION	08/07/2010
+VERSION	21/07/2010
  ***/
-void	prof_moments(profitstruct *profit, profstruct *prof)
+int	prof_moments(profitstruct *profit, profstruct *prof, double *jac)
   {
-   double	m20,m02, ct,st, bn, n;
+   double	*dmx2,*dmy2,*dmxy,
+		m20, a2, ct,st, mx2fac, my2fac,mxyfac, dc2,ds2,dcs,
+		bn,bn2, n,n2, nfac,nfac2, hscale2, dmdn;
+   int		nparam, index;
 
+  if (jac)
+/*-- Clear output Jacobian */
+    {
+    nparam = profit->nparam;
+    memset(jac, 0, nparam*3*sizeof(double));
+    dmx2 = jac;
+    dmy2 = jac + nparam;
+    dmxy = jac + 2*nparam;
+    }
   if (prof->posangle)
     {
+    a2 = *prof->aspect**prof->aspect;
+    ct = cos(*prof->posangle*DEG);
+    st = sin(*prof->posangle*DEG);
+    mx2fac = ct*ct + st*st*a2;
+    my2fac = st*st + ct*ct*a2;
+    mxyfac = ct*st * (1.0 - a2);
+    if (jac)
+      {
+      dc2 = -2.0*ct*st*DEG;
+      ds2 =  2.0*ct*st*DEG;
+      dcs = (ct*ct - st*st)*DEG;
+      }
     switch(prof->code)
       {
       case PROF_SERSIC:
         n = fabs(*prof->extra[0]);
         bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
 		+ 131.0/(1148175*n*n*n);	/* Ciotti & Bertin 1999 */
-        m20 = 0.5 * *prof->scale**prof->scale
-		* prof_gamma(4.0*n) / (prof_gamma(2.0*n) * pow(bn, 2.0*n));
+        nfac  = prof_gamma(4.0*n) / (prof_gamma(2.0*n)*pow(bn, 2.0*n));
+        hscale2 = 0.5 * *prof->scale**prof->scale;
+        m20 = hscale2 * nfac;
+        if (jac)
+          {
+          dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= *prof->scale * nfac * mx2fac;
+          dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= *prof->scale * nfac * my2fac;
+          dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= *prof->scale * nfac * mxyfac;
+          n2 = n+0.01;
+          bn2 = 2.0*n2 - 1.0/3.0 + 4.0/(405.0*n2) + 46.0/(25515.0*n2*n2)
+		+ 131.0/(1148175*n2*n2*n2);	/* Ciotti & Bertin 1999 */
+          nfac2 = prof_gamma(4.0*n2) / (prof_gamma(2.0*n2)*pow(bn2, 2.0*n2));
+          dmdn = 100.0 * hscale2 * (nfac2-nfac);
+          dmx2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mx2fac;
+          dmy2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * my2fac;
+          dmxy[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mxyfac;
+          dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= 2.0 * m20 * st*st * *prof->aspect;
+          dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= 2.0 * m20 * ct*ct * *prof->aspect;
+          dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= -2.0 * m20 * ct*st * *prof->aspect;
+          dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
+          dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
+          dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
+          }
+        index = profit->paramindex[PARAM_SPHEROID_FLUX];
         break; 
       case PROF_DEVAUCOULEURS:
         m20 = 10.83995 * *prof->scale**prof->scale;
+        if (jac)
+          {
+          dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= 21.680 * *prof->scale * mx2fac;
+          dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= 21.680 * *prof->scale * my2fac;
+          dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
+			= 21.680 * *prof->scale * mxyfac;
+          dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= 2.0 * m20 * st*st * *prof->aspect;
+          dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= 2.0 * m20 * ct*ct * *prof->aspect;
+          dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			= -2.0 * m20 * ct*st * *prof->aspect;
+          dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
+          dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
+          dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
+          }
+        index = profit->paramindex[PARAM_SPHEROID_FLUX];
         break;
       case PROF_EXPONENTIAL:
         m20 = 3.0 * *prof->scale**prof->scale;
+        if (jac)
+          {
+          dmx2[profit->paramindex[PARAM_DISK_SCALE]]
+			= 6.0 * *prof->scale * mx2fac;
+          dmy2[profit->paramindex[PARAM_DISK_SCALE]]
+			= 6.0 * *prof->scale * my2fac;
+          dmxy[profit->paramindex[PARAM_DISK_SCALE]]
+			= 6.0 * *prof->scale * mxyfac;
+          dmx2[profit->paramindex[PARAM_DISK_ASPECT]]
+			= 2.0 * m20 * st*st * *prof->aspect;
+          dmy2[profit->paramindex[PARAM_DISK_ASPECT]]
+			= 2.0 * m20 * ct*ct * *prof->aspect;
+          dmxy[profit->paramindex[PARAM_DISK_ASPECT]]
+			= -2.0 * m20 * ct*st * *prof->aspect;
+          dmx2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (dc2 + ds2*a2);
+          dmy2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (ds2 + dc2*a2);
+          dmxy[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (1.0 - a2) * dcs;
+          }
+        index = profit->paramindex[PARAM_DISK_FLUX];
         break;
       case PROF_ARMS:
+        m20 = 1.0;
+        index = profit->paramindex[PARAM_ARMS_FLUX];
+        break;
       case PROF_BAR:
+        m20 = 1.0;
+        index = profit->paramindex[PARAM_BAR_FLUX];
+        break;
       case PROF_INRING:
+        m20 = 1.0;
+        index = profit->paramindex[PARAM_INRING_FLUX];
+        break;
       case PROF_OUTRING:
-        m20 = 1.0 / 1.0;
+        m20 = 1.0;
+        index = profit->paramindex[PARAM_OUTRING_FLUX];
         break;
       default:
         m20 = 0.0;		/* to avoid gcc -Wall warnings */
+        index = 0;		/* to avoid gcc -Wall warnings */
         error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ",
 		"prof_moments()");
       break;
       }
 
-    m02 = m20**prof->aspect**prof->aspect;
-    ct = cos(*prof->posangle*DEG);
-    st = sin(*prof->posangle*DEG);
-    prof->mx2 = ct*ct*m20 + st*st*m02;
-    prof->my2 = st*st*m20 + ct*ct*m02;
-    prof->mxy = ct*st*(m20 - m02);
+    prof->mx2 = m20*mx2fac;
+    prof->my2 = m20*my2fac;
+    prof->mxy = m20*mxyfac;
+    
     }
   else
     prof->mx2 = prof->my2 = prof->mxy = 0.0;
 
-  return;
+  return index;
   }
 
 

src/profit.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	Include file for profit.c.
 *
-*	Last modify:	08/07/2010
+*	Last modify:	21/07/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -178,13 +178,14 @@ float		*profit_compresi(profitstruct *profit, float dynparam,
 int		profit_copyobjpix(profitstruct *profit, picstruct *field,
 			picstruct *wfield),
 		profit_minimize(profitstruct *profit, int niter),
+		prof_moments(profitstruct *profit, profstruct *prof,
+				double *jac),
 		profit_resample(profitstruct *profit, float *inpix,
 			PIXTYPE *outpix, float factor),
 		profit_setparam(profitstruct *profit, paramenum paramtype,
 			float param, float parammin, float parammax);
 
 void		prof_end(profstruct *prof),
-		prof_moments(profitstruct *profit, profstruct *prof),
 		profit_addparam(profitstruct *profit, paramenum paramindex,
 			float **param),
 		profit_boundtounbound(profitstruct *profit,

src/types.h

@@ -9,7 +9,7 @@
 *
 *	Contents:	global type definitions.
 *
-*	Last modify:	23/03/2010
+*	Last modify:	21/07/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -355,7 +355,11 @@ typedef struct
   float		prof_cxx, prof_cyy,
 		prof_cxy;			/* Model-fitting ellip. params*/
   float		prof_pol1, prof_pol2;		/* Model-fitting pol. vector*/
+  float		proferr_pol1, proferr_pol2,
+		profcorr_pol12;			/* Model-fitting pol. errors */
   float		prof_e1, prof_e2;		/* Model-fitting ellip.vector*/
+  float		proferr_e1, proferr_e2,
+		profcorr_e12;			/* Model-fitting ellip. errors */
   double	prof_mx2w, prof_my2w,
 		prof_mxyw;			/* WORLD model-fitting moments*/
   float		prof_aw, prof_bw,

src/winpos.h

@@ -19,7 +19,7 @@
 #define	WINPOS_NITERMAX	16	/* Maximum number of steps */
 #define	WINPOS_NSIG	4	/* Measurement radius */
 #define	WINPOS_OVERSAMP	11	/* oversampling in each dimension */
-#define	WINPOS_STEPMIN	0.001	/* Minimum change in position for continueing*/
+#define	WINPOS_STEPMIN	0.0001	/* Minimum change in position for continueing*/
 #define	WINPOS_FAC	2.0	/* Centroid offset factor (2 for a Gaussian) */
 
 /* NOTES:

src/xml.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	XML logging.
 *
-*	Last modify:	25/05/2010
+*	Last modify:	03/08/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -219,7 +219,7 @@ INPUT	Pointer to the output file (or stream),
 OUTPUT	RETURN_OK if everything went fine, RETURN_ERROR otherwise.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	05/02/2010
+VERSION	03/08/2010
  ***/
 int	write_xml_meta(FILE *file, char *error)
   {
@@ -388,8 +388,8 @@ int	write_xml_meta(FILE *file, char *error)
 	xmlstack[n].ident[0], xmlstack[n].ident[1],
 	xmlstack[n].backmean[0], xmlstack[n].backmean[1],
 	xmlstack[n].backsig[0], xmlstack[n].backsig[1],
-	xmlstack[n].sigfac[0], xmlstack[n].sigfac[1],
 	xmlstack[n].thresh[0], xmlstack[n].thresh[1],
+	xmlstack[n].sigfac[0], xmlstack[n].sigfac[1],
 	xmlstack[n].pixscale[0], xmlstack[n].pixscale[1],
 	xmlstack[n].epoch[0], xmlstack[n].epoch[1],
 	xmlstack[n].gain[0], xmlstack[n].gain[1],
@@ -412,8 +412,8 @@ int	write_xml_meta(FILE *file, char *error)
 	xmlstack[n].ident[0],
 	xmlstack[n].backmean[0],
 	xmlstack[n].backsig[0],
-	xmlstack[n].sigfac[0],
 	xmlstack[n].thresh[0],
+	xmlstack[n].sigfac[0],
 	xmlstack[n].pixscale[0],
         xmlstack[n].epoch[0],
 	xmlstack[n].gain[0],