Added new ASSOCCOORD_TYPE configuration parameter for ASSOCiating detections...

Added new ASSOCCOORD_TYPE configuration parameter for ASSOCiating detections either in PIXEL or WORLD coordinates.
Fixed sampling issues for Sersic, de Vaucouleurs and exponential model-fitting.
Version number pushed to 2.13.2.

configure

 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for sextractor 2.13.1.
+# Generated by GNU Autoconf 2.63 for sextractor 2.13.2.
 #
 # 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.13.1'
-PACKAGE_STRING='sextractor 2.13.1'
+PACKAGE_VERSION='2.13.2'
+PACKAGE_STRING='sextractor 2.13.2'
 PACKAGE_BUGREPORT='bertin@iap.fr'
 ac_unique_file="src/makeit.c"
@@ -1508,7 +1508,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.13.1 to adapt to many kinds of systems.
+\`configure' configures sextractor 2.13.2 to adapt to many kinds of systems.
 Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1578,7 +1578,7 @@ fi
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of sextractor 2.13.1:";;
+     short | recursive ) echo "Configuration of sextractor 2.13.2:";;
    esac
   cat <<\_ACEOF
@@ -1711,7 +1711,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-sextractor configure 2.13.1
+sextractor configure 2.13.2
 generated by GNU Autoconf 2.63
 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1725,7 +1725,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.13.1, which was
+It was created by sextractor $as_me 2.13.2, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   $ $0 $@
@@ -2428,7 +2428,7 @@ fi
 # Define the identity of the package.
  PACKAGE='sextractor'
- VERSION='2.13.1'
+ VERSION='2.13.2'
 cat >>confdefs.h <<_ACEOF
@@ -28593,7 +28593,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.13.1, which was
+This file was extended by sextractor $as_me 2.13.2, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   CONFIG_FILES    = $CONFIG_FILES
@@ -28656,7 +28656,7 @@ Report bugs to <bug-autoconf@gnu.org>."
 _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_version="\\
-sextractor config.status 2.13.1
+sextractor config.status 2.13.2
 configured by $0, generated by GNU Autoconf 2.63,
   with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"

configure.ac

@@ -7,7 +7,7 @@
 #
 #	This file part of:	SExtractor
 #
-#	Copyright:		(C) 2002-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+#	Copyright:		(C) 2002-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 #
 #	License:		GNU General Public License
 #
@@ -22,7 +22,7 @@
 #	You should have received a copy of the GNU General Public License
 #	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 #
-#	Last modified:		14/10/2010
+#	Last modified:		24/01/2011
 #
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -31,7 +31,7 @@ define([AC_CACHE_LOAD],)
 define([AC_CACHE_SAVE],)
 
 # This is your standard Bertin source code...
-AC_INIT(sextractor, 2.13.1, [bertin@iap.fr])
+AC_INIT(sextractor, 2.13.2, [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" "December 2010" "SExtractor 2.13.1" "User Commands"
+.TH SEXTRACTOR "1" "January 2011" "SExtractor 2.13.2" "User Commands"
 .SH NAME
 sex \- extract a source catalogue from an astronomical FITS image
 .SH SYNOPSIS

src/analyse.c

@@ -438,7 +438,7 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
 
 /* Association */
   if (prefs.assoc_flag)
-    obj2->assoc_number = do_assoc(field, obj2->sposx, obj2->sposy);
+    obj2->assoc_number = do_assoc(field, obj2->posx, obj2->posy);
 
   if (prefs.assoc_flag && prefs.assocselec_type!=ASSOCSELEC_ALL)
     selecflag = (prefs.assocselec_type==ASSOCSELEC_MATCHED)?

src/assoc.c

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1997-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1997-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		11/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -39,6 +39,7 @@
 #include	"globals.h"
 #include	"prefs.h"
 #include	"assoc.h"
+#include	"fitswcs.h"
 
 /********************************* comp_assoc ********************************/
 /*
@@ -46,10 +47,10 @@ Comparison function for sort_assoc().
 */
 int	comp_assoc(const void *i1, const void *i2)
   {
-   float	*f1,*f2;
+   double	*f1,*f2;
 
-  f1 = (float *)i1 + 1;
-  f2 = (float *)i2 + 1;
+  f1 = (double *)i1 + 1;
+  f2 = (double *)i2 + 1;
   if (*f1<*f2)
     return -1;
   else return (*f1==*f2)?0:1;
@@ -64,13 +65,13 @@ void  sort_assoc(picstruct *field, assocstruct *assoc)
 
   {
    int		comp_assoc(const void *i1, const void *i2);
-   float	*list, rad;
+   double	*list, rad;
    int		i,j, step,nobj, *hash;
 
   step = assoc->ncol;
   nobj = assoc->nobj;
   list = assoc->list;
-  qsort(assoc->list, assoc->nobj, step*sizeof(float), comp_assoc);
+  qsort(assoc->list, assoc->nobj, step*sizeof(double), comp_assoc);
 /* Build the hash table that contains the first object in the sorted list */
 /* which may be in reach from the current scanline */
   QMALLOC(assoc->hash, int, field->height);
@@ -94,13 +95,13 @@ void  sort_assoc(picstruct *field, assocstruct *assoc)
 Read an assoc-list, and returns a pointer to the new assoc struct (or NULL if
 no list was found).
 */
-assocstruct  *load_assoc(char *filename)
+assocstruct  *load_assoc(char *filename, wcsstruct *wcs)
 
   {
    assocstruct	*assoc;
    FILE		*file;
-   float	*list, val;
-   char		str[MAXCHAR], str2[MAXCHAR], *sstr;
+   double	*list, val;
+   char		str[MAXCHARL], str2[MAXCHARL], *sstr;
    int		*data,
 		i,ispoon,j,k,l, ncol, ndata, nlist, size,spoonsize,
 		xindex,yindex,mindex;
@@ -114,7 +115,7 @@ assocstruct  *load_assoc(char *filename)
   ispoon = ncol = ndata = nlist = size = spoonsize = xindex = yindex
 	= mindex = 0;
   NFPRINTF(OUTPUT, "Reading ASSOC input-list...");
-  for (i=0; fgets(str, MAXCHAR, file);)
+  for (i=0; fgets(str, MAXCHARL, file);)
     {
 /*-- Examine current input line (discard empty and comment lines) */
     if (!*str || strchr("#\t\n",*str))
@@ -171,14 +172,14 @@ assocstruct  *load_assoc(char *filename)
 
 /*---- Allocate memory for the ASSOC struct and the filtered list */
       QMALLOC(assoc, assocstruct, 1);
-      ispoon = ASSOC_BUFINC/(nlist*sizeof(float));
+      ispoon = ASSOC_BUFINC/(nlist*sizeof(double));
       spoonsize = ispoon*nlist;
-      QMALLOC(assoc->list, float, size = spoonsize);
+      QMALLOC(assoc->list, double, size = spoonsize);
       list = assoc->list;
       }
     else  if (!(i%ispoon))
       {
-      QREALLOC(assoc->list, float, size += spoonsize);
+      QREALLOC(assoc->list, double, size += spoonsize);
       list = assoc->list + i*nlist;
       }
 
@@ -192,7 +193,7 @@ assocstruct  *load_assoc(char *filename)
     *(list+2) = 0.0;
     for (sstr = str, j=0; j<ncol; j++)
       {
-      val = (float)strtod(sstr, &sstr);
+      val = (double)strtod(sstr, &sstr);
       if (j==xindex)
         *list = val;
       else if (j==yindex)
@@ -202,12 +203,14 @@ assocstruct  *load_assoc(char *filename)
       if ((k=data[j]))
         *(list+2+k) = val;
       }
+    if (wcs)
+      wcs_to_raw(wcs, list, list);
     list += nlist;
     }
 
   fclose(file);
   free(data);
-  QREALLOC(assoc->list, float, i*nlist);
+  QREALLOC(assoc->list, double, i*nlist);
   assoc->nobj = i;
   assoc->radius = prefs.assoc_radius;
   assoc->ndata = ndata;
@@ -227,7 +230,9 @@ void	init_assoc(picstruct *field)
    assocstruct	*assoc;
 
 /* Load the assoc-list */
-  if (!(assoc = field->assoc = load_assoc(prefs.assoc_name)))
+  if (!(assoc = field->assoc = load_assoc(prefs.assoc_name,
+				prefs.assoccoord_type==ASSOCCOORD_WORLD?
+					field->wcs : NULL)))
     error(EXIT_FAILURE, "*Error*: Assoc-list file not found: ",
 	prefs.assoc_name);
 
@@ -264,17 +269,16 @@ void	end_assoc(picstruct *field)
 /*
 Perform the association task for a source and return the number of IDs.
 */
-int	do_assoc(picstruct *field, float x, float y)
+int	do_assoc(picstruct *field, double x, double y)
   {
    assocstruct	*assoc;
-   double	aver;
-   float	dx,dy, dist, rad, rad2, comp, wparam,
+   double	aver, dx,dy, dist, rad, rad2, comp, wparam,
 		*list, *input, *data;
    int		h, step, i, flag, iy, nobj;
 
   assoc = field->assoc;
 /* Need to initialize the array */
-  memset(assoc->data, 0, prefs.assoc_size*sizeof(float));
+  memset(assoc->data, 0, prefs.assoc_size*sizeof(double));
   aver = 0.0;
 
   if (prefs.assoc_type == ASSOC_MIN || prefs.assoc_type == ASSOC_NEAREST)
@@ -304,7 +308,7 @@ int	do_assoc(picstruct *field, float x, float y)
       input = list+3;
       if (prefs.assoc_type == ASSOC_FIRST)
         {
-        memcpy(assoc->data, input, assoc->ndata*sizeof(float));
+        memcpy(assoc->data, input, assoc->ndata*sizeof(double));
         return 1;
         }
       wparam = *(list+2);
@@ -314,7 +318,7 @@ int	do_assoc(picstruct *field, float x, float y)
         case ASSOC_NEAREST:
           if (dist<comp)
             {
-            memcpy(data, input, assoc->ndata*sizeof(float));
+            memcpy(data, input, assoc->ndata*sizeof(double));
             comp = dist;
             }
           break;
@@ -340,14 +344,14 @@ int	do_assoc(picstruct *field, float x, float y)
         case ASSOC_MIN:
           if (wparam<comp)
             {
-            memcpy(data, input, assoc->ndata*sizeof(float));
+            memcpy(data, input, assoc->ndata*sizeof(double));
             comp = wparam;
             }
           break;
         case ASSOC_MAX:
           if (wparam>comp)
             {
-            memcpy(data, input, assoc->ndata*sizeof(float));
+            memcpy(data, input, assoc->ndata*sizeof(double));
             comp = wparam;
             }
           break;

src/assoc.h

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1997-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1997-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,30 +22,34 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		11/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
+#ifndef _FITSWCS_H_
+#include        "fitswcs.h"
+#endif
+
 #define		ASSOC_BUFINC	131072	/* Assoc buffer increment (bytes) */
 
 /*--------------------------------- typedefs --------------------------------*/
 
 typedef struct structassoc
   {
-  float		*list;			/* Pointer to the list of data */
+  double	*list;			/* Pointer to the list of data */
   int		nobj;			/* Number of data rows */
   int		ncol;			/* Total number of columns per row */
   int		ndata;			/* Number of retained cols per row */
   int		*hash;			/* Pointer to the hash table */
-  float		*data;			/* Copy of current parameters */
-  float		radius;			/* Radius of search for association */
+  double	*data;			/* Copy of current parameters */
+  double	radius;			/* Radius of search for association */
   }             assocstruct;
 
 /*------------------------------ Prototypes ---------------------------------*/
 
-assocstruct	*load_assoc(char *filename);
+assocstruct	*load_assoc(char *filename, wcsstruct *wcs);
 
-int		do_assoc(picstruct *field, float x, float y);
+int		do_assoc(picstruct *field, double x, double y);
 
 void		init_assoc(picstruct *field),
 		end_assoc(picstruct *field),

src/define.h

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		11/10/2010
+*	Last modified:		09/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -38,7 +38,7 @@
 #define		BANNER		"SExtractor"
 #define		MYVERSION	VERSION
 #define		EXECUTABLE	"sex"
-#define         COPYRIGHT       "2010 IAP/CNRS/UPMC"
+#define         COPYRIGHT       "2011 IAP/CNRS/UPMC"
 #define		DISCLAIMER	BANNER " comes with ABSOLUTELY NO WARRANTY\n" \
 		"You may redistribute copies of " BANNER "\n" \
 		"under the terms of the GNU General Public License."
@@ -140,6 +140,7 @@
 /*------------------------------- Other Macros -----------------------------*/
 
 #define	DEXP(x)		exp(2.30258509299*(x))		/* 10^x */
+#define	DEXPF(x)	expf(2.30258509299f*(x))	/* 10^x */
 
 #define QFREAD(ptr, size, afile, fname) \
 		if (fread(ptr, (size_t)(size), (size_t)1, afile)!=1) \

src/param.h

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		11/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -723,7 +723,7 @@ keystruct	objkey[] = {
 	&outobj2.vigshift, H_FLOAT, T_FLOAT, "%12.7g", "count",
 	"obs.image", "ct", 2, prefs.vigshiftsize},
   {"VECTOR_ASSOC", "ASSOCiated parameter vector",
-	&outobj2.assoc, H_FLOAT, T_FLOAT, "%12.7g", "",
+	&outobj2.assoc, H_FLOAT, T_DOUBLE, "%12.7g", "",
 	"src", "", 1, &prefs.assoc_size},
   {"NUMBER_ASSOC", "Number of ASSOCiated IDs",
 	&outobj2.assoc_number, H_INT, T_LONG, "%10d", "",

src/preflist.h

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		14/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -54,6 +54,8 @@
   {"ASSOC_TYPE", P_KEY, &prefs.assoc_type, 0,0, 0.0,0.0,
    {"FIRST", "NEAREST", "MEAN", "MAG_MEAN", "SUM", "MAG_SUM",
    "MIN", "MAX", ""}},
+  {"ASSOCCOORD_TYPE", P_KEY, &prefs.assoccoord_type, 0,0, 0.0,0.0,
+   {"PIXEL","WORLD",""}},
   {"ASSOCSELEC_TYPE", P_KEY, &prefs.assocselec_type, 0,0, 0.0,0.0,
    {"ALL","MATCHED","-MATCHED",""}},
   {"BACK_FILTERSIZE", P_INTLIST, prefs.backfsize, 1,11, 0.0,0.0,
@@ -275,6 +277,7 @@ char *default_prefs[] =
 "*ASSOC_NAME       sky.list       # name of the ASCII file to ASSOCiate",
 "*ASSOC_DATA       2,3,4          # columns of the data to replicate (0=all)",
 "*ASSOC_PARAMS     2,3,4          # columns of xpos,ypos[,mag]",
+"*ASSOCCOORD_TYPE  PIXEL          # ASSOC coordinates: PIXEL or WORLD",
 "*ASSOC_RADIUS     2.0            # cross-matching radius (pixels)",
 "*ASSOC_TYPE       NEAREST        # ASSOCiation method: FIRST, NEAREST, MEAN,",
 "*                                # MAG_MEAN, SUM, MAG_SUM, MIN or MAX",

src/prefs.h

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
+*	Copyright:		(C) 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		14/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -167,6 +167,8 @@ typedef struct
   enum	{ASSOC_FIRST, ASSOC_NEAREST, ASSOC_MEAN, ASSOC_MAGMEAN,
 	 ASSOC_SUM, ASSOC_MAGSUM, ASSOC_MIN, ASSOC_MAX}
 		assoc_type;				/* type of assoc. */
+  enum	{ASSOCCOORD_PIXEL, ASSOCCOORD_WORLD}
+		assoccoord_type;		       	/* type of coords */
   enum	{ASSOCSELEC_ALL, ASSOCSELEC_MATCHED, ASSOCSELEC_NOMATCHED}
 		assocselec_type;		       	/* type of assoc. */
   double	assoc_radius;				/* ASSOC range */

src/profit.c

@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		19/10/2010
+*	Last modified:		24/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -160,7 +160,6 @@ profitstruct	*profit_init(psfstruct *psf)
   QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
   QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
   profit->nprof = nprof;
-  profit->oversamp = PROFIT_OVERSAMP;
   profit->fluxfac = 1.0;	/* Default */
 
   return profit;
@@ -1140,7 +1139,7 @@ INPUT	Pointer to the profit structure involved in the fit,
 OUTPUT	Number of iterations used.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	11/10/2010
+VERSION	24/01/2010
  ***/
 int	profit_minimize(profitstruct *profit, int niter)
   {
@@ -1152,11 +1151,11 @@ int	profit_minimize(profitstruct *profit, int niter)
   memset(dcovar, 0, profit->nparam*profit->nparam*sizeof(double));
 
 /* Perform fit */
-  lm_opts[0] = 1.0e-2;
-  lm_opts[1] = 1.0e-12;
-  lm_opts[2] = 1.0e-12;
-  lm_opts[3] = 1.0e-12;
-  lm_opts[4] = 1.0e-4;
+  lm_opts[0] = 1.0e-3;		/* Initial mu */
+  lm_opts[1] = 1.0e-8;		/* ||J^T e||_inf stopping factor */
+  lm_opts[2] = 1.0e-8;		/* |Dp||_2 stopping factor */
+  lm_opts[3] = 1.0e-8;		/* ||e||_2 stopping factor */
+  lm_opts[4] = 1.0e-4;		/* Jacobian step */
 
   profit_boundtounbound(profit, profit->paraminit, dparam, PARAM_ALLPARAMS);
 
@@ -1231,7 +1230,7 @@ INPUT	Pointer to the vector of parameters,
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	07/10/2010
+VERSION	24/01/2011
  ***/
 void	profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
   {
@@ -1266,7 +1265,7 @@ void	profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
     if (f>0 && q==1)
       jflag = 1;
     }
-
+jflag = 0;	/* Temporarily deactivated (until problems are fixed) */
   if (jflag && !(profit->nprof==1 && profit->prof[0]->code == PROF_DIRAC))
     {
     prof = profit->prof;
@@ -2454,7 +2453,7 @@ INPUT	Pointer to the profile-fitting structure,
 OUTPUT	-.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	08/07/2010
+VERSION	24/01/2011
  ***/
 void	 profit_surface(profitstruct *profit, obj2struct *obj2)
   {
@@ -2466,7 +2465,6 @@ void	 profit_surface(profitstruct *profit, obj2struct *obj2)
    int		i,p, imax, npix, neff;
 
 /* Allocate "high-definition" raster only to make measurements */
-  hdprofit.oversamp = PROFIT_OVERSAMP;
   hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES;
   npix = hdprofit.nmodpix = hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1];
 /* Find best image size factor from fitting results */
@@ -2603,7 +2601,6 @@ endcheck(check);
     imax = npix-1 - imax;
 /*-- Recompute hi-def model raster without oversampling */
 /*-- and with the same flux correction factor */
-    hdprofit.oversamp = 0;
     memset(hdprofit.modpix,0, npix*sizeof(float));
     for (p=0; p<profit->nprof; p++)
       prof_add(&hdprofit, profit->prof[p], 1);
@@ -3323,11 +3320,12 @@ INPUT	Profile-fitting structure,
 OUTPUT	Total (asymptotic) flux contribution.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	08/10/2010
+VERSION	24/01/2011
  ***/
 float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
   {
-   double	xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n;
+   double	xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n,
+		dx1cout,dx2cout, ddx1[36],ddx2[36];
    float	posin[PROFIT_MAXEXTRA], posout[2], dnaxisn[2],
 		*pixin, *pixin2, *pixout,
 		fluxfac, amp,cd11,cd12,cd21,cd22, dcd11,dcd21, dx1,dx2,
@@ -3338,9 +3336,13 @@ float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
 		width, invwidth2,
 		r,r2,rmin, r2minxin,r2minxout, rmax, r2max,
 		r2max1, r2max2, r2min, invr2xdif,
-		val, theta, thresh, ra,rb,rao, num,num2,den;
-   int		npix, noversamp, threshflag,
-		d,e,i, ix1,ix2, idx1,idx2, nx2, npix2;
+		val, theta, thresh, ra,rb,rao, num,num2,den, ang,angstep,
+		invn, smoothfac, dr,deltar, krpinvn,dkrpinvn, rs,rs2,
+		a11,a12,a21,a22, invdet, dca,dsa, a0,a2,a3, p1,p2,
+		krspinvn, ekrspinvn, selem;
+   int		npix, threshflag,
+		a,d,e,i, ix1,ix2, ix1max,ix2max, ir,nang, idx1,idx2, nx2,
+		npix2;
 
   npix = profit->nmodpix;
 
@@ -3389,27 +3391,55 @@ float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
     num2 = x2max*x2max*num;
     r2max2 = num2<PROFIT_MAXR2MAX*den? num2 / den : PROFIT_MAXR2MAX;
     r2max = (r2max1 < r2max2? r2max1 : r2max2);
+    rmax = sqrtf(r2max);
     }
 
   switch(prof->code)
     {
     case PROF_DIRAC:
-      memset(prof->pix, 0, npix*sizeof(float));
       prof->pix[profit->modnaxisn[0]/2
 		+ (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0;
       prof->lostfluxfrac = 0.0;
       threshflag = 0;
       break;
     case PROF_SERSIC:
+    case PROF_DEVAUCOULEURS:
+    case PROF_EXPONENTIAL:
+/*---- Compute sharp/smooth transition radius */
+      rs = PROFIT_SMOOTHR*(xscale>yscale?xscale:yscale);
+      if (rs<=0)
+        rs = 1.0;
+      rs2 = rs*rs;
+/*---- The consequence of sampling on flux is compensated by PSF normalisation*/
+      if (prof->code==PROF_EXPONENTIAL)
+        bn = n = 1.0;
+      else if (prof->code==PROF_DEVAUCOULEURS)
+        {
+        n = 4.0;
+        bn = 7.66924944;
+        }
+      else
+        {
         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 */
-      k = -bn;
+        }
+      invn = 1.0/n;
       hinvn = 0.5/n;
-/*---- The consequence of sampling on flux is compensated by PSF normalisation*/
+      k = -bn;
+/*---- Compute central polynomial terms */
+      krspinvn = prof->code==PROF_EXPONENTIAL? -rs : k*expf(logf(rs)*invn);
+      ekrspinvn = expf(krspinvn);
+      p2 = krspinvn*invn*invn;
+      p1 = krspinvn*p2;
+      a0 = (1+(1.0/6.0)*(p1+(1.0-5.0*n)*p2))*ekrspinvn;
+      a2 = (-1.0/2.0)*(p1+(1.0-3.0*n)*p2)/rs2*ekrspinvn;
+      a3 = (1.0/3.0)*(p1+(1.0-2.0*n)*p2)/(rs2*rs)*ekrspinvn;
+/*---- Compute the smooth part of the profile */
       x10 = -x1cout - dx1;
       x2 = -x2cout - dx2;
       pixin = prof->pix;
+      if (prof->code==PROF_EXPONENTIAL)
         for (ix2=nx2; ix2--; x2+=1.0)
           {
           x1 = x10;
@@ -3423,54 +3453,11 @@ float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
               *(pixin++) = 0.0;
               continue;
               }
-          val = expf(k*expf(logf(ra)*hinvn));
-          noversamp  = (int)(val*profit->oversamp+0.1);
-          if (noversamp < 2)
+            val = ra<rs2? a0+ra*(a2+a3*sqrtf(ra)) : expf(-sqrtf(ra));
             *(pixin++) = val;
-          else
-            {
-            ostep = 1.0/noversamp;
-            dcd11 = cd11*ostep;
-            dcd21 = cd21*ostep;
-            odx = 0.5*(ostep-1.0);
-            x1t = x1+odx;
-            val = 0.0;
-            for (idx2=noversamp; idx2--; odx+=ostep)
-              {
-              x1in = cd12*(x2+odx) + cd11*x1t;
-              x2in = cd22*(x2+odx) + cd21*x1t;
-              for (idx1=noversamp; idx1--;)
-                {
-                rao = x1in*x1in+x2in*x2in;
-                val += expf(k*PROFIT_POWF(rao,hinvn));
-                x1in += dcd11;
-                x2in += dcd21;
-                }
-              }
-            *(pixin++) = val*ostep*ostep;
-            }
             }
           }
-/*---- Copy the symmetric part */
-      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
-        {
-        pixin2 = pixin - profit->modnaxisn[0] - 1;
-        if (!(profit->modnaxisn[0]&1))
-          {
-          *(pixin++) = 0.0;
-          npix2--;
-          }
-        for (i=npix2; i--;)
-          *(pixin++) = *(pixin2--);
-        }
-      prof->lostfluxfrac = 1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn));
-      threshflag = 0;
-      break;
-    case PROF_DEVAUCOULEURS:
-/*---- The consequence of sampling on flux is compensated by PSF normalisation*/
-      x10 = -x1cout - dx1;
-      x2 = -x2cout - dx2;
-      pixin = prof->pix;
+      else
         for (ix2=nx2; ix2--; x2+=1.0)
           {
           x1 = x10;
@@ -3484,32 +3471,8 @@ float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
               *(pixin++) = 0.0;
               continue;
               }
-          val = expf(-7.66924944f*PROFIT_POWF(ra,0.125));
-          noversamp  = (int)(sqrt(val)*profit->oversamp+0.1);
-          if (noversamp < 2)
+            val = ra<rs2? a0+ra*(a2+a3*sqrtf(ra)) : expf(k*expf(logf(ra)*hinvn));
             *(pixin++) = val;
-          else
-            {
-            ostep = 1.0/noversamp;
-            dcd11 = cd11*ostep;
-            dcd21 = cd21*ostep;
-            odx = 0.5*(ostep-1.0);
-            x1t = x1+odx;
-            val = 0.0;
-            for (idx2=noversamp; idx2--; odx+=ostep)
-              {
-              x1in = cd12*(x2+odx) + cd11*x1t;
-              x2in = cd22*(x2+odx) + cd21*x1t;
-              for (idx1=noversamp; idx1--;)
-                {
-                ra = x1in*x1in+x2in*x2in;
-                val += expf(-7.66924944f*PROFIT_POWF(ra,0.125));
-                x1in += dcd11;
-                x2in += dcd21;
-                }
-              }
-            *(pixin++) = val*ostep*ostep;
-            }
             }
           }
 /*---- Copy the symmetric part */
@@ -3524,68 +3487,54 @@ float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
         for (i=npix2; i--;)
           *(pixin++) = *(pixin2--);
         }
-      prof->lostfluxfrac = 1.0-prof_gammainc(8.0, 7.66924944*pow(r2max, 0.125));
-      threshflag = 0;
-      break;
-    case PROF_EXPONENTIAL:
-      x10 = -x1cout - dx1;
-      x2 = -x2cout - dx2;
+
+/*---- Compute the sharp part of the profile */
+      ix1max = profit->modnaxisn[0];
+      ix2max = profit->modnaxisn[1];
+      dx1cout = x1cout + 0.4999999;
+      dx2cout = x2cout + 0.4999999;
+      invdet = 1.0/fabsf(cd11*cd22 - cd12*cd21);
+      a11 = cd22*invdet;
+      a12 = -cd12*invdet;
+      a21 = -cd21*invdet;
+      a22 = cd11*invdet;
+      nang = 72 / 2;		/* 36 angles; only half of them are computed*/
+      angstep = PI/nang;
+      ang = 0.0;
+      for (a=0; a<nang; a++)
+        {
+        sincosf(ang, &dca, &dsa);
+        ddx1[a] = a11*dca+a12*dsa;
+        ddx2[a] = a21*dca+a22*dsa;
+        ang += angstep;
+        }
+      r = DEXPF(-4.0);
+      dr = DEXPF(0.05);
+      selem = 0.5*angstep*(dr - 1.0/dr)/(xscale*yscale);
+      krpinvn = k*DEXPF(-4.0*invn);
+      dkrpinvn = DEXPF(0.05*invn);
       pixin = prof->pix;
-      for (ix2=nx2; ix2--; x2+=1.0)
-        {
-        x1 = x10;
-        for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
-          {
-          x1in = cd12*x2 + cd11*x1;
-          x2in = cd22*x2 + cd21*x1;
-          ra = x1in*x1in+x2in*x2in;
-          if (ra>r2max)
-            {
-            *(pixin++) = 0.0;
-            continue;
-            }
-          val = expf(-sqrtf(ra));
-          noversamp  = (int)(val*sqrt(profit->oversamp)+0.1);
-          if (noversamp < 2)
-            *(pixin++) = val;
-          else
-            {
-            ostep = 1.0/noversamp;
-            dcd11 = cd11*ostep;
-            dcd21 = cd21*ostep;
-            odx = 0.5*(ostep-1.0);
-            x1t = x1+odx;
-            val = 0.0;
-            for (idx2=noversamp; idx2--; odx+=ostep)
+      for (; r<rs; r *= dr)
         {
-              x1in = cd12*(x2+odx) + cd11*x1t;
-              x2in = cd22*(x2+odx) + cd21*x1t;
-              for (idx1=noversamp; idx1--;)
+        r2 = r*r;
+        val = (expf(krpinvn) - (a0 + r2*(a2+a3*r)))*r2*selem;
+        for (a=0; a<nang; a++)
           {
-                ra = x1in*x1in+x2in*x2in;
-                val += expf(-sqrtf(ra));
-                x1in += dcd11;
-                x2in += dcd21;
-                }
-              }
-            *(pixin++) = val*ostep*ostep;
+          ix1 = (int)(dx1cout + r*ddx1[a]);
+          ix2 = (int)(dx2cout + r*ddx2[a]);
+          if (ix1>=0 && ix1<ix1max && ix2>=0 && ix2<ix2max)
+            pixin[ix2*ix1max+ix1] += val;
+          ix1 = (int)(dx1cout - r*ddx1[a]);
+          ix2 = (int)(dx2cout - r*ddx2[a]);
+          if (ix1>=0 && ix1<ix1max && ix2>=0 && ix2<ix2max)
+            pixin[ix2*ix1max+ix1] += val;
           }
+        krpinvn *= dkrpinvn;
         }
-        }
-/*---- Copy the symmetric part */
-      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
-        {
-        pixin2 = pixin - profit->modnaxisn[0] - 1;
-        if (!(profit->modnaxisn[0]&1))
-          {
-          *(pixin++) = 0.0;
-          npix2--;
-          }
-        for (i=npix2; i--;)
-          *(pixin++) = *(pixin2--);
-        }
-      rmax = sqrt(r2max);
-      prof->lostfluxfrac = (1.0 + rmax)*exp(-rmax);
+
+      prof->lostfluxfrac = prof->code==PROF_EXPONENTIAL?
+		(1.0 + rmax)*exp(-rmax)
+		:1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn));
       threshflag = 0;
       break;
     case PROF_ARMS:
@@ -3848,7 +3797,7 @@ width = 3.0;
     if (prof->lostfluxfrac < 1.0)
       flux /= (1.0 - prof->lostfluxfrac);
 
-    prof->fluxfac = fluxfac = fabs(flux)>0.0? profit->fluxfac/flux : 0.0;
+    prof->fluxfac = fluxfac = fabs(flux)>0.0? profit->fluxfac/fabs(flux) : 0.0;
     }
 
   pixin = prof->pix;

src/profit.h

@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		19/10/2010
+*	Last modified:		24/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -45,11 +45,11 @@
 #define	PARAM_ALLPARAMS	(-1)	/* All parameters */
 #define	PROFIT_MAXITER	1000	/* Max. nb of iterations in profile fitting */
 #define	PROFIT_MAXPROF	8	/* Max. nb of profile components */
-#define	PROFIT_OVERSAMP	5	/* Max. profile oversamp. factor on each axis */
 #define	PROFIT_HIDEFRES	201	/* Hi. def. model resol. (must be <MAXMODSIZE)*/
 #define	PROFIT_REFFFAC	3.0	/* Factor in r_eff for measurement radius*/
 #define	PROFIT_MAXR2MAX	1e6	/* Maximum r2_max for truncating profiles */
 #define	PROFIT_DYNPARAM	10.0	/* Dynamic compression param. in sigma units */
+#define	PROFIT_SMOOTHR	4.0	/* Profile smoothing radius (pixels) */
 #define	PROFIT_MAXMODSIZE  512	/* Maximum size allowed for the model raster */
 #define	PROFIT_MAXOBJSIZE  512	/* Maximum size allowed for the object raster */
 #define	PROFIT_BARXFADE	0.1	/* Fract. of bar length crossfaded with arms */
@@ -147,7 +147,6 @@ typedef struct
   float		pixstep;	/* Model/PSF sampling step */
   float		fluxfac;	/* Model flux scaling factor */
   float		subsamp;	/* Subsampling factor */
-  int		oversamp;	/* Oversampling factor */
   float		*psfdft;	/* Compressed Fourier Transform of the PSF */
   float		*psfpix;	/* Full res. pixmap of the PSF */
   float		*modpix;	/* Full res. pixmap of the complete model */

src/types.h

@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		19/10/2010
+*	Last modified:		12/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -215,7 +215,7 @@ typedef struct
   float		polarw;				/* WORLD "polarization" */
   float		sprob;				/* Stellarity index */
   float		fwhmw;				/* WORLD FWHM */
-  float		*assoc;				/* ASSOCiated data */
+  double	*assoc;				/* ASSOCiated data */
   int		assoc_number;			/* nb of ASSOCiated objects */
   float		*vignet;			/* Pixel data */
   float		*vigshift;			/* (Shifted) pixel data */

src/xml.c

@@ -22,7 +22,7 @@
 *	You should have received a copy of the GNU General Public License
 *	along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
 *
-*	Last modified:		11/10/2010
+*	Last modified:		24/01/2011
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
@@ -231,7 +231,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	03/08/2010
+VERSION	24/01/2011
  ***/
 int	write_xml_meta(FILE *file, char *error)
   {
@@ -754,6 +754,11 @@ int	write_xml_meta(FILE *file, char *error)
         fprintf(file, " %d", prefs.assoc_param[n]);
       fprintf(file, "\"/>\n");
       }
+    fprintf(file,
+	"   <PARAM name=\"AssocCoord_Type\" datatype=\"char\" arraysize=\"*\""
+	" ucd=\"meta.code;obs.param\" value=\"%s\"/>\n",
+    	key[findkeys("ASSOCCOORD_TYPE", keylist,
+			FIND_STRICT)].keylist[prefs.assoccoord_type]);
     fprintf(file, "   <PARAM name=\"Assoc_Radius\" datatype=\"float\""
 	" ucd=\"phys.size.radius;obs.param\" value=\"%g\" unit=\"pix\"/>\n",
     	prefs.assoc_radius);