Astromatic - public software.sextractor - Blame - Rev 309

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

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-2013 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:          05/04/2013
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
#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,*thedprofit, *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->presi, float, profit->nparam);
  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 12/07/2012
 ***/
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->presi);
  free(profit->prof);
  free(profit->covar);
  QFFTWF_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 05/04/2013
 ***/
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, err;
    PIXTYPE             valp,valq,sig2;
    float               param0[PARAM_NPARAM], param1[PARAM_NPARAM],
                        param[PARAM_NPARAM],
                        **list,
                        *cov,
                        psf_fwhm, dchi2, aspect, chi2;
    int                 *index,
                        i,j,p, nparam, nparam2, ncomp, nprof;
 
  nparam = profit->nparam;
  nparam2 = nparam*nparam;
  nprof = profit->nprof;
 
  if (profit->psfdft)
    QFFTWF_FREE(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];
 
/* Create pixmap at model 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];
 
/* Use (dirty) global variables to interface with lmfit */
  the_field = field;
  the_wfield = wfield;
  theprofit = profit;
  profit->obj = obj;
  profit->obj2 = obj2;
 
/* Compute the local PSF */
  profit_psf(profit);
 
  profit->nresi = profit_copyobjpix(profit, field, wfield);
  profit->npresi = 0;
/* 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;
    }
 
/* Set initial guesses and boundaries */
  profit->guesssigbkg = profit->sigma = obj->sigbkg;
  profit->guessdx = obj->mx - (int)(obj->mx+0.49999);
  profit->guessdy = obj->my - (int)(obj->my+0.49999);
  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;
  if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)
    profit->guessradius = obj2->hl_radius;
  profit->guessaspect = obj->b/obj->a;
  profit->guessposang = obj->theta;
 
  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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                1.0/(profit->subsamp*profit->subsamp));
    else
      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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                -1.0/(profit->subsamp*profit->subsamp));
    else
      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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                1.0/(profit->subsamp*profit->subsamp));
    else
      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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                -1.0/(profit->subsamp*profit->subsamp));
    else
      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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                1.0/(profit->subsamp*profit->subsamp));
    else
      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);
    if (profit->subsamp>1.0)
      addcheck_resample(check, profit->lmodpix,
                profit->objnaxisn[0],profit->objnaxisn[1],
                profit->ix,profit->iy, 1.0/profit->subsamp,
                -1.0/(profit->subsamp*profit->subsamp));
    else
      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>0.0? pmx2 : 0.0);
      obj2->poserrb_prof = (float)sqrt(pmy2>0.0? pmy2 : 0.0);
      obj2->poserrtheta_prof = (float)(theta/DEG);
      }
 
    if (FLAG(obj2.poserrcxx_prof))
      {
       double   temp;
 
      obj2->poserrcxx_prof = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
      obj2->poserrcyy_prof = (float)(emx2/temp);
      obj2->poserrcxy_prof = (float)(-2*emxy/temp);
      }
    }
 
/* Equivalent noise area */
  if (FLAG(obj2.prof_noisearea))
    obj2->prof_noisearea = profit_noisearea(profit);
 
/* Second order moments and ellipticities */
  if (FLAG(obj2.prof_mx2))
    profit_moments(profit, obj2);
 
/* Second order moments of the convolved model (used by other parameters) */
  if (FLAG(obj2.prof_convmx2))
    profit_convmoments(profit, obj2);
 
/* "Hybrid" magnitudes */
  if (FLAG(obj2.fluxcor_prof))
    {
    profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);
    profit_fluxcor(profit, obj, obj2);
    }
 
/* Do measurements on the rasterised model (surface brightnesses) */
  if (FLAG(obj2.fluxeff_prof))
    profit_surface(profit, obj2);
 
/* Background offset */
  if (FLAG(obj2.prof_offset_flux))
    {
    obj2->prof_offset_flux = *profit->paramlist[PARAM_BACK];
    obj2->prof_offset_fluxerr=profit->paramerr[profit->paramindex[PARAM_BACK]];
    }
 
/* Point source */
  if (FLAG(obj2.prof_dirac_flux))
    {
    obj2->prof_dirac_flux = *profit->paramlist[PARAM_DIRAC_FLUX];
    obj2->prof_dirac_fluxerr =
                profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX]];
    }
 
/* Spheroid */
  if (FLAG(obj2.prof_spheroid_flux))
    {
    if ((aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT]) > 1.0)
      {
      *profit->paramlist[PARAM_SPHEROID_REFF] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]]
                        /= (aspect*aspect);
      *profit->paramlist[PARAM_SPHEROID_ASPECT] = 1.0 / aspect;
      *profit->paramlist[PARAM_SPHEROID_POSANG] += 90.0;
      }
    obj2->prof_spheroid_flux = *profit->paramlist[PARAM_SPHEROID_FLUX];
    obj2->prof_spheroid_fluxerr =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX]];
    obj2->prof_spheroid_reff = *profit->paramlist[PARAM_SPHEROID_REFF];
    obj2->prof_spheroid_refferr =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]];
    obj2->prof_spheroid_aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT];
    obj2->prof_spheroid_aspecterr =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]];
    obj2->prof_spheroid_theta =
                        fmod_m90_p90(*profit->paramlist[PARAM_SPHEROID_POSANG]);
    obj2->prof_spheroid_thetaerr =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_POSANG]];
    if (FLAG(obj2.prof_spheroid_sersicn))
      {
      obj2->prof_spheroid_sersicn = *profit->paramlist[PARAM_SPHEROID_SERSICN];
      obj2->prof_spheroid_sersicnerr =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_SERSICN]];
      }
    else
      obj2->prof_spheroid_sersicn = 4.0;
    if (FLAG(obj2.prof_spheroid_peak))
      {
      n = obj2->prof_spheroid_sersicn;
      bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
                + 131.0/(1148175*n*n*n);        /* Ciotti & Bertin 1999 */
      cn = n * prof_gamma(2.0*n) * pow(bn, -2.0*n);
      obj2->prof_spheroid_peak = obj2->prof_spheroid_reff>0.0?
        obj2->prof_spheroid_flux
                / (2.0 * PI * cn
                * obj2->prof_spheroid_reff*obj2->prof_spheroid_reff
                * obj2->prof_spheroid_aspect)
        : 0.0;
      if (FLAG(obj2.prof_spheroid_fluxeff))
        obj2->prof_spheroid_fluxeff = obj2->prof_spheroid_peak * exp(-bn);
      if (FLAG(obj2.prof_spheroid_fluxmean))
        obj2->prof_spheroid_fluxmean = obj2->prof_spheroid_peak * cn;
      }
    }
 
