profit.c

/*
*				profit.c
*
* Fit a range of galaxy models to an image.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	SExtractor
*
*	Copyright:		(C) 2006-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
*
*	License:		GNU General Public License
*
*	SExtractor is free software: you can redistribute it and/or modify
*	it under the terms of the GNU General Public License as published by
*	the Free Software Foundation, either version 3 of the License, or
*	(at your option) any later version.
*	SExtractor is distributed in the hope that it will be useful,
*	but WITHOUT ANY WARRANTY; without even the implied warranty of
*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*	GNU General Public License for more details.
*	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:		06/09/2011
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#ifdef HAVE_CONFIG_H
#include        "config.h"
#endif

#ifndef HAVE_MATHIMF_H
#define _GNU_SOURCE
#endif

#include	<math.h>
#include	<stdio.h>
#include	<stdlib.h>
#include	<string.h>

#include	"define.h"
#include	"globals.h"
#include	"prefs.h"
#include	"fits/fitscat.h"
#include	"levmar/levmar.h"
#include	"fft.h"
#include	"fitswcs.h"
#include	"check.h"
#include	"image.h"
#include	"pattern.h"
#include	"psf.h"
#include	"profit.h"

static double	prof_gammainc(double x, double a),
		prof_gamma(double x);
static float	prof_interpolate(profstruct *prof, float *posin);
static float	interpolate_pix(float *posin, float *pix, int *naxisn,
		interpenum interptype);

static void	make_kernel(float pos, float *kernel, interpenum interptype);

/*------------------------------- variables ---------------------------------*/

const int	interp_kernwidth[5]={1,2,4,6,8};

const int	flux_flag[PARAM_NPARAM] = {0,
					1,0,0,
					1,0,0,0,0,
					1,0,0,0,
					1,0,0,0,0,0,0,0,
					1,0,0,
					1,0,0,
					1,0,0
					};

/* "Local" global variables for debugging purposes */
int theniter, the_gal;
static picstruct	*the_field, *the_wfield;
profitstruct		*theprofit, *thepprofit, *theqprofit;

/****** profit_init ***********************************************************
PROTO	profitstruct profit_init(psfstruct *psf, unsigned int modeltype)
PURPOSE	Allocate and initialize a new profile-fitting structure.
INPUT	Pointer to PSF structure,
	Model type.
OUTPUT	A pointer to an allocated profit structure.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	22/04/2011
 ***/
profitstruct	*profit_init(psfstruct *psf, unsigned int modeltype)
  {
   profitstruct		*profit;
   int			t, nmodels;

  QCALLOC(profit, profitstruct, 1);
  profit->psf = psf;
  QMALLOC(profit->prof, profstruct *, MODEL_NMAX);
  nmodels = 0;
  for (t=1; t<(1<<MODEL_NMAX); t<<=1)
    if (modeltype&t)
      profit->prof[nmodels++] = prof_init(profit, t);
/* Allocate memory for the complete model */
  QMALLOC16(profit->modpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->cmodpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->psfpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);
  QMALLOC16(profit->objpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->objweight, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->lmodpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);
  QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
  profit->nprof = nmodels;
  profit->fluxfac = 1.0;	/* Default */

  return profit;
  }  


/****** profit_end ************************************************************
PROTO	void prof_end(profstruct *prof)
PURPOSE	End (deallocate) a profile-fitting structure.
INPUT	Prof structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	06/04/2010
 ***/
void	profit_end(profitstruct *profit)
  {
   int	p;

  for (p=0; p<profit->nprof; p++)
    prof_end(profit->prof[p]);
  free(profit->modpix);
  free(profit->modpix2);
  free(profit->cmodpix);
  free(profit->psfpix);
  free(profit->lmodpix);
  free(profit->lmodpix2);
  free(profit->objpix);
  free(profit->objweight);
  free(profit->resi);
  free(profit->prof);
  free(profit->covar);
  free(profit->psfdft);
  free(profit);

  return;
  }


