profit.c

 /*
 				profit.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	Part of:	SExtractor
*
*	Authors:	E.BERTIN (IAP)
*
*	Contents:	Fit an arbitrary profile combination to a detection.
*
*	Last modify:	06/10/2009
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

#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/lm.h"
#include	"fft.h"
#include	"fitswcs.h"
#include	"check.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 char	profname[][32]={"background offset", "Sersic spheroid",
		"De Vaucouleurs spheroid", "exponential disk", "spiral arms",
		"bar", "inner ring", "outer ring", "tabulated model",
		""};

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

const int	flux_flag[PARAM_NPARAM] = {0,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
					};

int theniter, the_gal;
/* "Local" global variables; it seems dirty but it simplifies a lot */
/* interfacing to the LM routines */
static picstruct	*the_field, *the_wfield;
profitstruct		*theprofit;

/****** profit_init ***********************************************************
PROTO	profitstruct profit_init(psfstruct *psf)
PURPOSE	Allocate and initialize a new profile-fitting structure.
INPUT	Pointer to PSF structure.
OUTPUT	A pointer to an allocated profit structure.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	07/09/2009
 ***/
profitstruct	*profit_init(psfstruct *psf)
  {
   profitstruct		*profit;
   int			p, nprof,
			backflag, spheroidflag, diskflag, barflag, armsflag;

  QCALLOC(profit, profitstruct, 1);
  profit->psf = psf;
  profit->psfdft = NULL;

  profit->nparam = 0;
  QMALLOC(profit->prof, profstruct *, PROF_NPROF);
  backflag = spheroidflag = diskflag = barflag = armsflag = 0;
  nprof = 0;
  for (p=0; p<PROF_NPROF; p++)
    if (!backflag && FLAG(obj2.prof_offset_flux))
      {
      profit->prof[p] = prof_init(profit, PROF_BACK);
      backflag = 1;
      nprof++;
      }
    else if (!spheroidflag && FLAG(obj2.prof_spheroid_flux))
      {
      profit->prof[p] = prof_init(profit,
	FLAG(obj2.prof_spheroid_sersicn)? PROF_SERSIC : PROF_DEVAUCOULEURS);
      spheroidflag = 1;
      nprof++;
      }
    else if (!diskflag && FLAG(obj2.prof_disk_flux))
      {
      profit->prof[p] = prof_init(profit, PROF_EXPONENTIAL);
      diskflag = 1;
      nprof++;
      }
    else if (diskflag && !barflag && FLAG(obj2.prof_bar_flux))
      {
      profit->prof[p] = prof_init(profit, PROF_BAR);
      barflag = 1;
      nprof++;
      }
    else if (barflag && !armsflag && FLAG(obj2.prof_arms_flux))
      {
      profit->prof[p] = prof_init(profit, PROF_ARMS);
      armsflag = 1;
      nprof++;
      }

  QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
  profit->nprof = nprof;
  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	26/04/2008
 ***/
void	profit_end(profitstruct *profit)
  {
   int	p;

  for (p=0; p<profit->nprof; p++)
    prof_end(profit->prof[p]);
  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/10/2009
 ***/
void	profit_fit(profitstruct *profit,
		picstruct *field, picstruct *wfield,
		objstruct *obj, obj2struct *obj2)
  {
    profitstruct	pprofit;
    profitstruct	hdprofit;
    patternstruct *pattern;
    psfstruct		*psf;
    checkstruct		*check;
    double		emx2,emy2,emxy, a , cp,sp, cn, bn, n;
    float		**list,
			*cov,
			psf_fwhm, dchi2, err;
    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->nsubsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
    profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->nsubsamp);
    profit->pixstep *= profit->nsubsamp;
    obj2->prof_flag |= PROFLAG_OBJSUB;
    }
  else
    profit->nsubsamp = 1.0;

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

  QMALLOC16(profit->objpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]);
  QMALLOC16(profit->objweight, PIXTYPE,profit->objnaxisn[0]*profit->objnaxisn[1]);
  QMALLOC16(profit->lmodpix, PIXTYPE, profit->objnaxisn[0]*profit->objnaxisn[1]);
  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;
    }

  QMALLOC16(profit->resi, float, profit->nresi);

