Astromatic - public software.sextractor - Blame - Rev 285

[/] [branches/] [multi/] [src/] [profit.c] - Blame information for rev 285

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


/*
*                               profit.c
*
* Fit a range of galaxy models to an image.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*       This file part of:      SExtractor
*
*       Copyright:              (C) 2006-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
*
*       License:                GNU General Public License
*
*       SExtractor is free software: you can redistribute it and/or modify
*       it under the terms of the GNU General Public License as published by
*       the Free Software Foundation, either version 3 of the License, or
*       (at your option) any later version.
*       SExtractor is distributed in the hope that it will be useful,
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*       GNU General Public License for more details.
*       You should have received a copy of the GNU General Public License
*       along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
*       Last modified:          06/05/2012
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
#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        "graph.h"
#include        "image.h"
#include        "misc.h"
#include        "pattern.h"
#include        "psf.h"
#include        "profit.h"
#include        "subimage.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,2,3,4,5,0,0,
                                        1,2,3,4,5,0,0,0,0,
                                        1,2,3,4,5,0,0,0,
                                        1,0,0,0,0,0,0,0,
                                        1,0,0,
                                        1,0,0,
                                        1,0,0
                                        };
 
int     the_gal;
 
/****** profit_init ***********************************************************
PROTO   profitstruct profit_init(obj2struct *obj2, unsigned int modeltype)
PURPOSE Allocate and initialize a new profile-fitting structure.
INPUT   pointer to the obj2,
        model type.
OUTPUT  A pointer to an allocated profit structure.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
profitstruct    *profit_init(obj2struct *obj2, unsigned int modeltype)
  {
   profitstruct         *profit;
   subprofitstruct      *subprofit;
   subimagestruct       *subimage;
   double               psf_fwhm;
   int                  s,t, nmodels, npix, nsub;
 
  QCALLOC(profit, profitstruct, 1);
 
  profit->obj2 = obj2;
 
  QCALLOC(profit->subprofit, subprofitstruct, obj2->nsubimage);
  subprofit = profit->subprofit;
  subimage = obj2->subimage;
  nsub = 0;
  for (s=0; s<obj2->nsubimage; s++, subimage++)
    {
    subprofit->field = subimage->field;
    subprofit->wfield = subimage->wfield;
    subprofit->psf = subimage->field->psf;
    subprofit->pixstep = subprofit->psf->pixstep;
    subprofit->sigma = obj2->sigbkg[s];
/*-- Set initial guesses and boundaries */
    if ((subprofit->guessradius = 0.5*subprofit->psf->fwhm) < obj2->hl_radius)
      subprofit->guessradius = obj2->hl_radius;
    if ((subprofit->guessflux = obj2->flux_auto[s]) <= 0.0)
      subprofit->guessflux = 0.0;
    if ((subprofit->guessfluxmax = 10.0*obj2->fluxerr_auto[s])
        <= subprofit->guessflux)
      subprofit->guessfluxmax = subprofit->guessflux;
    if (subprofit->guessfluxmax <= 0.0)
      subprofit->guessfluxmax = 1.0;
    subprofit->fluxfac = 1.0;   /* Default */
/*-- Create pixmaps at image resolution */
    subprofit->ix = (int)(obj2->mx + 0.49999); /* 1st pix=0 */
    subprofit->iy = (int)(obj2->my + 0.49999); /* 1st pix=0 */
    psf_fwhm = subprofit->psf->masksize[0]*subprofit->psf->pixstep;
    subprofit->objnaxisn[0] = ((int)((subimage->imsize[0] + psf_fwhm + 0.499)
                *1.2)/2)*2 + 1;
    subprofit->objnaxisn[1] = ((int)((subimage->imsize[1] + psf_fwhm + 0.499)
                *1.2)/2)*2 + 1;
    if (subprofit->objnaxisn[1] < subprofit->objnaxisn[0])
      subprofit->objnaxisn[1] = subprofit->objnaxisn[0];
    else
      subprofit->objnaxisn[0] = subprofit->objnaxisn[1];
 
    if (subprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
      {
      subprofit->subsamp
                = ceil((float)subprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
      subprofit->objnaxisn[1]
                = (subprofit->objnaxisn[0] /= (int)subprofit->subsamp);
      profit->flag |= PROFLAG_OBJSUB;   /* !CHECK */
      }
    else
      subprofit->subsamp = 1.0;
    subprofit->nobjpix = subprofit->objnaxisn[0]*subprofit->objnaxisn[1];
 
/* Allocate memory for the data and the resampled model */
    QMALLOC16(subprofit->objpix, PIXTYPE, subprofit->nobjpix);
    QMALLOC16(subprofit->objweight, PIXTYPE, subprofit->nobjpix);
    QMALLOC16(subprofit->lmodpix, PIXTYPE, subprofit->nobjpix);
    QMALLOC16(subprofit->lmodpix2, PIXTYPE, subprofit->nobjpix);
 
/*-- Create pixmap at PSF resolution */
    subprofit->modnaxisn[0] =
        ((int)(subprofit->objnaxisn[0]*subprofit->subsamp
                /subprofit->pixstep+0.4999)/2+1)*2;
    subprofit->modnaxisn[1] =
        ((int)(subprofit->objnaxisn[1]*subprofit->subsamp
                /subprofit->pixstep+0.4999)/2+1)*2;
    if (subprofit->modnaxisn[1] < subprofit->modnaxisn[0])
      subprofit->modnaxisn[1] = subprofit->modnaxisn[0];
    else
      subprofit->modnaxisn[0] = subprofit->modnaxisn[1];
    if (subprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)
      {
      subprofit->pixstep = (double)subprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;
      subprofit->modnaxisn[0] = subprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;
      profit->flag |= PROFLAG_MODSUB;
      }
    subprofit->nmodpix = subprofit->modnaxisn[0]*subprofit->modnaxisn[1];
 
/* Allocate memory for the complete model */
    QMALLOC16(subprofit->modpix, float, subprofit->nmodpix);
    QMALLOC16(subprofit->modpix2, float, subprofit->nmodpix);
    QMALLOC16(subprofit->cmodpix, float, subprofit->nmodpix);
    QMALLOC16(subprofit->psfpix, float, subprofit->nmodpix);
    if ((npix = subprofit_copyobjpix(subprofit, subimage)))
      {
      profit->nresi += npix;
/*---- Compute the local PSF */
      subprofit_psf(subprofit, obj2);
      subprofit->index = nsub++;
      subprofit++;
      }
    else
      subprofit_end(subprofit);
    }
 
  profit->nsubprofit = nsub;
 
  QMALLOC16(profit->resi, float, profit->nresi);
 
  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);
  profit->nprof = nmodels;
 
  QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
 
  profit_resetparams(profit);
 
  return profit;
  }
 
 
/****** profit_end ************************************************************
PROTO   void profit_end(profitstruct *profit)
PURPOSE End (deallocate) a profile-fitting structure.
INPUT   Prof structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    profit_end(profitstruct *profit)
  {
   subprofitstruct      *subprofit;
   int                  p,s;
 
  for (p=0; p<profit->nprof; p++)
    prof_end(profit->prof[p]);
  free(profit->prof);
  subprofit = profit->subprofit;
  for (s=profit->nsubprofit; s--;)
    subprofit_end(subprofit++);
 
  free(profit->subprofit);
  free(profit->resi);
  free(profit->covar);
 
  free(profit);
 
  return;
  }
 
 
/****** subprofit_end *********************************************************
PROTO   void subprofit_end(subprofitstruct *profit)
PURPOSE End (free content) a subprofile-fitting structure.
INPUT   SubProfit structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 03/05/2012
 ***/
void    subprofit_end(subprofitstruct *subprofit)
  {
  free(subprofit->objpix);
  free(subprofit->objweight);
  free(subprofit->lmodpix);
  free(subprofit->lmodpix2);
  free(subprofit->modpix);
  free(subprofit->modpix2);
  free(subprofit->cmodpix);
  free(subprofit->psfpix);
  free(subprofit->psfdft);
 
  return;
  }
 
 
/****** profit_fit ************************************************************
PROTO   int profit_fit(profitstruct *profit, obj2struct *obj2)
PURPOSE Fit profile(s) convolved with the PSF to a detected object.
INPUT   Pointer to the model structure,
        pointer to the obj2.
OUTPUT  Number of minimization iterations (0 in case of error).
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/10/2011
 ***/
int     profit_fit(profitstruct *profit, obj2struct *obj2)
  {
   int                  p, nparam, nparam2;
 
  nparam = profit->nparam;
  nparam2 = nparam*nparam;
 
/* reset all flags except those set in profit_init() */
  profit->flag &= PROFLAG_OBJSUB|PROFLAG_MODSUB;
 
/* 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;
    profit->flag |= PROFLAG_NOTCONST;
    return 0;
    }
 
/* Actual minimisation */
  fft_reset();
  return (profit->niter = profit_minimize(profit, PROFIT_MAXITER));
  }
 
 