/****** profit_fit ************************************************************
PROTO	void profit_fit(profitstruct *profit, picstruct *field,
		picstruct *wfield, objstruct *obj, obj2struct *obj2)
PURPOSE	Fit profile(s) convolved with the PSF to a detected object.
INPUT	Array of profile structures,
	Number of profiles,
	Pointer to the profile-fitting structure,
	Pointer to the field,
	Pointer to the field weight,
	Pointer to the obj.
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	06/09/2011
 ***/
void	profit_fit(profitstruct *profit,
		picstruct *field, picstruct *wfield,
		objstruct *obj, obj2struct *obj2)
  {
    profitstruct	*pprofit, *qprofit;
    patternstruct	*pattern;
    psfstruct		*psf;
    checkstruct		*check;
    double		emx2,emy2,emxy, a , cp,sp, cn, bn, n,
			sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0;
    PIXTYPE		valp,valq,sig2;
    float		param0[PARAM_NPARAM], param1[PARAM_NPARAM],
			param[PARAM_NPARAM],
			**list,
			*cov,
			psf_fwhm, dchi2, err, aspect, chi2;
    int			*index,
			i,j,p, nparam, nparam2, ncomp, nprof;

  nparam = profit->nparam;
  nparam2 = nparam*nparam;
  nprof = profit->nprof;
  if (profit->psfdft)
    {
    QFREE(profit->psfdft);
    }

  psf = profit->psf;
  profit->pixstep = psf->pixstep;
  obj2->prof_flag = 0;

/* Create pixmaps at image resolution */
  profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
  profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
  psf_fwhm = psf->masksize[0]*psf->pixstep;
  profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
		*1.2)/2)*2 + 1;
  if (profit->objnaxisn[1]<profit->objnaxisn[0])
    profit->objnaxisn[1] = profit->objnaxisn[0];
  else
    profit->objnaxisn[0] = profit->objnaxisn[1];
  if (profit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
    {
    profit->subsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
    profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->subsamp);
    obj2->prof_flag |= PROFLAG_OBJSUB;
    }
  else
    profit->subsamp = 1.0;
  profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1];

/* Use (dirty) global variables to interface with lmfit */
  the_field = field;
  the_wfield = wfield;
  theprofit = profit;
  profit->obj = obj;
  profit->obj2 = obj2;

  profit->nresi = profit_copyobjpix(profit, field, wfield);
/* Check if the number of constraints exceeds the number of free parameters */
  if (profit->nresi < nparam)
    {
    if (FLAG(obj2.prof_vector))
      for (p=0; p<nparam; p++)
        obj2->prof_vector[p] = 0.0;
    if (FLAG(obj2.prof_errvector))
      for (p=0; p<nparam; p++)
        obj2->prof_errvector[p] = 0.0;
    if (FLAG(obj2.prof_errmatrix))
      for (p=0; p<nparam2; p++)
        obj2->prof_errmatrix[p] = 0.0;
    obj2->prof_niter = 0;
    obj2->prof_flag |= PROFLAG_NOTCONST;
    return;
    }

/* Create pixmap at PSF resolution */
  profit->modnaxisn[0] =
	((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep
		+0.4999)/2+1)*2; 
  profit->modnaxisn[1] =
	((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep
		+0.4999)/2+1)*2; 
  if (profit->modnaxisn[1] < profit->modnaxisn[0])
    profit->modnaxisn[1] = profit->modnaxisn[0];
  else
    profit->modnaxisn[0] = profit->modnaxisn[1];
  if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)
    {
    profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;
    profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;
    obj2->prof_flag |= PROFLAG_MODSUB;
    }
  profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];

/* Compute the local PSF */
  profit_psf(profit);

/* Set initial guesses and boundaries */
  profit->sigma = obj->sigbkg;
  if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)
    profit->guessradius = obj2->hl_radius;
  if ((profit->guessflux = obj2->flux_auto) <= 0.0)
    profit->guessflux = 0.0;
  if ((profit->guessfluxmax = 10.0*obj2->fluxerr_auto) <= profit->guessflux)
    profit->guessfluxmax = profit->guessflux;
  if (profit->guessfluxmax <= 0.0)
    profit->guessfluxmax = 1.0;

  profit_resetparams(profit);

/* Actual minimisation */
  fft_reset();
the_gal++;