/* Disk */
  if (FLAG(obj2.prof_disk_flux))
    {
    if ((aspect = *profit->paramlist[PARAM_DISK_ASPECT]) > 1.0)
      {
      *profit->paramlist[PARAM_DISK_SCALE] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]] *= aspect;
      profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]]
                        /= (aspect*aspect);
      *profit->paramlist[PARAM_DISK_ASPECT] = 1.0 / aspect;
      *profit->paramlist[PARAM_DISK_POSANG] += 90.0;
      }
    obj2->prof_disk_flux = *profit->paramlist[PARAM_DISK_FLUX];
    obj2->prof_disk_fluxerr =
                profit->paramerr[profit->paramindex[PARAM_DISK_FLUX]];
    obj2->prof_disk_scale = *profit->paramlist[PARAM_DISK_SCALE];
    obj2->prof_disk_scaleerr =
                profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]];
    obj2->prof_disk_aspect = *profit->paramlist[PARAM_DISK_ASPECT];
    obj2->prof_disk_aspecterr =
                profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]];
    obj2->prof_disk_theta = fmod_m90_p90(*profit->paramlist[PARAM_DISK_POSANG]);
    obj2->prof_disk_thetaerr =
                profit->paramerr[profit->paramindex[PARAM_DISK_POSANG]];
    if (FLAG(obj2.prof_disk_inclination))
      {
      obj2->prof_disk_inclination = acos(obj2->prof_disk_aspect) / DEG;
      if (FLAG(obj2.prof_disk_inclinationerr))
        {
        a = sqrt(1.0-obj2->prof_disk_aspect*obj2->prof_disk_aspect);
        obj2->prof_disk_inclinationerr = obj2->prof_disk_aspecterr
                                        /(a>0.1? a : 0.1)/DEG;
        }
      }
 
    if (FLAG(obj2.prof_disk_peak))
      {
      obj2->prof_disk_peak = obj2->prof_disk_scale>0.0?
        obj2->prof_disk_flux
        / (2.0 * PI * obj2->prof_disk_scale*obj2->prof_disk_scale
                * obj2->prof_disk_aspect)
        : 0.0;
      if (FLAG(obj2.prof_disk_fluxeff))
        obj2->prof_disk_fluxeff = obj2->prof_disk_peak * 0.186682; /* e^-(b_n)*/
      if (FLAG(obj2.prof_disk_fluxmean))
        obj2->prof_disk_fluxmean = obj2->prof_disk_peak * 0.355007;/* b_n^(-2)*/
      }
 
/* Disk pattern */
    if (prefs.pattern_flag)
      {
      profit_residuals(profit,field,wfield, PROFIT_DYNPARAM,
                        profit->paraminit,profit->resi);
      pattern = pattern_init(profit, prefs.pattern_type,
                prefs.prof_disk_patternncomp);
      pattern_fit(pattern, profit);
      if (FLAG(obj2.prof_disk_patternspiral))
        obj2->prof_disk_patternspiral = pattern_spiral(pattern);
      if (FLAG(obj2.prof_disk_patternvector))
        {
        ncomp = pattern->size[2];
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternvector[p] = (float)pattern->coeff[p];
        }
      if (FLAG(obj2.prof_disk_patternmodvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternmodvector[p] = (float)pattern->mcoeff[p];
        }
      if (FLAG(obj2.prof_disk_patternargvector))
        {
        ncomp = pattern->ncomp*pattern->nfreq;
        for (p=0; p<ncomp; p++)
          obj2->prof_disk_patternargvector[p] = (float)pattern->acoeff[p];
        }
      pattern_end(pattern);
      }
 
/* Bar */
    if (FLAG(obj2.prof_bar_flux))
      {
      obj2->prof_bar_flux = *profit->paramlist[PARAM_BAR_FLUX];
      obj2->prof_bar_fluxerr =
                profit->paramerr[profit->paramindex[PARAM_BAR_FLUX]];
      obj2->prof_bar_length = *profit->paramlist[PARAM_ARMS_START]
                                **profit->paramlist[PARAM_DISK_SCALE];
      obj2->prof_bar_lengtherr = *profit->paramlist[PARAM_ARMS_START]
                  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
                + *profit->paramlist[PARAM_DISK_SCALE]
                  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
      obj2->prof_bar_aspect = *profit->paramlist[PARAM_BAR_ASPECT];
      obj2->prof_bar_aspecterr =
                profit->paramerr[profit->paramindex[PARAM_BAR_ASPECT]];
      obj2->prof_bar_posang =
                        fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
      obj2->prof_bar_posangerr =
                profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
      if (FLAG(obj2.prof_bar_theta))
        {
        cp = cos(obj2->prof_bar_posang*DEG);
        sp = sin(obj2->prof_bar_posang*DEG);
        a = obj2->prof_disk_aspect;
        obj2->prof_bar_theta = fmod_m90_p90(atan2(a*sp,cp)/DEG
                                + obj2->prof_disk_theta);
        obj2->prof_bar_thetaerr = obj2->prof_bar_posangerr*a/(cp*cp+a*a*sp*sp);
        }
 
/* Arms */
      if (FLAG(obj2.prof_arms_flux))
        {
        obj2->prof_arms_flux = *profit->paramlist[PARAM_ARMS_FLUX];
        obj2->prof_arms_fluxerr =
                profit->paramerr[profit->paramindex[PARAM_ARMS_FLUX]];
        obj2->prof_arms_pitch =
                fmod_m90_p90(*profit->paramlist[PARAM_ARMS_PITCH]);
        obj2->prof_arms_pitcherr =
                profit->paramerr[profit->paramindex[PARAM_ARMS_PITCH]];
        obj2->prof_arms_start = *profit->paramlist[PARAM_ARMS_START]
                                **profit->paramlist[PARAM_DISK_SCALE];
        obj2->prof_arms_starterr = *profit->paramlist[PARAM_ARMS_START]
                  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
                + *profit->paramlist[PARAM_DISK_SCALE]
                  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
        obj2->prof_arms_quadfrac = *profit->paramlist[PARAM_ARMS_QUADFRAC];
        obj2->prof_arms_quadfracerr =
                profit->paramerr[profit->paramindex[PARAM_ARMS_QUADFRAC]];
        obj2->prof_arms_posang =
                        fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
        obj2->prof_arms_posangerr =
                profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
        }
      }
    }
 
