Astromatic - public software.sextractor - Blame - Rev 257

[/] [trunk/] [src/] [profit.c] - Blame information for rev 257

Go to most recent revision | Details | Compare with Previous | View Log


/*
*                               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:          25/07/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 31/05/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.peak_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 * 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))
    {
    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 * 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->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;
          var = 1.0/(weight*weight);
          }
        tvobj += pix;
        sigtvobj += var;
        tvmin += *(lmodpix++);
//        *(lmodpix++) = pix;
        }
      else
        tvmout += *(lmodpix++);
      }
    }
 
//  tv -= area*bkg;
 
  tvm = tvmin + tvmout;
  if (tvm != 0.0)
    {
    obj2->fluxcor_prof = tvobj+obj2->flux_prof*tvmout/tvm;
    obj2->fluxcorerr_prof = sqrt(sigtvobj
                        +obj2->fluxerr_prof*obj2->fluxerr_prof*tvmout/tvm);
    }
  else
    {
    obj2->fluxcor_prof = tvobj;
    obj2->fluxcorerr_prof = sqrt(sigtvobj);
    }
 
/*
  if ((check = prefs.check[CHECK_OTHER]))
    addcheck(check, profit->lmodpix, w, h, profit->ix,profit->iy, 1.0);
*/
  return;
  }
 
 
/****i* prof_gammainc *********************************************************
PROTO   double prof_gammainc(double x, double a)
PURPOSE Returns the incomplete Gamma function (based on algorithm described in
        Numerical Recipes in C, chap. 6.1).
INPUT   A double,
        upper integration limit.
OUTPUT  Incomplete Gamma function.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 08/10/2010
*/
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 (based on algorithm described in Numerical
        Recipes in C, chap 6.1).
INPUT   A double.
OUTPUT  Gamma function.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 11/09/2009
*/
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 08/10/2010
*/
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 MODEL_BACK:
      case MODEL_DIRAC:
        reff = 0.0;
      break;
      case MODEL_SERSIC:
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
        break;
      case MODEL_DEVAUCOULEURS:
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
        break;
      case MODEL_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 07/07/2010
 ***/
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)
    error(EXIT_FAILURE, "*Error*: PSF model is empty or negative: ", psf->name);
 
  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 20/05/2011
 ***/
int     profit_minimize(profitstruct *profit, int niter)
  {
   double       lm_opts[5], info[LM_INFO_SZ],
                dcovar[PARAM_NPARAM*PARAM_NPARAM], dparam[PARAM_NPARAM];
   int          nfree;
 
  profit->iter = 0;
  memset(dcovar, 0, profit->nparam*profit->nparam*sizeof(double));
 
/* Perform fit */
  lm_opts[0] = 1.0e-3;           /* Initial mu */
  lm_opts[1] = 1.0e-6;          /* ||J^T e||_inf stopping factor */
  lm_opts[2] = 1.0e-6;          /* |Dp||_2 stopping factor */
  lm_opts[3] = 1.0e-6;          /* ||e||_2 stopping factor */
  lm_opts[4] = 1.0e-4;          /* Jacobian step */
 
  nfree = profit_boundtounbound(profit, profit->paraminit, dparam,
                                PARAM_ALLPARAMS);
 
  niter = dlevmar_dif(profit_evaluate, dparam, NULL, nfree, profit->nresi,
                        niter, lm_opts, info, NULL, dcovar, profit);
 
  profit_unboundtobound(profit, dparam, profit->paraminit, PARAM_ALLPARAMS);
 
/* Convert covariance matrix to bounded space */
  profit_covarunboundtobound(profit, dcovar, profit->covar);
 
  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),
 
        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(double *par, double *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 (we use it for the profit structure here).
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 20/05/2011
 ***/
void    profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
  {
   profitstruct         *profit;
   profstruct           **prof;
   double               *dpar0, *dresi;
   float                *modpixt, *profpixt, *resi,
                        tparam, val;
   PIXTYPE              *lmodpixt,*lmodpix2t, *objpix,*weight,
                        wval;
   int                  c,f,i,p,q, fd,pd, jflag,sflag, nprof;
 
  profit = (profitstruct *)adata;
 
/* Detect "Jacobian-related" calls */
  jflag = pd = fd = 0;
  dpar0 = profit->dparam;
  if (profit->iter)
    {
    f = q = 0;
    for (p=0; p<profit->nparam; p++)
      {
      if (dpar[f] - dpar0[f] != 0.0)
        {
        pd = p;
        fd = f;
        q++;
        }
      if (profit->parfittype[p]!=PARFIT_FIXED)
        f++;
      }
    if (f>0 && q==1)
      jflag = 1;
    }
jflag = 0;       /* Temporarily deactivated (until problems are fixed) */
  if (jflag && !(profit->nprof==1 && profit->prof[0]->code == MODEL_DIRAC))
    {
    prof = profit->prof;
    nprof = profit->nprof;
 
/*-- "Jacobian call" */
    tparam = profit->param[pd];
    profit_unboundtobound(profit, &dpar[fd], &profit->param[pd], pd);
    sflag = 1;
    switch(profit->paramrevindex[pd])
      {
      case PARAM_BACK:
        lmodpixt = profit->lmodpix;
        lmodpix2t = profit->lmodpix2;
        val = (profit->param[pd] - tparam);
        for (i=profit->nobjpix;i--;)
          *(lmodpix2t++) = val;
        break;
      case PARAM_X:
      case PARAM_Y:
        profit_resample(profit, profit->cmodpix, profit->lmodpix2, 1.0);
        lmodpixt = profit->lmodpix;
        lmodpix2t = profit->lmodpix2;
        for (i=profit->nobjpix;i--;)
          *(lmodpix2t++) -= *(lmodpixt++);
        break;
      case PARAM_DIRAC_FLUX:
      case PARAM_SPHEROID_FLUX:
      case PARAM_DISK_FLUX:
      case PARAM_ARMS_FLUX:
      case PARAM_BAR_FLUX:
        if (nprof==1 && tparam != 0.0)
          {
          lmodpixt = profit->lmodpix;
          lmodpix2t = profit->lmodpix2;
          val = (profit->param[pd] - tparam) / tparam;
          for (i=profit->nobjpix;i--;)
            *(lmodpix2t++) = val**(lmodpixt++);
          }
        else
          {
          for (c=0; c<nprof; c++)
            if (prof[c]->flux == &profit->param[pd])
              break;
          memcpy(profit->modpix, prof[c]->pix, profit->nmodpix*sizeof(float));
          profit_convolve(profit, profit->modpix);
          profit_resample(profit, profit->modpix, profit->lmodpix2,
                profit->param[pd] - tparam);
          }
        break;
      case PARAM_SPHEROID_REFF:
      case PARAM_SPHEROID_ASPECT:
      case PARAM_SPHEROID_POSANG:
      case PARAM_SPHEROID_SERSICN:
        sflag = 0;                       /* We are in the same switch */
        for (c=0; c<nprof; c++)
          if (prof[c]->code == MODEL_SERSIC
                || prof[c]->code == MODEL_DEVAUCOULEURS)
            break;
      case PARAM_DISK_SCALE:
      case PARAM_DISK_ASPECT:
      case PARAM_DISK_POSANG:
        if (sflag)
          for (c=0; c<nprof; c++)
            if (prof[c]->code == MODEL_EXPONENTIAL)
              break;
        sflag = 0;
      case PARAM_ARMS_QUADFRAC:
      case PARAM_ARMS_SCALE:
      case PARAM_ARMS_START:
      case PARAM_ARMS_POSANG:
      case PARAM_ARMS_PITCH:
      case PARAM_ARMS_PITCHVAR:
      case PARAM_ARMS_WIDTH:
        if (sflag)
          for (c=0; c<nprof; c++)
            if (prof[c]->code == MODEL_ARMS)
              break;
        sflag = 0;
      case PARAM_BAR_ASPECT:
      case PARAM_BAR_POSANG:
        if (sflag)
          for (c=0; c<nprof; c++)
            if (prof[c]->code == MODEL_ARMS)
              break;
        modpixt = profit->modpix;
        profpixt = prof[c]->pix;
        val = -*prof[c]->flux;
        for (i=profit->nmodpix;i--;)
          *(modpixt++) = val**(profpixt++);
        memcpy(profit->modpix2, prof[c]->pix, profit->nmodpix*sizeof(float));
        prof_add(profit, prof[c], 0);
        memcpy(prof[c]->pix, profit->modpix2, profit->nmodpix*sizeof(float));
        profit_convolve(profit, profit->modpix);
        profit_resample(profit, profit->modpix, profit->lmodpix2, 1.0);
        break;
      default:
        error(EXIT_FAILURE, "*Internal Error*: ",
                        "unknown parameter index in profit_jacobian()");
        break;
      }
    objpix = profit->objpix;
    weight = profit->objweight;
    lmodpixt = profit->lmodpix;
    lmodpix2t = profit->lmodpix2;
    resi = profit->resi;
    dresi = fvec;
    if (PROFIT_DYNPARAM > 0.0)
      for (i=profit->nobjpix;i--; lmodpixt++, lmodpix2t++)
        {
        val = *(objpix++);
        if ((wval=*(weight++))>0.0)
          *(dresi++) = *(resi++) + *lmodpix2t
                * wval/(1.0+wval*fabs(*lmodpixt - val)/PROFIT_DYNPARAM);
        }
    else
      for (i=profit->nobjpix;i--; lmodpix2t++)
        if ((wval=*(weight++))>0.0)
          *(dresi++) = *(resi++) + *lmodpix2t * wval;
    }
  else
    {
/*-- "Regular call" */
    for (p=0; p<profit->nparam; p++)
      dpar0[p] = dpar[p];
    profit_unboundtobound(profit, dpar, profit->param, PARAM_ALLPARAMS);
 
    profit_residuals(profit, the_field, the_wfield, PROFIT_DYNPARAM,
        profit->param, profit->resi);
 
    for (p=0; p<profit->nresi; p++)
      fvec[p] = profit->resi[p];
    }
 
//  profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
  profit->iter++;
 
  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)
VERSION 08/07/2010
 ***/
float   *profit_residuals(profitstruct *profit, picstruct *field,
                picstruct *wfield, float dynparam, float *param, float *resi)
  {
   int          p, nmodpix;
 
  nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float);
  memset(profit->modpix, 0, nmodpix);
  for (p=0; p<profit->nparam; p++)
    profit->param[p] = param[p];
/* Simple PSF shortcut */
  if (profit->nprof == 1 && profit->prof[0]->code == MODEL_DIRAC)
    {
    profit_resample(profit, profit->psfpix, profit->lmodpix,
                *profit->prof[0]->flux);
    profit->flux = *profit->prof[0]->flux;
    }
  else
    {
    profit->flux = 0.0;
    for (p=0; p<profit->nprof; p++)
      profit->flux += prof_add(profit, profit->prof[p], 0);
    memcpy(profit->cmodpix, profit->modpix, profit->nmodpix*sizeof(float));
    profit_convolve(profit, profit->cmodpix);
    profit_resample(profit, profit->cmodpix, profit->lmodpix, 1.0);
    }
 
  if (resi)
    profit_compresi(profit, dynparam, resi);
 
  return resi;
  }
 
 