/****** profit_measure ********************************************************
PROTO   void profit_measure(profitstruct *profit, obj2struct *obj2)
PURPOSE Perform measurements on the fitted model.
INPUT   Pointer to the model structure,
        pointer to the obj2.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    profit_measure(profitstruct *profit, obj2struct *obj2)
  {
    profitstruct        *pprofit, *qprofit;
    subprofitstruct     *subprofit;
    patternstruct       *pattern;
    checkstruct         *check;
    double              emx2,emy2,emxy, a , cp,sp, cn, bn, n;
    float               param0[PARAM_NPARAM], param1[PARAM_NPARAM],
                        param[PARAM_NPARAM],
                        **list,
                        *cov,
                        psf_fwhm, err, aspect, chi2, flux;
    int                 *index,
                        i,j,p,s, nparam, nparam2, ncomp, nsub;
 
  nparam = profit->nparam;
  nparam2 = nparam*nparam;
  nsub = profit->nsubprofit;
  subprofit = profit->subprofit;
 
  if (profit->nlimmin)
    profit->flag |= PROFLAG_MINLIM;
  if (profit->nlimmax)
    profit->flag |= PROFLAG_MAXLIM;
 
  for (p=0; p<nparam; p++)
    profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);
 
/* Clean up FFTW plans */
  fft_reset();
 
/* CHECK-Images */
  if ((check = prefs.check[CHECK_PROFILES]))
    {
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->iy, 1.0);
    }
 
  if ((check = prefs.check[CHECK_SUBPROFILES]))
    {
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->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, 0.0, param, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->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, 0.0, param, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->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, 0.0, param, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->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, 0.0, param, NULL);
    check_add(check, subprofit->lmodpix,
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
                subprofit->ix,subprofit->iy, -1.0);
    }
 
/* Compute compressed residuals */
  profit_residuals(profit, 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;
  subprofit = profit->subprofit;
  for (s=0; s<nsub; s++, subprofit++)
    {
    if (FLAG(obj2.flux_prof))
      obj2->flux_prof[s] = subprofit->flux;
    if (FLAG(obj2.mag_prof))
      obj2->mag_prof[s] = subprofit->flux>0.0?
                -FDMAG * log(subprofit->flux) + prefs.mag_zeropoint[s] : 99.0f;
    if (FLAG(obj2.fluxerr_prof) || FLAG(obj2.magerr_prof))
      {
      err = 0.0;
      cov = profit->covar;
      index = profit->paramindex;
      list = profit->paramlist;
      for (i=0; i<PARAM_NPARAM; i++)
        if (flux_flag[i]==s+1 && list[i])
          {
          cov = profit->covar + nparam*index[i];
          for (j=0; j<PARAM_NPARAM; j++)
            if (flux_flag[j]==s+1 && list[j])
              err += cov[index[j]];
          }
      if (err<0.0)
        err = 0.0;
      if (FLAG(obj2.fluxerr_prof))
        obj2->fluxerr_prof[s] = sqrt(err);
      if (FLAG(obj2.magerr_prof))
        obj2->magerr_prof[s] = subprofit->flux>0.0?
                FDMAG * sqrt(err) / subprofit->flux : 99.0f;
      }
    }
 
  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->subprofit->ix + *profit->paramlist[PARAM_X]
                        + 1.0;  /* !CHECK */    /* FITS convention */
      obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];
      }
    else
      emx2 = 0.0;
    if (profit->paramlist[PARAM_Y])
      {
      obj2->y_prof = (double)profit->subprofit->iy + *profit->paramlist[PARAM_Y]
                        + 1.0;  /* !CHECK */    /* FITS convention */
      obj2->poserrmy2_prof = emy2 = profit->covar[j*(nparam+1)];
      }
    else
      emy2 = 0.0;
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
      obj2->poserrmxy_prof = emxy = profit->covar[i+j*nparam];
    else
      emxy = 0.0;
 
/*-- Error ellipse parameters */
    if (FLAG(obj2.poserra_prof))
      {
       double   pmx2,pmy2,temp,theta;
 
      if (fabs(temp=emx2-emy2) > 0.0)
        theta = atan2(2.0 * emxy,temp) / 2.0;
      else
        theta = PI/4.0;
 
      temp = sqrt(0.25*temp*temp+ emxy*emxy);
      pmy2 = pmx2 = 0.5*(emx2+emy2);
      pmx2+=temp;
      pmy2-=temp;
 
      obj2->poserra_prof = (float)sqrt(pmx2);
      obj2->poserrb_prof = (float)sqrt(pmy2);
      obj2->poserrtheta_prof = (float)(theta/DEG);
      }
 
    if (FLAG(obj2.poserrcxx_prof))
      {
       double   temp;
 
      obj2->poserrcxx_prof = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
      obj2->poserrcyy_prof = (float)(emx2/temp);
      obj2->poserrcxy_prof = (float)(-2*emxy/temp);
      }
    }
 
/* Equivalent noise area */
  if (FLAG(obj2.prof_noisearea))
    obj2->prof_noisearea = subprofit_noisearea(profit->subprofit); /* !CHECK */
 
/* 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))
    subprofit_convmoments(profit->subprofit, obj2); /* !CHECK */
 
/* "Hybrid" magnitudes */
  if (FLAG(obj2.fluxcor_prof))
    {
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
    subprofit_fluxcor(profit->subprofit, obj2); /* !CHECK */
    }
 
/* Do measurements on the rasterised model (surface brightnesses) */
  if (FLAG(obj2.fluxeff_prof))
    subprofit_surface(profit, profit->subprofit, obj2);  /* !CHECK */
 
/* 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))
    {
    for (s=0; s<nsub; s++)
      {
      flux = obj2->prof_dirac_flux[s] = *profit->paramlist[PARAM_DIRAC_FLUX+s];
      if (FLAG(obj2.prof_dirac_mag))
        obj2->prof_dirac_mag[s] = flux>0.0f?
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
      if (FLAG(obj2.prof_dirac_fluxerr))
        obj2->prof_dirac_fluxerr[s] =
                profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX+s]];
      if (FLAG(obj2.prof_dirac_magerr))
        obj2->prof_dirac_magerr[s] = flux>0.0f?
                 FDMAG*profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX+s]]
                / flux : 99.0f;
      }
    }
 
/* 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;
      }
    for (s=0; s<nsub; s++)
      {
      flux = obj2->prof_spheroid_flux[s]
                = *profit->paramlist[PARAM_SPHEROID_FLUX+s];
      if (FLAG(obj2.prof_spheroid_mag))
        obj2->prof_spheroid_mag[s] = flux>0.0f?
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
      if (FLAG(obj2.prof_spheroid_fluxerr))
        obj2->prof_spheroid_fluxerr[s] =
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX+s]];
      if (FLAG(obj2.prof_spheroid_magerr))
        obj2->prof_spheroid_magerr[s] = flux>0.0f?
              FDMAG*profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX+s]]
                / flux : 99.0f;
      }
    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[0]
                / (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;
      }
    for (s=0; s<nsub; s++)
      {
      flux = obj2->prof_disk_flux[s] = *profit->paramlist[PARAM_DISK_FLUX+s];
      if (FLAG(obj2.prof_disk_mag))
        obj2->prof_disk_mag[s] = flux>0.0f?
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
      if (FLAG(obj2.prof_disk_fluxerr))
        obj2->prof_disk_fluxerr[s] =
                profit->paramerr[profit->paramindex[PARAM_DISK_FLUX+s]];
      if (FLAG(obj2.prof_disk_magerr))
        obj2->prof_disk_magerr[s] = flux>0.0f?
                 FDMAG*profit->paramerr[profit->paramindex[PARAM_DISK_FLUX+s]]
                / flux : 99.0f;
      }
    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[0]
        / (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, 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]];
        }
      }
    }
 
  obj2->prof_flag = profit->flag;
 
  return;
  }
 
 
/****** profit_spread ********************************************************
PROTO   void profit_spread(profitstruct *profit, fieldstruct *field,
                fieldstruct *wfield, obj2struct *obj2)
PURPOSE Perform star/galaxy separation by comparing the best-fitting local PSF
        and a compact exponential profile to the actual data using linear
        discriminant analysis.
INPUT   Pointer to the profile-fitting structure,
        pointer to the field,
        pointer to the field weight,
        pointer to the obj2.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    profit_spread(profitstruct *profit,  fieldstruct *field,
                fieldstruct *wfield, obj2struct *obj2)
  {
   profitstruct         *pprofit,*qprofit;
   subprofitstruct      *subprofit, *psubprofit,*qsubprofit;
   PIXTYPE              valp,valq, sig2;
   double               sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0;
   float                dchi2;
   int                  p,s, nsub;
 
  nsub = profit->nsubprofit;
 
  pprofit = profit_init(obj2, MODEL_DIRAC);
  qprofit = profit_init(obj2, MODEL_EXPONENTIAL);
 
  profit_residuals(profit, PROFIT_DYNPARAM, profit->paraminit, profit->resi);
  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];
    }
 
  subprofit = profit->subprofit;
  for (s=0; s<nsub; s++)
    pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX+s]]
        = (subprofit++)->flux;
 
  pprofit->niter = profit_minimize(pprofit, PROFIT_MAXITER);
 
  profit_residuals(pprofit, PROFIT_DYNPARAM,pprofit->paraminit,pprofit->resi);
 
  qprofit->paraminit[qprofit->paramindex[PARAM_X]]
                = pprofit->paraminit[pprofit->paramindex[PARAM_X]];
  qprofit->paraminit[qprofit->paramindex[PARAM_Y]]
                = pprofit->paraminit[pprofit->paramindex[PARAM_Y]];
  for (s=0; s<nsub; s++)
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_FLUX+s]]
                = pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX+s]];
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_SCALE]]
        = profit->subprofit->psf->fwhm/16.0;    /* !CHECK */
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_ASPECT]] = 1.0;
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_POSANG]] = 0.0;
 
  profit_residuals(qprofit,PROFIT_DYNPARAM, qprofit->paraminit,qprofit->resi);
 
  if (FLAG(obj2.prof_class_star))
    {
    dchi2 = 0.5*(pprofit->chi2 - profit->chi2);
    obj2->prof_class_star = dchi2 < 50.0?
        (dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
    }
 
  psubprofit = pprofit->subprofit;
  qsubprofit = qprofit->subprofit;
  for (s=0; s<profit->nsubprofit; s++)
    {
    sump = sumq = sumpw2 = sumqw2 = sumpqw = sump0 = sumq0 = 0.0;
    for (p=0; p<psubprofit->nobjpix; p++)
      if (psubprofit->objweight[p]>0 && psubprofit->objpix[p]>-BIG)
        {
        valp = psubprofit->lmodpix[p];
        sump += (double)(valp*psubprofit->objpix[p]);
        valq = qsubprofit->lmodpix[p];
        sumq += (double)(valq*psubprofit->objpix[p]);
        sump0 += (double)(valp*valp);
        sumq0 += (double)(valp*valq);
        sig2 = 1.0f/(psubprofit->objweight[p]*psubprofit->objweight[p]);
        sumpw2 += valp*valp*sig2;
        sumqw2 += valq*valq*sig2;
        sumpqw += valp*valq*sig2;
        }
 
    if (FLAG(obj2.prof_concentration))
      {
      obj2->prof_concentration[s] = sump>0.0? (sumq/sump - sumq0/sump0) : 1.0;
      if (FLAG(obj2.prof_concentrationerr))
        obj2->prof_concentrationerr[s] = sump>0.0?
                sqrt(sumqw2*sump*sump+sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)
                        / (sump*sump) : 0.0;
      }
    }
 
  profit_end(pprofit);
  profit_end(qprofit);
 
  return;
  }
 
