profit.c

NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	06/09/2011
 ***/
void	 profit_surface(profitstruct *profit, obj2struct *obj2)
  {
   profitstruct	hdprofit;
   double	dsum,dhsum,dsumoff, dhval, frac, seff;
   float	*spix, *spixt,
		val,vmax,
		scalefac, imsizefac, flux, lost, sum, lostfluxfrac;
   int		i,p, imax, npix, neff;

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

  for (p=0; p<profit->nparam; p++)
    profit->param[p] = profit->paraminit[p];
  lost = sum = 0.0;

  for (p=0; p<profit->nprof; p++)
    {
    sum += (flux = prof_add(&hdprofit, profit->prof[p],0));
    lost += flux*profit->prof[p]->lostfluxfrac;
    }
  lostfluxfrac = sum > 0.0? lost / sum : 0.0;
/*
char filename[256];
checkstruct *check;
sprintf(filename, "raster_%02d.fits", the_gal);
check=initcheck(filename, CHECK_OTHER, 0);
check->width = hdprofit.modnaxisn[0];
check->height = hdprofit.modnaxisn[1];
reinitcheck(the_field, check);
memcpy(check->pix,hdprofit.modpix,check->npix*sizeof(float));
reendcheck(the_field, check);
endcheck(check);
*/
  if (FLAG(obj2.fluxeff_prof))
    {
/*-- Sort model pixel values */
    spix = NULL;			/* to avoid gcc -Wall warnings */
    QMEMCPY(hdprofit.modpix, spix, float, npix);
    fqmedian(spix, npix);
/*-- Build a cumulative distribution */
    dsum = 0.0;
    spixt = spix;
    for (i=npix; i--;)
      dsum += (double)*(spixt++);
/*-- Find matching surface brightness */
    if (lostfluxfrac > 1.0)
      lostfluxfrac = 0.0;
    dhsum = 0.5 * dsum / (1.0-lostfluxfrac);
    dsum = lostfluxfrac * dsum / (1.0-lostfluxfrac);
    neff = 0;
    spixt = spix;
    for (i=npix; i--;)
      if ((dsum += (double)*(spixt++)) >= dhsum)
        {
        neff = i;
        break;
        }
    dhval = (double)*(spixt-1);
    seff = neff;
    dsumoff = 0.0;
    if (spixt>=spix+2)
      if (dhval > 0.0 && (frac = (dsum - dhsum) / dhval) < 1.0)
        {
        seff += frac;
        dsumoff = frac*dhval;
        dhval = dsumoff + (1.0 - frac)*(double)*(spixt-2);
        }
    obj2->fluxeff_prof = dhval;
    if (FLAG(obj2.fluxmean_prof))
      {
      dsum = dsumoff;
      for (i=neff; i--;)
        dsum += (double)*(spixt++);
      obj2->fluxmean_prof = seff > 0.0? dsum / seff : 0.0;
      }
    free(spix);
    }

/* Compute model peak */
  if (FLAG(obj2.peak_prof))
    {
/*-- Find position of maximum pixel in current hi-def raster */
    imax = 0;
    vmax = -BIG;
    spixt = hdprofit.modpix;
    for (i=npix; i--;)
      if ((val=*(spixt++))>vmax)
        {
        vmax = val;
        imax = i;
        }
    imax = npix-1 - imax;
    obj2->peak_prof = hdprofit.modpix[imax];
    }

/* Free hi-def model raster */
  free(hdprofit.modpix);

  return;
  }


/****** profit_addparam *******************************************************
PROTO	void profit_addparam(profitstruct *profit, paramenum paramindex,
		float **param)
PURPOSE	Add a profile parameter to the list of fitted items.
INPUT	Pointer to the profit structure,
	Parameter index,
	Pointer to the parameter pointer.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	29/03/2010
 ***/
void	profit_addparam(profitstruct *profit, paramenum paramindex,
		float **param)
  {
/* Check whether the parameter has already be registered */
  if (profit->paramlist[paramindex])
/*-- Yes */
    *param = profit->paramlist[paramindex];
  else
/*-- No */
    {
    *param = profit->paramlist[paramindex] = &profit->param[profit->nparam];
    profit->paramindex[paramindex] = profit->nparam;
    profit->paramrevindex[profit->nparam++] = paramindex;
    }

  return;
  }