/* Create pixmap at PSF resolution */
  profit->modnaxisn[0] =
	((int)(profit->objnaxisn[0]/profit->pixstep +0.4999)/2+1)*2; 
  profit->modnaxisn[1] =
	((int)(profit->objnaxisn[1]/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;
    }

/* Allocate memory for the complete model */
  QMALLOC16(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
  memset(profit->modpix, 0, profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float));
  QMALLOC16(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
/* Allocate memory for the partial model */
  QMALLOC16(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
/* Compute the local PSF */
  profit_psf(profit);

/* Set initial guesses and boundaries */
  profit->sigma = obj->sigbkg;

  profit_resetparams(profit);

//the_gal++;

/* Actual minimisation */
  fft_reset();
the_gal++;
  profit->niter = profit_minimize(profit, PROFIT_MAXITER);
  profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi);

/* Convert covariance matrix to bound space */
  profit_covarunboundtobound(profit);
  for (p=0; p<nparam; p++)
    profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);

/* Equate param and paraminit vectors to avoid confusion later on */
  for (p=0; p<profit->nparam; p++)
    profit->param[p] = profit->paraminit[p];

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

/* Fill measurement parameters */
  if (FLAG(obj2.prof_vector))
    {
    for (p=0; p<nparam; p++)
      obj2->prof_vector[p]= profit->param[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;

  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 = 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 = 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);
      }
    }

/* Do measurements on the rasterised model (shear and surface brightnesses) */
  if (FLAG(obj2.prof_mx2) || FLAG(obj2.peak_prof))
    {
     float	scalefac, imsizefac, flux, lost, sum, lostfluxfrac;

/*-- Allocate "high-definition" rasters only to make measurements */
    hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES;
/*-- Find best image size factor from fitting results */
    imsizefac = 2.0*profit_minradius(profit, PROFIT_REFFFAC)/profit->pixstep
	/ (float)profit->modnaxisn[0];
    if (imsizefac<0.01)
      imsizefac = 0.01;
    else if (imsizefac>100.0)
      imsizefac = 100.0;
    scalefac = (float)hdprofit.modnaxisn[0] / (float)profit->modnaxisn[0]
		/ imsizefac;
    hdprofit.pixstep = profit->pixstep / scalefac;
    hdprofit.fluxfac = scalefac*scalefac;
    QCALLOC(hdprofit.modpix, float,
	hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1]*sizeof(float));
    QCALLOC(hdprofit.pmodpix, float,
	hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1]*sizeof(float));

    lost = sum = 0.0;
    for (p=0; p<profit->nprof; p++)
      {
      sum += (flux = prof_add(profit->prof[p], &hdprofit));
      lost += flux*profit->prof[p]->lostfluxfrac;
      }
    lostfluxfrac = sum > 0.0? lost / sum : 0.0;

    if (FLAG(obj2.prof_mx2))
      profit_moments(&hdprofit, obj2);

    if (FLAG(obj2.peak_prof))
      profit_surface(&hdprofit, obj2, lostfluxfrac);

/*-- Free rasters */
    free(hdprofit.modpix);
    free(hdprofit.pmodpix);
    }

/* Spheroid */
  if (FLAG(obj2.prof_spheroid_flux))
    {
    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 * profit->pixstep*profit->pixstep
		/ (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))
    {
    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 * profit->pixstep*profit->pixstep
	/ (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->param,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))
    {
    pprofit = *profit;
    memset(pprofit.paramindex, 0, PARAM_NPARAM*sizeof(int));
    memset(pprofit.paramlist, 0, PARAM_NPARAM*sizeof(float *));
    pprofit.nparam = 0;
    QMALLOC(pprofit.prof, profstruct *, 1);
    pprofit.prof[0] = prof_init(&pprofit, PROF_DIRAC);
    QMALLOC16(pprofit.covar, float, pprofit.nparam*pprofit.nparam);
    pprofit.nprof = 1;
    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];
      }
    pprofit.paraminit[pprofit.paramindex[PARAM_DISK_FLUX]] = profit->flux;
    pprofit.niter = profit_minimize(&pprofit, PROFIT_MAXITER);
    profit_residuals(&pprofit,field,wfield, 10.0, pprofit.param,pprofit.resi);
    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))
      {
      if (profit->flux > 0.0 && pprofit.flux > 0.0)
        obj2->prof_concentration = -2.5*log10(pprofit.flux / profit->flux);
      else  if (profit->flux > 0.0)
        obj2->prof_concentration = 99.0;
      else  if (pprofit.flux > 0.0)
        obj2->prof_concentration = -99.0;
      if (FLAG(obj2.prof_concentrationerr))
        obj2->prof_concentrationerr = (obj2->flux_prof > 0.0?
		1.086*(obj2->fluxerr_prof / obj2->flux_prof) : 99.0);
      }
    prof_end(pprofit.prof[0]);
    free(pprofit.prof);
    free(pprofit.covar);
    }

/* clean up. */
  fft_reset();
  free(profit->modpix);
  free(profit->psfpix);
  free(profit->pmodpix);
  free(profit->lmodpix);
  free(profit->objpix);
  free(profit->objweight);
  free(profit->resi);

  return;
  }