/*
char str[1024];
sprintf(str, "obj_%04d.fits", the_gal);
catstruct *bcat;
float *bpix, *opix,*lmodpix,*objpix;
bcat=read_cat("base.fits");
QMALLOC(bpix, float, field->npix);
QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);
read_body(bcat->tab, bpix, field->npix); 
free_cat(&bcat,1);
bcat=read_cat(str);
QMALLOC(opix, float, profit->nobjpix);
QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);
read_body(bcat->tab, opix, profit->nobjpix); 
free_cat(&bcat,1);
addfrombig(bpix, field->width, field->height,
		profit->objpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
objpix = profit->objpix;
lmodpix = opix;
for (i=profit->nobjpix; i--;)
*(objpix++) += *(lmodpix++);
free(bpix);
free(opix);
*/
  profit->niter = profit_minimize(profit, PROFIT_MAXITER);
/*
profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
check=initcheck(str, CHECK_OTHER,1);
check->width = profit->objnaxisn[0];
check->height = profit->objnaxisn[1];
reinitcheck(field,check);
memcpy(check->pix, profit->lmodpix, profit->nobjpix*sizeof(float));
reendcheck(field,check);
endcheck(check);

  chi2 = profit->chi2;
  for (p=0; p<nparam; p++)
    param1[p] = profit->paraminit[p];
  profit_resetparams(profit);
  for (p=0; p<nparam; p++)
    profit->paraminit[p] = param1[p] + (profit->paraminit[p]<param1[p]?1.0:-1.0)
			* sqrt(profit->covar[p*(nparam+1)]);

  profit->niter = profit_minimize(profit, PROFIT_MAXITER);
  if (chi2<profit->chi2)
    for (p=0; p<nparam; p++)
      profit->paraminit[p] = param1[p];

list = profit->paramlist;
index = profit->paramindex;
for (i=0; i<PARAM_NPARAM; i++)
if (list[i] && i!= PARAM_SPHEROID_ASPECT && i!=PARAM_SPHEROID_POSANG)
profit->freeparam_flag[index[i]] = 0;
profit->niter = profit_minimize(profit, PROFIT_MAXITER);
*/
  if (profit->nlimmin)
    obj2->prof_flag |= PROFLAG_MINLIM;
  if (profit->nlimmax)
    obj2->prof_flag |= PROFLAG_MAXLIM;

  for (p=0; p<nparam; p++)
    profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);

/* CHECK-Images */
  if ((check = prefs.check[CHECK_PROFILES]))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }

  if ((check = prefs.check[CHECK_SUBPROFILES]))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
  if ((check = prefs.check[CHECK_SPHEROIDS]))
    {
/*-- Set to 0 flux components that do not belong to spheroids */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }
  if ((check = prefs.check[CHECK_SUBSPHEROIDS]))
    {
/*-- Set to 0 flux components that do not belong to spheroids */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
  if ((check = prefs.check[CHECK_DISKS]))
    {
/*-- Set to 0 flux components that do not belong to disks */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);
    }
  if ((check = prefs.check[CHECK_SUBDISKS]))
    {
/*-- Set to 0 flux components that do not belong to disks */
    for (p=0; p<profit->nparam; p++)
      param[p] = profit->paraminit[p];
    list = profit->paramlist;
    index = profit->paramindex;
    for (i=0; i<PARAM_NPARAM; i++)
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
        param[index[i]] = 0.0;
    profit_residuals(profit,field,wfield, 0.0, param, NULL);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, -1.0);
    }
 
/* Compute compressed residuals */
  profit_residuals(profit,field,wfield, 10.0, profit->paraminit,profit->resi);

