Fixed world coordinates in PSF-fitting.

Fixed error ellipse measurements in PSF-fitting (but works only for PSF_NMAX=1 currently). Simplified PSF parameter handling (removed PSFDISPLAY_TYPE configuration option). Simplified triggering of local north axis computation. Allowed model axis ratios to exceed 1 to facilitate convergence. Fixed alignment bug for symmetric models. Restricted maximum allowed range for model coordinates. Pushed version number to 2.10.0.

Fixed world coordinates in PSF-fitting.
Fixed error ellipse measurements in PSF-fitting (but works only for PSF_NMAX=1 currently). Simplified PSF parameter handling (removed PSFDISPLAY_TYPE configuration option). Simplified triggering of local north axis computation. Allowed model axis ratios to exceed 1 to facilitate convergence. Fixed alignment bug for symmetric models. Restricted maximum allowed range for model coordinates. Pushed version number to 2.10.0.
4d8f6615 · Emmanuel Bertin · 98adcc78 · 4d8f6615 · 4d8f6615 · 4d8f6615
Commit 4d8f6615 authored Dec 16, 2009 by Emmanuel Bertin
--- a/configure
+++ b/configure
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for sextractor 2.9.9.
+# Generated by GNU Autoconf 2.63 for sextractor 2.10.0.
 #
 # 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.9.9'
+PACKAGE_VERSION='2.10.0'
-PACKAGE_STRING='sextractor 2.9.9'
+PACKAGE_STRING='sextractor 2.10.0'
 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.9.9 to adapt to many kinds of systems.