/* Star/galaxy classification */
  if (FLAG(obj2.prof_class_star) || FLAG(obj2.prof_concentration))
    {
    profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, profit->paraminit,
        FLAG(obj2.prof_class_star) ? profit->resi : NULL);
    pprofit = thepprofit;
    nparam = pprofit->nparam;
    if (pprofit->psfdft)
      QFFTWF_FREE(pprofit->psfdft);
    psf = pprofit->psf;
    pprofit->pixstep = profit->pixstep;
    pprofit->guesssigbkg = profit->guesssigbkg;
    pprofit->guessdx = profit->guessdx;
    pprofit->guessdy = profit->guessdy;
    pprofit->guessflux = profit->guessflux;
    pprofit->guessfluxmax = profit->guessfluxmax;
    pprofit->guessradius = profit->guessradius;
    pprofit->guessaspect = profit->guessaspect;
    pprofit->guessposang = profit->guessposang;
    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,
                        FLAG(obj2.prof_class_star)? pprofit->resi : NULL);
    qprofit = theqprofit;
    nparam = qprofit->nparam;
    if (qprofit->psfdft)
      QFFTWF_FREE(qprofit->psfdft);
    qprofit->pixstep = profit->pixstep;
    qprofit->guesssigbkg = profit->guesssigbkg;
    qprofit->guessdx = profit->guessdx;
    qprofit->guessdy = profit->guessdy;
    qprofit->guessflux = profit->guessflux;
    qprofit->guessfluxmax = profit->guessfluxmax;
    qprofit->guessradius = profit->guessradius;
    qprofit->guessaspect = profit->guessaspect;
    qprofit->guessposang = profit->guessposang;
    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,
                        FLAG(obj2.prof_class_star)? qprofit->resi : NULL);
    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 && (err = sumqw2*sump*sump
                +sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)>0.0)?
                        sqrt(err) / (sump*sump) : 1.0;
      }
    }
 
/* clean up. */
  fft_reset();
 
  return;
  }
 
 
/****** profit_dfit ***********************************************************
PROTO   void profit_dfit(profitstruct *profit, profitstruct *dprofit,
                picstruct *field, picstruct *dfield,
                picstruct *wfield, picstruct *dwfield,
                objstruct *obj, obj2struct *obj2)
PURPOSE Fit profile(s) convolved with the PSF to a detected object on the
        detection image, and use the measurement image to scale the flux.
INPUT   Pointer to the measurement profile-fitting structure,
        pointer to the detection profile-fitting structure,
        pointer to the measurement field,
        pointer to the detection field,
        pointer to the measurement field weight,
        pointer to the detection 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 16/02/2013
 ***/