/****** profit_resetparam ****************************************************
PROTO	void profit_resetparam(profitstruct *profit, paramenum paramtype)
PURPOSE	Set the initial, lower and upper boundary values of a profile parameter.
INPUT	Pointer to the profit structure,
	Parameter index.
OUTPUT	.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	24/04/2013
 ***/
void	profit_resetparam(profitstruct *profit, paramenum paramtype)
  {
   obj2struct	*obj2;
   float	param, parammin,parammax, range, priorcen,priorsig;
   parfitenum	fittype;

  obj2 = profit->obj2;
  param = parammin = parammax = priorcen = priorsig = 0.0;

  switch(paramtype)
    {
    case PARAM_BACK:
      fittype = PARFIT_LINBOUND;
      param = 0.0;
      parammin = -6.0*profit->guesssigbkg;
      parammax =  6.0*profit->guesssigbkg;
      break;
    case PARAM_X:
      fittype = PARFIT_LINBOUND;
      param = profit->guessdx;
      range = profit->guessradius*4.0;
      if (range>profit->objnaxisn[0]*profit->subsamp*2.0)
        range = profit->objnaxisn[0]*profit->subsamp*2.0;
      parammin = -range;
      parammax =  range;
      break;
    case PARAM_Y:
      fittype = PARFIT_LINBOUND;
      param = profit->guessdy;
      range = profit->guessradius*4.0;
      if (range>profit->objnaxisn[1]*profit->subsamp*2.0)
        range = profit->objnaxisn[1]*profit->subsamp*2.0;
      parammin = -range;
      parammax =  range;
      break;
    case PARAM_DIRAC_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 10.0*profit->guessfluxmax;
      break;
    case PARAM_SPHEROID_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 10.0*profit->guessfluxmax;
      break;
    case PARAM_SPHEROID_REFF:
      fittype = PARFIT_LOGBOUND;
      param = FLAG(obj2.prof_disk_flux)? profit->guessradius
			: profit->guessradius/sqrtf(profit->guessaspect);
      parammin = 0.01;
      parammax = param * 10.0;
      break;
    case PARAM_SPHEROID_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = FLAG(obj2.prof_disk_flux)? 1.0 : profit->guessaspect;
      parammin = FLAG(obj2.prof_disk_flux)? 0.5 : 0.01;
      parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0;
      priorcen = 0.3;
      priorsig = 0.0 /*0.4*/;
      break;
    case PARAM_SPHEROID_POSANG:
      fittype = PARFIT_UNBOUND;
      param = profit->guessposang;
      parammin = 90.0;
      parammax =  90.0;
      break;
    case PARAM_SPHEROID_SERSICN:
      fittype = PARFIT_LINBOUND;
      param = 4.0;
      parammin = FLAG(obj2.prof_disk_flux)? 1.0 : 0.3;
      parammax = 10.0;
      priorcen = 1.0;
      priorsig = 0.0 /*2.0*/;
      break;
    case PARAM_DISK_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 10.0*profit->guessfluxmax;
      break;
    case PARAM_DISK_SCALE:	/* From scalelength to Re */
      fittype = PARFIT_LOGBOUND;
      param = profit->guessradius/(1.67835*sqrtf(profit->guessaspect));
      parammin = 0.01/1.67835;
      parammax = param * 10.0;
      break;
    case PARAM_DISK_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessaspect;
      parammin = 0.01;
      parammax = 100.0;
      break;
    case PARAM_DISK_POSANG:
      fittype = PARFIT_UNBOUND;
      param = profit->guessposang;
      parammin = 90.0;
      parammax =  90.0;
      break;
    case PARAM_ARMS_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 10.0*profit->guessfluxmax;
      break;
    case PARAM_ARMS_QUADFRAC:
      fittype = PARFIT_LINBOUND;
      param = 0.5;
      parammin = 0.0;
      parammax = 1.0;
      break;
    case PARAM_ARMS_SCALE:
      fittype = PARFIT_LINBOUND;
      param = 1.0;
      parammin = 0.5;
      parammax = 10.0;
      break;
    case PARAM_ARMS_START:
      fittype = PARFIT_LINBOUND;
      param = 0.5;
      parammin = 0.0;
      parammax = 3.0;
      break;
    case PARAM_ARMS_PITCH:
      fittype = PARFIT_LINBOUND;
      param = 20.0;
      parammin = 5.0;
      parammax = 50.0;
      break;
    case PARAM_ARMS_PITCHVAR:
      fittype = PARFIT_LINBOUND;
      param = 0.0;
      parammin = -1.0;
      parammax = 1.0;
      break;
//      if ((profit->spirindex=profit_spiralindex(profit, obj, obj2)) > 0.0)
//        {
//        param = -param;
//        parammin = -parammax;
//        parammax = -parammin;
//        }
//      printf("spiral index: %g  \n", profit->spirindex);
//      break;
    case PARAM_ARMS_POSANG:
      fittype = PARFIT_UNBOUND;
      param = 0.0;
      parammin = 0.0;
      parammax = 0.0;
      break;
    case PARAM_ARMS_WIDTH:
      fittype = PARFIT_LINBOUND;
      param = 3.0;
      parammin = 1.5;
      parammax = 11.0;
      break;
    case PARAM_BAR_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 4.0*profit->guessfluxmax;
      break;
    case PARAM_BAR_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = 0.3;
      parammin = 0.2;
      parammax = 0.5;
      break;
    case PARAM_BAR_POSANG:
      fittype = PARFIT_UNBOUND;
      param = 0.0;
      parammin = 90.0;
      parammax = 90.0;
      break;
    case PARAM_INRING_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 4.0*profit->guessfluxmax;
      break;
    case PARAM_INRING_WIDTH:
      fittype = PARFIT_LINBOUND;
      param = 0.3;
      parammin = 0.0;
      parammax = 0.5;
      break;
    case PARAM_INRING_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = 0.8;
      parammin = 0.4;
      parammax = 1.0;
      break;
    case PARAM_OUTRING_FLUX:
      fittype = PARFIT_LOGBOUND;
      param = profit->guessflux/profit->nprof;
      parammin = 0.00001*profit->guessfluxmax;
      parammax = 4.0*profit->guessfluxmax;
      break;
    case PARAM_OUTRING_START:
      fittype = PARFIT_LINBOUND;
      param = 4.0;
      parammin = 3.5;
      parammax = 6.0;
      break;
    case PARAM_OUTRING_WIDTH:
      fittype = PARFIT_LINBOUND;
      param = 0.3;
      parammin = 0.0;
      parammax = 0.5;
      break;
    default:
      error(EXIT_FAILURE, "*Internal Error*: Unknown profile parameter in ",
		"profit_resetparam()");
      break;
   }

  if (parammin!=parammax && (param<=parammin || param>=parammax))
    param = (parammin+parammax)/2.0;
  else if (parammin==0.0 && parammax==0.0)
    parammax = 1.0;
  profit_setparam(profit, paramtype, param, parammin, parammax, fittype,
	priorcen, priorsig);

  return;
  }