/* Fill measurement parameters */
  if (FLAG(obj2.prof_vector))
    {
    for (p=0; p<nparam; p++)
      obj2->prof_vector[p]= profit->paraminit[p];
    }
  if (FLAG(obj2.prof_errvector))
    {
    for (p=0; p<nparam; p++)
      obj2->prof_errvector[p]= profit->paramerr[p];
    }
  if (FLAG(obj2.prof_errmatrix))
    {
    for (p=0; p<nparam2; p++)
      obj2->prof_errmatrix[p]= profit->covar[p];
    }

  obj2->prof_niter = profit->niter;
  obj2->flux_prof = profit->flux;
  if (FLAG(obj2.fluxerr_prof))
    {
    err = 0.0;
    cov = profit->covar;
    index = profit->paramindex;
    list = profit->paramlist;
    for (i=0; i<PARAM_NPARAM; i++)
      if (flux_flag[i] && list[i])
        {
        cov = profit->covar + nparam*index[i];
        for (j=0; j<PARAM_NPARAM; j++)
          if (flux_flag[j] && list[j])
            err += cov[index[j]];
        }
    obj2->fluxerr_prof = err>0.0? sqrt(err): 0.0;
    }

  obj2->prof_chi2 = (profit->nresi > profit->nparam)?
		profit->chi2 / (profit->nresi - profit->nparam) : 0.0;

/* Position */
  if (FLAG(obj2.x_prof))
    {
    i = profit->paramindex[PARAM_X];
    j = profit->paramindex[PARAM_Y];
/*-- Model coordinates follow the FITS convention (first pixel at 1,1) */
    if (profit->paramlist[PARAM_X])
      {
      obj2->x_prof = (double)profit->ix + *profit->paramlist[PARAM_X] + 1.0;
      obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];
      }
    else
      emx2 = 0.0;
    if (profit->paramlist[PARAM_Y])
      {
      obj2->y_prof = (double)profit->iy + *profit->paramlist[PARAM_Y] + 1.0;
      obj2->poserrmy2_prof = emy2 = profit->covar[j*(nparam+1)];
      }
    else
      emy2 = 0.0;
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
      obj2->poserrmxy_prof = emxy = profit->covar[i+j*nparam];
    else
      emxy = 0.0;

/*-- Error ellipse parameters */
    if (FLAG(obj2.poserra_prof))
      {
       double	pmx2,pmy2,temp,theta;

      if (fabs(temp=emx2-emy2) > 0.0)
        theta = atan2(2.0 * emxy,temp) / 2.0;
      else
        theta = PI/4.0;

      temp = sqrt(0.25*temp*temp+ emxy*emxy);
      pmy2 = pmx2 = 0.5*(emx2+emy2);
      pmx2+=temp;
      pmy2-=temp;

      obj2->poserra_prof = (float)sqrt(pmx2);
      obj2->poserrb_prof = (float)sqrt(pmy2);
      obj2->poserrtheta_prof = (float)(theta/DEG);
      }

    if (FLAG(obj2.poserrcxx_prof))
      {
       double	temp;

      obj2->poserrcxx_prof = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
      obj2->poserrcyy_prof = (float)(emx2/temp);
      obj2->poserrcxy_prof = (float)(-2*emxy/temp);
      }
    }

/* Equivalent noise area */
  if (FLAG(obj2.prof_noisearea))
    obj2->prof_noisearea = profit_noisearea(profit);

/* Second order moments and ellipticities */
  if (FLAG(obj2.prof_mx2))
    profit_moments(profit, obj2);

/* Second order moments of the convolved model (used by other parameters) */
  if (FLAG(obj2.prof_convmx2))
    profit_convmoments(profit, obj2);

/* "Hybrid" magnitudes */
  if (FLAG(obj2.fluxcor_prof))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    profit_fluxcor(profit, obj, obj2);
    }

/* Do measurements on the rasterised model (surface brightnesses) */
  if (FLAG(obj2.fluxeff_prof))
    profit_surface(profit, obj2); 

/* Background offset */
  if (FLAG(obj2.prof_offset_flux))
    {
    obj2->prof_offset_flux = *profit->paramlist[PARAM_BACK];
    obj2->prof_offset_fluxerr=profit->paramerr[profit->paramindex[PARAM_BACK]];
    }

/* Point source */
  if (FLAG(obj2.prof_dirac_flux))
    {
    obj2->prof_dirac_flux = *profit->paramlist[PARAM_DIRAC_FLUX];
    obj2->prof_dirac_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX]];
    }