/****** profit_compresi ******************************************************
PROTO   float *prof_compresi(profitstruct *profit, float dynparam,
                        float *resi)
PURPOSE Compute the vector of residuals between the data and the galaxy
        profile model.
INPUT   Profile-fitting structure,
        dynamic-compression parameter (0=no compression),
        vector of residuals (output).
OUTPUT  Vector of residuals.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 07/09/2009
 ***/
float   *profit_compresi(profitstruct *profit, float dynparam, float *resi)
  {
   double       error;
   float        *resit;
   PIXTYPE      *objpix, *objweight, *lmodpix,
                val,val2,wval, invsig;
   int          npix, i;
 
/* Compute vector of residuals */
  resit = resi;
  objpix = profit->objpix;
  objweight = profit->objweight;
  lmodpix = profit->lmodpix;
  error = 0.0;
  npix = profit->objnaxisn[0]*profit->objnaxisn[1];
  if (dynparam > 0.0)
    {
    invsig = (PIXTYPE)(1.0/dynparam);
    for (i=npix; i--; lmodpix++)
      {
      val = *(objpix++);
      if ((wval=*(objweight++))>0.0)
        {
        val2 = (*lmodpix - val)*wval*invsig;
        val2 = val2>0.0? logf(1.0+val2) : -logf(1.0-val2);
        *(resit++) = val2*dynparam;
        error += val2*val2;
        }
      }
    profit->chi2 = dynparam*dynparam*error;
    }
  else
    {
    for (i=npix; i--; lmodpix++)
      {
      val = *(objpix++);
      if ((wval=*(objweight++))>0.0)
        {
        val2 = (*lmodpix - val)*wval;
        *(resit++) = val2;
        error += val2*val2;
        }
      }
    profit->chi2 = error;
    }
 
  return resi;
  }
 
 
/****** profit_resample ******************************************************
PROTO   int     prof_resample(profitstruct *profit, float *inpix,
                PIXTYPE *outpix, float factor)
PURPOSE Resample the current full resolution model to image resolution.
INPUT   Profile-fitting structure,
        pointer to input raster,
        pointer to output raster,
        multiplicating factor.
OUTPUT  RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 12/09/2010
 ***/
int     profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
        float factor)
  {
   PIXTYPE      *pixout,*pixout0;
   float        *pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0,
                *dpixout,*dpixout0, *dx,*dy,
                xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val,
                invpixstep, norm;
   int          *start,*startt, *nmask,*nmaskt,
                i,j,k,n,t,
                ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout,
                iysina, nyin, hmw,hmh, ix,iy, ixin,iyin;
 
  invpixstep = profit->subsamp/profit->pixstep;
 
  xcin = (profit->modnaxisn[0]/2);
  xcout = (float)(profit->objnaxisn[0]/2);
  if ((dx=(profit->paramlist[PARAM_X])))
    xcout += *dx;
 
  xsin = xcin - xcout*invpixstep;                       /* Input start x-coord*/
  if ((int)xsin >= profit->modnaxisn[0])
    return RETURN_ERROR;
  ixsout = 0;                            /* Int. part of output start x-coord */
  if (xsin<0.0)
    {
    dixout = (int)(1.0-xsin/invpixstep);
/*-- Simply leave here if the images do not overlap in x */
    if (dixout >= profit->objnaxisn[0])
      return RETURN_ERROR;
    ixsout += dixout;
    xsin += dixout*invpixstep;
    }
  nxout = (int)((profit->modnaxisn[0]-xsin)/invpixstep);/* nb of interpolated
                                                        input pixels along x */
  if (nxout>(ixout=profit->objnaxisn[0]-ixsout))
    nxout = ixout;
  if (!nxout)
    return RETURN_ERROR;
 
  ycin = (profit->modnaxisn[1]/2);
  ycout = (float)(profit->objnaxisn[1]/2);
  if ((dy=(profit->paramlist[PARAM_Y])))
    ycout += *dy;
 
  ysin = ycin - ycout*invpixstep;               /* Input start y-coord*/
  if ((int)ysin >= profit->modnaxisn[1])
    return RETURN_ERROR;
  iysout = 0;                            /* Int. part of output start y-coord */
  if (ysin<0.0)
    {
    diyout = (int)(1.0-ysin/invpixstep);
/*-- Simply leave here if the images do not overlap in y */
    if (diyout >= profit->objnaxisn[1])
      return RETURN_ERROR;
    iysout += diyout;
    ysin += diyout*invpixstep;
    }
  nyout = (int)((profit->modnaxisn[1]-ysin)/invpixstep);/* nb of interpolated
                                                        input pixels along y */
  if (nyout>(iyout=profit->objnaxisn[1]-iysout))
    nyout = iyout;
  if (!nyout)
    return RETURN_ERROR;
 
/* Set the yrange for the x-resampling with some margin for interpolation */
  iysina = (int)ysin;   /* Int. part of Input start y-coord with margin */
  hmh = INTERPW/2 - 1;  /* Interpolant start */
  if (iysina<0 || ((iysina -= hmh)< 0))
    iysina = 0;
  nyin = (int)(ysin+nyout*invpixstep)+INTERPW-hmh;/* Interpolated Input y size*/
  if (nyin>profit->modnaxisn[1])                                                /* with margin */
    nyin = profit->modnaxisn[1];
/* Express everything relative to the effective Input start (with margin) */
  nyin -= iysina;
  ysin -= (float)iysina;
 
/* Allocate interpolant stuff for the x direction */
  QMALLOC(mask, float, nxout*INTERPW);  /* Interpolation masks */
  QMALLOC(nmask, int, nxout);           /* Interpolation mask sizes */
  QMALLOC(start, int, nxout);           /* Int. part of Input conv starts */
/* Compute the local interpolant and data starting points in x */
  hmw = INTERPW/2 - 1;
  xin = xsin;
  maskt = mask;
  nmaskt = nmask;
  startt = start;
  for (j=nxout; j--; xin+=invpixstep)
    {
    ix = (ixin=(int)xin) - hmw;
    dxm = ixin - xin - hmw;     /* starting point in the interpolation func */
    if (ix < 0)
      {
      n = INTERPW+ix;
      dxm -= (float)ix;
      ix = 0;
      }
    else
      n = INTERPW;
    if (n>(t=profit->modnaxisn[0]-ix))
      n=t;
    *(startt++) = ix;
    *(nmaskt++) = n;
    norm = 0.0;
    for (x=dxm, i=n; i--; x+=1.0)
      norm += (*(maskt++) = INTERPF(x));
    norm = norm>0.0? 1.0/norm : 1.0;
    maskt -= n;
    for (i=n; i--;)
      *(maskt++) *= norm;
    }
 
  QCALLOC(pixinout, float, nxout*nyin); /* Intermediary frame-buffer */
 
/* Make the interpolation in x (this includes transposition) */
  pixin0 = inpix + iysina*profit->modnaxisn[0];
  dpixout0 = pixinout;
  for (k=nyin; k--; pixin0+=profit->modnaxisn[0], dpixout0++)
    {
    maskt = mask;
    nmaskt = nmask;
    startt = start;
    dpixout = dpixout0;
    for (j=nxout; j--; dpixout+=nyin)
      {
      pixin = pixin0+*(startt++);
      val = 0.0;
      for (i=*(nmaskt++); i--;)
        val += *(maskt++)**(pixin++);
      *dpixout = val;
      }
    }
 
/* Reallocate interpolant stuff for the y direction */
  QREALLOC(mask, float, nyout*INTERPW); /* Interpolation masks */
  QREALLOC(nmask, int, nyout);                  /* Interpolation mask sizes */
  QREALLOC(start, int, nyout);          /* Int. part of Input conv starts */
 
/* Compute the local interpolant and data starting points in y */
  hmh = INTERPW/2 - 1;
  yin = ysin;
  maskt = mask;
  nmaskt = nmask;
  startt = start;
  for (j=nyout; j--; yin+=invpixstep)
    {
    iy = (iyin=(int)yin) - hmh;
    dym = iyin - yin - hmh;     /* starting point in the interpolation func */
    if (iy < 0)
      {
      n = INTERPW+iy;
      dym -= (float)iy;
      iy = 0;
      }
    else
      n = INTERPW;
    if (n>(t=nyin-iy))
      n=t;
    *(startt++) = iy;
    *(nmaskt++) = n;
    norm = 0.0;
    for (y=dym, i=n; i--; y+=1.0)
      norm += (*(maskt++) = INTERPF(y));
    norm = norm>0.0? 1.0/norm : 1.0;
    maskt -= n;
    for (i=n; i--;)
      *(maskt++) *= norm;
    }
 
/* Initialize destination buffer to zero */
  memset(outpix, 0, (size_t)profit->nobjpix*sizeof(PIXTYPE));
 
/* Make the interpolation in y and transpose once again */
  dpixin0 = pixinout;
  pixout0 = outpix+ixsout+iysout*profit->objnaxisn[0];
  for (k=nxout; k--; dpixin0+=nyin, pixout0++)
    {
    maskt = mask;
    nmaskt = nmask;
    startt = start;
    pixout = pixout0;
    for (j=nyout; j--; pixout+=profit->objnaxisn[0])
      {
      dpixin = dpixin0+*(startt++);
      val = 0.0;
      for (i=*(nmaskt++); i--;)
        val += *(maskt++)**(dpixin++);
       *pixout = (PIXTYPE)(factor*val);
      }
    }
 
/* Free memory */
  free(pixinout);
  free(mask);
  free(nmask);
  free(start);
 
  return RETURN_OK;
  }
 
 