/****** profit_resetparams ****************************************************
PROTO	void profit_resetparams(profitstruct *profit)
PURPOSE	Set the initial, lower and upper boundary values of profile parameters.
INPUT	Pointer to the profit structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	18/09/2008
 ***/
void	profit_resetparams(profitstruct *profit)
  {
   int		p;


  for (p=0; p<PARAM_NPARAM; p++)
    profit_resetparam(profit, (paramenum)p);

  return;
  }


/****** profit_setparam ****************************************************
PROTO	void profit_setparam(profitstruct *profit, paramenum paramtype,
		float param, float parammin, float parammax,
		parfitenum parfittype,
		float priorcen, float priorsig)
PURPOSE	Set the actual, lower and upper boundary values of a profile parameter.
INPUT	Pointer to the profit structure,
	parameter index,
	actual value,
	lower boundary to the parameter,
	upper boundary to the parameter,
	parameter fitting type.
	prior central value,
	prior standard deviation.
OUTPUT	RETURN_OK if the parameter is registered, RETURN_ERROR otherwise.
AUTHOR	E. Bertin (IAP)
VERSION	16/01/2013
 ***/
int	profit_setparam(profitstruct *profit, paramenum paramtype,
		float param, float parammin,float parammax,
		parfitenum parfittype,
		float priorcen, float priorsig)
  {
   double	dtemp;
   float	*paramptr;
   int		index;

/* Check whether the parameter has already be registered */
  if ((paramptr=profit->paramlist[(int)paramtype]))
    {
    index = profit->paramindex[(int)paramtype];
    profit->paraminit[index] = param;
    profit->parammin[index] = parammin;
    profit->parammax[index] = parammax;
    profit->parfittype[index] = parfittype;
    profit_boundtounbound(profit, &priorcen, &profit->dparampcen[index], index);
    profit->dparampsig[index] = priorsig;
    return RETURN_OK;
    }
  else
    return RETURN_ERROR;
  }

  