/* 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]];
    obj2->prof_spheroid_reff = *profit->paramlist[PARAM_SPHEROID_REFF];
    obj2->prof_spheroid_refferr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]];
    obj2->prof_spheroid_aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT];
    obj2->prof_spheroid_aspecterr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]];
    obj2->prof_spheroid_theta =
			fmod_m90_p90(*profit->paramlist[PARAM_SPHEROID_POSANG]);
    obj2->prof_spheroid_thetaerr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_POSANG]];
    if (FLAG(obj2.prof_spheroid_sersicn))
      {
      obj2->prof_spheroid_sersicn = *profit->paramlist[PARAM_SPHEROID_SERSICN];
      obj2->prof_spheroid_sersicnerr = 
		profit->paramerr[profit->paramindex[PARAM_SPHEROID_SERSICN]];
      }
    else
      obj2->prof_spheroid_sersicn = 4.0;
    if (FLAG(obj2.prof_spheroid_peak))
      {
      n = obj2->prof_spheroid_sersicn;
      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 */
      cn = n * prof_gamma(2.0*n) * pow(bn, -2.0*n);
      obj2->prof_spheroid_peak = obj2->prof_spheroid_reff>0.0?
	obj2->prof_spheroid_flux
		/ (2.0 * PI * cn
		* obj2->prof_spheroid_reff*obj2->prof_spheroid_reff
		* obj2->prof_spheroid_aspect)
	: 0.0;
      if (FLAG(obj2.prof_spheroid_fluxeff))
        obj2->prof_spheroid_fluxeff = obj2->prof_spheroid_peak * exp(-bn);
      if (FLAG(obj2.prof_spheroid_fluxmean))
        obj2->prof_spheroid_fluxmean = obj2->prof_spheroid_peak * cn;
      }
    }

/* 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]];
    obj2->prof_disk_scale = *profit->paramlist[PARAM_DISK_SCALE];
    obj2->prof_disk_scaleerr =
		profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]];
    obj2->prof_disk_aspect = *profit->paramlist[PARAM_DISK_ASPECT];
    obj2->prof_disk_aspecterr =
		profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]];
    obj2->prof_disk_theta = fmod_m90_p90(*profit->paramlist[PARAM_DISK_POSANG]);
    obj2->prof_disk_thetaerr =
		profit->paramerr[profit->paramindex[PARAM_DISK_POSANG]];
    if (FLAG(obj2.prof_disk_inclination))
      {
      obj2->prof_disk_inclination = acos(obj2->prof_disk_aspect) / DEG;
      if (FLAG(obj2.prof_disk_inclinationerr))
        {
        a = sqrt(1.0-obj2->prof_disk_aspect*obj2->prof_disk_aspect);
        obj2->prof_disk_inclinationerr = obj2->prof_disk_aspecterr
					/(a>0.1? a : 0.1)/DEG;
        }
      }

    if (FLAG(obj2.prof_disk_peak))
      {
      obj2->prof_disk_peak = obj2->prof_disk_scale>0.0?
	obj2->prof_disk_flux
	/ (2.0 * PI * obj2->prof_disk_scale*obj2->prof_disk_scale
		* obj2->prof_disk_aspect)
	: 0.0;
      if (FLAG(obj2.prof_disk_fluxeff))
        obj2->prof_disk_fluxeff = obj2->prof_disk_peak * 0.186682; /* e^-(b_n)*/
      if (FLAG(obj2.prof_disk_fluxmean))
        obj2->prof_disk_fluxmean = obj2->prof_disk_peak * 0.355007;/* b_n^(-2)*/
      }