/****** subprofit_noisearea ***************************************************
PROTO   float subprofit_noisearea(profitstruct *profit)
PURPOSE Return the equivalent noise area (see King 1983) of a model.
INPUT   Sub-profile-fitting structure,
OUTPUT  Equivalent noise area, in pixels.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
float   subprofit_noisearea(subprofitstruct *subprofit)
  {
   double       dval, flux,flux2;
   PIXTYPE      *pix;
   int          p;
 
  flux = flux2 = 0.0;
  pix = subprofit->lmodpix;
  for (p=subprofit->nobjpix; p--;)
    {
    dval = (double)*(pix++);
    flux += dval;
    flux2 += dval*dval;
    }
 
  return (float)(flux2>0.0? flux*flux / flux2 : 0.0);
  }
 
 
/****** subprofit_fluxcor ****************************************************
PROTO   void subprofit_fluxcor(subprofitstruct *subprofit, obj2struct *obj2)
PURPOSE Integrate the flux within an ellipse and complete it with the wings of
        the fitted model.
INPUT   Sub-profile-fitting structure,
        pointer to the obj2 structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 16/03/2012
 ***/
void    subprofit_fluxcor(subprofitstruct *subprofit, 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 = subprofit->objnaxisn[0];
  h = subprofit->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->auto_kronfactor>0.0)
    {
    cx2 = obj2->cxx;
    cy2 = obj2->cyy;
    cxy = obj2->cxy;
    klim2 = 0.64*obj2->auto_kronfactor*obj2->auto_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]))
sexellipse(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 = subprofit->lmodpix;
  objpixt = objpix = subprofit->objpix;
  objweightt = objweight = subprofit->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[0]*tvmout/tvm;
    obj2->fluxcorerr_prof = sqrt(sigtvobj
                +obj2->fluxerr_prof[0]*obj2->fluxerr_prof[0]*tvmout/tvm);
    }
  else
    {
    obj2->fluxcor_prof = tvobj;
    obj2->fluxcorerr_prof = sqrt(sigtvobj);
    }
 