/****** profit_convolve *******************************************************
PROTO   void profit_convolve(profitstruct *profit, float *modpix)
PURPOSE Convolve a model image with the local PSF.
INPUT   Pointer to the profit structure,
        Pointer to the image raster.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 15/09/2008
 ***/
void    profit_convolve(profitstruct *profit, float *modpix)
  {
  if (!profit->psfdft)
    profit_makedft(profit);
 
  fft_conv(modpix, profit->psfdft, profit->modnaxisn);
 
  return;
  }
 
 
/****** profit_makedft *******************************************************
PROTO   void profit_makedft(profitstruct *profit)
PURPOSE Create the Fourier transform of the descrambled PSF component.
INPUT   Pointer to the profit structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 22/04/2008
 ***/
void    profit_makedft(profitstruct *profit)
  {
   psfstruct    *psf;
   float      *mask,*maskt, *ppix;
   float       dx,dy, r,r2,rmin,rmin2,rmax,rmax2,rsig,invrsig2;
   int          width,height,npix,offset, psfwidth,psfheight,psfnpix,
                cpwidth, cpheight,hcpwidth,hcpheight, i,j,x,y;
 
  if (!(psf=profit->psf))
    return;
 
  psfwidth = profit->modnaxisn[0];
  psfheight = profit->modnaxisn[1];
  psfnpix = psfwidth*psfheight;
  width = profit->modnaxisn[0];
  height = profit->modnaxisn[1];
  npix = width*height;
  QCALLOC(mask, float, npix);
  cpwidth = (width>psfwidth)?psfwidth:width;
  hcpwidth = cpwidth>>1;
  cpwidth = hcpwidth<<1;
  offset = width - cpwidth;
  cpheight = (height>psfheight)?psfheight:height;
  hcpheight = cpheight>>1;
  cpheight = hcpheight<<1;
 
/* Frame and descramble the PSF data */
  ppix = profit->psfpix + (psfheight/2)*psfwidth + psfwidth/2;
  maskt = mask;
  for (j=hcpheight; j--; ppix+=psfwidth)
    {
    for (i=hcpwidth; i--;)
      *(maskt++) = *(ppix++);
    ppix -= cpwidth;
    maskt += offset;
    for (i=hcpwidth; i--;)
      *(maskt++) = *(ppix++);
    }
 
  ppix = profit->psfpix + ((psfheight/2)-hcpheight)*psfwidth + psfwidth/2;
  maskt += width*(height-cpheight);
  for (j=hcpheight; j--; ppix+=psfwidth)
    {
    for (i=hcpwidth; i--;)
      *(maskt++) = *(ppix++);
    ppix -= cpwidth;
    maskt += offset;
    for (i=hcpwidth; i--;)
      *(maskt++) = *(ppix++);
    }
 
/* Truncate to a disk that has diameter = (box width) */
  rmax = cpwidth - 1.0 - hcpwidth;
  if (rmax > (r=hcpwidth))
    rmax = r;
  if (rmax > (r=cpheight-1.0-hcpheight))
    rmax = r;
  if (rmax > (r=hcpheight))
    rmax = r;
  if (rmax<1.0)
    rmax = 1.0;
  rmax2 = rmax*rmax;
  rsig = psf->fwhm/profit->pixstep;
  invrsig2 = 1/(2*rsig*rsig);
  rmin = rmax - (3*rsig);     /* 3 sigma annulus (almost no aliasing) */
  rmin2 = rmin*rmin;
 
  maskt = mask;
  dy = 0.0;
  for (y=hcpheight; y--; dy+=1.0)
    {
    dx = 0.0;
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
      if ((r2=dx*dx+dy*dy)>rmin2)
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
    dx = -hcpwidth;
    maskt += offset;
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
      if ((r2=dx*dx+dy*dy)>rmin2)
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
    }
  dy = -hcpheight;
  maskt += width*(height-cpheight);
  for (y=hcpheight; y--; dy+=1.0)
    {
    dx = 0.0;
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
      if ((r2=dx*dx+dy*dy)>rmin2)
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
    dx = -hcpwidth;
    maskt += offset;
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
      if ((r2=dx*dx+dy*dy)>rmin2)
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
    }
 
/* Finally move to Fourier space */
  profit->psfdft = fft_rtf(mask, profit->modnaxisn);
 
  free(mask);
 
  return;
  }
 
 
/****** profit_copyobjpix *****************************************************
PROTO   int profit_copyobjpix(profitstruct *profit, picstruct *field,
                        picstruct *wfield)
PURPOSE Copy a piece of the input field image to a profit structure.
INPUT   Pointer to the profit structure,
        Pointer to the field structure,
        Pointer to the field weight structure.
OUTPUT  The number of valid pixels copied.
NOTES   Global preferences are used.
AUTHOR  E. Bertin (IAP)
VERSION 01/12/2009
 ***/