/****** profit_boundtounbound *************************************************
PROTO	int profit_boundtounbound(profitstruct *profit,
				float *param, double *dparam, int index)
PURPOSE	Convert parameters from bounded to unbounded space.
INPUT	Pointer to the profit structure,
	input array of single precision parameters,
	output (incomplete) array of double precision parameters,
	parameter selection index (<0 = all parameters)
OUTPUT	Number of free parameters.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	20/05/2011
 ***/
int	profit_boundtounbound(profitstruct *profit,
				float *param, double *dparam, int index)
  {
   double	num,den, tparam;
   int		f,p, pstart,np;

  if (index<0)
    {
    pstart = 0;
    np = profit->nparam;
    }
  else
    {
    pstart = index;
    np = pstart+1;
    }

  f = 0;
  for (p=pstart ; p<np; p++)
    switch(profit->parfittype[p])
      {
      case PARFIT_FIXED:
        break;
      case PARFIT_UNBOUND:
        if (profit->parammax[p] > 0.0 || profit->parammax[p] < 0.0)
          dparam[f] = param[p-pstart] / profit->parammax[p];
        f++;
        break;
      case PARFIT_LINBOUND:
        tparam = param[p-pstart];
        num = tparam - profit->parammin[p];
        den = profit->parammax[p] - tparam;
        dparam[f] = num>1e-50? (den>1e-50? log(num/den): 50.0) : -50.0;
        f++;
        break;
      case PARFIT_LOGBOUND:
        tparam = log(param[p-pstart]);
        num = tparam - log(profit->parammin[p]);
        den = log(profit->parammax[p]) - tparam;
        dparam[f] = num>1e-50? (den>1e-50? log(num/den): 50.0) : -50.0;
        f++;
        break;
      default:
        error(EXIT_FAILURE,
		"*Internal Error*: Unknown parameter fitting type in ",
		"profit_boundtounbound()");
      }

  return f;
  }


/****** profit_unboundtobound *************************************************
PROTO	int profit_unboundtobound(profitstruct *profit,
				double *dparam, float *param, int index)
PURPOSE	Convert parameters from unbounded to bounded space.
INPUT	Pointer to the profit structure,
	input (incomplete) array of double precision parameters,
	output array of single precision parameters.
	parameter selection index (<0 = all parameters)
OUTPUT	Number of free parameters.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	20/05/2011
 ***/
int	profit_unboundtobound(profitstruct *profit,
				double *dparam, float *param, int index)
  {
   int		f,p, pstart,np;

  if (index<0)
    {
    pstart = 0;
    np = profit->nparam;
    }
  else
    {
    pstart = index;
    np = pstart+1;
    }

  f = 0;
  for (p=pstart; p<np; p++)
    switch(profit->parfittype[p])
      {
      case PARFIT_FIXED:
        param[p-pstart] = profit->paraminit[p];
        break;
      case PARFIT_UNBOUND:
        param[p-pstart] = dparam[f]*profit->parammax[p];
        f++;
        break;
      case PARFIT_LINBOUND:
        param[p-pstart] = (profit->parammax[p] - profit->parammin[p])
		/ (1.0 + exp(-(dparam[f]>50.0? 50.0
				: (dparam[f]<-50.0? -50.0: dparam[f]))))
		+ profit->parammin[p];
        f++;
        break;
      case PARFIT_LOGBOUND:
        param[p-pstart] = exp(log(profit->parammax[p]/profit->parammin[p])
		/ (1.0 + exp(-(dparam[f]>50.0? 50.0
				: (dparam[f]<-50.0? -50.0: dparam[f])))))
		*profit->parammin[p];
        f++;
        break;
      default:
        error(EXIT_FAILURE,
		"*Internal Error*: Unknown parameter fitting type in ",
		"profit_unboundtobound()");
      }
 
  return f;
  }