/****i* prof_gammainc *********************************************************
PROTO	double prof_gammainc(double x, double a)
PURPOSE	Returns the incomplete Gamma function (from Num. Recipes in C, p.216).
INPUT	A double,
	upper integration limit.
OUTPUT	Incomplete Gamma function.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	18/09/009
*/
static double	prof_gammainc (double x, double a)

  {
   double	b,c,d,h, xn,xp, del,sum;
   int		i;

  if (a < 0.0 || x <= 0.0)
    return 0.0;

  if (a < (x+1.0))
    {
/*-- Use the series representation */
    xp = x;
    del = sum = 1.0/x;
    for (i=100;i--;)	/* Iterate to convergence */
      {
      sum += (del *= a/(++xp));
      if (fabs(del) < fabs(sum)*3e-7)
        return sum*exp(-a+x*log(a)) / prof_gamma(x);
      }
    }
  else
    {
/*-- Use the continued fraction representation and take its complement */
    b = a + 1.0 - x;
    c = 1e30;
    h = d = 1.0/b;
    for (i=1; i<=100; i++)	/* Iterate to convergence */
      {
      xn = -i*(i-x);
      b += 2.0;
      if (fabs(d=xn*d+b) < 1e-30)
        d = 1e-30;
      if (fabs(c=b+xn/c) < 1e-30)
        c = 1e-30;
      del= c * (d = 1.0/d);
      h *= del;
      if (fabs(del-1.0) < 3e-7)
        return 1.0 - exp(-a+x*log(a))*h / prof_gamma(x);
      }
    }
  error(EXIT_FAILURE, "*Error*: out of bounds in ",
		"prof_gammainc()");
  return 0.0;
  }


/****i* prof_gamma ************************************************************
PROTO	double prof_gamma(double xx)
PURPOSE	Returns the Gamma function (from Num. Recipes in C, p.213).
INPUT	A double.
OUTPUT	Gamma function.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	11/09/009
*/
static double	prof_gamma(double xx)

  {
   double		x,tmp,ser;
   static double	cof[6]={76.18009173,-86.50532033,24.01409822,
			-1.231739516,0.120858003e-2,-0.536382e-5};
   int			j;

  tmp=(x=xx-1.0)+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<6;j++)
    ser += cof[j]/(x+=1.0);

  return 2.50662827465*ser*exp(-tmp);
  }


/****** profit_minradius ******************************************************
PROTO	float profit_minradius(profitstruct *profit, float refffac)
PURPOSE	Returns the minimum disk radius that guarantees that each and
	every model component fits within some margin in that disk.
INPUT	Profit structure pointer,
	margin in units of (r/r_eff)^(1/n)).
OUTPUT	Radius (in pixels).
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	21/09/009
*/
float	profit_minradius(profitstruct *profit, float refffac)

  {
   double	r,reff,rmax;
   int		p;

  rmax = reff = 0.0;
  for (p=0; p<profit->nprof; p++)
    {
    switch (profit->prof[p]->code)
      {
      case PROF_SERSIC:
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
      break;
      case PROF_DEVAUCOULEURS:
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
       break;
      case PROF_EXPONENTIAL:
        reff = *profit->paramlist[PARAM_DISK_SCALE]*1.67835;
      break;
      default:
        error(EXIT_FAILURE, "*Internal Error*: Unknown profile parameter in ",
		"profit_minradius()");
      break;
      }
    r = reff*(double)refffac;
    if (r>rmax)
      rmax = r;
    }

  return (float)rmax;
  }


/****** profit_psf ************************************************************
PROTO	void	profit_psf(profitstruct *profit)
PURPOSE	Build the local PSF at a given resolution.
INPUT	Profile-fitting structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	05/10/2009
 ***/
void	profit_psf(profitstruct *profit)
  {
   double	flux;
   float	posin[2], posout[2], dnaxisn[2],
		*pixout,
		xcout,ycout, xcin,ycin, invpixstep, norm;
   int		d,i;

  psf = profit->psf;
  psf_build(psf);

  xcout = (float)(profit->modnaxisn[0]/2) + 1.0;	/* FITS convention */
  ycout = (float)(profit->modnaxisn[1]/2) + 1.0;	/* FITS convention */
  xcin = (psf->masksize[0]/2) + 1.0;			/* FITS convention */
  ycin = (psf->masksize[1]/2) + 1.0;			/* FITS convention */
  invpixstep = profit->pixstep / psf->pixstep;

/* Initialize multi-dimensional counters */
  for (d=0; d<2; d++)
    {
    posout[d] = 1.0;					/* FITS convention */
    dnaxisn[d] = profit->modnaxisn[d]+0.5;
    }

/* Remap each pixel */
  pixout = profit->psfpix;
  flux = 0.0;
  for (i=profit->modnaxisn[0]*profit->modnaxisn[1]; i--;)
    {
    posin[0] = (posout[0] - xcout)*invpixstep + xcin;
    posin[1] = (posout[1] - ycout)*invpixstep + ycin;
    flux += ((*(pixout++) = interpolate_pix(posin, psf->maskloc,
		psf->masksize, INTERP_LANCZOS3)));
    for (d=0; d<2; d++)
      if ((posout[d]+=1.0) < dnaxisn[d])
        break;
      else
        posout[d] = 1.0;
    }

/* Normalize PSF flux (just in case...) */
  flux *= profit->pixstep*profit->pixstep;
  if (fabs(flux) > 0.0)
    {
    norm = 1.0/flux;
    pixout = profit->psfpix;
    for (i=profit->modnaxisn[0]*profit->modnaxisn[1]; i--;)
      *(pixout++) *= norm;
    }

  return;
  }