int     profit_copyobjpix(profitstruct *profit, picstruct *field,
                        picstruct *wfield)
  {
   float        dx, dy2, dr2, rad2;
   PIXTYPE      *pixin,*spixin, *wpixin,*swpixin, *pixout,*wpixout,
                backnoise2, invgain, satlevel, wthresh, pix,spix, wpix,swpix;
   int          i,x,y, xmin,xmax,ymin,ymax, w,h,dw, npix, off, gainflag,
                badflag, sflag, sx,sy,sn,sw, ix,iy;
 
/* First put the image background to -BIG */
  pixout = profit->objpix;
  wpixout = profit->objweight;
  for (i=profit->objnaxisn[0]*profit->objnaxisn[1]; i--;)
    {
    *(pixout++) = -BIG;
    *(wpixout++) = 0.0;
    }
 
/* Don't go further if out of frame!! */
  ix = profit->ix;
  iy = profit->iy;
  if (ix<0 || ix>=field->width || iy<field->ymin || iy>=field->ymax)
    return 0;
 
  backnoise2 = field->backsig*field->backsig;
  sn = (int)profit->subsamp;
  sflag = (sn>1);
  w = profit->objnaxisn[0]*sn;
  h = profit->objnaxisn[1]*sn;
  if (sflag)
    backnoise2 *= (PIXTYPE)sn;
  invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0;
  satlevel = field->satur_level - profit->obj->bkg;
  rad2 = h/2.0;
  if (rad2 > w/2.0)
    rad2 = w/2.0;
  rad2 *= rad2;
 
/* Set the image boundaries */
  pixout = profit->objpix;
  wpixout = profit->objweight;
  ymin = iy-h/2;
  ymax = ymin + h;
  if (ymin<field->ymin)
    {
    off = (field->ymin-ymin-1)/sn + 1;
    pixout += off*profit->objnaxisn[0];
    wpixout += off*profit->objnaxisn[0];
    ymin += off*sn;
    }
  if (ymax>field->ymax)
    ymax -= ((ymax-field->ymax-1)/sn + 1)*sn;
 
  xmin = ix-w/2;
  xmax = xmin + w;
  dw = 0;
  if (xmax>field->width)
    {
    off = (xmax-field->width-1)/sn + 1;
    dw += off;
    xmax -= off*sn;
    }
  if (xmin<0)
    {
    off = (-xmin-1)/sn + 1;
    pixout += off;
    wpixout += off;
    dw += off;
    xmin += off*sn;
    }
/* Make sure the input frame size is a multiple of the subsampling step */
  if (sflag)
    {
/*
    if (((rem=ymax-ymin)%sn))
      {
      ymin += rem/2;
      ymax -= (rem-rem/2);
      }
    if (((rem=xmax-xmin)%sn))
      {
      xmin += rem/2;
      pixout += rem/2;
      wpixout += rem/2;
      dw += rem;
      xmax -= (rem-rem/2);
      }
*/
    sw = field->width;
    }
 
/* Copy the right pixels to the destination */
  npix = 0;
  if (wfield)
    {
    wthresh = wfield->weight_thresh;
    gainflag = prefs.weightgain_flag;
    if (sflag)
      {
/*---- Sub-sampling case */
      for (y=ymin; y<ymax; y+=sn, pixout+=dw,wpixout+=dw)
        {
        for (x=xmin; x<xmax; x+=sn)
          {
          pix = wpix = 0.0;
          badflag = 0;
          for (sy=0; sy<sn; sy++)
            {
            dy2 = (y+sy-iy);
            dy2 *= dy2;
            dx = (x-ix);
            spixin = &PIX(field, x, y+sy);
            swpixin = &PIX(wfield, x, y+sy);
            for (sx=sn; sx--;)
              {
              dr2 = dy2 + dx*dx;
              dx++;
              spix = *(spixin++);
              swpix = *(swpixin++);
              if (dr2<rad2 && spix>-BIG && spix<satlevel && swpix<wthresh)
                {
                pix += spix;
                wpix += swpix;
                }
              else
                badflag=1;
              }
            }
          *(pixout++) = pix;
          if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
            {
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
                (gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
            npix++;
            }
          else
            *(wpixout++) = 0.0;
          }
        }
      }
    else
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
        {
        dy2 = y-iy;
        dy2 *= dy2;
        pixin = &PIX(field, xmin, y);
        wpixin = &PIX(wfield, xmin, y);
        for (x=xmin; x<xmax; x++)
          {
          dx = x-ix;
          dr2 = dy2 + dx*dx;
          pix = *(pixin++);
          wpix = *(wpixin++);
          if (dr2<rad2 && pix>-BIG && pix<satlevel && wpix<wthresh)
            {
            *(pixout++) = pix;
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
                (gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
            npix++;
            }
          else
            *(pixout++) = *(wpixout++) = 0.0;
          }
        }
    }
  else
    {
    if (sflag)
      {
/*---- Sub-sampling case */
      for (y=ymin; y<ymax; y+=sn, pixout+=dw, wpixout+=dw)
        {
        for (x=xmin; x<xmax; x+=sn)
          {
          pix = 0.0;
          badflag = 0;
          for (sy=0; sy<sn; sy++)
            {
            dy2 = y+sy-iy;
            dy2 *= dy2;
            dx = x-ix;
            spixin = &PIX(field, x, y+sy);
            for (sx=sn; sx--;)
              {
              dr2 = dy2 + dx*dx;
              dx++;
              spix = *(spixin++);
              if (dr2<rad2 && spix>-BIG && spix<satlevel)
                pix += spix;
              else
                badflag=1;
              }
            }
          *(pixout++) = pix;
          if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
            {
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain:0.0));
            npix++;
            }
          else
            *(wpixout++) = 0.0;
          }
        }
      }
    else
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
        {
        dy2 = y-iy;
        dy2 *= dy2;
        pixin = &PIX(field, xmin, y);
        for (x=xmin; x<xmax; x++)
          {
          dx = x-ix;
          dr2 = dy2 + dx*dx;
          pix = *(pixin++);
          if (dr2<rad2 && pix>-BIG && pix<satlevel)
            {
            *(pixout++) = pix;
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain : 0.0));
            npix++;
            }
          else
            *(pixout++) = *(wpixout++) = 0.0;
          }
        }
    }
 
  return npix;
  }
 
 
/****** profit_spiralindex ****************************************************
PROTO   float profit_spiralindex(profitstruct *profit)
PURPOSE Compute the spiral index of a galaxy image (positive for arms
        extending counter-clockwise and negative for arms extending CW, 0 for
        no spiral pattern).
INPUT   Profile-fitting structure.
OUTPUT  Vector of residuals.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 18/05/2011
 ***/
float profit_spiralindex(profitstruct *profit)
  {
   objstruct    *obj;
   obj2struct   *obj2;
   float        *dx,*dy, *fdx,*fdy, *gdx,*gdy, *gdxt,*gdyt, *pix,
                fwhm, invtwosigma2, hw,hh, ohw,ohh, x,y,xstart, tx,ty,txstart,
                gx,gy, r2, spirindex, invsig, val, sep;
   PIXTYPE      *fpix;
   int          i,j, npix;
 
  npix = profit->objnaxisn[0]*profit->objnaxisn[1];
 
  obj = profit->obj;
  obj2 = profit->obj2;
/* Compute simple derivative vectors at a fraction of the object scale */
  fwhm = profit->guessradius * 2.0 / 4.0;
  if (fwhm < 2.0)
    fwhm = 2.0;
  sep = 2.0;
 
  invtwosigma2 = -(2.35*2.35/(2.0*fwhm*fwhm));
  hw = (float)(profit->objnaxisn[0]/2);
  ohw = profit->objnaxisn[0] - hw;
  hh = (float)(profit->objnaxisn[1]/2);
  ohh = profit->objnaxisn[1] - hh;
  txstart = -hw;
  ty = -hh;
  QMALLOC(dx, float, npix);
  pix = dx;
  for (j=profit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=profit->objnaxisn[0]; i--; tx+=1.0)
      {
      x = tx < -0.5? tx + hw : tx - ohw;
      *(pix++) = exp(invtwosigma2*((x+sep)*(x+sep)+y*y))
                - exp(invtwosigma2*((x-sep)*(x-sep)+y*y));
      }
    }
  QMALLOC(dy, float, npix);
  pix = dy;
  ty = -hh;
  for (j=profit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=profit->objnaxisn[0]; i--; tx+=1.0)
      {
      x = tx < -0.5? tx + hw : tx - ohw;
      *(pix++) = exp(invtwosigma2*(x*x+(y+sep)*(y+sep)))
                - exp(invtwosigma2*(x*x+(y-sep)*(y-sep)));
      }
    }
 
  QMALLOC(gdx, float, npix);
  gdxt = gdx;
  fpix = profit->objpix;
  invsig = npix/profit->sigma;
  for (i=npix; i--; fpix++)
    {
    val = *fpix > -1e29? *fpix*invsig : 0.0;
    *(gdxt++) = (val>0.0? log(1.0+val) : -log(1.0-val));
    }
  gdy = NULL;                   /* to avoid gcc -Wall warnings */
  QMEMCPY(gdx, gdy, float, npix);
  fdx = fft_rtf(dx, profit->objnaxisn);
  fft_conv(gdx, fdx, profit->objnaxisn);
  fdy = fft_rtf(dy, profit->objnaxisn);
  fft_conv(gdy, fdy, profit->objnaxisn);
 
/* Compute estimator */
  invtwosigma2 = -1.18*1.18 / (2.0*profit->guessradius*profit->guessradius);
  xstart = -hw - obj->mx + (int)(obj->mx+0.49999);
  y = -hh -  obj->my + (int)(obj->my+0.49999);;
  spirindex = 0.0;
  gdxt = gdx;
  gdyt = gdy;
  for (j=profit->objnaxisn[1]; j--; y+=1.0)
    {
    x = xstart;
    for (i=profit->objnaxisn[0]; i--; x+=1.0)
      {
      gx = *(gdxt++);
      gy = *(gdyt++);
      if ((r2=x*x+y*y)>0.0)
        spirindex += (x*y*(gx*gx-gy*gy)+gx*gy*(y*y-x*x))/r2
                        * exp(invtwosigma2*r2);
      }
    }
 
  free(dx);
  free(dy);
  free(fdx);
  free(fdy);
  free(gdx);
  free(gdy);
 
  return spirindex;
  }
 
 