+\`configure' configures sextractor 2.10.0 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.9.9:";;
+     short | recursive ) echo "Configuration of sextractor 2.10.0:";;
   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.9.9
+sextractor configure 2.10.0
 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.9.9, which was
+It was created by sextractor $as_me 2.10.0, 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.9.9'
+ VERSION='2.10.0'
 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.9.9, which was
+This file was extended by sextractor $as_me 2.10.0, 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.9.9
+sextractor config.status 2.10.0
 configured by $0, generated by GNU Autoconf 2.63,
  with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"

--- a/configure.ac
+++ b/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.9.9, [bertin@iap.fr])
+AC_INIT(sextractor, 2.10.0, [bertin@iap.fr])
 AC_CONFIG_SRCDIR(src/makeit.c)
 AC_CONFIG_AUX_DIR(autoconf)
 AM_CONFIG_HEADER(config.h)

--- a/man/sex.1
+++ b/man/sex.1
-.TH SEXTRACTOR "1" "December 2009" "SWarp 2.9.9" "User Commands"
+.TH SEXTRACTOR "1" "December 2009" "SWarp 2.10.0" "User Commands"
 .SH NAME
 sex \- extract a source catalog from an astronomical FITS image
 .SH SYNOPSIS

--- a/src/analyse.c
+++ b/src/analyse.c
@@ -9,7 +9,7 @@
 *
 *	Contents:	analyse(), endobject()...: measurements on detections.
 *
-*	Last modify:	01/12/2009
+*	Last modify:	15/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -488,6 +488,7 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
 	|| FLAG(obj2.win_mx2w)
 	|| FLAG(obj2.poserr_mx2w)
 	|| FLAG(obj2.winposerr_mx2w)
+	|| FLAG(obj2.poserrmx2w_psf)
 	|| FLAG(obj2.poserrmx2w_prof)
 	|| FLAG(obj2.prof_flagw)
 	|| ((!prefs.pixel_scale) && FLAG(obj2.area_flagw)))
@@ -671,30 +672,9 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
      else
        psf_fit(thepsf, field, wfield, obj);
      obj2->npsf = thepsfit->npsf;
-      if (prefs.psfdisplay_type == PSFDISPLAY_SPLIT)
+      nsub = thepsfit->npsf;
-        {
+      if (nsub<1)
-        nsub = thepsfit->npsf;
+        nsub = 1;
-        if (nsub<1)
-          nsub = 1;
-        }
-      else
-        for (j=0; j<thepsfit->npsf; j++)
-          {
-          if (FLAG(obj2.x_psf) && j<prefs.psf_xsize)
-            obj2->x_psf[j] = thepsfit->x[j];
-          if (FLAG(obj2.y_psf) && j<prefs.psf_ysize)
-            obj2->y_psf[j] = thepsfit->y[j];
-          if (FLAG(obj2.flux_psf) && j<prefs.psf_fluxsize)
-            obj2->flux_psf[j] = thepsfit->flux[j];
-          if (FLAG(obj2.magerr_psf) && j<prefs.psf_magerrsize)
-            obj2->magerr_psf[j] = obj2->fluxerr_psf[j]>0.0?
-		1.086*obj2->fluxerr_psf[j]/thepsfit->flux[j] : 99.0;
-          if (FLAG(obj2.fluxerr_psf) && j<prefs.psf_fluxerrsize)
-            obj2->fluxerr_psf[j] = obj2->fluxerr_psf[j];     
-          if (FLAG(obj2.mag_psf) && j<prefs.psf_magsize)
-            obj2->mag_psf[j] = thepsfit->flux[j]>0.0?
-		prefs.mag_zeropoint -2.5*log10(thepsfit->flux[j]) : 99.0;
-          }
      }
 /*----------------------------- Profile fitting -----------------------------*/
@@ -726,23 +706,26 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
 /*-- Go through each newly identified component */
    for (j=0; j<nsub; j++)
      {
-      if (prefs.psf_flag && prefs.psfdisplay_type == PSFDISPLAY_SPLIT)
+      if (prefs.psf_flag)
        {
-        if (FLAG(obj2.x_psf))
+        obj2->x_psf = thepsfit->x[j];
-          obj2->x_psf[0] = thepsfit->x[j];
+        obj2->y_psf = thepsfit->y[j];
-        if (FLAG(obj2.y_psf))
+        if (FLAG(obj2.xf_psf) || FLAG(obj2.xw_psf))
-          obj2->y_psf[0] = thepsfit->y[j];
+          astrom_psfpos(field, obj);
+/*------ Express position error parameters in the FOCAL or WORLD frame */
+        if (FLAG(obj2.poserrmx2w_psf))
+          astrom_psferrparam(field, obj);
        if (FLAG(obj2.flux_psf))
-          obj2->flux_psf[0] = thepsfit->flux[j]>0.0? thepsfit->flux[j]:0.0; /*?*/
+          obj2->flux_psf = thepsfit->flux[j]>0.0? thepsfit->flux[j]:0.0; /*?*/
        if (FLAG(obj2.mag_psf))
-          obj2->mag_psf[0] = thepsfit->flux[j]>0.0?
+          obj2->mag_psf = thepsfit->flux[j]>0.0?
 		prefs.mag_zeropoint -2.5*log10(thepsfit->flux[j]) : 99.0;
-        if (FLAG(obj2.magerr_psf))
-          obj2->magerr_psf[0]=
-		(thepsfit->flux[j]>0.0 && obj2->fluxerr_psf[j]>0.0) ? /*?*/
-			1.086*obj2->fluxerr_psf[j]/thepsfit->flux[j] : 99.0;
        if (FLAG(obj2.fluxerr_psf))
-          obj2->fluxerr_psf[0]= obj2->fluxerr_psf[j];
+          obj2->fluxerr_psf= thepsfit->fluxerr[j];
+        if (FLAG(obj2.magerr_psf))
+          obj2->magerr_psf =
+		(thepsfit->flux[j]>0.0 && thepsfit->fluxerr[j]>0.0) ? /*?*/
+			1.086*thepsfit->fluxerr[j]/thepsfit->flux[j] : 99.0;
        if (j)
          obj->number = ++thecat.ntotal;
        }

--- a/src/astrom.c
+++ b/src/astrom.c
@@ -9,7 +9,7 @@
 *
 *	Contents:	Astrometrical computations.
 *
-*	Last modify:	27/08/2009
+*	Last modify:	15/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -46,15 +46,7 @@ void	initastrom(picstruct *field)
 /* Test if the WCS is in use */
  if (wcs->lng != wcs->lat)
    {
-    if (FLAG(obj2.theta2000) || FLAG(obj2.theta1950)
+    if (FLAG(obj2.dtheta2000) || FLAG(obj2.dtheta1950))
-	|| FLAG(obj2.poserr_theta2000) || FLAG(obj2.poserr_theta1950)
-	|| FLAG(obj2.win_theta2000) || FLAG(obj2.win_theta1950)
-	|| FLAG(obj2.winposerr_theta2000) || FLAG(obj2.winposerr_theta1950)
-	|| FLAG(obj2.prof_spheroid_theta2000)
-	|| FLAG(obj2.prof_spheroid_theta1950)
-	|| FLAG(obj2.prof_disk_theta2000) || FLAG(obj2.prof_disk_theta1950)
-//	|| FLAG(obj2.prof_arms_theta2000) || FLAG(obj2.prof_arms_theta1950)
-	|| FLAG(obj2.prof_bar_theta2000) || FLAG(obj2.prof_bar_theta1950))
      {
      if (fabs(wcs->equinox-2000.0)>0.003)
        precess(wcs->equinox, 0.0, 90.0, 2000.0, &wcs->ap2000, &wcs->dp2000);
@@ -278,6 +270,61 @@ void	astrom_winpos(picstruct *field, objstruct *obj)
  }
+/****************************** astrom_psfpos *******************************/
+/*
+Compute real FOCAL and WORLD PSF coordinates according to FITS info.
+*/
+void	astrom_psfpos(picstruct *field, objstruct *obj)
+  {
+   wcsstruct	*wcs;
+   double	rawpos[NAXIS], wcspos[NAXIS];
+   int		lng,lat;
+  wcs = field->wcs;
+  lng = wcs->lng;
+  lat = wcs->lat;
+  if (FLAG(obj2.xf_psf))
+    {
+    rawpos[0] = obj2->x_psf;
+    rawpos[1] = obj2->y_psf;
+    raw_to_red(wcs, rawpos, wcspos);
+    obj2->xf_psf = wcspos[0];
+    obj2->yf_psf = wcspos[1];
+    }
+  if (FLAG(obj2.xw_psf))
+    {
+    rawpos[0] = obj2->x_psf;
+    rawpos[1] = obj2->y_psf;
+    raw_to_wcs(wcs, rawpos, wcspos);
+    obj2->xw_psf = wcspos[0];
+    obj2->yw_psf = wcspos[1];
+    if (lng != lat)
+      {
+      obj2->alphas_psf = lng<lat? obj2->xw_psf : obj2->yw_psf;
+      obj2->deltas_psf = lng<lat? obj2->yw_psf : obj2->xw_psf;
+      if (FLAG(obj2.alpha2000_psf))
+        {
+        if (fabs(wcs->equinox-2000.0)>0.003)
+          precess(wcs->equinox, wcspos[0], wcspos[1],
+		2000.0, &obj2->alpha2000_psf, &obj2->delta2000_psf);
+        else
+          {
+          obj2->alpha2000_psf = lng<lat? obj2->xw_psf : obj2->yw_psf;
+          obj2->delta2000_psf = lng<lat? obj2->yw_psf : obj2->xw_psf;
+          }
+        if (FLAG(obj2.alpha1950_psf))
+          j2b(wcs->equinox, obj2->alpha2000_psf, obj2->delta2000_psf,
+		&obj2->alpha1950_psf, &obj2->delta1950_psf);
+        }
+      }
+    }
+  return;
+  }
 /****************************** astrom_profpos *******************************/
 /*
 Compute real FOCAL and WORLD profit coordinates according to FITS info.
@@ -644,6 +691,84 @@ void	astrom_winerrparam(picstruct *field, objstruct *obj)
  }
+/***************************** astrom_psferrparam ***************************/
+/*
+Compute error ellipse parameters in WORLD and SKY coordinates.
+*/
+void	astrom_psferrparam(picstruct *field, objstruct *obj)
+  {
+   wcsstruct	*wcs;
+   double	dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
+   int		lng,lat, naxis;
+  wcs = field->wcs;
+  naxis = wcs->naxis;
+  lng = wcs->lng;
+  lat = wcs->lat;
+  if (lng == lat)
+    {
+    lng = 0;
+    lat = 1;
+    }
+  lm0 = obj2->jacob[lng+naxis*lng];
+  lm1 = obj2->jacob[lat+naxis*lng];
+  lm2 = obj2->jacob[lng+naxis*lat];
+  lm3 = obj2->jacob[lat+naxis*lat];
+/* All WORLD params based on 2nd order moments have to pass through here */
+  dx2 = obj2->poserrmx2_psf;
+  dy2 = obj2->poserrmy2_psf;
+  dxy = obj2->poserrmxy_psf;
+  obj2->poserrmx2w_psf = xm2 = lm0*lm0*dx2+lm1*lm1*dy2+lm0*lm1*dxy;
+  obj2->poserrmy2w_psf = ym2 = lm2*lm2*dx2+lm3*lm3*dy2+lm2*lm3*dxy;
+  obj2->poserrmxyw_psf = xym = lm0*lm2*dx2+lm1*lm3*dy2+(lm0*lm3+lm1*lm2)*dxy;
+  temp=xm2-ym2;
+  if (FLAG(obj2.poserrthetaw_psf))
+    {
+    obj2->poserrthetaw_psf = (fmod_m90_p90(temp==0.0)?
+				(45.0):(0.5*atan2(2.0*xym,temp)/DEG));
+/*-- Compute position angles in J2000 or B1950 reference frame */
+    if (wcs->lng != wcs->lat)
+      {
+      if (FLAG(obj2.poserrthetas_psf))
+        obj2->poserrthetas_psf = fmod_m90_p90(lng<lat?
+		((obj2->poserrthetaw_psf>0.0?90:-90.0)-obj2->poserrthetaw_psf)
+		: obj2->poserrthetaw_psf);
+      if (FLAG(obj2.poserrtheta2000_psf))
+        obj2->poserrtheta2000_psf = fmod_m90_p90(obj2->poserrthetas_psf
+				+ obj2->dtheta2000);
+      if (FLAG(obj2.poserrtheta1950_psf))
+        obj2->poserrtheta1950_psf = fmod_m90_p90(obj2->poserrthetas_psf
+				+ obj2->dtheta1950);
+      }
+    }
+  if (FLAG(obj2.poserraw_psf))
+    {
+    temp = sqrt(0.25*temp*temp+xym*xym);
+    pm2 = 0.5*(xm2+ym2);
+    obj2->poserraw_psf = (float)sqrt(pm2+temp);
+    obj2->poserrbw_psf = (float)sqrt(pm2-temp);
+    }
+  if (FLAG(obj2.poserrcxxw_psf))
+    {
+/*-- Handle large, fully correlated profiles (can cause a singularity...) */
+    if ((temp=xm2*ym2-xym*xym)<1e-6)
+      {
+      temp = 1e-6;
+      xym *= 0.99999;
+      }
+    obj2->poserrcxxw_psf = (float)(ym2/temp);
+    obj2->poserrcyyw_psf = (float)(xm2/temp);
+    obj2->poserrcxyw_psf = (float)(-2*xym/temp);
+    }
+  return;
+  }
 /***************************** astrom_proferrparam ***************************/
 /*
 Compute error ellipse parameters in WORLD and SKY coordinates.

--- a/src/astrom.h
+++ b/src/astrom.h
@@ -9,7 +9,7 @@
 *
 *	Contents:	Astrometrical stuff.
 *
-*	Last modify:	20/07/2009
+*	Last modify:	15/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -35,6 +35,8 @@ extern void		astrom_errparam(picstruct *, objstruct *),
 			astrom_proferrparam(picstruct *, objstruct *),
 			astrom_profpos(picstruct *, objstruct *),
 			astrom_profshapeparam(picstruct *, objstruct *),
+			astrom_psferrparam(picstruct *, objstruct *),
+			astrom_psfpos(picstruct *, objstruct *),
 			astrom_shapeparam(picstruct *, objstruct *),
 			astrom_winerrparam(picstruct *, objstruct *),
 			astrom_winpos(picstruct *, objstruct *),

--- a/src/param.h
+++ b/src/param.h
@@ -9,7 +9,7 @@
 *
 *	Contents:	parameter list for catalog data.
 *
-*	Last modify:	01/10/2009
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -738,51 +738,51 @@ keystruct	objkey[] = {
 	"phys.size.radius;instr.det", "pix",  1, &prefs.flux_radiussize},
  {"XPSF_IMAGE", "X coordinate from PSF-fitting",
-	&outobj2.x_psf, H_FLOAT, T_FLOAT, "%10.3f", "pixel",
+	&outobj2.x_psf, H_FLOAT, T_DOUBLE, "%10.3f", "pixel",
-	"pos.cartesian.x;stat.fit.param;instr.det", "pix", 1, &prefs.psf_xsize},
+	"pos.cartesian.x;stat.fit.param;instr.det", "pix"},
  {"YPSF_IMAGE", "Y coordinate from PSF-fitting",
-	&outobj2.y_psf, H_FLOAT, T_FLOAT, "%10.3f", "pixel",
+	&outobj2.y_psf, H_FLOAT, T_DOUBLE, "%10.3f", "pixel",
-	"pos.cartesian.y;stat.fit.param;instr.det", "pix", 1, &prefs.psf_ysize},
+	"pos.cartesian.y;stat.fit.param;instr.det", "pix"},
  {"XPSF_WORLD", "PSF position along world x axis",
 	&outobj2.xw_psf, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
-	"pos.eq.ra;stat.fit.param", "deg", 1, &prefs.psf_xwsize},
+	"pos.eq.ra;stat.fit.param", "deg"},
  {"YPSF_WORLD", "PSF position along world y axis",
 	&outobj2.yw_psf, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
-	"pos.eq.dec;stat.fit.param", "deg", 1, &prefs.psf_ywsize},
+	"pos.eq.dec;stat.fit.param", "deg"},
  {"ALPHAPSF_SKY", "Right ascension of the fitted PSF (native)",
 	&outobj2.alphas_psf, H_FLOAT, T_DOUBLE, "%11.7f", "deg",
-	"pos.eq.ra;stat.fit.param", "deg", 1, &prefs.psf_alphassize},
+	"pos.eq.ra;stat.fit.param", "deg"},
  {"DELTAPSF_SKY", "Declination of the fitted PSF (native)",
 	&outobj2.deltas_psf, H_FLOAT, T_DOUBLE, "%+11.7f", "deg",
-	"pos.eq.dec;stat.fit.param", "deg", 1, &prefs.psf_deltassize},
+	"pos.eq.dec;stat.fit.param", "deg"},
  {"ALPHAPSF_J2000", "Right ascension of the fitted PSF (J2000)",
 	&outobj2.alpha2000_psf, H_FLOAT, T_DOUBLE, "%11.7f", "deg",
-	"pos.eq.ra;stat.fit.param", "deg", 1, &prefs.psf_alpha2000size},
+	"pos.eq.ra;stat.fit.param", "deg"},
  {"DELTAPSF_J2000", "Declination of the fitted PSF (J2000)",
 	&outobj2.delta2000_psf, H_FLOAT, T_DOUBLE, "%+11.7f", "deg",
-	"pos.eq.dec;stat.fit.param", "deg", 1, &prefs.psf_delta2000size},
+	"pos.eq.dec;stat.fit.param", "deg"},
  {"ALPHAPSF_B1950", "Right ascension of the fitted PSF (B1950)",
 	&outobj2.alpha1950_psf, H_FLOAT, T_DOUBLE, "%11.7f", "deg",
-	"pos.eq.ra;stat.fit.param", "deg", 1, &prefs.psf_alpha1950size},
+	"pos.eq.ra;stat.fit.param", "deg"},
  {"DELTAPSF_B1950", "Declination of the fitted PSF (B1950)",
 	&outobj2.delta1950_psf, H_FLOAT, T_DOUBLE, "%+11.7f", "deg",
-	"pos.eq.dec;stat.fit.param", "deg", 1, &prefs.psf_delta1950size},
+	"pos.eq.dec;stat.fit.param", "deg"},
  {"FLUX_PSF", "Flux from PSF-fitting",
 	&outobj2.flux_psf, H_FLOAT, T_FLOAT, "%12.7g", "count",
-	"phot.count;stat.fit.param", "ct", 1, &prefs.psf_fluxsize},
+	"phot.count;stat.fit.param", "ct"},
  {"FLUXERR_PSF", "RMS flux error for PSF-fitting",
 	&outobj2.fluxerr_psf, H_FLOAT, T_FLOAT, "%12.7g", "count",
-	"stat.stdev;phot.count", "ct", 1, &prefs.psf_fluxerrsize},
+	"stat.stdev;phot.count", "ct"},
  {"MAG_PSF", "Magnitude from PSF-fitting",
 	&outobj2.mag_psf, H_FLOAT, T_FLOAT, "%8.4f", "mag",
-	"phot.mag", "mag", 1, &prefs.psf_magsize},
+	"phot.mag", "mag"},
  {"MAGERR_PSF", "RMS magnitude error from PSF-fitting",
 	&outobj2.magerr_psf, H_FLOAT, T_FLOAT, "%8.4f", "mag",
-	"stat.stdev;phot.mag", "mag", 1, &prefs.psf_magsize},
+	"stat.stdev;phot.mag", "mag"},
  {"NITER_PSF", "Number of iterations for PSF-fitting",
 	&outobj2.niter_psf, H_INT, T_SHORT, "%3d", "",

--- a/src/preflist.h
+++ b/src/preflist.h
@@ -9,7 +9,7 @@
 *
 *	Contents:	Keywords for the configuration file.
 *
-*	Last modify:	24/09/2009
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -125,8 +125,6 @@
  {"PSF_NAME", P_STRINGLIST, prefs.psf_name, 0,0, 0.0,0.0,
   {""}, 1, 2, &prefs.npsf_name},	/*?*/
  {"PSF_NMAX", P_INT, &prefs.psf_npsfmax, 1, PSF_NPSFMAX},
-  {"PSFDISPLAY_TYPE", P_KEY, &prefs.psfdisplay_type, 0,0, 0.0,0.0,
-   {"SPLIT","VECTOR",""}},
  {"SATUR_KEY", P_STRING, prefs.satur_key},
  {"SATUR_LEVEL", P_FLOAT, &prefs.satur_level, 0,0, -1e+30, 1e+30},
  {"SEEING_FWHM", P_FLOAT, &prefs.seeing_fwhm, 0,0, 1e-10, 1e+10},
@@ -288,7 +286,6 @@ char *default_prefs[] =
 "*",
 "*PSF_NAME         default.psf    # File containing the PSF model",
 "*PSF_NMAX         2              # Max.number of PSFs fitted simultaneously",
-"*PSFDISPLAY_TYPE  SPLIT          # Catalog type for PSF-fitting: SPLIT or VECTOR",
 "*PATTERN_TYPE     RINGS-HARMONIC # can RINGS-QUADPOLE, RINGS-OCTOPOLE,",
 "*                                # RINGS-HARMONICS or GAUSS-LAGUERRE",
 "*SOM_NAME         default.som    # File containing Self-Organizing Map weights",

--- a/src/prefs.h
+++ b/src/prefs.h
@@ -9,7 +9,7 @@
 *
 *	Contents:	Keywords for the configuration file.
 *
-*	Last modify:	24/09/2009
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -197,8 +197,6 @@ typedef struct
  char		*(psf_name[2]);				/* PSF filename */
  int		npsf_name;				/* nb of params */
  int		psf_npsfmax;				/* Max # of PSFs */
-  enum	{PSFDISPLAY_SPLIT, PSFDISPLAY_VECTOR}
-		psfdisplay_type;			/* PSF display type */
  int		psf_xsize,psf_ysize;			/* nb of params */
  int		psf_xwsize,psf_ywsize;			/* nb of params */
  int		psf_alphassize,psf_deltassize;		/* nb of params */

--- a/src/profit.c
+++ b/src/profit.c
@@ -9,7 +9,7 @@
 *
 *	Contents:	Fit an arbitrary profile combination to a detection.
 *
-*	Last modify:	01/12/2009
+*	Last modify:	15/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -173,7 +173,7 @@ OUTPUT	Pointer to an allocated fit structure (containing details about the
 	fit).
 NOTES	It is a modified version of the lm_minimize() of lmfit.
 AUTHOR	E. Bertin (IAP)
-VERSION	01/12/2009
+VERSION	14/12/2009
 ***/
 void	profit_fit(profitstruct *profit,
 		picstruct *field, picstruct *wfield,
@@ -187,7 +187,7 @@ void	profit_fit(profitstruct *profit,
    double		emx2,emy2,emxy, a , cp,sp, cn, bn, n;
    float		**list,
 			*cov,
-			psf_fwhm, dchi2, err;
+			psf_fwhm, dchi2, err, aspect;
    int			*index,
 			i,j,p, nparam, nparam2, ncomp, nprof;
@@ -455,6 +455,15 @@ the_gal++;
 /* Spheroid */
  if (FLAG(obj2.prof_spheroid_flux))
    {
+    if ((aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT]) > 1.0)
+      {
+      *profit->paramlist[PARAM_SPHEROID_REFF] *= aspect;
+      profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]] *= aspect;
+      profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]]
+			/= (aspect*aspect);
+      *profit->paramlist[PARAM_SPHEROID_ASPECT] = 1.0 / aspect;
+      *profit->paramlist[PARAM_SPHEROID_POSANG] += 90.0;
+      }
    obj2->prof_spheroid_flux = *profit->paramlist[PARAM_SPHEROID_FLUX];
    obj2->prof_spheroid_fluxerr =
 		profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX]];
@@ -498,6 +507,15 @@ the_gal++;
 /* Disk */
  if (FLAG(obj2.prof_disk_flux))
    {
+    if ((aspect = *profit->paramlist[PARAM_DISK_ASPECT]) > 1.0)
+      {
+      *profit->paramlist[PARAM_DISK_SCALE] *= aspect;
+      profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]] *= aspect;
+      profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]]
+			/= (aspect*aspect);
+      *profit->paramlist[PARAM_DISK_ASPECT] = 1.0 / aspect;
+      *profit->paramlist[PARAM_DISK_POSANG] += 90.0;
+      }
    obj2->prof_disk_flux = *profit->paramlist[PARAM_DISK_FLUX];
    obj2->prof_disk_fluxerr =
 		profit->paramerr[profit->paramindex[PARAM_DISK_FLUX]];
@@ -1876,13 +1894,13 @@ INPUT	Pointer to the profit structure,
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	25/09/2008
+VERSION	15/12/2009
 ***/
 void	profit_resetparam(profitstruct *profit, paramenum paramtype)
  {
   objstruct	*obj;
   obj2struct	*obj2;
-   float	param, parammin,parammax;
+   float	param, parammin,parammax, range;
  obj = profit->obj;
  obj2 = profit->obj2;
@@ -1897,13 +1915,19 @@ void	profit_resetparam(profitstruct *profit, paramenum paramtype)
      break;
    case PARAM_X:
      param = obj->mx - (int)(obj->mx+0.49999);
-      parammin = -obj2->hl_radius*4;
+      range = fabs(obj2->hl_radius*4.0);
-      parammax =  obj2->hl_radius*4;
+      if (range>profit->objnaxisn[0]*2.0)
+        range = profit->objnaxisn[0]*2.0;
+      parammin = -range;
+      parammax =  range;
      break;
    case PARAM_Y:
      param = obj->my - (int)(obj->my+0.49999);
-      parammin = -obj2->hl_radius*4;
+      range = fabs(obj2->hl_radius*4);
-      parammax =  obj2->hl_radius*4;
+      if (range>profit->objnaxisn[1]*2)
+        range = profit->objnaxisn[1]*2;
+      parammin = -range;
+      parammax =  range;
      break;
    case PARAM_SPHEROID_FLUX:
      param = obj2->flux_auto/2.0;
@@ -1911,7 +1935,7 @@ void	profit_resetparam(profitstruct *profit, paramenum paramtype)
      parammax = 2*obj2->flux_auto;
      break;
    case PARAM_SPHEROID_REFF:
-      param = obj2->hl_radius;
+      param = fabs(obj2->hl_radius);
      parammin = 0.1;
      parammax = param * 4.0;
      break;
@@ -1919,6 +1943,7 @@ void	profit_resetparam(profitstruct *profit, paramenum paramtype)
      param = FLAG(obj2.prof_disk_flux)? 1.0 : obj->b/obj->a;
      parammin = FLAG(obj2.prof_disk_flux)? 0.5 : 0.01;
      parammax = 1.0;
+parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0;
      break;
    case PARAM_SPHEROID_POSANG:
      param = obj->theta;
@@ -1936,14 +1961,14 @@ void	profit_resetparam(profitstruct *profit, paramenum paramtype)
      parammax = 2*obj2->flux_auto;
      break;
    case PARAM_DISK_SCALE:	/* From scalelength to Re */
-      param = obj2->hl_radius/1.67835*sqrt(obj->a/obj->b);
+      param = fabs(obj2->hl_radius)/1.67835*sqrt(obj->a/obj->b);
      parammin = param/4.0;
      parammax = param * 4.0;
      break;
    case PARAM_DISK_ASPECT:
      param = obj->b/obj->a;
      parammin = 0.01;
-      parammax = 1.0;
+      parammax = 100.0;
      break;
    case PARAM_DISK_POSANG:
      param = obj->theta;
@@ -2425,7 +2450,7 @@ INPUT	Profile structure,
 OUTPUT	Corrected flux contribution.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	09/10/2009
+VERSION	14/12/2009
 ***/
 float	prof_add(profstruct *prof, profitstruct *profit)
  {
@@ -2544,7 +2569,12 @@ float	prof_add(profstruct *prof, profitstruct *profit)
 /*---- Copy the symmetric part */
      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
        {
-        pixin2 = pixin - profit->modnaxisn[0];
+        pixin2 = pixin - profit->modnaxisn[0] - 1;
+        if (!(profit->modnaxisn[0]&1))
+          {
+          *(pixin++) = 0.0;
+          npix2--;
+          }
        for (i=npix2; i--;)
          *(pixin++) = *(pixin2--);
        }
@@ -2600,7 +2630,12 @@ float	prof_add(profstruct *prof, profitstruct *profit)
 /*---- Copy the symmetric part */
      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
        {
-        pixin2 = pixin - profit->modnaxisn[0];
+        pixin2 = pixin - profit->modnaxisn[0] - 1;
+        if (!(profit->modnaxisn[0]&1))
+          {
+          *(pixin++) = 0.0;
+          npix2--;
+          }
        for (i=npix2; i--;)
          *(pixin++) = *(pixin2--);
        }
@@ -2655,7 +2690,12 @@ float	prof_add(profstruct *prof, profitstruct *profit)
 /*---- Copy the symmetric part */
      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
        {
-        pixin2 = pixin - profit->modnaxisn[0];
+        pixin2 = pixin - profit->modnaxisn[0] - 1;
+        if (!(profit->modnaxisn[0]&1))
+          {
+          *(pixin++) = 0.0;
+          npix2--;
+          }
        for (i=npix2; i--;)
          *(pixin++) = *(pixin2--);
        }

--- a/src/psf.c
+++ b/src/psf.c
@@ -10,7 +10,7 @@
 *
 *	Contents:	Fit the PSF to a detection.
 *
-*	Last modify:	01/12/2009
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -47,13 +47,15 @@ Allocate memory and stuff for the PSF-fitting.
 void	psf_init(psfstruct *psf)
  {
  QMALLOC(thepsfit, psfitstruct, 1);
-  QMALLOC(thepsfit->x, float, prefs.psf_npsfmax);
+  QMALLOC(thepsfit->x, double, prefs.psf_npsfmax);
-  QMALLOC(thepsfit->y, float, prefs.psf_npsfmax);
+  QMALLOC(thepsfit->y, double, prefs.psf_npsfmax);
  QMALLOC(thepsfit->flux, float, prefs.psf_npsfmax);
+  QMALLOC(thepsfit->fluxerr, float, prefs.psf_npsfmax);
  QMALLOC(ppsfit, psfitstruct, 1); /*?*/
-  QMALLOC(ppsfit->x, float, prefs.psf_npsfmax);
+  QMALLOC(ppsfit->x, double, prefs.psf_npsfmax);
-  QMALLOC(ppsfit->y, float, prefs.psf_npsfmax);
+  QMALLOC(ppsfit->y, double, prefs.psf_npsfmax);
  QMALLOC(ppsfit->flux, float, prefs.psf_npsfmax);
+  QMALLOC(ppsfit->fluxerr, float, prefs.psf_npsfmax);
  return;
  }  
@@ -89,6 +91,7 @@ void	psf_end(psfstruct *psf, psfitstruct *psfit)
    free(psfit->x);
    free(psfit->y);
    free(psfit->flux);
+    free(psfit->fluxerr);
    free(psfit);
    }
@@ -315,9 +318,10 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
  static obj2struct     *obj2 = &outobj2;
  static double		x2[PSF_NPSFMAX],y2[PSF_NPSFMAX],xy[PSF_NPSFMAX],
 			deltax[PSF_NPSFMAX],
-			deltay[PSF_NPSFMAX],flux[PSF_NPSFMAX],
+			deltay[PSF_NPSFMAX],
+			flux[PSF_NPSFMAX],fluxerr[PSF_NPSFMAX],
 			deltaxb[PSF_NPSFMAX],deltayb[PSF_NPSFMAX],
-			fluxb[PSF_NPSFMAX],
+			fluxb[PSF_NPSFMAX],fluxerrb[PSF_NPSFMAX],
 			sol[PSF_NTOT], covmat[PSF_NTOT*PSF_NTOT], 
 			vmat[PSF_NTOT*PSF_NTOT], wmat[PSF_NTOT];
  float			*data, *data2, *data3, *weight, *d, *w;
@@ -328,7 +332,7 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
 			backnoise2, gain, radmin2,radmax2,satlevel,chi2,
 			r2, valmax, psf_fwhm;
  float			**psfmasks, **psfmaskx,**psfmasky,
-			*ps, *dh, *wh, pixstep,fluxerr;
+			*ps, *dh, *wh, pixstep;
  PIXTYPE		*datah, *weighth;
  int			i,j,p, npsf,npsfmax, npix, nppix, ix,iy,niter,
 			width, height, pwidth,pheight, x,y,
@@ -355,6 +359,7 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
      thepsfit->x[j] = obj2->posx;
      thepsfit->y[j] = obj2->posy;
      thepsfit->flux[j] = 0.0;
+      thepsfit->fluxerr[j] = 0.0;
    }
  /* Scale data area with object "size" */
@@ -369,19 +374,12 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
  npix = width*height;
  radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
  radmax2 = npix/2.0;
-  fluxerr = 0.0;
  /* Scale total area with PSF FWHM */
  pwidth = (int)(psf->masksize[0]*psf->pixstep)+width;;
  pheight = (int)(psf->masksize[1]*psf->pixstep)+height;
  nppix = pwidth*pheight;
-  /* Allocate working space */
-   if (prefs.psf_flag==1)
-      if (prefs.dpsf_flag!=1)
-        if(!FLAG(obj2.fluxerr_psf))
-            QMALLOC(obj2->fluxerr_psf, float, prefs.psf_npsfmax);  
  QMALLOC(weighth, PIXTYPE, npix);
  QMALLOC(weight, float, npix);
  QMALLOC(datah, PIXTYPE, npix);
@@ -453,7 +451,7 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
  npsfflag = 1;
  r2 = psf_fwhm*psf_fwhm/2.0;
-  fluxb[0] = deltaxb[0] = deltayb[0] = 0.0;
+  fluxb[0] = fluxerrb[0] = deltaxb[0] = deltayb[0] = 0.0;
  for (npsf=1; npsf<=npsfmax && npsfflag; npsf++)
    {
@@ -534,16 +532,17 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
          compute_pos( &npsf, &convflag, &npsfflag,radmin2,radmax2,
                       r2, sol,flux, deltax, deltay,&dx,&dy);
        }
-      for (j=0; j<npsf; j++)
-        {
 /*-- Compute variances and covariances */
-          svdvar(vmat, wmat, npsf*PSF_NA, covmat);
+      svdvar(vmat, wmat, npsf*PSF_NA, covmat);
-          var = covmat;
+      var = covmat;
+      for (j=0; j<npsf; j++, var += (npsf*PSF_NA+1)*PSF_NA)
+        {
 /*---- First, the error on the flux estimate */      
-          fluxerr = sqrt(*var)>0.0?  sqrt(*var):999999.0;
+          fluxerr[j] = sqrt(*var)>0.0?  sqrt(*var):999999.0;
          //if (flux[j]<12*fluxerr && j>0)
          //  npsfmax--,flux[j]=0;
-          if (flux[j]<12*fluxerr && j>0)
+          if (flux[j]<12*fluxerr[j] && j>0)
                 {
                   flux[j]=0,kill++,npsfmax--;
                   //if(j==npsfmax-1)
@@ -558,7 +557,7 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
              deltaxb[j] = deltax[j];
              deltayb[j] = deltay[j];
              fluxb[j] = flux[j];
-              obj2->fluxerr_psf[j]=fluxerr;
+              fluxerrb[j]=fluxerr[j];
            }
        }
    }
@@ -580,17 +579,14 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
        continue; 
      if (FLAG(obj2.poserrmx2_psf))
-        {
+        compute_poserr(j,covmat,sol,obj2,x2,y2,xy, npsf);
-          compute_poserr(j,var,sol,obj2,x2,y2,xy);
-        }
-      else
-        var += 3*PSF_NA+3;  
      deltax[i] = deltaxb[j];
      deltay[i] = deltayb[j];
-      flux[i++] = fluxb[j];
+      flux[i] = fluxb[j];
+      fluxerr[i++] = fluxerrb[j];
    }
  npsf = i;
  /* Compute chi2 if asked to 
@@ -651,6 +647,7 @@ void	psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
      thepsfit->x[j] = ix+deltax[j]+1.0;
      thepsfit->y[j] = iy+deltay[j]+1.0;
      thepsfit->flux[j] = flux[j];
+      thepsfit->fluxerr[j] = fluxerr[j];
    }
@@ -732,7 +729,8 @@ void    double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
 {
  static double      /* sum[PSF_NPSFMAX]*/ pdeltax[PSF_NPSFMAX],
    pdeltay[PSF_NPSFMAX],psol[PSF_NPSFMAX], pcovmat[PSF_NPSFMAX*PSF_NPSFMAX], 
-    pvmat[PSF_NPSFMAX*PSF_NPSFMAX], pwmat[PSF_NPSFMAX],pflux[PSF_NPSFMAX];
+    pvmat[PSF_NPSFMAX*PSF_NPSFMAX], pwmat[PSF_NPSFMAX],pflux[PSF_NPSFMAX],
+    pfluxerr[PSF_NPSFMAX];
    double *pmat,
     *pm, /* *pps,  *px, *py,*/
@@ -769,15 +767,13 @@ void    double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
  /* Initialize outputs */
  ppsfit->niter = 0;
  ppsfit->npsf = 0;
-  if(!FLAG(obj2.fluxerr_psf))
-    QMALLOC(obj2->fluxerr_psf, float, npsfmax);
  for (j=0; j<npsfmax; j++) 
    {
      ppsfit->x[j] = 999999.0;
      ppsfit->y[j] = 999999.0;
      ppsfit->flux[j] = 0.0;
-      obj2->fluxerr_psf[j] = 0.0;
+      ppsfit->fluxerr[j] = 0.0;
-      pdeltax[j]= pdeltay[j]=psol[j]=  pwmat[j]=pflux[j]=0.0;
+      pdeltax[j]= pdeltay[j]=psol[j]= pwmat[j]=pflux[j]=pfluxerr[j]=0.0;
    }
@@ -819,6 +815,7 @@ void    double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
      pdeltax[j] =thepsfit->x[j]-ix-1 ;
      pdeltay[j] =thepsfit->y[j]-iy-1 ;
      ppsfit->flux[j] = 0;
+      ppsfit->fluxerr[j] = 0;
    }
 /*-------------------  Now the photometry fit ---------------------*/
@@ -920,7 +917,7 @@ void    double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
        {
          svdvar(pvmat, pwmat, j, pcovmat);
          pvar = pcovmat;
-          obj2->fluxerr_psf[k]= sqrt(*pvar)>0.0 && sqrt(*pvar)<99? 
+          pfluxerr[k]= sqrt(*pvar)>0.0 && sqrt(*pvar)<99? 
            sqrt(*pvar):99;
        }
    }
@@ -971,9 +968,11 @@ void    double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
      thepsfit->x[j] = ix+pdeltax[j]+1.0;
      thepsfit->y[j] = iy+pdeltay[j]+1.0;
      thepsfit->flux[j] = pflux[j];
+      thepsfit->fluxerr[j] = pfluxerr[j];
      ppsfit->x[j] = ix+pdeltax[j]+1.0;
      ppsfit->y[j] = iy+pdeltay[j]+1.0;
      ppsfit->flux[j] = pflux[j];
+      ppsfit->fluxerr[j] = pfluxerr[j];
    }
@@ -1179,23 +1178,20 @@ void compute_pos_phot(int *pnpsf,double *sol,double *flux)
 *********/
 void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
-                    double *y2,double *xy)
+                    double *y2,double *xy, int npsf)
 {
-  double vara,covab,varb;
+  double vara,covab,varb, f2;
  /*------ Variances and covariance along x and y */
-  vara = *(var += PSF_NA+1);
+  vara = *(var += (PSF_NA*npsf+1)*(j*PSF_NA+1));
  covab = *(++var);
-  varb = *(var += PSF_NA);
+  varb = *(var += PSF_NA*npsf);
-  var += PSF_NA+1;
+  f2 = sol[PSF_NA*j];
-  obj2->poserrmx2_psf = (vara*x2[j]*x2[j]+varb*xy[j]*xy[j]
+  f2 *= f2;
-                         +2*covab*x2[j]*xy[j])/(sol[0]*sol[0]);
+  obj2->poserrmx2_psf = vara/f2;
-  obj2->poserrmy2_psf = (varb*y2[j]*y2[j]+vara*xy[j]*xy[j]
+  obj2->poserrmy2_psf = varb/f2;
-                         +2*covab*y2[j]*xy[j])/(sol[0]*sol[0]);
+  obj2->poserrmxy_psf = covab/f2;
-  obj2->poserrmxy_psf = (vara*x2[j]*xy[j]+varb*y2[j]*xy[j]
-                         +covab*(x2[j]*y2[j]+xy[j]*xy[j]))
-    /(sol[0]*sol[0]);
  /*------ If requested, translate variances to major and minor error axes... */
  if (FLAG(obj2.poserra_psf))
    {
@@ -1210,7 +1206,7 @@ void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
      pmy2 = pmx2 = 0.5*(obj2->poserrmx2_psf+obj2->poserrmy2_psf);
      pmx2+=temp;
      pmy2-=temp;
      obj2->poserra_psf = (float)sqrt(pmx2);
      obj2->poserrb_psf = (float)sqrt(pmy2);
      obj2->poserrtheta_psf = theta*180.0/PI;

--- a/src/psf.h
+++ b/src/psf.h
@@ -10,7 +10,7 @@
 *
 *	Contents:	Include file for psffit.c.
 *
-*	Last modify:	13/09/2009
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -84,8 +84,9 @@ typedef struct
  {
  int		niter;		/* Number of iterations required */
  int		npsf;		/* Number of fitted stars for this detection */
-  float		*x,*y;		/* Position derived from the PSF-fitting */
+  double	*x,*y;		/* Position derived from the PSF-fitting */
  float		*flux;		/* Flux derived from the PSF-fitting */
+  float		*fluxerr;	/* Flux error estimated from the PSF-fitting */
  }	psfitstruct;
 /*----------------------------- Global variables ----------------------------*/
@@ -100,7 +101,7 @@ extern void	compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,
 			double *pdx,double *pdy),
 		compute_pos_phot(int *pnpsf,double *sol,double *flux),
 		compute_poserr(int j,double *var,double *sol,obj2struct *obj2,
-			double *x2, double *y2,double *xy),
+			double *x2, double *y2,double *xy, int npsf),
 		psf_build(psfstruct *psf),
 		psf_end(psfstruct *psf, psfitstruct *psfit),
 		psf_init(psfstruct *psf),

--- a/src/types.h
+++ b/src/types.h
@@ -263,11 +263,12 @@ typedef struct
  float		*flux_radius;			/* f-light-radii */
  float		hl_radius;			/* Scalar half-light radius */
 /* ---- PSF-fitting */
-  float		*flux_psf;			/* Flux from PSF-fitting */
+  float		flux_psf;			/* Flux from PSF-fitting */
-  float		*fluxerr_psf;			/* RMS error on PSF flux */
+  float		fluxerr_psf;			/* RMS error on PSF flux */
-  float		*mag_psf;			/* Mag from PSF-fitting */
+  float		mag_psf;			/* Mag from PSF-fitting */
-  float		*magerr_psf;			/* RMS mag from PSF-fitting */
+  float		magerr_psf;			/* RMS mag from PSF-fitting */
-  float		*x_psf, *y_psf;			/* Coords from PSF-fitting */
+  double	x_psf, y_psf;			/* Coords from PSF-fitting */
+  double	xf_psf, yf_psf;			/* FOCAL coordinates */
  short		niter_psf;			/* # of PSF-fitting iterat. */
  short		npsf;				/* # of fitted PSFs */
  float		chi2_psf;			/* Red. chi2 of PSF-fitting */

--- a/src/xml.c
+++ b/src/xml.c
@@ -9,7 +9,7 @@
 *
 *	Contents:	XML logging.
 *
-*	Last modify:	19/12/2007
+*	Last modify:	14/12/2009
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -217,7 +217,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	19/12/2007
+VERSION	14/12/2009
 ***/
 int	write_xml_meta(FILE *file, char *error)
  {
@@ -763,11 +763,6 @@ int	write_xml_meta(FILE *file, char *error)
    fprintf(file, "   <PARAM name=\"PSF_NMax\" datatype=\"int\""
 	" ucd=\"meta.number;obs.param\" value=\"%d\"/>\n",
    	prefs.psf_npsfmax);
-    fprintf(file,
-	"   <PARAM name=\"PSFDisplay_Type\" datatype=\"char\" arraysize=\"*\""
-	" ucd=\"meta.code;obs.param\" value=\"%s\"/>\n",
-    	key[findkeys("PSFDISPLAY_TYPE", keylist,
-			FIND_STRICT)].keylist[prefs.psfdisplay_type]);
    fprintf(file,
 	"   <PARAM name=\"SOM_Name\" datatype=\"char\" arraysize=\"*\""