/* Disk pattern */
    if (prefs.pattern_flag)
      {
      profit_residuals(profit,field,wfield, PROFIT_DYNPARAM,
			profit->paraminit,profit->resi);
      pattern = pattern_init(profit, prefs.pattern_type,
		prefs.prof_disk_patternncomp);
      pattern_fit(pattern, profit);
      if (FLAG(obj2.prof_disk_patternspiral))
        obj2->prof_disk_patternspiral = pattern_spiral(pattern);
      if (FLAG(obj2.prof_disk_patternvector))
        {
        ncomp = pattern->size[2];
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternvector[p] = (float)pattern->coeff[p];
        }
      if (FLAG(obj2.prof_disk_patternmodvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternmodvector[p] = (float)pattern->mcoeff[p];
        }
      if (FLAG(obj2.prof_disk_patternargvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternargvector[p] = (float)pattern->acoeff[p];
        }
      pattern_end(pattern);
      }

/* Bar */
    if (FLAG(obj2.prof_bar_flux))
      {
      obj2->prof_bar_flux = *profit->paramlist[PARAM_BAR_FLUX];
      obj2->prof_bar_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_BAR_FLUX]];
      obj2->prof_bar_length = *profit->paramlist[PARAM_ARMS_START]
				**profit->paramlist[PARAM_DISK_SCALE];
      obj2->prof_bar_lengtherr = *profit->paramlist[PARAM_ARMS_START]
		  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
		+ *profit->paramlist[PARAM_DISK_SCALE]
		  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
      obj2->prof_bar_aspect = *profit->paramlist[PARAM_BAR_ASPECT];
      obj2->prof_bar_aspecterr =
		profit->paramerr[profit->paramindex[PARAM_BAR_ASPECT]];
      obj2->prof_bar_posang = 
			fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
      obj2->prof_bar_posangerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
      if (FLAG(obj2.prof_bar_theta))
        {
        cp = cos(obj2->prof_bar_posang*DEG);
        sp = sin(obj2->prof_bar_posang*DEG);
        a = obj2->prof_disk_aspect;
        obj2->prof_bar_theta = fmod_m90_p90(atan2(a*sp,cp)/DEG
				+ obj2->prof_disk_theta);
        obj2->prof_bar_thetaerr = obj2->prof_bar_posangerr*a/(cp*cp+a*a*sp*sp);
        }

/* Arms */
      if (FLAG(obj2.prof_arms_flux))
        {
        obj2->prof_arms_flux = *profit->paramlist[PARAM_ARMS_FLUX];
        obj2->prof_arms_fluxerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_FLUX]];
        obj2->prof_arms_pitch =
		fmod_m90_p90(*profit->paramlist[PARAM_ARMS_PITCH]);
        obj2->prof_arms_pitcherr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_PITCH]];
        obj2->prof_arms_start = *profit->paramlist[PARAM_ARMS_START]
				**profit->paramlist[PARAM_DISK_SCALE];
        obj2->prof_arms_starterr = *profit->paramlist[PARAM_ARMS_START]
		  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
		+ *profit->paramlist[PARAM_DISK_SCALE]
		  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
        obj2->prof_arms_quadfrac = *profit->paramlist[PARAM_ARMS_QUADFRAC];
        obj2->prof_arms_quadfracerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_QUADFRAC]];
        obj2->prof_arms_posang =
			fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
        obj2->prof_arms_posangerr =
		profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
        }
      }
    }