/****** profit_moments ****************************************************
PROTO   void profit_moments(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute the 2nd order moments from the unconvolved object model.
INPUT   Profile-fitting structure,
        Pointer to obj2 structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 22/04/2011
 ***/
void     profit_moments(profitstruct *profit, obj2struct *obj2)
  {
   profstruct   *prof;
   double       dpdmx2[6], cov[4],
                *jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
                *dmx2,*dmy2,*dmxy,
                m0,invm0, mx2,my2,mxy, den,invden,
                temp, temp2,invtemp2,invstemp2,
                pmx2,theta, flux, dval;
   float         *covart;
   int          findex[MODEL_NMAX],
                i,j,p, nparam;
 
/*  hw = (float)(profit->modnaxisn[0]/2);*/
/*  hh = (float)(profit->modnaxisn[1]/2);*/
/*  r2max = hw<hh? hw*hw : hh*hh;*/
/*  xstart = -hw;*/
/*  y = -hh;*/
/*  pix = profit->modpix;*/
/*  mx2 = my2 = mxy = mx = my = sum = 0.0;*/
/*  for (iy=profit->modnaxisn[1]; iy--; y+=1.0)*/
/*    {*/
/*    x = xstart;*/
/*    for (ix=profit->modnaxisn[0]; ix--; x+=1.0)*/
/*      if (y*y+x*x <= r2max)*/
/*        {*/
/*        val = *(pix++);*/
/*        sum += val;*/
/*        mx  += val*x;*/
/*        my  += val*y;*/
/*        mx2 += val*x*x;*/
/*        mxy += val*x*y;*/
/*        my2 += val*y*y;*/
/*        }*/
/*      else*/
/*        pix++;*/
/*    }*/
 
/*  if (sum <= 1.0/BIG)*/
/*    sum = 1.0;*/
/*  mx /= sum;*/
/*  my /= sum;*/
/*  obj2->prof_mx2 = mx2 = mx2/sum - mx*mx;*/
/*  obj2->prof_my2 = my2 = my2/sum - my*my;*/
/*  obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
 
  nparam = profit->nparam;
  if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
    {
/*-- Set up Jacobian matrices */
    QCALLOC(jac, double, nparam*3);
    QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
    QMALLOC(dcovar, double, nparam*nparam);
    dcovart = dcovar;
    covart = profit->covar;
    for (i=nparam*nparam; i--;)
      *(dcovart++) = (double)(*(covart++));
    dmx2 = jac;
    dmy2 = jac+nparam;
    dmxy = jac+2*nparam;
    }
  else
    jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;
 
  m0 = mx2 = my2 = mxy = 0.0;
  for (p=0; p<profit->nprof; p++)
    {
    prof = profit->prof[p];
    findex[p] = prof_moments(profit, prof, pjac);
    flux = *prof->flux;
    m0 += flux;
    mx2 += prof->mx2*flux;
    my2 += prof->my2*flux;
    mxy += prof->mxy*flux;
    if (jac)
      {
      jact = jac;
      pjact = pjac;
      for (j=nparam*3; j--;)
        *(jact++) += flux * *(pjact++);
      }
    }
  invm0 = 1.0 / m0;
  obj2->prof_mx2 = (mx2 *= invm0);
  obj2->prof_my2 = (my2 *= invm0);
  obj2->prof_mxy = (mxy *= invm0);
/* Complete the flux derivative of moments */
  if (jac)
    {
    for (p=0; p<profit->nprof; p++)
      {
      prof = profit->prof[p];
      dmx2[findex[p]] = prof->mx2 - mx2;
      dmy2[findex[p]] = prof->my2 - my2;
      dmxy[findex[p]] = prof->mxy - mxy;
      }
    jact = jac;
    for (j=nparam*3; j--;)
      *(jact++) *= invm0;
    }
 
/* Handle fully correlated profiles (which cause a singularity...) */
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
    {
    mx2 += 0.0833333;
    my2 += 0.0833333;
    temp2 = mx2*my2-mxy*mxy;
    }
 
/* Use the Jacobians to compute the moment covariance matrix */
  if (jac)
    propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
                                                pjac);  /* We re-use pjac */
 
  if (FLAG(obj2.prof_pol1))
    {
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
    if (mx2+my2 > 1.0/BIG)
      {
      obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
      obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
      if (FLAG(obj2.prof_pol1err))
        {
/*------ Compute the Jacobian of polarisation */
        invden = 1.0/(mx2+my2);
        dpdmx2[0] =  2.0*my2*invden*invden;
        dpdmx2[1] = -2.0*mx2*invden*invden;
        dpdmx2[2] =  0.0;
        dpdmx2[3] = -2.0*mxy*invden*invden;
        dpdmx2[4] = -2.0*mxy*invden*invden;
        dpdmx2[5] =  2.0*invden;
 
/*------ Use the Jacobian to compute the polarisation covariance matrix */
        propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
                                                pjac);  /* We re-use pjac */
        obj2->prof_pol1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
        obj2->prof_pol2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
        obj2->prof_pol12corr = (dval=cov[0]*cov[3]) > 0.0?
                                        (float)(cov[1]/sqrt(dval)) : 0.0;
        }
      }
    else
      obj2->prof_pol1 = obj2->prof_pol2
        = obj2->prof_pol1err = obj2->prof_pol2err = obj2->prof_pol12corr = 0.0;
    }
 
  if (FLAG(obj2.prof_e1))
    {
/*--- "Ellipticity", i.e. module = (a-b)/(a+b) */
    if (mx2+my2 > 1.0/BIG)
      {
      den = (temp2>=0.0) ? mx2+my2+2.0*sqrt(temp2) : mx2+my2;
      invden = 1.0/den;
      obj2->prof_e1 = (float)(invden * (mx2 - my2));
      obj2->prof_e2 = (float)(2.0 * invden * mxy);
      if (FLAG(obj2.prof_e1err))
        {
/*------ Compute the Jacobian of ellipticity */
        invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0;
        dpdmx2[0] = ( den - (1.0+my2*invstemp2)*(mx2-my2))*invden*invden;
        dpdmx2[1] = (-den - (1.0+mx2*invstemp2)*(mx2-my2))*invden*invden;
        dpdmx2[2] = 2.0*mxy*invstemp2*(mx2-my2)*invden*invden;
        dpdmx2[3] = -2.0*mxy*(1.0+my2*invstemp2)*invden*invden;
        dpdmx2[4] = -2.0*mxy*(1.0+mx2*invstemp2)*invden*invden;
        dpdmx2[5] =  (2.0*den+4.0*mxy*mxy*invstemp2)*invden*invden;
 
/*------ Use the Jacobian to compute the ellipticity covariance matrix */
        propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
                                        pjac);  /* We re-use pjac */
        obj2->prof_e1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
        obj2->prof_e2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
        obj2->prof_e12corr = (dval=cov[0]*cov[3]) > 0.0?
                                        (float)(cov[1]/sqrt(dval)) : 0.0;
        }
      }
    else
      obj2->prof_e1 = obj2->prof_e2
        = obj2->prof_e1err = obj2->prof_e2err = obj2->prof_e12corr = 0.0;
    }
 
  if (FLAG(obj2.prof_cxx))
    {
    invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
    obj2->prof_cxx = (float)(my2*invtemp2);
    obj2->prof_cyy = (float)(mx2*invtemp2);
    obj2->prof_cxy = (float)(-2*mxy*invtemp2);
    }
 
  if (FLAG(obj2.prof_a))
    {
    if ((fabs(temp=mx2-my2)) > 0.0)
      theta = atan2(2.0 * mxy,temp) / 2.0;
    else
      theta = PI/4.0;
 
    temp = sqrt(0.25*temp*temp+mxy*mxy);
    pmx2 = 0.5*(mx2+my2);
    obj2->prof_a = (float)sqrt(pmx2 + temp);
    obj2->prof_b = (float)sqrt(pmx2 - temp);
    obj2->prof_theta = theta*180.0/PI;
    }
 
/* Free memory used by Jacobians */
  free(jac);
  free(pjac);
  free(dcovar);
 
  return;
  }
 
 
/****** profit_convmoments ****************************************************
PROTO   void profit_convmoments(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute the 2nd order moments of the convolved object model.
INPUT   Profile-fitting structure,
        Pointer to obj2 structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 12/04/2011
 ***/