/****** profit_covarunboundtobound ********************************************
PROTO	int profit_covarunboundtobound(profitstruct *profit,
				double *dcovar, float *covar)
PURPOSE	Convert covariance matrix from unbounded to bounded space.
INPUT	Pointer to the profit structure,
	input (incomplete) matrix of double precision covariances,
	output matrix of single precision covariances.
OUTPUT	Number of parameters that hit a boundary.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	20/05/2011
 ***/
int	profit_covarunboundtobound(profitstruct *profit,
				double *dcovar, float *covar)
  {
   double	*dxdy,
		xmmin,maxmx, maxmmin;
   float	*x,*xmin,*xmax;
   parfitenum	*fittype;
   int		*fflag,
		f,f1,f2, p,p1,p2, nfree, nparam, nmin,nmax;

  nparam = profit->nparam;
  fittype = profit->parfittype;
  QMALLOC16(dxdy, double, PARAM_NPARAM);
  x = profit->paraminit;
  xmin = profit->parammin;
  xmax = profit->parammax;
  nmin = nmax = 0;
  f = 0;
  for (p=0; p<profit->nparam; p++)
    switch(fittype[p])
      {
      case PARFIT_FIXED:
        break;
      case PARFIT_UNBOUND:
        dxdy[f++] = xmax[p];
        break;
      case PARFIT_LINBOUND:
        xmmin   = x[p] - xmin[p];
        maxmx   = xmax[p] - x[p];
        maxmmin = xmax[p] - xmin[p];
        dxdy[f++] = xmmin * maxmx / maxmmin;
        if (xmmin<0.001*maxmmin)
          nmin++;
        else if (maxmx<0.001*maxmmin)
          nmax++;
        break;
      case PARFIT_LOGBOUND:
        xmmin   = log(x[p]/xmin[p]);
        maxmx   = log(xmax[p]/x[p]);
        maxmmin = log(xmax[p]/xmin[p]);
        dxdy[f++] = x[p] * xmmin * maxmx / maxmmin;
        if (xmmin<0.001*maxmmin)
          nmin++;
        else if (maxmx<0.001*maxmmin)
          nmax++;
        break;
      default:
        error(EXIT_FAILURE,
		"*Internal Error*: Unknown parameter fitting type in ",
		"profit_boundtounbound()");
      }

  nfree = f;

  memset(profit->covar, 0, nparam*nparam*sizeof(float));
  f2 = 0;
  for (p2=0; p2<nparam; p2++)
    {
    if (fittype[p2]!=PARFIT_FIXED)
      {
      f1 = 0;
      for (p1=0; p1<nparam; p1++)
        if (fittype[p1]!=PARFIT_FIXED)
          {
          covar[p2*nparam+p1] = (float)(dcovar[f2*nfree+f1]*dxdy[f1]*dxdy[f2]);
          f1++;
          }
      f2++;
      }
    }

  free(dxdy);

  profit->nlimmin = nmin;
  profit->nlimmax = nmax;

  return nmin+nmax;
  }


/****** prof_init *************************************************************
PROTO	profstruct prof_init(profitstruct *profit, unsigned int modeltype)
PURPOSE	Allocate and initialize a new profile structure.
INPUT	Pointer to the profile-fitting structure,
	model type.
OUTPUT	A pointer to an allocated prof structure.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	08/12/2011
 ***/
profstruct	*prof_init(profitstruct *profit, unsigned int modeltype)
  {
   profstruct	*prof;
   float	*pix,
		rmax2, re2, dy2,r2, scale, zero, k,n, hinvn;
   int		width,height, ixc,iyc, ix,iy, nsub,
		d,s;

  QCALLOC(prof, profstruct, 1);
  prof->code = modeltype;
  switch(modeltype)
    {
    case MODEL_BACK:
      prof->name = "background offset";
      prof->naxis = 2;
      prof->pix = NULL;
      profit_addparam(profit, PARAM_BACK, &prof->flux);
      prof->typscale = 1.0;
      break;
    case MODEL_DIRAC:
      prof->name = "point source";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DIRAC_FLUX, &prof->flux);
      break;
    case MODEL_SERSIC:
      prof->name = "Sersic spheroid";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_SPHEROID_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_SPHEROID_REFF, &prof->scale);
      profit_addparam(profit, PARAM_SPHEROID_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_SPHEROID_SERSICN, &prof->extra[0]);
      break;
    case MODEL_DEVAUCOULEURS:
      prof->name = "de Vaucouleurs spheroid";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_SPHEROID_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_SPHEROID_REFF, &prof->scale);
      profit_addparam(profit, PARAM_SPHEROID_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle);
      break;
    case MODEL_EXPONENTIAL:
      prof->name = "exponential disk";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DISK_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_DISK_SCALE, &prof->scale);
      profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
      break;
    case MODEL_ARMS:
      prof->name = "spiral arms";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DISK_SCALE, &prof->scale);
      profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_ARMS_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_ARMS_QUADFRAC, &prof->featfrac);