/* Star/galaxy classification */
  if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration))
    {
profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, profit->paraminit,
profit->resi);
    pprofit = thepprofit;
    nparam = pprofit->nparam;
    if (pprofit->psfdft)
      {
      QFREE(pprofit->psfdft);
      }
    psf = pprofit->psf;
    pprofit->pixstep = profit->pixstep;
    pprofit->guessradius = profit->guessradius;
    pprofit->guessflux = profit->guessflux;
    pprofit->guessfluxmax = profit->guessfluxmax;
    pprofit->ix = profit->ix;
    pprofit->iy = profit->iy;
    pprofit->objnaxisn[0] = profit->objnaxisn[0];
    pprofit->objnaxisn[1] = profit->objnaxisn[1];
    pprofit->subsamp = profit->subsamp;
    pprofit->nobjpix = profit->nobjpix;
    pprofit->obj = obj;
    pprofit->obj2 = obj2;
    pprofit->nresi = profit_copyobjpix(pprofit, field, wfield);
    pprofit->modnaxisn[0] = profit->modnaxisn[0];
    pprofit->modnaxisn[1] = profit->modnaxisn[1];
    pprofit->nmodpix = profit->nmodpix;
    profit_psf(pprofit);
    pprofit->sigma = obj->sigbkg;
    profit_resetparams(pprofit);
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
      {
      pprofit->paraminit[pprofit->paramindex[PARAM_X]] = *profit->paramlist[PARAM_X];
      pprofit->paraminit[pprofit->paramindex[PARAM_Y]] = *profit->paramlist[PARAM_Y];
      }
    fft_reset();
    pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]] = profit->flux;
    pprofit->niter = profit_minimize(pprofit, PROFIT_MAXITER);
    profit_residuals(pprofit,field,wfield, PROFIT_DYNPARAM, pprofit->paraminit,
			pprofit->resi);
    qprofit = theqprofit;
    nparam = qprofit->nparam;
    if (qprofit->psfdft)
      {
      QFREE(qprofit->psfdft);
      }
    qprofit->pixstep = profit->pixstep;
    qprofit->guessradius = profit->guessradius;
    qprofit->guessflux = profit->guessflux;
    qprofit->guessfluxmax = profit->guessfluxmax;
    qprofit->ix = profit->ix;
    qprofit->iy = profit->iy;
    qprofit->objnaxisn[0] = profit->objnaxisn[0];
    qprofit->objnaxisn[1] = profit->objnaxisn[1];
    qprofit->subsamp = profit->subsamp;
    qprofit->nobjpix = profit->nobjpix;
    qprofit->obj = obj;
    qprofit->obj2 = obj2;
    qprofit->nresi = profit_copyobjpix(qprofit, field, wfield);
    qprofit->modnaxisn[0] = profit->modnaxisn[0];
    qprofit->modnaxisn[1] = profit->modnaxisn[1];
    qprofit->nmodpix = profit->nmodpix;
    profit_psf(qprofit);
    qprofit->sigma = obj->sigbkg;
    profit_resetparams(qprofit);
    fft_reset();
    qprofit->paraminit[qprofit->paramindex[PARAM_X]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_X]];
    qprofit->paraminit[qprofit->paramindex[PARAM_Y]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_Y]];
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_FLUX]]
		= pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX]];
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_SCALE]] = psf->fwhm/16.0;
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_ASPECT]] = 1.0;
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_POSANG]] = 0.0;
    profit_residuals(qprofit,field,wfield, PROFIT_DYNPARAM, qprofit->paraminit,
			qprofit->resi);
    sump = sumq = sumpw2 = sumqw2 = sumpqw = sump0 = sumq0 = 0.0;
    for (p=0; p<pprofit->nobjpix; p++)
      if (pprofit->objweight[p]>0 && pprofit->objpix[p]>-BIG)
        {
        valp = pprofit->lmodpix[p];
        sump += (double)(valp*pprofit->objpix[p]);
	valq = qprofit->lmodpix[p];
        sumq += (double)(valq*pprofit->objpix[p]);
	sump0 += (double)(valp*valp);
	sumq0 += (double)(valp*valq);
        sig2 = 1.0f/(pprofit->objweight[p]*pprofit->objweight[p]);
        sumpw2 += valp*valp*sig2;
        sumqw2 += valq*valq*sig2;
        sumpqw += valp*valq*sig2;
        }

    if (FLAG(obj2.prof_class_star))
      {
      dchi2 = 0.5*(pprofit->chi2 - profit->chi2);
      obj2->prof_class_star = dchi2 < 50.0?
	(dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
      }
    if (FLAG(obj2.prof_concentration))
      {
      obj2->prof_concentration = sump>0.0? (sumq/sump - sumq0/sump0) : 1.0;
      if (FLAG(obj2.prof_concentrationerr))
        obj2->prof_concentrationerr = sump>0.0?
		sqrt(sumqw2*sump*sump+sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)
			/ (sump*sump) : 0.0;
      }
    }

/* clean up. */
  fft_reset();

  return;
  }


/****** profit_noisearea ******************************************************
PROTO	float profit_noisearea(profitstruct *profit)
PURPOSE	Return the equivalent noise area (see King 1983) of a model.
INPUT	Profile-fitting structure,
OUTPUT	Equivalent noise area, in pixels.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	19/10/2010
 ***/
float	profit_noisearea(profitstruct *profit)
  {
   double	dval, flux,flux2;
   PIXTYPE	*pix;
   int		p;

  flux = flux2 = 0.0;
  pix = profit->lmodpix;
  for (p=profit->nobjpix; p--;)
    {
    dval = (double)*(pix++);
    flux += dval;
    flux2 += dval*dval;
    }

  return (float)(flux2>0.0? flux*flux / flux2 : 0.0);
  }