void     profit_convmoments(profitstruct *profit, obj2struct *obj2)
  {
   double       hw,hh, r2max, x,xstart,y, mx2,my2,mxy,mx,my,sum, dval,
                temp,temp2,invtemp2, pmx2, theta;
   PIXTYPE      *pix;
   int          ix,iy, w,h;
 
  w = profit->modnaxisn[0];
  h = profit->modnaxisn[1];
  hw = (double)(w/2);
  hh = (double)(h/2);
 
  r2max = hw<hh? hw*hw : hh*hh;
  xstart = -hw;
  y = -hh;
  pix = profit->cmodpix;
  mx2 = my2 = mxy = mx = my = sum = 0.0;
  for (iy=h; iy--; y+=1.0)
    {
    x = xstart;
    for (ix=w; ix--; x+=1.0)
      if (y*y+x*x <= r2max)
        {
        dval = *(pix++);
        sum += dval;
        mx  += dval*x;
        my  += dval*y;
        mx2 += dval*x*x;
        mxy += dval*x*y;
        my2 += dval*y*y;
        }
      else
        pix++;
    }
 
  if (sum <= 1.0/BIG)
    sum = 1.0;
  mx /= sum;
  my /= sum;
  obj2->prof_convmx2 = (mx2 = mx2/sum - mx*mx)*profit->pixstep*profit->pixstep;
  obj2->prof_convmy2 = (my2 = my2/sum - my*my)*profit->pixstep*profit->pixstep;
  obj2->prof_convmxy = (mxy = mxy/sum - mx*my)*profit->pixstep*profit->pixstep;
 
/* Handle fully correlated profiles (which cause a singularity...) */
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
    {
    mx2 += 0.0833333;
    my2 += 0.0833333;
    temp2 = mx2*my2-mxy*mxy;
    }
 
  temp2 *= profit->pixstep*profit->pixstep;
 
  if (FLAG(obj2.prof_convcxx))
    {
    invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
    obj2->prof_convcxx = (float)(my2*invtemp2);
    obj2->prof_convcyy = (float)(mx2*invtemp2);
    obj2->prof_convcxy = (float)(-2*mxy*invtemp2);
    }
 
  if (1 /*FLAG(obj2.prof_conva)*/)
    {
    if ((fabs(temp=mx2-my2)) > 0.0)
      theta = atan2(2.0 * mxy,temp) / 2.0;
    else
      theta = PI/4.0;
 
    temp = sqrt(0.25*temp*temp+mxy*mxy);
    pmx2 = 0.5*(mx2+my2);
    obj2->prof_conva = (float)sqrt(pmx2 + temp)*profit->pixstep;
    obj2->prof_convb = (float)sqrt(pmx2 - temp)*profit->pixstep;
    obj2->prof_convtheta = theta/DEG;
    }
 
  return;
  }
 
 
/****** profit_surface ****************************************************
PROTO   void profit_surface(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute surface brightnesses from the unconvolved object model.
INPUT   Pointer to the profile-fitting structure,
        Pointer to obj2 structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 24/01/2011
 ***/
void     profit_surface(profitstruct *profit, obj2struct *obj2)
  {
   profitstruct hdprofit;
   double       dsum,dhsum,dsumoff, dhval, frac, seff;
   float        *spix, *spixt,
                val,vmax,
                scalefac, imsizefac, flux, lost, sum, lostfluxfrac;
   int          i,p, imax, npix, neff;
 
/* Allocate "high-definition" raster only to make measurements */
  hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES;
  npix = hdprofit.nmodpix = hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1];