/*
  if ((check = prefs.check[CHECK_OTHER]))
    check_add(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;
  }
 
 
/****** subprofit_psf *********************************************************
PROTO   void    subprofit_psf(subprofitstruct *subprofit, obj2struct *obj2)
PURPOSE Build the local PSF at a given resolution.
INPUT   Pointer to sub-profile-fitting structure,
        pointer to obj2.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    subprofit_psf(subprofitstruct *subprofit, obj2struct *obj2)
  {
   psfstruct    *psf;
   double       flux;
   float        posin[2], posout[2], dnaxisn[2],
                *pixout,
                xcout,ycout, xcin,ycin, invpixstep, norm;
   int          d,i;
 
  psf = subprofit->psf;
  psf_build(psf, obj2);
 
  xcout = (float)(subprofit->modnaxisn[0]/2) + 1.0;      /* FITS convention */
  ycout = (float)(subprofit->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 = subprofit->pixstep / psf->pixstep;
 
/* Initialize multi-dimensional counters */
  for (d=0; d<2; d++)
    {
    posout[d] = 1.0;                                    /* FITS convention */
    dnaxisn[d] = subprofit->modnaxisn[d]+0.5;
    }
 
/* Remap each pixel */
  pixout = subprofit->psfpix;
  flux = 0.0;
  for (i=subprofit->modnaxisn[0]*subprofit->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 *= subprofit->pixstep*subprofit->pixstep;
  if (fabs(flux) <= 0.0)
    error(EXIT_FAILURE, "*Error*: PSF model is empty or negative: ", psf->name);
 
  norm = 1.0/flux;
  pixout = subprofit->psfpix;
  for (i=subprofit->modnaxisn[0]*subprofit->modnaxisn[1]; i--;)
    *(pixout++) *= norm;
 
  return;
  }
 
 
/****** profit_minimize *******************************************************
PROTO   void profit_minimize(profitstruct *profit)
PURPOSE Provide a function returning residuals to lmfit.
INPUT   Pointer to the profit structure involved in the fit,
        maximum number of iterations.
OUTPUT  Number of iterations used.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 20/05/2011
 ***/
int     profit_minimize(profitstruct *profit, int niter)
  {
   double       lm_opts[5], info[LM_INFO_SZ],
                dcovar[PARAM_NPARAM*PARAM_NPARAM], dparam[PARAM_NPARAM];
   int          nfree;
 
  profit->iter = 0;
  memset(dcovar, 0, profit->nparam*profit->nparam*sizeof(double));
 
/* Perform fit */
  lm_opts[0] = 1.0e-3;           /* Initial mu */
  lm_opts[1] = 1.0e-8;          /* ||J^T e||_inf stopping factor */
  lm_opts[2] = 1.0e-8;          /* |Dp||_2 stopping factor */
  lm_opts[3] = 1.0e-8;          /* ||e||_2 stopping factor */
  lm_opts[4] = 1.0e-4;          /* Jacobian step */
 
  nfree = profit_boundtounbound(profit, profit->paraminit, dparam,
                                PARAM_ALLPARAMS);
 
  niter = dlevmar_dif(profit_evaluate, dparam, NULL, nfree, profit->nresi,
                        niter, lm_opts, info, NULL, dcovar, profit);
 
  profit_unboundtobound(profit, dparam, profit->paraminit, PARAM_ALLPARAMS);
 
/* Convert covariance matrix to bounded space */
  profit_covarunboundtobound(profit, dcovar, profit->covar);
 
  return niter;
  }
 
 
/****** profit_printout *******************************************************
PROTO   void profit_printout(int n_par, float* par, int m_dat, float* fvec,
                void *data, int iflag, int iter, int nfev )
PURPOSE Provide a function to print out results to lmfit.
INPUT   Number of fitted parameters,
        pointer to the vector of parameters,
        number of data points,
        pointer to the vector of residuals (output),
        pointer to the data structure (unused),
 
        outer loop counter,
        number of calls to evaluate().
OUTPUT  -.
NOTES   Input arguments are there only for compatibility purposes (unused)
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    profit_printout(int n_par, float* par, int m_dat, float* fvec,
                void *data, int iflag, int iter, int nfev )
  {
   fieldstruct  *field;
   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;
    field = NULL;       /*! Should be replaced with a global variable to work!*/
    sprintf(filename, "check_%d_%04d.fits", the_gal, itero);
    check=check_init(filename, CHECK_PROFILES, 0, 1);
    check_reinit(field, check);
    check_add(check, profit->subprofit->lmodpix,
                profit->subprofit->objnaxisn[0],profit->subprofit->objnaxisn[1],
                profit->subprofit->ix,profit->subprofit->iy, 1.0);
 
    check_reend(field, check);
    check_end(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 26/07/2011
 ***/
void    profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
  {
   profitstruct         *profit;
   profstruct           **prof;
   double               *dpar0, *dresi;
   float                *modpixt, *profpixt, *resi,
                        tparam, val;
   PIXTYPE              *lmodpixt,*lmodpix2t, *objpix,*weight,
                        wval;
   int                  c,f,i,p,q, fd,pd, jflag,sflag, nprof;
 
  profit = (profitstruct *)adata;
 
/* Detect "Jacobian-related" calls */
  jflag = pd = fd = 0;
  dpar0 = profit->dparam;
  if (profit->iter)
    {
    f = q = 0;
    for (p=0; p<profit->nparam; p++)
      {
      if (dpar[f] - dpar0[f] != 0.0)
        {
        pd = p;
        fd = f;
        q++;
        }
      if (profit->parfittype[p]!=PARFIT_FIXED)
        f++;
      }
    if (f>0 && q==1)
      jflag = 1;
    }
jflag = 0;       /* Temporarily deactivated (until problems are fixed) */
/*
  if (jflag && !(profit->nprof==1 && profit->prof[0]->code == MODEL_DIRAC))
    {
    prof = profit->prof;
    nprof = profit->nprof;
 
*-- "Jacobian call" *
    tparam = profit->param[pd];
    profit_unboundtobound(profit, &dpar[fd], &profit->param[pd], pd);
    sflag = 1;
    switch(profit->paramrevindex[pd])
      {
      case PARAM_BACK:
        lmodpixt = profit->lmodpix;
        lmodpix2t = profit->lmodpix2;
        val = (profit->param[pd] - tparam);
        for (i=profit->nobjpix;i--;)
          *(lmodpix2t++) = val;
        break;
      case PARAM_X:
      case PARAM_Y:
        profit_resample(profit, profit->cmodpix, profit->lmodpix2, 1.0);
        lmodpixt = profit->lmodpix;
        lmodpix2t = profit->lmodpix2;
        for (i=profit->nobjpix;i--;)
          *(lmodpix2t++) -= *(lmodpixt++);
        break;
      case PARAM_DIRAC_FLUX:
      case PARAM_SPHEROID_FLUX:
      case PARAM_DISK_FLUX:
      case PARAM_ARMS_FLUX:
      case PARAM_BAR_FLUX:
        if (nprof==1 && tparam != 0.0)
          {
          lmodpixt = profit->lmodpix;
          lmodpix2t = profit->lmodpix2;
          val = (profit->param[pd] - tparam) / tparam;
          for (i=profit->nobjpix;i--;)
            *(lmodpix2t++) = val**(lmodpixt++);
          }
        else
          {
          for (c=0; c<nprof; c++)
            if (prof[c]->flux[0] == &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[0];
        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, PROFIT_DYNPARAM, profit->param, profit->resi);
 
    for (p=0; p<profit->nresi; p++)
      fvec[p] = profit->resi[p];
    }
 
//  profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
  profit->iter++;
 
  return;
  }
 
 
/****** profit_residuals ******************************************************
PROTO   float *profit_residuals(profitstruct *profit, float dynparam,
                        float *param, 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),
        pointer to the model parameters (output),
        pointer to the computed residuals (output).
OUTPUT  Vector of residuals.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
float   *profit_residuals(profitstruct *profit, float dynparam, float *param,
                float *resi)
  {
   subprofitstruct      *subprofit;
   float                *pixin, *pixout,
                        pflux;
   int                  i,p,s;
 
  subprofit = profit->subprofit;
  for (s=0; s<profit->nsubprofit; s++, subprofit++)
    memset(subprofit->modpix, 0, subprofit->nmodpix*sizeof(float));
  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)
    {
    subprofit = profit->subprofit;
    for (s=0; s<profit->nsubprofit; s++, subprofit++)
      {
      profit_resample(profit, subprofit, subprofit->psfpix,
                subprofit->lmodpix, *profit->prof[0]->flux[s]);
      subprofit->flux = *profit->prof[0]->flux[s];
      }
    }
  else
    {
    subprofit = profit->subprofit;
    for (s=0; s<profit->nsubprofit; s++, subprofit++)
      {
      subprofit->flux = 0.0;
      for (p=0; p<profit->nprof; p++)
if(!s)
        subprofit->flux += prof_add(subprofit, profit->prof[p], 0);
else
{
pflux = *profit->prof[p]->flux[s];
pixin = profit->prof[p]->pix;
pixout = subprofit->modpix;
for (i=subprofit->nmodpix; i--;)
  *(pixout++) += pflux**(pixin++);
subprofit->flux += pflux;
        }
      memcpy(subprofit->cmodpix, subprofit->modpix,
                subprofit->nmodpix*sizeof(float));
      subprofit_convolve(subprofit, subprofit->cmodpix);
      profit_resample(profit, subprofit,
        subprofit->cmodpix, subprofit->lmodpix, 1.0);
      }
    }
 
  if (resi)
    profit_compresi(profit, dynparam, resi);
 
  return resi;
  }
 
 
/****** 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 06/05/2012
 ***/
float   *profit_compresi(profitstruct *profit, float dynparam, float *resi)
  {
   subprofitstruct      *subprofit;
   double               error;
   float                *resit;
   PIXTYPE              *objpix, *objweight, *lmodpix,
                        val,val2,wval, invsig;
   int                  i,s;
 
/* Compute vector of residuals */
  resit = resi;
  error = 0.0;
  subprofit = profit->subprofit;
  if (dynparam > 0.0)
    {
    invsig = (PIXTYPE)(1.0/dynparam);
    for (s=profit->nsubprofit; s--; subprofit++)
      {
      objpix = subprofit->objpix;
      objweight = subprofit->objweight;
      lmodpix = subprofit->lmodpix;
      for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; 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 (s=profit->nsubprofit; s--; subprofit++)
      {
      objpix = subprofit->objpix;
      objweight = subprofit->objweight;
      lmodpix = subprofit->lmodpix;
      for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; 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, subprofitstruct *subprofit,
                float *inpix, PIXTYPE *outpix, float factor)
PURPOSE Resample the current full resolution model to image resolution.
INPUT   Profile-fitting structure,
        sub-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 05/05/2012
 ***/
int     profit_resample(profitstruct *profit, subprofitstruct *subprofit,
                float *inpix, PIXTYPE *outpix, float factor)
  {
   PIXTYPE      *pixout,*pixout0;
   float        *pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0,
                *dpixout,*dpixout0, *dx,*dy,
                xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val,
                invpixstep, norm;
   int          *start,*startt, *nmask,*nmaskt,
                i,j,k,n,t,
                ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout,
                iysina, nyin, hmw,hmh, ix,iy, ixin,iyin;
 
  invpixstep = subprofit->subsamp/subprofit->pixstep;
 
  xcin = (subprofit->modnaxisn[0]/2);
  xcout = (float)(subprofit->objnaxisn[0]/2);
  if ((dx=(profit->paramlist[PARAM_X])))
    xcout += *dx;
 
  xsin = xcin - xcout*invpixstep;                       /* Input start x-coord*/
  if ((int)xsin >= subprofit->modnaxisn[0])
    return RETURN_ERROR;
  ixsout = 0;                            /* Int. part of output start x-coord */
  if (xsin<0.0)
    {
    dixout = (int)(1.0-xsin/invpixstep);
/*-- Simply leave here if the images do not overlap in x */
    if (dixout >= subprofit->objnaxisn[0])
      return RETURN_ERROR;
    ixsout += dixout;
    xsin += dixout*invpixstep;
    }
  nxout = (int)((subprofit->modnaxisn[0]-xsin)/invpixstep);/* nb of interpolated
                                                        input pixels along x */
  if (nxout>(ixout=subprofit->objnaxisn[0]-ixsout))
    nxout = ixout;
  if (!nxout)
    return RETURN_ERROR;
 
  ycin = (subprofit->modnaxisn[1]/2);
  ycout = (float)(subprofit->objnaxisn[1]/2);
  if ((dy=(profit->paramlist[PARAM_Y])))
    ycout += *dy;
 
  ysin = ycin - ycout*invpixstep;               /* Input start y-coord*/
  if ((int)ysin >= subprofit->modnaxisn[1])
    return RETURN_ERROR;
  iysout = 0;                            /* Int. part of output start y-coord */
  if (ysin<0.0)
    {
    diyout = (int)(1.0-ysin/invpixstep);
/*-- Simply leave here if the images do not overlap in y */
    if (diyout >= subprofit->objnaxisn[1])
      return RETURN_ERROR;
    iysout += diyout;
    ysin += diyout*invpixstep;
    }
  nyout = (int)((subprofit->modnaxisn[1]-ysin)/invpixstep);/* nb of interpolated
                                                        input pixels along y */
  if (nyout>(iyout=subprofit->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>subprofit->modnaxisn[1])                                             /* with margin */
    nyin = subprofit->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=subprofit->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*subprofit->modnaxisn[0];
  dpixout0 = pixinout;
  for (k=nyin; k--; pixin0+=subprofit->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)subprofit->nobjpix*sizeof(PIXTYPE));
 
/* Make the interpolation in y and transpose once again */
  dpixin0 = pixinout;
  pixout0 = outpix+ixsout+iysout*subprofit->objnaxisn[0];
  for (k=nxout; k--; dpixin0+=nyin, pixout0++)
    {
    maskt = mask;
    nmaskt = nmask;
    startt = start;
    pixout = pixout0;
    for (j=nyout; j--; pixout+=subprofit->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;
  }
 
 
/****** subprofit_convolve ****************************************************
PROTO   void subprofit_convolve(subprofitstruct *subprofit, float *modpix)
PURPOSE Convolve a model image with the local PSF.
INPUT   Pointer to the subprofit structure,
        Pointer to the image raster.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 05/05/2012
 ***/
void    subprofit_convolve(subprofitstruct *subprofit, float *modpix)
  {
  if (!subprofit->psfdft)
    subprofit_makedft(subprofit);
 
  fft_conv(modpix, subprofit->psfdft, subprofit->modnaxisn);
 
  return;
  }
 
 
/****** subprofit_makedft ****************************************************
PROTO   void subprofit_makedft(subprofitstruct *subprofit)
PURPOSE Create the Fourier transform of the descrambled PSF component.
INPUT   Pointer to the subprofit structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
void    subprofit_makedft(subprofitstruct *subprofit)
  {
   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=subprofit->psf))
    return;
 
  psfwidth = subprofit->modnaxisn[0];
  psfheight = subprofit->modnaxisn[1];
  psfnpix = psfwidth*psfheight;
  width = subprofit->modnaxisn[0];
  height = subprofit->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 = subprofit->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 = subprofit->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/subprofit->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 */
  subprofit->psfdft = fft_rtf(mask, subprofit->modnaxisn);
 
  free(mask);
 
  return;
  }
 
 
/****** subprofit_copyobjpix **************************************************
PROTO   int subprofit_copyobjpix(subprofitstruct *subprofit,
                                subimagestruct *subimage)
PURPOSE Copy a piece of the input object image to a subprofit structure.
INPUT   Pointer to the subprofit structure,
        subimagestruct *subimage.
OUTPUT  The number of valid pixels copied.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 06/05/2012
 ***/
int     subprofit_copyobjpix(subprofitstruct *subprofit,
                                subimagestruct *subimage)
  {
   fieldstruct          *field, *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, ix,iy, win,hin;
 
  field = subimage->field;
  wfield = subimage->wfield;
 
/* First put the image background to -BIG */
  pixout = subprofit->objpix;
  wpixout = subprofit->objweight;
  for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; i--;)
    {
    *(pixout++) = -BIG;
    *(wpixout++) = 0.0;
    }
 
  ix = subprofit->ix - subimage->immin[0];
  iy = subprofit->iy - subimage->immin[1];
  win = subimage->imsize[0];
  hin = subimage->imsize[1];
 
  backnoise2 = field->backsig*field->backsig;
  sn = (int)subprofit->subsamp;
  sflag = (sn>1);
  w = subprofit->objnaxisn[0]*sn;
  h = subprofit->objnaxisn[1]*sn;
  if (sflag)
    backnoise2 *= (PIXTYPE)sn;
  invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0;
  satlevel = field->satur_level - subimage->bkg;
  rad2 = h/2.0;
  if (rad2 > w/2.0)
    rad2 = w/2.0;
  rad2 *= rad2;
 
/* Set the image boundaries */
  pixout = subprofit->objpix;
  wpixout = subprofit->objweight;
  ymin = iy-h/2;
  ymax = ymin + h;
  if (ymin<0)
    {
    off = (-ymin-1)/sn + 1;
    pixout += off*subprofit->objnaxisn[0];
    wpixout += off*subprofit->objnaxisn[0];
    ymin += off*sn;
    }
  if (ymax>hin)
    ymax -= ((ymax-hin-1)/sn + 1)*sn;
 
  xmin = ix-w/2;
  xmax = xmin + w;
  dw = 0;
  if (xmax>win)
    {
    off = (xmax-win-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;
    }
 
/* Copy the right pixels to the destination */
  npix = 0;
  if (wfield)
    {
    wthresh = wfield->weight_thresh;
    gainflag = field->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);
            off = x + (y+sy)*win;
            spixin = subimage->image + off;
            swpixin = subimage->weight + off;
            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 = subimage->image + xmin + y*win;
        wpixin = subimage->weight + xmin + y*win;
        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 = subimage->image + x + (y+sy)*win;
            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 = subimage->image + xmin + y*win;
        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;
  }
 
 
/****** subprofit_submodpix *****************************************************
PROTO   void    subprofit_submodpix(subprofitstruct *subprofitobj,
                        subprofitstruct *subprofitmod, float fac)
PURPOSE Subtract a rasterized model from the objpix of another model structure.
INPUT   Pointer to the sub-profile structure containing the object pixels,
        Pointer to the sub-profile structure containing the model raster,
        multiplicative factor to apply to the model pixels.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 05/05/2012
 ***/
void    subprofit_submodpix(subprofitstruct *subprofitobj,
                        subprofitstruct *subprofitmod, float fac)
  {
   float        dx, dy2, dr2, rad2;
   PIXTYPE      *pixin,*spixin, *pixout,
                pix,spix;
   int          i,x,y, xmin,xmax,ymin,ymax, w,h,dw, npix, off, gainflag,
                badflag, sflag, sx,sy,sn, ix,iy, win,hin;
 
/* Don't go further if out of frame!! */
  win = subprofitmod->objnaxisn[0];
  hin = subprofitmod->objnaxisn[1];
  ix = subprofitobj->ix - (subprofitmod->ix - win/2);
  iy = subprofitobj->iy - (subprofitmod->iy - hin/2);
 
  sn = (int)subprofitobj->subsamp;
  sflag = (sn>1);
  w = subprofitobj->objnaxisn[0]*sn;
  h = subprofitobj->objnaxisn[1]*sn;
  rad2 = h/2.0;
  if (rad2 > w/2.0)
    rad2 = w/2.0;
  rad2 *= rad2;
 
/* Set the image boundaries */
  pixout = subprofitobj->objpix;
  ymin = iy-h/2;
  ymax = ymin + h;
  if (ymin<0)
    {
    off = (-ymin-1)/sn + 1;
    pixout += off*subprofitobj->objnaxisn[0];
    ymin += off*sn;
    }
  if (ymax>hin)
    ymax -= ((ymax-hin-1)/sn + 1)*sn;
 
  xmin = ix-w/2;
  xmax = xmin + w;
  dw = 0;
  if (xmax>win)
    {
    off = (xmax-win-1)/sn + 1;
    dw += off;
    xmax -= off*sn;
    }
  if (xmin<0)
    {
    off = (-xmin-1)/sn + 1;
    pixout += off;
    dw += off;
    xmin += off*sn;
    }
 
/* Copy the right pixels to the destination */
  npix = 0;
  if (sflag)
    {
/*-- Sub-sampling case */
    for (y=ymin; y<ymax; y+=sn, pixout+=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 = subprofitmod->lmodpix + x + (y+sy)*win;
          for (sx=sn; sx--;)
            {
            dr2 = dy2 + dx*dx;
            dx++;
            spix = *(spixin++);
            pix += spix;
            }
          }
        *(pixout++) -= fac*pix;
        }
      }
    }
  else
    for (y=ymin; y<ymax; y++, pixout+=dw)
      {
      dy2 = y-iy;
      dy2 *= dy2;
      pixin = subprofitmod->lmodpix + xmin + y*win;
      for (x=xmin; x<xmax; x++)
        {
        dx = x-ix;
        dr2 = dy2 + dx*dx;
        *(pixout++) -= fac**(pixin++);
        }
      }
 
  return;
  }
 
 
/****** subprofit_spiralindex *************************************************
PROTO   float subprofit_spiralindex(subprofitstruct *subprofit)
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 06/05/2012
 ***/
float subprofit_spiralindex(subprofitstruct *subprofit, 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 = subprofit->objnaxisn[0]*subprofit->objnaxisn[1];
 
/* Compute simple derivative vectors at a fraction of the object scale */
  fwhm = subprofit->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)(subprofit->objnaxisn[0]/2);
  ohw = subprofit->objnaxisn[0] - hw;
  hh = (float)(subprofit->objnaxisn[1]/2);
  ohh = subprofit->objnaxisn[1] - hh;
  txstart = -hw;
  ty = -hh;
  QMALLOC(dx, float, npix);
  pix = dx;
  for (j=subprofit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=subprofit->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=subprofit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=subprofit->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 = subprofit->objpix;
  invsig = npix/subprofit->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, subprofit->objnaxisn);
  fft_conv(gdx, fdx, subprofit->objnaxisn);
  fdy = fft_rtf(dy, subprofit->objnaxisn);
  fft_conv(gdy, fdy, subprofit->objnaxisn);
 
/* Compute estimator */
  invtwosigma2 = -1.18*1.18/(2.0*subprofit->guessradius*subprofit->guessradius);
  xstart = -hw - obj2->mx + (int)(obj2->mx+0.49999);
  y = -hh -  obj2->my + (int)(obj2->my+0.49999);;
  spirindex = 0.0;
  gdxt = gdx;
  gdyt = gdy;
  for (j=subprofit->objnaxisn[1]; j--; y+=1.0)
    {
    x = xstart;
    for (i=subprofit->objnaxisn[0]; i--; x+=1.0)
      {
      gx = *(gdxt++);
      gy = *(gdyt++);
      if ((r2=x*x+y*y)>0.0)
        spirindex += (x*y*(gx*gx-gy*gy)+gx*gy*(y*y-x*x))/r2
                        * exp(invtwosigma2*r2);
      }
    }
 
  free(dx);
  free(dy);
  free(fdx);
  free(fdy);
  free(gdx);
  free(gdy);
 
  return spirindex;
  }
 
 
/****** profit_moments *******************************************************
PROTO   void profit_moments(profitstruct *profit, obj2struct *obj2)
PURPOSE Compute the 2nd order moments from the unconvolved object model.
INPUT   Profile-fitting structure,
        Pointer to obj2 structure.
OUTPUT  -.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 22/04/2011
 ***/
void     profit_moments(profitstruct *profit, obj2struct *obj2)
  {
   profstruct   *prof;
   double       dpdmx2[6], cov[4],
                *jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
                *dmx2,*dmy2,*dmxy,
                m0,invm0, mx2,my2,mxy, den,invden,
                temp, temp2,invtemp2,invstemp2,
                pmx2,theta, flux, dval;
   float         *covart;
   int          findex[MODEL_NMAX],
                i,j,p, nparam;
 
/*  hw = (float)(profit->modnaxisn[0]/2);*/
/*  hh = (float)(profit->modnaxisn[1]/2);*/
/*  r2max = hw<hh? hw*hw : hh*hh;*/
/*  xstart = -hw;*/
/*  y = -hh;*/
/*  pix = profit->modpix;*/
/*  mx2 = my2 = mxy = mx = my = sum = 0.0;*/
/*  for (iy=profit->modnaxisn[1]; iy--; y+=1.0)*/
/*    {*/
/*    x = xstart;*/
/*    for (ix=profit->modnaxisn[0]; ix--; x+=1.0)*/
/*      if (y*y+x*x <= r2max)*/
/*        {*/
/*        val = *(pix++);*/
/*        sum += val;*/
/*        mx  += val*x;*/
/*        my  += val*y;*/
/*        mx2 += val*x*x;*/
/*        mxy += val*x*y;*/
/*        my2 += val*y*y;*/
/*        }*/
/*      else*/
/*        pix++;*/
/*    }*/
 
/*  if (sum <= 1.0/BIG)*/
/*    sum = 1.0;*/
/*  mx /= sum;*/
/*  my /= sum;*/
/*  obj2->prof_mx2 = mx2 = mx2/sum - mx*mx;*/
/*  obj2->prof_my2 = my2 = my2/sum - my*my;*/
/*  obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
 
  nparam = profit->nparam;
  if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
    {
/*-- Set up Jacobian matrices */
    QCALLOC(jac, double, nparam*3);
    QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
    QMALLOC(dcovar, double, nparam*nparam);
    dcovart = dcovar;
    covart = profit->covar;
    for (i=nparam*nparam; i--;)
      *(dcovart++) = (double)(*(covart++));
    dmx2 = jac;
    dmy2 = jac+nparam;
    dmxy = jac+2*nparam;
    }
  else
    jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;
 
  m0 = mx2 = my2 = mxy = 0.0;
  for (p=0; p<profit->nprof; p++)
    {
    prof = profit->prof[p];
    findex[p] = prof_moments(profit, prof, pjac);
    flux = *prof->flux[0];
    m0 += flux;
    mx2 += prof->mx2*flux;
    my2 += prof->my2*flux;
    mxy += prof->mxy*flux;
    if (jac)
      {
      jact = jac;
      pjact = pjac;
      for (j=nparam*3; j--;)
        *(jact++) += flux * *(pjact++);
      }
    }
  invm0 = 1.0 / m0;
  obj2->prof_mx2 = (mx2 *= invm0);
  obj2->prof_my2 = (my2 *= invm0);
  obj2->prof_mxy = (mxy *= invm0);
/* Complete the flux derivative of moments */
  if (jac)
    {
    for (p=0; p<profit->nprof; p++)
      {
      prof = profit->prof[p];
      dmx2[findex[p]] = prof->mx2 - mx2;
      dmy2[findex[p]] = prof->my2 - my2;
      dmxy[findex[p]] = prof->mxy - mxy;
      }
    jact = jac;
    for (j=nparam*3; j--;)
      *(jact++) *= invm0;
    }
 
/* Handle fully correlated profiles (which cause a singularity...) */
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
    {
    mx2 += 0.0833333;
    my2 += 0.0833333;
    temp2 = mx2*my2-mxy*mxy;
    }
 
/* Use the Jacobians to compute the moment covariance matrix */
  if (jac)
    propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
                                                pjac);  /* We re-use pjac */
 
  if (FLAG(obj2.prof_pol1))
    {
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
    if (mx2+my2 > 1.0/BIG)
      {
      obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
      obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
      if (FLAG(obj2.prof_pol1err))
        {
/*------ Compute the Jacobian of polarisation */
        invden = 1.0/(mx2+my2);
        dpdmx2[0] =  2.0*my2*invden*invden;
        dpdmx2

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	278	bertin	`* Copyright: (C) 2006-2012 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	285	bertin	`* Last modified: 06/05/2012`
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	274	bertin	`#include "graph.h"`
51	229	bertin	`#include "image.h"`
52	274	bertin	`#include "misc.h"`
53	141	bertin	`#include "pattern.h"`
54	50	bertin	`#include "psf.h"`
55	42	bertin	`#include "profit.h"`
56	283	bertin	`#include "subimage.h"`
57	42	bertin
58	207	bertin	`static double prof_gammainc(double x, double a),`
59			`prof_gamma(double x);`
60	201	bertin	`static float prof_interpolate(profstruct prof, float posin);`
61			`static float interpolate_pix(float posin, float pix, int *naxisn,`
62	52	bertin	`interpenum interptype);`
63	50	bertin
64	201	bertin	`static void make_kernel(float pos, float *kernel, interpenum interptype);`
65	50	bertin
66	42	bertin	`/------------------------------- variables ---------------------------------/`
67
68	209	bertin	`const int interp_kernwidth[5]={1,2,4,6,8};`
69
70	233	bertin	`const int flux_flag[PARAM_NPARAM] = {0,`
71	285	bertin	`1,2,3,4,5,0,0,`
72			`1,2,3,4,5,0,0,0,0,`
73			`1,2,3,4,5,0,0,0,`
74	209	bertin	`1,0,0,0,0,0,0,0,`
75			`1,0,0,`
76			`1,0,0,`
77			`1,0,0`
78			`};`
79
80	267	bertin	`int the_gal;`
81
82	46	bertin	`/**** profit_init *********************************************************`
83	285	bertin	`PROTO profitstruct profit_init(obj2struct *obj2, unsigned int modeltype)`
84	46	bertin	`PURPOSE Allocate and initialize a new profile-fitting structure.`
85	285	bertin	`INPUT pointer to the obj2,`
86	258	bertin	`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	285	bertin	`VERSION 06/05/2012`
91	42	bertin	`***/`
92	285	bertin	`profitstruct profit_init(obj2struct obj2, unsigned int modeltype)`
93	42	bertin	`{`
94	46	bertin	`profitstruct *profit;`
95	285	bertin	`subprofitstruct *subprofit;`
96			`subimagestruct *subimage;`
97	267	bertin	`double psf_fwhm;`
98	285	bertin	`int s,t, nmodels, npix, nsub;`
99	42	bertin
100	49	bertin	`QCALLOC(profit, profitstruct, 1);`
101	285	bertin
102	258	bertin	`profit->obj2 = obj2;`
103
104	285	bertin	`QCALLOC(profit->subprofit, subprofitstruct, obj2->nsubimage);`
105			`subprofit = profit->subprofit;`
106			`subimage = obj2->subimage;`
107			`nsub = 0;`
108			`for (s=0; s<obj2->nsubimage; s++, subimage++)`
109			`{`
110			`subprofit->field = subimage->field;`
111			`subprofit->wfield = subimage->wfield;`
112			`subprofit->psf = subimage->field->psf;`
113			`subprofit->pixstep = subprofit->psf->pixstep;`
114			`subprofit->sigma = obj2->sigbkg[s];`
115			`/-- Set initial guesses and boundaries /`
116			`if ((subprofit->guessradius = 0.5*subprofit->psf->fwhm) < obj2->hl_radius)`
117			`subprofit->guessradius = obj2->hl_radius;`
118			`if ((subprofit->guessflux = obj2->flux_auto[s]) <= 0.0)`
119			`subprofit->guessflux = 0.0;`
120			`if ((subprofit->guessfluxmax = 10.0*obj2->fluxerr_auto[s])`
121			`<= subprofit->guessflux)`
122			`subprofit->guessfluxmax = subprofit->guessflux;`
123			`if (subprofit->guessfluxmax <= 0.0)`
124			`subprofit->guessfluxmax = 1.0;`
125			`subprofit->fluxfac = 1.0; /* Default */`
126			`/-- Create pixmaps at image resolution /`
127			`subprofit->ix = (int)(obj2->mx + 0.49999); /* 1st pix=0 */`
128			`subprofit->iy = (int)(obj2->my + 0.49999); /* 1st pix=0 */`
129			`psf_fwhm = subprofit->psf->masksize[0]*subprofit->psf->pixstep;`
130			`subprofit->objnaxisn[0] = ((int)((subimage->imsize[0] + psf_fwhm + 0.499)`
131	258	bertin	`1.2)/2)2 + 1;`
132	285	bertin	`subprofit->objnaxisn[1] = ((int)((subimage->imsize[1] + psf_fwhm + 0.499)`
133	258	bertin	`1.2)/2)2 + 1;`
134	285	bertin	`if (subprofit->objnaxisn[1] < subprofit->objnaxisn[0])`
135			`subprofit->objnaxisn[1] = subprofit->objnaxisn[0];`
136			`else`
137			`subprofit->objnaxisn[0] = subprofit->objnaxisn[1];`
138	258	bertin
139	285	bertin	`if (subprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)`
140			`{`
141			`subprofit->subsamp`
142			`= ceil((float)subprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);`
143			`subprofit->objnaxisn[1]`
144			`= (subprofit->objnaxisn[0] /= (int)subprofit->subsamp);`
145			`profit->flag \|= PROFLAG_OBJSUB; /* !CHECK */`
146			`}`
147			`else`
148			`subprofit->subsamp = 1.0;`
149			`subprofit->nobjpix = subprofit->objnaxisn[0]*subprofit->objnaxisn[1];`
150	258	bertin
151			`/* Allocate memory for the data and the resampled model */`
152	285	bertin	`QMALLOC16(subprofit->objpix, PIXTYPE, subprofit->nobjpix);`
153			`QMALLOC16(subprofit->objweight, PIXTYPE, subprofit->nobjpix);`
154			`QMALLOC16(subprofit->lmodpix, PIXTYPE, subprofit->nobjpix);`
155			`QMALLOC16(subprofit->lmodpix2, PIXTYPE, subprofit->nobjpix);`
156	258	bertin
157	285	bertin	`/-- Create pixmap at PSF resolution /`
158			`subprofit->modnaxisn[0] =`
159			`((int)(subprofit->objnaxisn[0]*subprofit->subsamp`
160			`/subprofit->pixstep+0.4999)/2+1)*2;`
161			`subprofit->modnaxisn[1] =`
162			`((int)(subprofit->objnaxisn[1]*subprofit->subsamp`
163			`/subprofit->pixstep+0.4999)/2+1)*2;`
164			`if (subprofit->modnaxisn[1] < subprofit->modnaxisn[0])`
165			`subprofit->modnaxisn[1] = subprofit->modnaxisn[0];`
166			`else`
167			`subprofit->modnaxisn[0] = subprofit->modnaxisn[1];`
168			`if (subprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)`
169			`{`
170			`subprofit->pixstep = (double)subprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;`
171			`subprofit->modnaxisn[0] = subprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;`
172			`profit->flag \|= PROFLAG_MODSUB;`
173			`}`
174			`subprofit->nmodpix = subprofit->modnaxisn[0]*subprofit->modnaxisn[1];`
175	258	bertin
176	285	bertin	`/* Allocate memory for the complete model */`
177			`QMALLOC16(subprofit->modpix, float, subprofit->nmodpix);`
178			`QMALLOC16(subprofit->modpix2, float, subprofit->nmodpix);`
179			`QMALLOC16(subprofit->cmodpix, float, subprofit->nmodpix);`
180			`QMALLOC16(subprofit->psfpix, float, subprofit->nmodpix);`
181			`if ((npix = subprofit_copyobjpix(subprofit, subimage)))`
182			`{`
183			`profit->nresi += npix;`
184			`/---- Compute the local PSF /`
185			`subprofit_psf(subprofit, obj2);`
186			`subprofit->index = nsub++;`
187			`subprofit++;`
188			`}`
189			`else`
190			`subprofit_end(subprofit);`
191	258	bertin	`}`
192
193	285	bertin	`profit->nsubprofit = nsub;`
194	258	bertin
195	285	bertin	`QMALLOC16(profit->resi, float, profit->nresi);`
196	258	bertin
197	245	bertin	`QMALLOC(profit->prof, profstruct *, MODEL_NMAX);`
198			`nmodels = 0;`
199			`for (t=1; t<(1<<MODEL_NMAX); t<<=1)`
200			`if (modeltype&t)`
201			`profit->prof[nmodels++] = prof_init(profit, t);`
202			`profit->nprof = nmodels;`
203	50	bertin
204	258	bertin	`QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);`
205
206			`profit_resetparams(profit);`
207
208	46	bertin	`return profit;`
209	42	bertin	`}`
210
211
212	46	bertin	`/**** profit_end **********************************************************`
213	285	bertin	`PROTO void profit_end(profitstruct *profit)`
214	46	bertin	`PURPOSE End (deallocate) a profile-fitting structure.`
215	42	bertin	`INPUT Prof structure.`
216			`OUTPUT -.`
217	117	bertin	`NOTES -.`
218	42	bertin	`AUTHOR E. Bertin (IAP)`
219	285	bertin	`VERSION 06/05/2012`
220	42	bertin	`***/`
221	46	bertin	`void profit_end(profitstruct *profit)`
222	42	bertin	`{`
223	285	bertin	`subprofitstruct *subprofit;`
224			`int p,s;`
225	49	bertin
226			`for (p=0; p<profit->nprof; p++)`
227			`prof_end(profit->prof[p]);`
228	270	bertin	`free(profit->prof);`
229	285	bertin	`subprofit = profit->subprofit;`
230			`for (s=profit->nsubprofit; s--;)`
231			`subprofit_end(subprofit++);`
232
233			`free(profit->subprofit);`
234	221	bertin	`free(profit->resi);`
235	123	bertin	`free(profit->covar);`
236	285	bertin
237	46	bertin	`free(profit);`
238	42	bertin
239			`return;`
240			`}`
241
242
243	285	bertin	`/**** subprofit_end *******************************************************`
244			`PROTO void subprofit_end(subprofitstruct *profit)`
245			`PURPOSE End (free content) a subprofile-fitting structure.`
246			`INPUT SubProfit structure.`
247			`OUTPUT -.`
248			`NOTES -.`
249			`AUTHOR E. Bertin (IAP)`
250			`VERSION 03/05/2012`
251			`***/`
252			`void subprofit_end(subprofitstruct *subprofit)`
253			`{`
254			`free(subprofit->objpix);`
255			`free(subprofit->objweight);`
256			`free(subprofit->lmodpix);`
257			`free(subprofit->lmodpix2);`
258			`free(subprofit->modpix);`
259			`free(subprofit->modpix2);`
260			`free(subprofit->cmodpix);`
261			`free(subprofit->psfpix);`
262			`free(subprofit->psfdft);`
263
264			`return;`
265			`}`
266
267
268	65	bertin	`/**** profit_fit **********************************************************`
269	259	bertin	`PROTO int profit_fit(profitstruct profit, obj2struct obj2)`
270	46	bertin	`PURPOSE Fit profile(s) convolved with the PSF to a detected object.`
271	259	bertin	`INPUT Pointer to the model structure,`
272			`pointer to the obj2.`
273			`OUTPUT Number of minimization iterations (0 in case of error).`
274			`NOTES -.`
275	46	bertin	`AUTHOR E. Bertin (IAP)`
276	267	bertin	`VERSION 06/10/2011`
277	46	bertin	`***/`
278	259	bertin	`int profit_fit(profitstruct profit, obj2struct obj2)`
279	46	bertin	`{`
280	259	bertin	`int p, nparam, nparam2;`
281	46	bertin
282	128	bertin	`nparam = profit->nparam;`
283	208	bertin	`nparam2 = nparam*nparam;`
284	52	bertin
285	258	bertin	`/* reset all flags except those set in profit_init() */`
286	259	bertin	`profit->flag &= PROFLAG_OBJSUB\|PROFLAG_MODSUB;`
287	52	bertin
288	211	bertin	`/* Check if the number of constraints exceeds the number of free parameters */`
289	128	bertin	`if (profit->nresi < nparam)`
290	49	bertin	`{`
291	176	bertin	`if (FLAG(obj2.prof_vector))`
292			`for (p=0; p<nparam; p++)`
293			`obj2->prof_vector[p] = 0.0;`
294	208	bertin	`if (FLAG(obj2.prof_errvector))`
295			`for (p=0; p<nparam; p++)`
296			`obj2->prof_errvector[p] = 0.0;`
297			`if (FLAG(obj2.prof_errmatrix))`
298			`for (p=0; p<nparam2; p++)`
299			`obj2->prof_errmatrix[p] = 0.0;`
300	50	bertin	`obj2->prof_niter = 0;`
301	259	bertin	`profit->flag \|= PROFLAG_NOTCONST;`
302			`return 0;`
303	49	bertin	`}`
304	48	bertin
305	123	bertin	`/* Actual minimisation */`
306	209	bertin	`fft_reset();`
307	259	bertin	`return (profit->niter = profit_minimize(profit, PROFIT_MAXITER));`
308			`}`
309	235	bertin
310	259	bertin
311			`/**** profit_measure ******************************************************`
312			`PROTO void profit_measure(profitstruct profit, obj2struct obj2)`
313			`PURPOSE Perform measurements on the fitted model.`
314			`INPUT Pointer to the model structure,`
315			`pointer to the obj2.`
316			`OUTPUT -.`
317			`NOTES -.`
318			`AUTHOR E. Bertin (IAP)`
319	285	bertin	`VERSION 06/05/2012`
320	259	bertin	`***/`
321	267	bertin	`void profit_measure(profitstruct profit, obj2struct obj2)`
322	259	bertin	`{`
323			`profitstruct pprofit, qprofit;`
324	285	bertin	`subprofitstruct *subprofit;`
325	259	bertin	`patternstruct *pattern;`
326			`checkstruct *check;`
327			`double emx2,emy2,emxy, a , cp,sp, cn, bn, n;`
328			`float param0[PARAM_NPARAM], param1[PARAM_NPARAM],`
329			`param[PARAM_NPARAM],`
330			`**list,`
331			`*cov,`
332	285	bertin	`psf_fwhm, err, aspect, chi2, flux;`
333	259	bertin	`int *index,`
334	285	bertin	`i,j,p,s, nparam, nparam2, ncomp, nsub;`
335	269	bertin
336			`nparam = profit->nparam;`
337			`nparam2 = nparam*nparam;`
338	285	bertin	`nsub = profit->nsubprofit;`
339			`subprofit = profit->subprofit;`
340	269	bertin
341	247	bertin	`if (profit->nlimmin)`
342	259	bertin	`profit->flag \|= PROFLAG_MINLIM;`
343	247	bertin	`if (profit->nlimmax)`
344	259	bertin	`profit->flag \|= PROFLAG_MAXLIM;`
345	247	bertin
346	128	bertin	`for (p=0; p<nparam; p++)`
347			`profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);`
348	126	bertin
349	268	bertin	`/* Clean up FFTW plans */`
350			`fft_reset();`
351
352	50	bertin	`/* CHECK-Images */`
353	221	bertin	`if ((check = prefs.check[CHECK_PROFILES]))`
354			`{`
355	258	bertin	`profit_residuals(profit, 0.0, profit->paraminit, NULL);`
356	285	bertin	`check_add(check, subprofit->lmodpix,`
357			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
358			`subprofit->ix,subprofit->iy, 1.0);`
359	221	bertin	`}`
360	235	bertin
361	50	bertin	`if ((check = prefs.check[CHECK_SUBPROFILES]))`
362	51	bertin	`{`
363	258	bertin	`profit_residuals(profit, 0.0, profit->paraminit, NULL);`
364	285	bertin	`check_add(check, subprofit->lmodpix,`
365			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
366			`subprofit->ix,subprofit->iy, -1.0);`
367	51	bertin	`}`
368	221	bertin	`if ((check = prefs.check[CHECK_SPHEROIDS]))`
369	51	bertin	`{`
370	221	bertin	`/-- Set to 0 flux components that do not belong to spheroids /`
371			`for (p=0; p<profit->nparam; p++)`
372			`param[p] = profit->paraminit[p];`
373			`list = profit->paramlist;`
374			`index = profit->paramindex;`
375			`for (i=0; i<PARAM_NPARAM; i++)`
376			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
377			`param[index[i]] = 0.0;`
378	258	bertin	`profit_residuals(profit, 0.0, param, NULL);`
379	285	bertin	`check_add(check, subprofit->lmodpix,`
380			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
381			`subprofit->ix,subprofit->iy, 1.0);`
382	51	bertin	`}`
383	221	bertin	`if ((check = prefs.check[CHECK_SUBSPHEROIDS]))`
384			`{`
385			`/-- Set to 0 flux components that do not belong to spheroids /`
386			`for (p=0; p<profit->nparam; p++)`
387			`param[p] = profit->paraminit[p];`
388			`list = profit->paramlist;`
389			`index = profit->paramindex;`
390			`for (i=0; i<PARAM_NPARAM; i++)`
391			`if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)`
392			`param[index[i]] = 0.0;`
393	258	bertin	`profit_residuals(profit, 0.0, param, NULL);`
394	285	bertin	`check_add(check, subprofit->lmodpix,`
395			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
396			`subprofit->ix,subprofit->iy, -1.0);`
397	221	bertin	`}`
398			`if ((check = prefs.check[CHECK_DISKS]))`
399			`{`
400			`/-- Set to 0 flux components that do not belong to disks /`
401			`for (p=0; p<profit->nparam; p++)`
402			`param[p] = profit->paraminit[p];`
403			`list = profit->paramlist;`
404			`index = profit->paramindex;`
405			`for (i=0; i<PARAM_NPARAM; i++)`
406			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
407			`param[index[i]] = 0.0;`
408	258	bertin	`profit_residuals(profit, 0.0, param, NULL);`
409	285	bertin	`check_add(check, subprofit->lmodpix,`
410			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
411			`subprofit->ix,subprofit->iy, 1.0);`
412	221	bertin	`}`
413			`if ((check = prefs.check[CHECK_SUBDISKS]))`
414			`{`
415			`/-- Set to 0 flux components that do not belong to disks /`
416			`for (p=0; p<profit->nparam; p++)`
417			`param[p] = profit->paraminit[p];`
418			`list = profit->paramlist;`
419			`index = profit->paramindex;`
420			`for (i=0; i<PARAM_NPARAM; i++)`
421			`if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)`
422			`param[index[i]] = 0.0;`
423	258	bertin	`profit_residuals(profit, 0.0, param, NULL);`
424	285	bertin	`check_add(check, subprofit->lmodpix,`
425			`subprofit->objnaxisn[0],subprofit->objnaxisn[1],`
426			`subprofit->ix,subprofit->iy, -1.0);`
427	221	bertin	`}`
428	235	bertin
429	221	bertin	`/* Compute compressed residuals */`
430	258	bertin	`profit_residuals(profit, 10.0, profit->paraminit,profit->resi);`
431	221	bertin
432	73	bertin	`/* Fill measurement parameters */`
433	69	bertin	`if (FLAG(obj2.prof_vector))`
434			`{`
435	128	bertin	`for (p=0; p<nparam; p++)`
436	221	bertin	`obj2->prof_vector[p]= profit->paraminit[p];`
437	69	bertin	`}`
438	127	bertin	`if (FLAG(obj2.prof_errvector))`
439			`{`
440	128	bertin	`for (p=0; p<nparam; p++)`
441			`obj2->prof_errvector[p]= profit->paramerr[p];`
442	127	bertin	`}`
443	208	bertin	`if (FLAG(obj2.prof_errmatrix))`
444			`{`
445			`for (p=0; p<nparam2; p++)`
446			`obj2->prof_errmatrix[p]= profit->covar[p];`
447			`}`
448	89	bertin
449	79	bertin	`obj2->prof_niter = profit->niter;`
450	285	bertin	`subprofit = profit->subprofit;`
451			`for (s=0; s<nsub; s++, subprofit++)`
452	194	bertin	`{`
453	285	bertin	`if (FLAG(obj2.flux_prof))`
454			`obj2->flux_prof[s] = subprofit->flux;`
455			`if (FLAG(obj2.mag_prof))`
456			`obj2->mag_prof[s] = subprofit->flux>0.0?`
457			`-FDMAG * log(subprofit->flux) + prefs.mag_zeropoint[s] : 99.0f;`
458			`if (FLAG(obj2.fluxerr_prof) \|\| FLAG(obj2.magerr_prof))`
459			`{`
460			`err = 0.0;`
461			`cov = profit->covar;`
462			`index = profit->paramindex;`
463			`list = profit->paramlist;`
464			`for (i=0; i<PARAM_NPARAM; i++)`
465			`if (flux_flag[i]==s+1 && list[i])`
466			`{`
467			`cov = profit->covar + nparam*index[i];`
468			`for (j=0; j<PARAM_NPARAM; j++)`
469			`if (flux_flag[j]==s+1 && list[j])`
470			`err += cov[index[j]];`
471			`}`
472			`if (err<0.0)`
473			`err = 0.0;`
474			`if (FLAG(obj2.fluxerr_prof))`
475			`obj2->fluxerr_prof[s] = sqrt(err);`
476			`if (FLAG(obj2.magerr_prof))`
477			`obj2->magerr_prof[s] = subprofit->flux>0.0?`
478			`FDMAG * sqrt(err) / subprofit->flux : 99.0f;`
479			`}`
480	194	bertin	`}`
481
482	180	bertin	`obj2->prof_chi2 = (profit->nresi > profit->nparam)?`
483			`profit->chi2 / (profit->nresi - profit->nparam) : 0.0;`
484
485	245	bertin	`/* Position */`
486	131	bertin	`if (FLAG(obj2.x_prof))`
487			`{`
488			`i = profit->paramindex[PARAM_X];`
489			`j = profit->paramindex[PARAM_Y];`
490	145	bertin	`/-- Model coordinates follow the FITS convention (first pixel at 1,1) /`
491	179	bertin	`if (profit->paramlist[PARAM_X])`
492			`{`
493	285	bertin	`obj2->x_prof = (double)profit->subprofit->ix + *profit->paramlist[PARAM_X]`
494			`+ 1.0; /* !CHECK / / FITS convention */`
495	179	bertin	`obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];`
496			`}`
497			`else`
498			`emx2 = 0.0;`
499			`if (profit->paramlist[PARAM_Y])`
500			`{`