/****** profit_fluxcor ******************************************************
PROTO	void profit_fluxcor(profitstruct *profit, objstruct *obj,
			obj2struct *obj2)
PURPOSE	Integrate the flux within an ellipse and complete it with the wings of
		the fitted model.
INPUT		Profile-fitting structure,
		pointer to the obj structure,
		pointer to the obj2 structure.
OUTPUT	Model-corrected flux.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	12/04/2011
 ***/
void	profit_fluxcor(profitstruct *profit, objstruct *obj, obj2struct *obj2)
  {
    checkstruct		*check;
    double		mx,my, dx,dy, cx2,cy2,cxy, klim,klim2, tvobj,sigtvobj,
			tvm,tvmin,tvmout, r1,v1;
    PIXTYPE		*objpix,*objpixt,*objweight,*objweightt, *lmodpix,
			pix, weight,var;
    int			x,y, x2,y2, pos, w,h, area, corrflag;


  corrflag = (prefs.mask_type==MASK_CORRECT);
  w = profit->objnaxisn[0];
  h = profit->objnaxisn[1];
  mx = (float)(w/2);
  my = (float)(h/2);
  if (FLAG(obj2.x_prof))
    {
    if (profit->paramlist[PARAM_X])
      mx += *profit->paramlist[PARAM_X];
    if (profit->paramlist[PARAM_Y])
      my += *profit->paramlist[PARAM_Y];
    }

  if (obj2->kronfactor>0.0)
    {
    cx2 = obj->cxx;
    cy2 = obj->cyy;
    cxy = obj->cxy;
    klim2 = 0.64*obj2->kronfactor*obj2->kronfactor;
    }
  else
/*-- ...if not, use the circular aperture provided by the user */
    {
    cx2 = cy2 = 1.0;
    cxy = 0.0;
    klim2 = (prefs.autoaper[1]/2.0)*(prefs.autoaper[1]/2.0);
    }
/*
  cx2 = obj2->prof_convcxx;
  cy2 = obj2->prof_convcyy;
  cxy = obj2->prof_convcxy;

  lmodpix = profit->lmodpix;
  r1 = v1 = 0.0;
  for (y=0; y<h; y++)
    {
    dy = y - my;
    for (x=0; x<w; x++)
      {
      dx = x - mx;
      pix = *(lmodpix++);
      r1 += sqrt(cx2*dx*dx + cy2*dy*dy + cxy*dx*dy)*pix;
      v1 += pix;
      }
    }

  klim = r1/v1*2.0;
  klim2 = klim*klim;

if ((check = prefs.check[CHECK_APERTURES]))
sexellips(check->pix, check->width, check->height,
obj2->x_prof-1.0, obj2->y_prof-1.0, klim*obj2->prof_conva,klim*obj2->prof_convb,
obj2->prof_convtheta, check->overlay, 0);
*/

  area = 0;
  tvmin = tvmout = tvobj = sigtvobj = 0.0;
  lmodpix = profit->lmodpix;
  objpixt = objpix = profit->objpix;
  objweightt = objweight = profit->objweight;
  for (y=0; y<h; y++)
    {
    for (x=0; x<w; x++, objpixt++,objweightt++)
      {
      dx = x - mx;
      dy = y - my;
      if ((cx2*dx*dx + cy2*dy*dy + cxy*dx*dy) <= klim2)
        {
        area++;
/*------ Here begin tests for pixel and/or weight overflows. Things are a */
/*------ bit intricated to have it running as fast as possible in the most */
/*------ common cases */
        if ((weight=*objweightt)<=0.0)
          {
          if (corrflag
		&& (x2=(int)(2*mx+0.49999-x))>=0 && x2<w
		&& (y2=(int)(2*my+0.49999-y))>=0 && y2<h
		&& (weight=objweight[pos = y2*w + x2])>0.0)
            {
            pix = objpix[pos];
            var = 1.0/(weight*weight);
            }
          else
            pix = var = 0.0;
          }
        else
          {
          pix = *objpixt;