Line No.	Rev	Author	Line
1	233	bertin	`/*`
2			`* profit.c`
3	42	bertin	`*`
4	233	bertin	`* Fit a range of galaxy models to an image.`
5	42	bertin	`*`
6	233	bertin	`*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%`
7	42	bertin	`*`
8	233	bertin	`* This file part of: SExtractor`
9	42	bertin	`*`
10	246	bertin	`* Copyright: (C) 2006-2011 Emmanuel Bertin -- IAP/CNRS/UPMC`
11	42	bertin	`*`
12	233	bertin	`* License: GNU General Public License`
13			`*`
14			`* SExtractor is free software: you can redistribute it and/or modify`
15			`* it under the terms of the GNU General Public License as published by`
16			`* the Free Software Foundation, either version 3 of the License, or`
17			`* (at your option) any later version.`
18			`* SExtractor is distributed in the hope that it will be useful,`
19			`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
20			`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
21			`* GNU General Public License for more details.`
22			`* You should have received a copy of the GNU General Public License`
23			`* along with SExtractor. If not, see <http://www.gnu.org/licenses/>.`
24			`*`
25	254	bertin	`* Last modified: 25/07/2011`
26	233	bertin	`*`
27			`%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/`
28	42	bertin
29			`#ifdef HAVE_CONFIG_H`
30			`#include "config.h"`
31			`#endif`
32
33	206	bertin	`#ifndef HAVE_MATHIMF_H`
34			`#define _GNU_SOURCE`
35			`#endif`
36
37	201	bertin	`#include <math.h>`
38	42	bertin	`#include <stdio.h>`
39			`#include <stdlib.h>`
40			`#include <string.h>`
41
42			`#include "define.h"`
43			`#include "globals.h"`
44			`#include "prefs.h"`
45			`#include "fits/fitscat.h"`
46	228	bertin	`#include "levmar/levmar.h"`
47	50	bertin	`#include "fft.h"`
48	123	bertin	`#include "fitswcs.h"`
49	42	bertin	`#include "check.h"`
50	229	bertin	`#include "image.h"`
51	141	bertin	`#include "pattern.h"`
52	50	bertin	`#include "psf.h"`
53	42	bertin	`#include "profit.h"`
54
55	207	bertin	`static double prof_gammainc(double x, double a),`
56			`prof_gamma(double x);`
57	201	bertin	`static float prof_interpolate(profstruct prof, float posin);`
58			`static float interpolate_pix(float posin, float pix, int *naxisn,`
59	52	bertin	`interpenum interptype);`
60	50	bertin
61	201	bertin	`static void make_kernel(float pos, float *kernel, interpenum interptype);`
62	50	bertin
63	42	bertin	`/------------------------------- variables ---------------------------------/`
64
65	209	bertin	`const int interp_kernwidth[5]={1,2,4,6,8};`
66
67	233	bertin	`const int flux_flag[PARAM_NPARAM] = {0,`
68			`1,0,0,`
69	209	bertin	`1,0,0,0,0,`
70			`1,0,0,0,`
71			`1,0,0,0,0,0,0,0,`
72			`1,0,0,`
73			`1,0,0,`
74			`1,0,0`
75			`};`
76
77	235	bertin	`/* "Local" global variables for debugging purposes */`
78	145	bertin	`int theniter, the_gal;`
79	48	bertin	`static picstruct the_field, the_wfield;`
80	244	bertin	`profitstruct theprofit, thepprofit, *theqprofit;`
81	48	bertin
82	46	bertin	`/**** profit_init *********************************************************`
83	244	bertin	`PROTO profitstruct profit_init(psfstruct *psf, unsigned int modeltype)`
84	46	bertin	`PURPOSE Allocate and initialize a new profile-fitting structure.`
85	244	bertin	`INPUT Pointer to PSF structure,`
86			`Model type.`
87	46	bertin	`OUTPUT A pointer to an allocated profit structure.`
88	117	bertin	`NOTES -.`
89	42	bertin	`AUTHOR E. Bertin (IAP)`
90	244	bertin	`VERSION 22/04/2011`
91	42	bertin	`***/`
92	244	bertin	`profitstruct profit_init(psfstruct psf, unsigned int modeltype)`
93	42	bertin	`{`
94	46	bertin	`profitstruct *profit;`
95	244	bertin	`int t, nmodels;`
96	42	bertin
97	49	bertin	`QCALLOC(profit, profitstruct, 1);`
98	58	bertin	`profit->psf = psf;`
99	244	bertin	`QMALLOC(profit->prof, profstruct *, MODEL_NMAX);`
100			`nmodels = 0;`
101			`for (t=1; t<(1<<MODEL_NMAX); t<<=1)`
102			`if (modeltype&t)`
103			`profit->prof[nmodels++] = prof_init(profit, t);`
104	221	bertin	`/* Allocate memory for the complete model */`
105			`QMALLOC16(profit->modpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);`
106			`QMALLOC16(profit->modpix2, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);`
107			`QMALLOC16(profit->cmodpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);`
108			`QMALLOC16(profit->psfpix, float, PROFIT_MAXMODSIZE*PROFIT_MAXMODSIZE);`
109			`QMALLOC16(profit->objpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);`
110			`QMALLOC16(profit->objweight, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);`
111			`QMALLOC16(profit->lmodpix, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);`
112			`QMALLOC16(profit->lmodpix2, PIXTYPE, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);`
113			`QMALLOC16(profit->resi, float, PROFIT_MAXOBJSIZE*PROFIT_MAXOBJSIZE);`
114	211	bertin	`QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);`
115	244	bertin	`profit->nprof = nmodels;`
116	207	bertin	`profit->fluxfac = 1.0; /* Default */`
117	50	bertin
118	46	bertin	`return profit;`
119	42	bertin	`}`
120
121
122	46	bertin	`/**** profit_end **********************************************************`
123	42	bertin	`PROTO void prof_end(profstruct *prof)`
124	46	bertin	`PURPOSE End (deallocate) a profile-fitting structure.`
125	42	bertin	`INPUT Prof structure.`
126			`OUTPUT -.`
127	117	bertin	`NOTES -.`
128	42	bertin	`AUTHOR E. Bertin (IAP)`
129	221	bertin	`VERSION 06/04/2010`
130	42	bertin	`***/`
131	46	bertin	`void profit_end(profitstruct *profit)`
132	42	bertin	`{`
133	49	bertin	`int p;`
134
135			`for (p=0; p<profit->nprof; p++)`
136			`prof_end(profit->prof[p]);`
137	221	bertin	`free(profit->modpix);`
138			`free(profit->modpix2);`
139			`free(profit->cmodpix);`
140			`free(profit->psfpix);`
141			`free(profit->lmodpix);`
142			`free(profit->lmodpix2);`
143			`free(profit->objpix);`
144			`free(profit->objweight);`
145			`free(profit->resi);`
146	49	bertin	`free(profit->prof);`
147	123	bertin	`free(profit->covar);`
148	60	bertin	`free(profit->psfdft);`
149	46	bertin	`free(profit);`
150	42	bertin
151			`return;`
152			`}`
153
154
155	65	bertin	`/**** profit_fit **********************************************************`
156			`PROTO void profit_fit(profitstruct profit, picstruct field,`
157			`picstruct wfield, objstruct obj, obj2struct *obj2)`
158	46	bertin	`PURPOSE Fit profile(s) convolved with the PSF to a detected object.`
159	42	bertin	`INPUT Array of profile structures,`
160			`Number of profiles,`
161	52	bertin	`Pointer to the profile-fitting structure,`
162	42	bertin	`Pointer to the field,`
163			`Pointer to the field weight,`
164			`Pointer to the obj.`
165	46	bertin	`OUTPUT Pointer to an allocated fit structure (containing details about the`
166			`fit).`
167	48	bertin	`NOTES It is a modified version of the lm_minimize() of lmfit.`
168	46	bertin	`AUTHOR E. Bertin (IAP)`
169	246	bertin	`VERSION 31/05/2011`
170	46	bertin	`***/`
171	65	bertin	`void profit_fit(profitstruct *profit,`
172	52	bertin	`picstruct field, picstruct wfield,`
173	50	bertin	`objstruct obj, obj2struct obj2)`
174	46	bertin	`{`
175	244	bertin	`profitstruct pprofit, qprofit;`
176	221	bertin	`patternstruct *pattern;`
177	52	bertin	`psfstruct *psf;`
178	50	bertin	`checkstruct *check;`
179	244	bertin	`double emx2,emy2,emxy, a , cp,sp, cn, bn, n,`
180			`sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0;`
181			`PIXTYPE valp,valq,sig2;`
182	221	bertin	`float param0[PARAM_NPARAM], param1[PARAM_NPARAM],`
183			`param[PARAM_NPARAM],`
184			`**list,`
185	209	bertin	`*cov,`
186	221	bertin	`psf_fwhm, dchi2, err, aspect, chi2;`
187	209	bertin	`int *index,`
188			`i,j,p, nparam, nparam2, ncomp, nprof;`
189	46	bertin
190	128	bertin	`nparam = profit->nparam;`
191	208	bertin	`nparam2 = nparam*nparam;`
192	194	bertin	`nprof = profit->nprof;`
193	52	bertin	`if (profit->psfdft)`
194			`{`
195			`QFREE(profit->psfdft);`
196			`}`
197
198			`psf = profit->psf;`
199	118	bertin	`profit->pixstep = psf->pixstep;`
200	211	bertin	`obj2->prof_flag = 0;`
201	52	bertin
202			`/* Create pixmaps at image resolution */`
203	211	bertin	`profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */`
204			`profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */`
205	51	bertin	`psf_fwhm = psf->masksize[0]*psf->pixstep;`
206	233	bertin	`profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)`
207	75	bertin	`1.2)/2)2 + 1;`
208	233	bertin	`profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)`
209	75	bertin	`1.2)/2)2 + 1;`
210	117	bertin	`if (profit->objnaxisn[1]<profit->objnaxisn[0])`
211			`profit->objnaxisn[1] = profit->objnaxisn[0];`
212			`else`
213			`profit->objnaxisn[0] = profit->objnaxisn[1];`
214	211	bertin	`if (profit->objnaxisn[0]>PROFIT_MAXOBJSIZE)`
215			`{`
216	212	bertin	`profit->subsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);`
217			`profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->subsamp);`
218	211	bertin	`obj2->prof_flag \|= PROFLAG_OBJSUB;`
219			`}`
220			`else`
221	212	bertin	`profit->subsamp = 1.0;`
222	221	bertin	`profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1];`
223	117	bertin
224	78	bertin	`/* Use (dirty) global variables to interface with lmfit */`
225			`the_field = field;`
226			`the_wfield = wfield;`
227			`theprofit = profit;`
228	146	bertin	`profit->obj = obj;`
229			`profit->obj2 = obj2;`
230	78	bertin
231	146	bertin	`profit->nresi = profit_copyobjpix(profit, field, wfield);`
232	211	bertin	`/* Check if the number of constraints exceeds the number of free parameters */`
233	128	bertin	`if (profit->nresi < nparam)`
234	49	bertin	`{`
235	176	bertin	`if (FLAG(obj2.prof_vector))`
236			`for (p=0; p<nparam; p++)`
237			`obj2->prof_vector[p] = 0.0;`
238	208	bertin	`if (FLAG(obj2.prof_errvector))`
239			`for (p=0; p<nparam; p++)`
240			`obj2->prof_errvector[p] = 0.0;`
241			`if (FLAG(obj2.prof_errmatrix))`
242			`for (p=0; p<nparam2; p++)`
243			`obj2->prof_errmatrix[p] = 0.0;`
244	50	bertin	`obj2->prof_niter = 0;`
245	211	bertin	`obj2->prof_flag \|= PROFLAG_NOTCONST;`
246	50	bertin	`return;`
247	49	bertin	`}`
248	48	bertin
249	52	bertin	`/* Create pixmap at PSF resolution */`
250	149	bertin	`profit->modnaxisn[0] =`
251	212	bertin	`((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep`
252			`+0.4999)/2+1)*2;`
253	149	bertin	`profit->modnaxisn[1] =`
254	212	bertin	`((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep`
255			`+0.4999)/2+1)*2;`
256	56	bertin	`if (profit->modnaxisn[1] < profit->modnaxisn[0])`
257			`profit->modnaxisn[1] = profit->modnaxisn[0];`
258			`else`
259			`profit->modnaxisn[0] = profit->modnaxisn[1];`
260	211	bertin	`if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)`
261			`{`
262			`profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;`
263			`profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;`
264			`obj2->prof_flag \|= PROFLAG_MODSUB;`
265			`}`
266	221	bertin	`profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];`
267	118	bertin
268			`/* Compute the local PSF */`
269			`profit_psf(profit);`
270
271	50	bertin	`/* Set initial guesses and boundaries */`
272	78	bertin	`profit->sigma = obj->sigbkg;`
273	246	bertin	`if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)`
274			`profit->guessradius = obj2->hl_radius;`
275			`if ((profit->guessflux = obj2->flux_auto) <= 0.0)`
276			`profit->guessflux = 0.0;`
277			`if ((profit->guessfluxmax = 10.0*obj2->fluxerr_auto) <= profit->guessflux)`
278			`profit->guessfluxmax = profit->guessflux;`
279			`if (profit->guessfluxmax <= 0.0)`
280			`profit->guessfluxmax = 1.0;`
281	115	bertin
282	146	bertin	`profit_resetparams(profit);`
283	50	bertin
284	123	bertin	`/* Actual minimisation */`
285	209	bertin	`fft_reset();`
286	211	bertin	`the_gal++;`
287	221	bertin
288	235	bertin	`/*`
289			`char str[1024];`
290			`sprintf(str, "obj_%04d.fits", the_gal);`
291			`catstruct *bcat;`
292			`float bpix, opix,lmodpix,objpix;`
293			`bcat=read_cat("base.fits");`
294			`QMALLOC(bpix, float, field->npix);`
295			`QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);`
296			`read_body(bcat->tab, bpix, field->npix);`
297			`free_cat(&bcat,1);`
298			`bcat=read_cat(str);`
299			`QMALLOC(opix, float, profit->nobjpix);`
300			`QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);`
301			`read_body(bcat->tab, opix, profit->nobjpix);`
302			`free_cat(&bcat,1);`
303			`addfrombig(bpix, field->width, field->height,`
304			`profit->objpix, profit->objnaxisn[0],profit->objnaxisn[1],`
305			`profit->ix,profit->iy, -1.0);`
306			`objpix = profit->objpix;`
307			`lmodpix = opix;`
308			`for (i=profit->nobjpix; i--;)`
309			`(objpix++) += (lmodpix++);`
310			`free(bpix);`
311			`free(opix);`
312			`*/`
313	79	bertin	`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
314	221	bertin	`/*`
315	235	bertin	`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
316			`check=initcheck(str, CHECK_OTHER,1);`
317			`check->width = profit->objnaxisn[0];`
318			`check->height = profit->objnaxisn[1];`
319			`reinitcheck(field,check);`
320			`memcpy(check->pix, profit->lmodpix, profit->nobjpix*sizeof(float));`
321			`reendcheck(field,check);`
322			`endcheck(check);`
323
324	221	bertin	`chi2 = profit->chi2;`
325			`for (p=0; p<nparam; p++)`
326			`param1[p] = profit->paraminit[p];`
327			`profit_resetparams(profit);`
328			`for (p=0; p<nparam; p++)`
329			`profit->paraminit[p] = param1[p] + (profit->paraminit[p]<param1[p]?1.0:-1.0)`
330			`* sqrt(profit->covar[p*(nparam+1)]);`
331	200	bertin
332	221	bertin	`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
333			`if (chi2<profit->chi2)`
334			`for (p=0; p<nparam; p++)`
335			`profit->paraminit[p] = param1[p];`
336
337			`list = profit->paramlist;`
338			`index = profit->paramindex;`
339			`for (i=0; i<PARAM_NPARAM; i++)`
340			`if (list[i] && i!= PARAM_SPHEROID_ASPECT && i!=PARAM_SPHEROID_POSANG)`
341			`profit->freeparam_flag[index[i]] = 0;`
342			`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
343			`*/`
344	246	bertin	`if (profit->nlimmin)`
345			`obj2->prof_flag \|= PROFLAG_MINLIM;`
346			`if (profit->nlimmax)`
347			`obj2->prof_flag \|= PROFLAG_MAXLIM;`
348
349	128	bertin	`for (p=0; p<nparam; p++)`
350			`profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);`
351	126	bertin
352	50	bertin	`/* CHECK-Images */`
353	221	bertin	`if ((check = prefs.check[CHECK_PROFILES]))`
354			`{`
355			`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
356			`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
357			`profit->ix,profit->iy, 1.0);`
358			`}`
359	235	bertin
360	50	bertin	`if ((check = prefs.check[CHECK_SUBPROFILES]))`
361	51	bertin	`{`
362	221	bertin	`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
363	50	bertin	`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
364	145	bertin	`profit->ix,profit->iy, -1.0);`
365	51	bertin	`}`
366	221	bertin	`if ((check = prefs.check[CHECK_SPHEROIDS]))`
367	51	bertin	`{`
368	221	bertin	`/-- Set to 0 flux components that do not belong to spheroids /`
369			`for (p=0; p<profit->nparam; p++)`
370			`param[p] = profit->paraminit[p];`
371			`list = profit->paramlist;`
372			`index = profit->paramindex;`
373			`for (i=0; i<PARAM_NPARAM; i++)`
374			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
375			`param[index[i]] = 0.0;`
376			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
377	50	bertin	`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
378	145	bertin	`profit->ix,profit->iy, 1.0);`
379	51	bertin	`}`
380	221	bertin	`if ((check = prefs.check[CHECK_SUBSPHEROIDS]))`
381			`{`
382			`/-- Set to 0 flux components that do not belong to spheroids /`
383			`for (p=0; p<profit->nparam; p++)`
384			`param[p] = profit->paraminit[p];`
385			`list = profit->paramlist;`
386			`index = profit->paramindex;`
387			`for (i=0; i<PARAM_NPARAM; i++)`
388			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
389			`param[index[i]] = 0.0;`
390			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
391			`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
392			`profit->ix,profit->iy, -1.0);`
393			`}`
394			`if ((check = prefs.check[CHECK_DISKS]))`
395			`{`
396			`/-- Set to 0 flux components that do not belong to disks /`
397			`for (p=0; p<profit->nparam; p++)`
398			`param[p] = profit->paraminit[p];`
399			`list = profit->paramlist;`
400			`index = profit->paramindex;`
401			`for (i=0; i<PARAM_NPARAM; i++)`
402			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
403			`param[index[i]] = 0.0;`
404			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
405			`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
406			`profit->ix,profit->iy, 1.0);`
407			`}`
408			`if ((check = prefs.check[CHECK_SUBDISKS]))`
409			`{`
410			`/-- Set to 0 flux components that do not belong to disks /`
411			`for (p=0; p<profit->nparam; p++)`
412			`param[p] = profit->paraminit[p];`
413			`list = profit->paramlist;`
414			`index = profit->paramindex;`
415			`for (i=0; i<PARAM_NPARAM; i++)`
416			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
417			`param[index[i]] = 0.0;`
418			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
419			`addcheck(check, profit->lmodpix, profit->objnaxisn[0],profit->objnaxisn[1],`
420			`profit->ix,profit->iy, -1.0);`
421			`}`
422	235	bertin
423	221	bertin	`/* Compute compressed residuals */`
424			`profit_residuals(profit,field,wfield, 10.0, profit->paraminit,profit->resi);`
425
426	73	bertin	`/* Fill measurement parameters */`
427	69	bertin	`if (FLAG(obj2.prof_vector))`
428			`{`
429	128	bertin	`for (p=0; p<nparam; p++)`
430	221	bertin	`obj2->prof_vector[p]= profit->paraminit[p];`
431	69	bertin	`}`
432	127	bertin	`if (FLAG(obj2.prof_errvector))`
433			`{`
434	128	bertin	`for (p=0; p<nparam; p++)`
435			`obj2->prof_errvector[p]= profit->paramerr[p];`
436	127	bertin	`}`
437	208	bertin	`if (FLAG(obj2.prof_errmatrix))`
438			`{`
439			`for (p=0; p<nparam2; p++)`
440			`obj2->prof_errmatrix[p]= profit->covar[p];`
441			`}`
442	89	bertin
443	79	bertin	`obj2->prof_niter = profit->niter;`
444	73	bertin	`obj2->flux_prof = profit->flux;`
445	194	bertin	`if (FLAG(obj2.fluxerr_prof))`
446			`{`
447			`err = 0.0;`
448	209	bertin	`cov = profit->covar;`
449			`index = profit->paramindex;`
450			`list = profit->paramlist;`
451			`for (i=0; i<PARAM_NPARAM; i++)`
452			`if (flux_flag[i] && list[i])`
453			`{`
454			`cov = profit->covar + nparam*index[i];`
455			`for (j=0; j<PARAM_NPARAM; j++)`
456			`if (flux_flag[j] && list[j])`
457			`err += cov[index[j]];`
458			`}`
459	194	bertin	`obj2->fluxerr_prof = err>0.0? sqrt(err): 0.0;`
460			`}`
461
462	180	bertin	`obj2->prof_chi2 = (profit->nresi > profit->nparam)?`
463			`profit->chi2 / (profit->nresi - profit->nparam) : 0.0;`
464
465	244	bertin	`/* Position */`
466	131	bertin	`if (FLAG(obj2.x_prof))`
467			`{`
468			`i = profit->paramindex[PARAM_X];`
469			`j = profit->paramindex[PARAM_Y];`
470	145	bertin	`/-- Model coordinates follow the FITS convention (first pixel at 1,1) /`
471	179	bertin	`if (profit->paramlist[PARAM_X])`
472			`{`
473	217	bertin	`obj2->x_prof = (double)profit->ix + *profit->paramlist[PARAM_X] + 1.0;`
474	179	bertin	`obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];`
475			`}`
476			`else`
477			`emx2 = 0.0;`
478			`if (profit->paramlist[PARAM_Y])`
479			`{`
480	217	bertin	`obj2->y_prof = (double)profit->iy + *profit->paramlist[PARAM_Y] + 1.0;`
481	179	bertin	`obj2->poserrmy2_prof = emy2 = profit->covar[j*(nparam+1)];`
482			`}`
483			`else`
484			`emy2 = 0.0;`
485			`if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])`
486			`obj2->poserrmxy_prof = emxy = profit->covar[i+j*nparam];`
487			`else`
488			`emxy = 0.0;`
489	131	bertin
490	135	bertin	`/-- Error ellipse parameters /`
491			`if (FLAG(obj2.poserra_prof))`
492			`{`
493			`double pmx2,pmy2,temp,theta;`
494	131	bertin
495	135	bertin	`if (fabs(temp=emx2-emy2) > 0.0)`
496			`theta = atan2(2.0 * emxy,temp) / 2.0;`
497			`else`
498			`theta = PI/4.0;`
499	131	bertin
500	135	bertin	`temp = sqrt(0.25temptemp+ emxy*emxy);`