//      profit_addparam(profit, PARAM_ARMS_SCALE, &prof->featscale);
      profit_addparam(profit, PARAM_ARMS_START, &prof->featstart);
      profit_addparam(profit, PARAM_ARMS_PITCH, &prof->featpitch);
//      profit_addparam(profit, PARAM_ARMS_PITCHVAR, &prof->featpitchvar);
      profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang);
//      profit_addparam(profit, PARAM_ARMS_WIDTH, &prof->featwidth);
      break;
    case MODEL_BAR:
      prof->name = "bar";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DISK_SCALE, &prof->scale);
      profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_ARMS_START, &prof->featstart);
      profit_addparam(profit, PARAM_BAR_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_BAR_ASPECT, &prof->feataspect);
      profit_addparam(profit, PARAM_ARMS_POSANG, &prof->featposang);
      break;
    case MODEL_INRING:
      prof->name = "inner ring";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DISK_SCALE, &prof->scale);
      profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_ARMS_START, &prof->featstart);
      profit_addparam(profit, PARAM_INRING_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_INRING_WIDTH, &prof->featwidth);
      profit_addparam(profit, PARAM_INRING_ASPECT, &prof->feataspect);
      break;
    case MODEL_OUTRING:
      prof->name = "outer ring";
      prof->naxis = 2;
      prof->naxisn[0] = PROFIT_MAXMODSIZE;
      prof->naxisn[1] = PROFIT_MAXMODSIZE;
      prof->naxisn[2] = 1;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      prof->typscale = 1.0;
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_DISK_SCALE, &prof->scale);
      profit_addparam(profit, PARAM_DISK_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_DISK_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_OUTRING_START, &prof->featstart);
      profit_addparam(profit, PARAM_OUTRING_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_OUTRING_WIDTH, &prof->featwidth);
      break;
    case MODEL_TABULATED:	/* An example of tabulated profile */
      prof->name = "tabulated model";
      prof->naxis = 3;
      width =  prof->naxisn[0] = PROFIT_MAXMODSIZE;
      height = prof->naxisn[1] = PROFIT_MAXMODSIZE;
      nsub = prof->naxisn[2] = PROFIT_MAXSMODSIZE;
      prof->npix = prof->naxisn[0]*prof->naxisn[1]*prof->naxisn[2];
      QMALLOC(prof->pix, float, prof->npix);
      ixc = width/2;
      iyc = height/2;
      rmax2 = (ixc - 1.0)*(ixc - 1.0);
      re2 = width/64.0;
      prof->typscale = re2;
      re2 *= re2;
      zero = prof->extrazero[0] = 0.2;
      scale = prof->extrascale[0]= 8.0/PROFIT_MAXSMODSIZE;
      pix = prof->pix;
      for (s=0; s<nsub; s++)
        {
        n = s*scale + zero;
        hinvn = 0.5/n;
        k = -1.0/3.0 + 2.0*n + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
		+ 131.0/(1148175*n*n*n);
        for (iy=0; iy<height; iy++)
          {
          dy2 = (iy-iyc)*(iy-iyc);
          for (ix=0; ix<width; ix++)
            {
            r2 = dy2 + (ix-ixc)*(ix-ixc);
            *(pix++) = (r2<rmax2)? exp(-k*pow(r2/re2,hinvn)) : 0.0;
            }
          }
        }
      profit_addparam(profit, PARAM_X, &prof->x[0]);
      profit_addparam(profit, PARAM_Y, &prof->x[1]);
      profit_addparam(profit, PARAM_SPHEROID_FLUX, &prof->flux);
      profit_addparam(profit, PARAM_SPHEROID_REFF, &prof->scale);
      profit_addparam(profit, PARAM_SPHEROID_ASPECT, &prof->aspect);
      profit_addparam(profit, PARAM_SPHEROID_POSANG, &prof->posangle);
      profit_addparam(profit, PARAM_SPHEROID_SERSICN, &prof->extra[0]);
      break;
    default:
      error(EXIT_FAILURE, "*Internal Error*: Unknown profile in ",
		"prof_init()");
      break;
    }

  if (prof->naxis>2)
    {
    prof->kernelnlines = 1;
    for (d=0; d<prof->naxis; d++)
      {
      prof->interptype[d] = INTERP_BILINEAR;
      prof->kernelnlines *= 
	(prof->kernelwidth[d] = interp_kernwidth[prof->interptype[d]]);
      }
    prof->kernelnlines /= prof->kernelwidth[0];
    QMALLOC16(prof->kernelbuf, float, prof->kernelnlines);
    }

  return prof;
  }  