void    profit_dfit(profitstruct *profit, profitstruct *dprofit,
                picstruct *field, picstruct *dfield,
                picstruct *wfield, picstruct *dwfield,
                objstruct *obj, obj2struct *obj2)
  {
    psfstruct           *dpsf;
    double              emx2,emy2,emxy, a , cp,sp, cn, bn, n,
                        sumn,sumd;
    PIXTYPE             valn,vald,w2;
    float               param0[PARAM_NPARAM], param1[PARAM_NPARAM],
                        param[PARAM_NPARAM],
                        **list,
                        *cov, *pix,
                        psf_fwhm, dchi2, err, aspect, chi2, ffac;
    int                 *index,
                        i,j,p, nparam, nparam2, ncomp, nprof;
 
  nparam = dprofit->nparam;
  nparam2 = nparam*nparam;
  nprof = dprofit->nprof;
  if (dprofit->psfdft)
    QFFTWF_FREE(dprofit->psfdft);
 
  dpsf = dprofit->psf;
  dprofit->pixstep = dpsf->pixstep;
  obj2->dprof_flag = 0;
 
/* Create pixmaps at image resolution */
  dprofit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */
  dprofit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */
  psf_fwhm = dpsf->masksize[0]*dpsf->pixstep;
  dprofit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)
                *1.2)/2)*2 + 1;
  dprofit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)
                *1.2)/2)*2 + 1;
  if (dprofit->objnaxisn[1]<dprofit->objnaxisn[0])
    dprofit->objnaxisn[1] = dprofit->objnaxisn[0];
  else
    dprofit->objnaxisn[0] = dprofit->objnaxisn[1];
  if (dprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
    {
    dprofit->subsamp = ceil((float)dprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
    dprofit->objnaxisn[1] = (dprofit->objnaxisn[0] /= (int)dprofit->subsamp);
    obj2->dprof_flag |= PROFLAG_OBJSUB;
    }
  else
    dprofit->subsamp = 1.0;
  dprofit->nobjpix = dprofit->objnaxisn[0]*dprofit->objnaxisn[1];
 
/* Use (dirty) global variables to interface with lmfit */
  the_field = dfield;
  the_wfield = dwfield;
  dprofit->obj = obj;
  dprofit->obj2 = obj2;
 
  dprofit->nresi = profit_copyobjpix(dprofit, dfield, dwfield);
 
/* Check if the number of constraints exceeds the number of free parameters */
  if (dprofit->nresi < nparam || profit->nresi < 1)
    {
    obj2->dprof_flag |= PROFLAG_NOTCONST;
    return;
    }
 
/* Create pixmap at PSF resolution */
  dprofit->modnaxisn[0] =
        ((int)(dprofit->objnaxisn[0]*dprofit->subsamp/dprofit->pixstep
                +0.4999)/2+1)*2;
  dprofit->modnaxisn[1] =
        ((int)(dprofit->objnaxisn[1]*dprofit->subsamp/dprofit->pixstep
                +0.4999)/2+1)*2;
  if (dprofit->modnaxisn[1] < dprofit->modnaxisn[0])
    dprofit->modnaxisn[1] = dprofit->modnaxisn[0];
  else
    dprofit->modnaxisn[0] = dprofit->modnaxisn[1];
  if (dprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)
    {
    dprofit->pixstep = (double)dprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;
    dprofit->modnaxisn[0] = dprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;
    obj2->dprof_flag |= PROFLAG_MODSUB;
    }
  dprofit->nmodpix = dprofit->modnaxisn[0]*dprofit->modnaxisn[1];
 
/* Compute the local PSF */
  profit_psf(dprofit);
 
/* Set initial guesses and boundaries */
  dprofit->guesssigbkg = dprofit->sigma = obj->dsigbkg;
  dprofit->guessdx = obj->mx - (int)(obj->mx+0.49999);
  dprofit->guessdy = obj->my - (int)(obj->my+0.49999);
  if ((dprofit->guessflux = obj->dflux) <= 0.0)
    dprofit->guessflux = 0.0;
  if ((dprofit->guessfluxmax = 10.0*obj->dnpix*obj->dsigbkg*obj->dsigbkg)
                         <= dprofit->guessflux)
    dprofit->guessfluxmax = dprofit->guessflux;
  if (dprofit->guessfluxmax <= 0.0)
    dprofit->guessfluxmax = 1.0;
  if ((dprofit->guessradius = sqrtf(obj->a*obj->b)*1.17)<0.5*dpsf->fwhm)
    dprofit->guessradius = 0.5*dpsf->fwhm;
  dprofit->guessaspect = obj->b/obj->a;
  dprofit->guessposang = obj->theta;
 
  profit_resetparams(dprofit);
 
/* Actual minimisation */
  fft_reset();
 
  dprofit->niter = profit_minimize(dprofit, PROFIT_MAXITER);
 
  if (dprofit->nlimmin)
    obj2->dprof_flag |= PROFLAG_MINLIM;
  if (dprofit->nlimmax)
    obj2->dprof_flag |= PROFLAG_MAXLIM;
 
  for (p=0; p<nparam; p++)
    dprofit->paramerr[p]= sqrt(dprofit->covar[p*(nparam+1)]);
 
  obj2->dprof_niter = dprofit->niter;
 
/* Now inject fitted parameters into the measurement model */
  fft_reset();
  profit_residuals(profit,field,wfield, 0.0, dprofit->paraminit, NULL);
 
/* Compute flux correction */
  sumn = sumd = 0.0;
  for (p=0; p<profit->nobjpix; p++)
    if (profit->objweight[p]>0 && profit->objpix[p]>-BIG)
      {
      w2 = profit->objweight[p]*profit->objweight[p] * profit->lmodpix[p];
      sumn += (double)(w2*profit->objpix[p]);
      sumd += (double)(w2*profit->lmodpix[p]);
      }
 
  ffac = (float)(sumn/sumd);
  obj2->flux_dprof = sumd!=0.0? dprofit->flux*ffac: 0.0f;
  obj2->fluxerr_dprof = sumd!=0.0? dprofit->flux/sqrtf((float)sumd): 0.0f;
 
  if (FLAG(obj2.dprof_chi2))
    {
/*-- Compute reduced chi2 on measurement image */
    pix = profit->lmodpix;
    for (p=profit->nobjpix; p--;)
      *(pix++) *= ffac;
    profit_compresi(profit, 0.0, profit->resi);
    obj2->dprof_chi2 = (profit->nresi > dprofit->nparam)?
                profit->chi2 / (profit->nresi - dprofit->nparam) : 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 13/02/2013
 ***/
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 / (profit->subsamp*profit->subsamp);
  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 15/01/2013
 ***/
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 + profit->npresi,
                        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 16/01/2013
 ***/
void    profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
  {
   profitstruct         *profit;
   profstruct           **prof;
   double               *dpar0, *dresi, *fvect;
   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);
 
    profit_presiduals(profit, dpar, profit->presi);
 
    fvect = fvec;
    for (p=0; p<profit->nresi; p++)
      *(fvect++) = profit->resi[p];
    for (p=0; p<profit->npresi; p++)
      *(fvect++) = profit->presi[p];
 
    }
 
//  profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
  profit->iter++;
 
  return;
  }
 
 
/****** profit_residuals ******************************************************
PROTO   float *profit_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
        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_presiduals *****************************************************
PROTO   float *profit_presiduals(profitstruct *profit, float *param,
                float *presi)
PURPOSE Compute the vector of prior "residuals" for the model parameters.
INPUT   Profile-fitting structure,
        pointer to the (unbound) model parameters,
        pointer to the computed prior "residuals" (output).
OUTPUT  Vector of residuals.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 15/01/2013
 ***/
float   *profit_presiduals(profitstruct *profit, double *dparam, float *presi)
  {
   float        *presit;
   int          p;
 
  presit = presi;
  for (p=0; p<profit->nparam; p++)
    if (profit->dparampsig[p]>0.0)
      *(presit++) = (float)((dparam[p] - profit->dparampcen[p])
                                / profit->dparampsig[p]);
 
  return presi;
  }
 
 
/****** profit_compresi ******************************************************
PROTO   float *profit_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 24/04/2013
 ***/
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, fluxnorm;
   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 = (float)(profit->modnaxisn[0]/2);
  xcout = ((int)(profit->subsamp*profit->objnaxisn[0])/2 + 0.5)
                / profit->subsamp - 0.5;
  if ((dx=profit->paramlist[PARAM_X]))
    xcout += *dx/profit->subsamp;
 
  xsin = xcin - xcout*invpixstep;                       /* Input start x-coord*/
 
  if ((int)xsin >= profit->modnaxisn[0] || !finitef(xsin))
    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 = (float)(profit->modnaxisn[1]/2);
  ycout = ((int)(profit->subsamp*profit->objnaxisn[1])/2 + 0.5)
                / profit->subsamp - 0.5;
  if ((dy=profit->paramlist[PARAM_Y]))
    ycout += *dy/profit->subsamp;
 
  ysin = ycin - ycout*invpixstep;               /* Input start y-coord*/
  if ((int)ysin >= profit->modnaxisn[1] || !finitef(ysin))
    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 12/07/2012
 ***/
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);
  QFFTWF_FREE(fdx);
  QFFTWF_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))


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	305	bertin	`* Copyright: (C) 2006-2013 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	307	bertin	`* Last modified: 05/04/2013`
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	284	bertin	`profitstruct theprofit,thedprofit, 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	305	bertin	`QMALLOC16(profit->presi, float, profit->nparam);`
115	211	bertin	`QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);`
116	244	bertin	`profit->nprof = nmodels;`
117	207	bertin	`profit->fluxfac = 1.0; /* Default */`
118	50	bertin
119	46	bertin	`return profit;`
120	42	bertin	`}`
121
122
123	46	bertin	`/**** profit_end **********************************************************`
124	42	bertin	`PROTO void prof_end(profstruct *prof)`
125	46	bertin	`PURPOSE End (deallocate) a profile-fitting structure.`
126	42	bertin	`INPUT Prof structure.`
127			`OUTPUT -.`
128	117	bertin	`NOTES -.`
129	42	bertin	`AUTHOR E. Bertin (IAP)`
130	297	bertin	`VERSION 12/07/2012`
131	42	bertin	`***/`
132	46	bertin	`void profit_end(profitstruct *profit)`
133	42	bertin	`{`
134	49	bertin	`int p;`
135
136			`for (p=0; p<profit->nprof; p++)`
137			`prof_end(profit->prof[p]);`
138	221	bertin	`free(profit->modpix);`
139			`free(profit->modpix2);`
140			`free(profit->cmodpix);`
141			`free(profit->psfpix);`
142			`free(profit->lmodpix);`
143			`free(profit->lmodpix2);`
144			`free(profit->objpix);`
145			`free(profit->objweight);`
146			`free(profit->resi);`
147	305	bertin	`free(profit->presi);`
148	49	bertin	`free(profit->prof);`
149	123	bertin	`free(profit->covar);`
150	297	bertin	`QFFTWF_FREE(profit->psfdft);`
151	46	bertin	`free(profit);`
152	42	bertin
153			`return;`
154			`}`
155
156
157	65	bertin	`/**** profit_fit **********************************************************`
158			`PROTO void profit_fit(profitstruct profit, picstruct field,`
159			`picstruct wfield, objstruct obj, obj2struct *obj2)`
160	46	bertin	`PURPOSE Fit profile(s) convolved with the PSF to a detected object.`
161	42	bertin	`INPUT Array of profile structures,`
162			`Number of profiles,`
163	52	bertin	`Pointer to the profile-fitting structure,`
164	42	bertin	`Pointer to the field,`
165			`Pointer to the field weight,`
166			`Pointer to the obj.`
167	46	bertin	`OUTPUT Pointer to an allocated fit structure (containing details about the`
168			`fit).`
169	48	bertin	`NOTES It is a modified version of the lm_minimize() of lmfit.`
170	46	bertin	`AUTHOR E. Bertin (IAP)`
171	307	bertin	`VERSION 05/04/2013`
172	46	bertin	`***/`
173	65	bertin	`void profit_fit(profitstruct *profit,`
174	52	bertin	`picstruct field, picstruct wfield,`
175	50	bertin	`objstruct obj, obj2struct obj2)`
176	46	bertin	`{`
177	244	bertin	`profitstruct pprofit, qprofit;`
178	221	bertin	`patternstruct *pattern;`
179	52	bertin	`psfstruct *psf;`
180	50	bertin	`checkstruct *check;`
181	244	bertin	`double emx2,emy2,emxy, a , cp,sp, cn, bn, n,`
182	292	bertin	`sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0, err;`
183	244	bertin	`PIXTYPE valp,valq,sig2;`
184	221	bertin	`float param0[PARAM_NPARAM], param1[PARAM_NPARAM],`
185			`param[PARAM_NPARAM],`
186			`**list,`
187	209	bertin	`*cov,`
188	292	bertin	`psf_fwhm, dchi2, aspect, chi2;`
189	209	bertin	`int *index,`
190			`i,j,p, nparam, nparam2, ncomp, nprof;`
191	46	bertin
192	128	bertin	`nparam = profit->nparam;`
193	208	bertin	`nparam2 = nparam*nparam;`
194	194	bertin	`nprof = profit->nprof;`
195	294	bertin
196	52	bertin	`if (profit->psfdft)`
197	297	bertin	`QFFTWF_FREE(profit->psfdft);`
198	52	bertin
199	297	bertin
200	52	bertin	`psf = profit->psf;`
201	118	bertin	`profit->pixstep = psf->pixstep;`
202	211	bertin	`obj2->prof_flag = 0;`
203	52	bertin
204			`/* Create pixmaps at image resolution */`
205	211	bertin	`profit->ix = (int)(obj->mx + 0.49999);/* internal convention: 1st pix = 0 */`
206			`profit->iy = (int)(obj->my + 0.49999);/* internal convention: 1st pix = 0 */`
207	51	bertin	`psf_fwhm = psf->masksize[0]*psf->pixstep;`
208	233	bertin	`profit->objnaxisn[0] = (((int)((obj->xmax-obj->xmin+1) + psf_fwhm + 0.499)`
209	75	bertin	`1.2)/2)2 + 1;`
210	233	bertin	`profit->objnaxisn[1] = (((int)((obj->ymax-obj->ymin+1) + psf_fwhm + 0.499)`
211	75	bertin	`1.2)/2)2 + 1;`
212	117	bertin	`if (profit->objnaxisn[1]<profit->objnaxisn[0])`
213			`profit->objnaxisn[1] = profit->objnaxisn[0];`
214			`else`
215			`profit->objnaxisn[0] = profit->objnaxisn[1];`
216	211	bertin	`if (profit->objnaxisn[0]>PROFIT_MAXOBJSIZE)`
217			`{`
218	212	bertin	`profit->subsamp = ceil((float)profit->objnaxisn[0]/PROFIT_MAXOBJSIZE);`
219			`profit->objnaxisn[1] = (profit->objnaxisn[0] /= (int)profit->subsamp);`
220	211	bertin	`obj2->prof_flag \|= PROFLAG_OBJSUB;`
221			`}`
222			`else`
223	212	bertin	`profit->subsamp = 1.0;`
224	221	bertin	`profit->nobjpix = profit->objnaxisn[0]*profit->objnaxisn[1];`
225	117	bertin
226	305	bertin	`/* Create pixmap at model resolution */`
227			`profit->modnaxisn[0] =`
228			`((int)(profit->objnaxisn[0]*profit->subsamp/profit->pixstep`
229			`+0.4999)/2+1)*2;`
230			`profit->modnaxisn[1] =`
231			`((int)(profit->objnaxisn[1]*profit->subsamp/profit->pixstep`
232			`+0.4999)/2+1)*2;`
233			`if (profit->modnaxisn[1] < profit->modnaxisn[0])`
234			`profit->modnaxisn[1] = profit->modnaxisn[0];`
235			`else`
236			`profit->modnaxisn[0] = profit->modnaxisn[1];`
237			`if (profit->modnaxisn[0]>PROFIT_MAXMODSIZE)`
238			`{`
239			`profit->pixstep = (double)profit->modnaxisn[0] / PROFIT_MAXMODSIZE;`
240			`profit->modnaxisn[0] = profit->modnaxisn[1] = PROFIT_MAXMODSIZE;`
241			`obj2->prof_flag \|= PROFLAG_MODSUB;`
242			`}`
243			`profit->nmodpix = profit->modnaxisn[0]*profit->modnaxisn[1];`
244
245	78	bertin	`/* Use (dirty) global variables to interface with lmfit */`
246			`the_field = field;`
247			`the_wfield = wfield;`
248			`theprofit = profit;`
249	146	bertin	`profit->obj = obj;`
250			`profit->obj2 = obj2;`
251	78	bertin
252	305	bertin	`/* Compute the local PSF */`
253			`profit_psf(profit);`
254
255	146	bertin	`profit->nresi = profit_copyobjpix(profit, field, wfield);`
256	305	bertin	`profit->npresi = 0;`
257	211	bertin	`/* Check if the number of constraints exceeds the number of free parameters */`
258	128	bertin	`if (profit->nresi < nparam)`
259	49	bertin	`{`
260	176	bertin	`if (FLAG(obj2.prof_vector))`
261			`for (p=0; p<nparam; p++)`
262			`obj2->prof_vector[p] = 0.0;`
263	208	bertin	`if (FLAG(obj2.prof_errvector))`
264			`for (p=0; p<nparam; p++)`
265			`obj2->prof_errvector[p] = 0.0;`
266			`if (FLAG(obj2.prof_errmatrix))`
267			`for (p=0; p<nparam2; p++)`
268			`obj2->prof_errmatrix[p] = 0.0;`
269	50	bertin	`obj2->prof_niter = 0;`
270	211	bertin	`obj2->prof_flag \|= PROFLAG_NOTCONST;`
271	50	bertin	`return;`
272	49	bertin	`}`
273	48	bertin
274	50	bertin	`/* Set initial guesses and boundaries */`
275	284	bertin	`profit->guesssigbkg = profit->sigma = obj->sigbkg;`
276			`profit->guessdx = obj->mx - (int)(obj->mx+0.49999);`
277			`profit->guessdy = obj->my - (int)(obj->my+0.49999);`
278	246	bertin	`if ((profit->guessflux = obj2->flux_auto) <= 0.0)`
279			`profit->guessflux = 0.0;`
280			`if ((profit->guessfluxmax = 10.0*obj2->fluxerr_auto) <= profit->guessflux)`
281			`profit->guessfluxmax = profit->guessflux;`
282			`if (profit->guessfluxmax <= 0.0)`
283			`profit->guessfluxmax = 1.0;`
284	284	bertin	`if ((profit->guessradius = 0.5*psf->fwhm) < obj2->hl_radius)`
285			`profit->guessradius = obj2->hl_radius;`
286			`profit->guessaspect = obj->b/obj->a;`
287			`profit->guessposang = obj->theta;`
288	115	bertin
289	146	bertin	`profit_resetparams(profit);`
290	50	bertin
291	123	bertin	`/* Actual minimisation */`
292	209	bertin	`fft_reset();`
293	211	bertin	`the_gal++;`
294	221	bertin
295	235	bertin	`/*`
296			`char str[1024];`
297			`sprintf(str, "obj_%04d.fits", the_gal);`
298			`catstruct *bcat;`
299			`float bpix, opix,lmodpix,objpix;`
300			`bcat=read_cat("base.fits");`
301			`QMALLOC(bpix, float, field->npix);`
302			`QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);`
303			`read_body(bcat->tab, bpix, field->npix);`
304			`free_cat(&bcat,1);`
305			`bcat=read_cat(str);`
306			`QMALLOC(opix, float, profit->nobjpix);`
307			`QFSEEK(bcat->file, bcat->tab->bodypos, SEEK_SET, bcat->filename);`
308			`read_body(bcat->tab, opix, profit->nobjpix);`
309			`free_cat(&bcat,1);`
310			`addfrombig(bpix, field->width, field->height,`
311			`profit->objpix, profit->objnaxisn[0],profit->objnaxisn[1],`
312			`profit->ix,profit->iy, -1.0);`
313			`objpix = profit->objpix;`
314			`lmodpix = opix;`
315			`for (i=profit->nobjpix; i--;)`
316			`(objpix++) += (lmodpix++);`
317			`free(bpix);`
318			`free(opix);`
319			`*/`
320	79	bertin	`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
321	221	bertin	`/*`
322	235	bertin	`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
323			`check=initcheck(str, CHECK_OTHER,1);`
324			`check->width = profit->objnaxisn[0];`
325			`check->height = profit->objnaxisn[1];`
326			`reinitcheck(field,check);`
327			`memcpy(check->pix, profit->lmodpix, profit->nobjpix*sizeof(float));`
328			`reendcheck(field,check);`
329			`endcheck(check);`
330
331	221	bertin	`chi2 = profit->chi2;`
332			`for (p=0; p<nparam; p++)`
333			`param1[p] = profit->paraminit[p];`
334			`profit_resetparams(profit);`
335			`for (p=0; p<nparam; p++)`
336			`profit->paraminit[p] = param1[p] + (profit->paraminit[p]<param1[p]?1.0:-1.0)`
337			`* sqrt(profit->covar[p*(nparam+1)]);`
338	200	bertin
339	221	bertin	`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
340			`if (chi2<profit->chi2)`
341			`for (p=0; p<nparam; p++)`
342			`profit->paraminit[p] = param1[p];`
343
344			`list = profit->paramlist;`
345			`index = profit->paramindex;`
346			`for (i=0; i<PARAM_NPARAM; i++)`
347			`if (list[i] && i!= PARAM_SPHEROID_ASPECT && i!=PARAM_SPHEROID_POSANG)`
348			`profit->freeparam_flag[index[i]] = 0;`
349			`profit->niter = profit_minimize(profit, PROFIT_MAXITER);`
350			`*/`
351	284	bertin
352	246	bertin	`if (profit->nlimmin)`
353			`obj2->prof_flag \|= PROFLAG_MINLIM;`
354			`if (profit->nlimmax)`
355			`obj2->prof_flag \|= PROFLAG_MAXLIM;`
356
357	128	bertin	`for (p=0; p<nparam; p++)`
358			`profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);`
359	126	bertin
360	50	bertin	`/* CHECK-Images */`
361	221	bertin	`if ((check = prefs.check[CHECK_PROFILES]))`
362			`{`
363			`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
364	305	bertin	`if (profit->subsamp>1.0)`
365			`addcheck_resample(check, profit->lmodpix,`
366			`profit->objnaxisn[0],profit->objnaxisn[1],`
367			`profit->ix,profit->iy, 1.0/profit->subsamp,`
368			`1.0/(profit->subsamp*profit->subsamp));`
369			`else`
370			`addcheck(check, profit->lmodpix,`
371			`profit->objnaxisn[0],profit->objnaxisn[1],`
372	221	bertin	`profit->ix,profit->iy, 1.0);`
373			`}`
374	235	bertin
375	50	bertin	`if ((check = prefs.check[CHECK_SUBPROFILES]))`
376	51	bertin	`{`
377	221	bertin	`profit_residuals(profit,field,wfield, 0.0, profit->paraminit, NULL);`
378	305	bertin	`if (profit->subsamp>1.0)`
379			`addcheck_resample(check, profit->lmodpix,`
380			`profit->objnaxisn[0],profit->objnaxisn[1],`
381			`profit->ix,profit->iy, 1.0/profit->subsamp,`
382			`-1.0/(profit->subsamp*profit->subsamp));`
383			`else`
384			`addcheck(check, profit->lmodpix,`
385			`profit->objnaxisn[0],profit->objnaxisn[1],`
386	145	bertin	`profit->ix,profit->iy, -1.0);`
387	51	bertin	`}`
388	221	bertin	`if ((check = prefs.check[CHECK_SPHEROIDS]))`
389	51	bertin	`{`
390	221	bertin	`/-- Set to 0 flux components that do not belong to spheroids /`
391			`for (p=0; p<profit->nparam; p++)`
392			`param[p] = profit->paraminit[p];`
393			`list = profit->paramlist;`
394			`index = profit->paramindex;`
395			`for (i=0; i<PARAM_NPARAM; i++)`
396			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
397			`param[index[i]] = 0.0;`
398			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
399	307	bertin	`if (profit->subsamp>1.0)`
400			`addcheck_resample(check, profit->lmodpix,`
401			`profit->objnaxisn[0],profit->objnaxisn[1],`
402			`profit->ix,profit->iy, 1.0/profit->subsamp,`
403			`1.0/(profit->subsamp*profit->subsamp));`
404			`else`
405			`addcheck(check, profit->lmodpix,`
406			`profit->objnaxisn[0],profit->objnaxisn[1],`
407	145	bertin	`profit->ix,profit->iy, 1.0);`
408	51	bertin	`}`
409	221	bertin	`if ((check = prefs.check[CHECK_SUBSPHEROIDS]))`
410			`{`
411			`/-- Set to 0 flux components that do not belong to spheroids /`
412			`for (p=0; p<profit->nparam; p++)`
413			`param[p] = profit->paraminit[p];`
414			`list = profit->paramlist;`
415			`index = profit->paramindex;`
416			`for (i=0; i<PARAM_NPARAM; i++)`
417			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
418			`param[index[i]] = 0.0;`
419			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
420	307	bertin	`if (profit->subsamp>1.0)`
421			`addcheck_resample(check, profit->lmodpix,`
422			`profit->objnaxisn[0],profit->objnaxisn[1],`
423			`profit->ix,profit->iy, 1.0/profit->subsamp,`
424			`-1.0/(profit->subsamp*profit->subsamp));`
425			`else`
426			`addcheck(check, profit->lmodpix,`
427			`profit->objnaxisn[0],profit->objnaxisn[1],`
428	221	bertin	`profit->ix,profit->iy, -1.0);`
429			`}`
430			`if ((check = prefs.check[CHECK_DISKS]))`
431			`{`
432			`/-- Set to 0 flux components that do not belong to disks /`
433			`for (p=0; p<profit->nparam; p++)`
434			`param[p] = profit->paraminit[p];`
435			`list = profit->paramlist;`
436			`index = profit->paramindex;`
437			`for (i=0; i<PARAM_NPARAM; i++)`
438			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
439			`param[index[i]] = 0.0;`
440			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
441	307	bertin	`if (profit->subsamp>1.0)`
442			`addcheck_resample(check, profit->lmodpix,`
443			`profit->objnaxisn[0],profit->objnaxisn[1],`
444			`profit->ix,profit->iy, 1.0/profit->subsamp,`
445			`1.0/(profit->subsamp*profit->subsamp));`
446			`else`
447			`addcheck(check, profit->lmodpix,`
448			`profit->objnaxisn[0],profit->objnaxisn[1],`
449	221	bertin	`profit->ix,profit->iy, 1.0);`
450			`}`
451			`if ((check = prefs.check[CHECK_SUBDISKS]))`
452			`{`
453			`/-- Set to 0 flux components that do not belong to disks /`
454			`for (p=0; p<profit->nparam; p++)`
455			`param[p] = profit->paraminit[p];`
456			`list = profit->paramlist;`
457			`index = profit->paramindex;`
458			`for (i=0; i<PARAM_NPARAM; i++)`
459			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
460			`param[index[i]] = 0.0;`
461			`profit_residuals(profit,field,wfield, 0.0, param, NULL);`
462	307	bertin	`if (profit->subsamp>1.0)`
463			`addcheck_resample(check, profit->lmodpix,`
464			`profit->objnaxisn[0],profit->objnaxisn[1],`
465			`profit->ix,profit->iy, 1.0/profit->subsamp,`
466			`-1.0/(profit->subsamp*profit->subsamp));`
467			`else`
468			`addcheck(check, profit->lmodpix,`
469			`profit->objnaxisn[0],profit->objnaxisn[1],`
470	221	bertin	`profit->ix,profit->iy, -1.0);`
471			`}`
472	235	bertin
473	221	bertin	`/* Compute compressed residuals */`
474			`profit_residuals(profit,field,wfield, 10.0, profit->paraminit,profit->resi);`
475
476	73	bertin	`/* Fill measurement parameters */`
477	69	bertin	`if (FLAG(obj2.prof_vector))`
478			`{`
479	128	bertin	`for (p=0; p<nparam; p++)`
480	221	bertin	`obj2->prof_vector[p]= profit->paraminit[p];`
481	69	bertin	`}`
482	127	bertin	`if (FLAG(obj2.prof_errvector))`
483			`{`
484	128	bertin	`for (p=0; p<nparam; p++)`
485			`obj2->prof_errvector[p]= profit->paramerr[p];`
486	127	bertin	`}`
487	208	bertin	`if (FLAG(obj2.prof_errmatrix))`
488			`{`
489			`for (p=0; p<nparam2; p++)`
490			`obj2->prof_errmatrix[p]= profit->covar[p];`
491			`}`
492	89	bertin
493	79	bertin	`obj2->prof_niter = profit->niter;`
494	73	bertin	`obj2->flux_prof = profit->flux;`
495	194	bertin	`if (FLAG(obj2.fluxerr_prof))`
496			`{`
497			`err = 0.0;`
498	209	bertin	`cov = profit->covar;`
499			`index = profit->paramindex;`
500			`list = profit->paramlist;`