/****** profit_minimize *******************************************************
PROTO	void profit_minimize(profitstruct *profit)
PURPOSE	Provide a function returning residuals to lmfit.
INPUT	Pointer to the profit structure involved in the fit,
	maximum number of iterations.
OUTPUT	Number of iterations used.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	23/05/2008
 ***/
int	profit_minimize(profitstruct *profit, int niter)
  {
   float		lm_opts[5], info[LM_INFO_SZ];
   int			m,n;

/* Allocate work space */
  n = profit->nparam;
  m = profit->nresi;

  memset(profit->resi, 0, profit->nresi*sizeof(float));
  memset(profit->covar, 0, profit->nparam*profit->nparam*sizeof(float));
  profit_boundtounbound(profit, profit->paraminit);

/* Perform fit */
  lm_opts[0] = 1.0e-3;
  lm_opts[1] = 1.0e-17;
  lm_opts[2] = 1.0e-17;
  lm_opts[3] = 1.0e-17;
  lm_opts[4] = 1.0e-6;

  niter = slevmar_dif(profit_evaluate, profit->paraminit, profit->resi,
	n, m, niter, lm_opts, info, NULL, profit->covar, profit);

  profit_unboundtobound(profit, profit->paraminit);


  return niter;
  }


/****** profit_printout *******************************************************
PROTO	void profit_printout(int n_par, float* par, int m_dat, float* fvec,
		void *data, int iflag, int iter, int nfev )
PURPOSE	Provide a function to print out results to lmfit.
INPUT	Number of fitted parameters,
	pointer to the vector of parameters,
	number of data points,
	pointer to the vector of residuals (output),
	pointer to the data structure (unused),
	0 (init) 1 (outer loop) 2(inner loop) -1(terminated),
	outer loop counter,
	number of calls to evaluate().
OUTPUT	-.
NOTES	Input arguments are there only for compatibility purposes (unused)
AUTHOR	E. Bertin (IAP)
VERSION	17/09/2008
 ***/
void	profit_printout(int n_par, float* par, int m_dat, float* fvec,
		void *data, int iflag, int iter, int nfev )
  {
   checkstruct	*check;
   profitstruct	*profit;
   char		filename[256];
   static int	itero;

  profit = (profitstruct *)data;

  if (0 && (iter!=itero || iter<0))
    {
    if (iter<0)
      itero++;
    else
      itero = iter;
    sprintf(filename, "check_%d_%04d.fits", the_gal, itero);
    check=initcheck(filename, CHECK_PROFILES, 0);
    reinitcheck(the_field, check);
    addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],
		profit->ix,profit->iy, 1.0);

    reendcheck(the_field, check);
    endcheck(check);
    }

  return;
  }


/****** profit_evaluate ******************************************************
PROTO	void profit_evaluate(float *par, float *fvec, int m, int n,
		void *adata)
PURPOSE	Provide a function returning residuals to levmar.
INPUT	Pointer to the vector of parameters,
	pointer to the vector of residuals (output),
	number of model parameters,
	number of data points,
	pointer to a data structure (unused).
OUTPUT	-.
NOTES	Input arguments are there only for compatibility purposes (unused)
AUTHOR	E. Bertin (IAP)
VERSION	18/09/2008
 ***/
void	profit_evaluate(float *par, float *fvec, int m, int n,
			void *adata)
  {
   profitstruct	*profit;

  profit = (profitstruct *)adata;
  profit_unboundtobound(profit, par);
  profit_residuals(profit, the_field, the_wfield, PROFIT_DYNPARAM, par, fvec);
  profit_boundtounbound(profit, par);
  profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
  return;
  }


/****** profit_residuals ******************************************************
PROTO	float *prof_residuals(profitstruct *profit, picstruct *field,
		picstruct *wfield, float dynparam, float *param, float *resi)
PURPOSE	Compute the vector of residuals between the data and the galaxy
	profile model.
INPUT	Profile-fitting structure,
	pointer to the field,
	pointer to the field weight,
	dynamic compression parameter (0=no compression),
	pointer to the model parameters (output),
	pointer to the computed residuals (output).
OUTPUT	Vector of residuals.
NOTES	-.
AUTHOR	E. Bertin (IAP)