/****** prof_end **************************************************************
PROTO	void prof_end(profstruct *prof)
PURPOSE	End (deallocate) a profile structure.
INPUT	Prof structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	09/04/2007
 ***/
void	prof_end(profstruct *prof)
  {
  if (prof->pix)
    {
    free(prof->pix);
    free(prof->kernelbuf);
    }
  free(prof);

  return;
  }


/****** prof_add **************************************************************
PROTO	float prof_add(profitstruct *profit, profstruct *prof,
		int extfluxfac_flag)
PURPOSE	Add a model profile to an image.
INPUT	Profile-fitting structure,
	profile structure,
	flag (0 if flux correction factor is to be computed internally) 
OUTPUT	Total (asymptotic) flux contribution.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	13/02/2013
 ***/
float	prof_add(profitstruct *profit, profstruct *prof, int extfluxfac_flag)
  {
   double	xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n,
		dx1cout,dx2cout, ddx1[36],ddx2[36];
   float	posin[PROFIT_MAXEXTRA], posout[2], dnaxisn[2],
		*pixin, *pixin2, *pixout,
		fluxfac, amp,cd11,cd12,cd21,cd22, dx1,dx2,
		x1,x10,x2, x1cin,x2cin, x1cout,x2cout, x1max,x2max, x1in,x2in,
		k, hinvn, x1t,x2t, ca,sa, u,umin,
		armamp,arm2amp, armrdphidr, armrdphidrvar, posang,
		width, invwidth2,
		r,r2,rmin, r2minxin,r2minxout, rmax, r2max,
		r2max1, r2max2, r2min, invr2xdif,
		val, theta, thresh, ra,rb, num,num2,den, ang,angstep,
		invn, dr, krpinvn,dkrpinvn, rs,rs2,
		a11,a12,a21,a22, invdet, dca,dsa, a0,a2,a3, p1,p2,
		krspinvn, ekrspinvn, selem;
   int		npix, threshflag,
		a,d,e,i, ix1,ix2, ix1max,ix2max, nang, nx2,
		npix2;

  npix = profit->nmodpix;

  if (prof->code==MODEL_BACK)
    {
    amp = fabs(*prof->flux);
    pixout = profit->modpix;
    for (i=npix; i--;)
      *(pixout++) += amp;
    prof->lostfluxfrac = 0.0;
    return 0.0;
    }

  scaling = profit->pixstep / prof->typscale;

  if (prof->code!=MODEL_DIRAC)
    {
/*-- Compute Profile CD matrix */
    ctheta = cos(*prof->posangle*DEG);
    stheta = sin(*prof->posangle*DEG);
    saspect = fabs(*prof->aspect);
    xscale = (*prof->scale==0.0)?
		0.0 : fabs(scaling / (*prof->scale*prof->typscale));
    yscale = (*prof->scale*saspect == 0.0)?
		0.0 : fabs(scaling / (*prof->scale*prof->typscale*saspect));
    cd11 = (float)(xscale*ctheta);
    cd12 = (float)(xscale*stheta);
    cd21 = (float)(-yscale*stheta);
    cd22 = (float)(yscale*ctheta);
    dx1 = 0.0;	/* Shifting operations have been moved to profit_resample() */
    dx2 = 0.0;	/* Shifting operations have been moved to profit_resample() */

    x1cout = (float)(profit->modnaxisn[0]/2);
    x2cout = (float)(profit->modnaxisn[1]/2);
    nx2 = profit->modnaxisn[1]/2 + 1;

/*-- Compute the largest r^2 that fits in the frame */
    num = cd11*cd22-cd12*cd21;