profit.c

    }

  return;
  }


/****** profit_covarunboundtobound ********************************************
PROTO	void 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	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	27/07/2010
 ***/
void	profit_covarunboundtobound(profitstruct *profit,
				double *dcovar, float *covar)
  {
   double	*dxdy;
   float	*x,*xmin,*xmax;
   int		*fflag,
		f,f1,f2, nfree, p,p1,p2, nparam;

  nparam = profit->nparam;
  fflag = profit->freeparam_flag;
  nfree = profit->nfreeparam;
  QMALLOC16(dxdy, double, nfree);
  x = profit->paraminit;
  xmin = profit->parammin;
  xmax = profit->parammax;
  f = 0;
  for (p=0; p<profit->nparam; p++)
    if (fflag[p])
      {
      if (xmin[p]!=xmax[p])
        dxdy[f++] = (x[p] - xmin[p]) * (xmax[p] - x[p]) / (xmax[p] - xmin[p]);
      else
        dxdy[f++] = xmax[p];
      }

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

  free(dxdy);

  return;
  }


/****** prof_init *************************************************************
PROTO	profstruct prof_init(profitstruct *profit, proftypenum profcode)
PURPOSE	Allocate and initialize a new profile structure.
INPUT	Pointer to the profile-fitting structure,
	profile type.
OUTPUT	A pointer to an allocated prof structure.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	08/10/2010
 ***/
profstruct	*prof_init(profitstruct *profit, proftypenum profcode)
  {
   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 = profcode;
  switch(profcode)
    {
    case PROF_BACK:
      prof->naxis = 2;
      prof->pix = NULL;
      profit_addparam(profit, PARAM_BACK, &prof->flux);
      prof->typscale = 1.0;
      break;
    case PROF_DIRAC:
      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];
      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 PROF_SERSIC:
      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 PROF_DEVAUCOULEURS:
      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, profit->nmodpix);
      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 PROF_EXPONENTIAL:
      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, profit->nmodpix);
      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 PROF_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, profit->nmodpix);
      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 PROF_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, profit->nmodpix);
      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 PROF_INRING:
      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, profit->nmodpix);
      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 PROF_OUTRING:
      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, profit->nmodpix);
      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 PROF_SERSIC_TABEX:	/* An example of tabulated profile */
      prof->naxis = 3;
      width =  prof->naxisn[0] = PROFIT_PROFRES;
      height = prof->naxisn[1] = PROFIT_PROFRES;
      nsub = prof->naxisn[2] = PROFIT_PROFSRES;
      QCALLOC(prof->pix, float, width*height*nsub);
      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_PROFSRES;
      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;
    }

  QMALLOC(prof->pix, float, prof->npix);

  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	11/03/2011
 ***/
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==PROF_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!=PROF_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;
    num *= num;
    x1max = x1cout - 1.0;
    x2max = x2cout - 1.0;
    den = fabs(cd12*cd12+cd22*cd22);
    num2 = x1max*x1max*num;
    r2max1 = num2<PROFIT_MAXR2MAX*den? num2 / den : PROFIT_MAXR2MAX;
    den = fabs(cd11*cd11+cd21*cd21);
    num2 = x2max*x2max*num;
    r2max2 = num2<PROFIT_MAXR2MAX*den? num2 / den : PROFIT_MAXR2MAX;
    r2max = (r2max1 < r2max2? r2max1 : r2max2);
    rmax = sqrtf(r2max);
    }

  switch(prof->code)
    {
    case PROF_DIRAC:
      prof->pix[profit->modnaxisn[0]/2
		+ (profit->modnaxisn[1]/2)*profit->modnaxisn[0]] = 1.0;
      prof->lostfluxfrac = 0.0;
      threshflag = 0;
      break;
    case PROF_SERSIC:
    case PROF_DEVAUCOULEURS:
    case PROF_EXPONENTIAL:
/*---- Compute sharp/smooth transition radius */
      rs = PROFIT_SMOOTHR*(xscale>yscale?xscale:yscale);
      if (rs<=0)
        rs = 1.0;
      rs2 = rs*rs;
/*---- The consequence of sampling on flux is compensated by PSF normalisation*/
      if (prof->code==PROF_EXPONENTIAL)
        bn = n = 1.0;
      else if (prof->code==PROF_DEVAUCOULEURS)
        {
        n = 4.0;
        bn = 7.66924944;
        }
      else
        {
        n = fabs(*prof->extra[0]);
        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 */
        }
      invn = 1.0/n;
      hinvn = 0.5/n;
      k = -bn;
/*---- Compute central polynomial terms */
      krspinvn = prof->code==PROF_EXPONENTIAL? -rs : k*expf(logf(rs)*invn);
      ekrspinvn = expf(krspinvn);
      p2 = krspinvn*invn*invn;
      p1 = krspinvn*p2;
      a0 = (1+(1.0/6.0)*(p1+(1.0-5.0*n)*p2))*ekrspinvn;
      a2 = (-1.0/2.0)*(p1+(1.0-3.0*n)*p2)/rs2*ekrspinvn;
      a3 = (1.0/3.0)*(p1+(1.0-2.0*n)*p2)/(rs2*rs)*ekrspinvn;
/*---- Compute the smooth part of the profile */
      x10 = -x1cout - dx1;
      x2 = -x2cout - dx2;
      pixin = prof->pix;
      if (prof->code==PROF_EXPONENTIAL)
        for (ix2=nx2; ix2--; x2+=1.0)
          {
          x1 = x10;
          for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
            {
            x1in = cd12*x2 + cd11*x1;
            x2in = cd22*x2 + cd21*x1;
            ra = x1in*x1in+x2in*x2in;
            if (ra>r2max)
              {
              *(pixin++) = 0.0;
              continue;
              }
            val = ra<rs2? a0+ra*(a2+a3*sqrtf(ra)) : expf(-sqrtf(ra));
            *(pixin++) = val;
            }
          }
      else
        for (ix2=nx2; ix2--; x2+=1.0)
          {
          x1 = x10;
          for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
            {
            x1in = cd12*x2 + cd11*x1;
            x2in = cd22*x2 + cd21*x1;
            ra = x1in*x1in+x2in*x2in;
            if (ra>r2max)
              {
              *(pixin++) = 0.0;
              continue;
              }
            val = ra<rs2? a0+ra*(a2+a3*sqrtf(ra)) : expf(k*expf(logf(ra)*hinvn));
            *(pixin++) = val;
            }
          }
/*---- Copy the symmetric part */
      if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
        {
        pixin2 = pixin - profit->modnaxisn[0] - 1;
        if (!(profit->modnaxisn[0]&1))
          {
          *(pixin++) = 0.0;
          npix2--;
          }
        for (i=npix2; i--;)
          *(pixin++) = *(pixin2--);
        }

/*---- Compute the sharp part of the profile */
      ix1max = profit->modnaxisn[0];
      ix2max = profit->modnaxisn[1];
      dx1cout = x1cout + 0.4999999;
      dx2cout = x2cout + 0.4999999;
      invdet = 1.0/fabsf(cd11*cd22 - cd12*cd21);
      a11 = cd22*invdet;
      a12 = -cd12*invdet;
      a21 = -cd21*invdet;
      a22 = cd11*invdet;
      nang = 72 / 2;		/* 36 angles; only half of them are computed*/
      angstep = PI/nang;
      ang = 0.0;
      for (a=0; a<nang; a++)
        {
        sincosf(ang, &dca, &dsa);
        ddx1[a] = a11*dca+a12*dsa;
        ddx2[a] = a21*dca+a22*dsa;
        ang += angstep;
        }
      r = DEXPF(-4.0);
      dr = DEXPF(0.05);
      selem = 0.5*angstep*(dr - 1.0/dr)/(xscale*yscale);
      krpinvn = k*DEXPF(-4.0*invn);
      dkrpinvn = DEXPF(0.05*invn);
      pixin = prof->pix;
      for (; r<rs; r *= dr)
        {
        r2 = r*r;
        val = (expf(krpinvn) - (a0 + r2*(a2+a3*r)))*r2*selem;
        for (a=0; a<nang; a++)
          {
          ix1 = (int)(dx1cout + r*ddx1[a]);
          ix2 = (int)(dx2cout + r*ddx2[a]);
          if (ix1>=0 && ix1<ix1max && ix2>=0 && ix2<ix2max)
            pixin[ix2*ix1max+ix1] += val;
          ix1 = (int)(dx1cout - r*ddx1[a]);
          ix2 = (int)(dx2cout - r*ddx2[a]);
          if (ix1>=0 && ix1<ix1max && ix2>=0 && ix2<ix2max)
            pixin[ix2*ix1max+ix1] += val;
          }
        krpinvn *= dkrpinvn;
        }

      prof->lostfluxfrac = prof->code==PROF_EXPONENTIAL?
		(1.0 + rmax)*exp(-rmax)
		:1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn));
      threshflag = 0;
      break;
    case PROF_ARMS:
      r2min = *prof->featstart**prof->featstart;
      r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE);
      r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE);
      if ((invr2xdif = (r2minxout - r2minxin)) > 0.0)
        invr2xdif = 1.0 / invr2xdif;
      else
        invr2xdif = 1.0;
      umin = 0.5*logf(r2minxin + 0.00001);
      arm2amp = *prof->featfrac;
      armamp = 1.0 - arm2amp;
      armrdphidr = 1.0/tan(*prof->featpitch*DEG);
      armrdphidrvar = 0.0 /**prof->featpitchvar*/;
      posang = *prof->featposang*DEG;
      width = fabs(*prof->featwidth);
width = 3.0;
      x1 = -x1cout - dx1;
      x2 = -x2cout - dx2;
      pixin = prof->pix;
      for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0)
        {
        x1t = cd12*x2 + cd11*x1;
        x2t = cd22*x2 + cd21*x1;
        for (ix1=profit->modnaxisn[0]; ix1--;)
          {
          r2 = x1t*x1t+x2t*x2t;
          if (r2>r2minxin)
            {
            u = 0.5*logf(r2 + 0.00001);
            theta = (armrdphidr+armrdphidrvar*(u-umin))*u+posang;
            ca = cosf(theta);
            sa = sinf(theta);
            x1in = (x1t*ca - x2t*sa);
            x2in = (x1t*sa + x2t*ca);
            amp = expf(-sqrtf(x1t*x1t+x2t*x2t));
            if (r2<r2minxout)
              amp *= (r2 - r2minxin)*invr2xdif;
            ra = x1in*x1in/r2;
            rb = x2in*x2in/r2;
            *(pixin++) = amp * (armamp*PROFIT_POWF(ra,width)
				+ arm2amp*PROFIT_POWF(rb,width));
            }
          else
            *(pixin++) = 0.0;
          x1t += cd11;
          x2t += cd21; 
         }
        }
      prof->lostfluxfrac = 0.0;
      threshflag = 1;
      break;
    case PROF_BAR:
      r2min = *prof->featstart**prof->featstart;
      r2minxin = r2min * (1.0 - PROFIT_BARXFADE) * (1.0 - PROFIT_BARXFADE);
      r2minxout = r2min * (1.0 + PROFIT_BARXFADE) * (1.0 + PROFIT_BARXFADE);
      if ((invr2xdif = (r2minxout - r2minxin)) > 0.0)
        invr2xdif = 1.0 / invr2xdif;
      else
        invr2xdif = 1.0;
      invwidth2 = fabs(1.0 / (*prof->featstart**prof->feataspect));
      posang = *prof->featposang*DEG;
      ca = cosf(posang);
      sa = sinf(posang);
      x1 = -x1cout - dx1;
      x2 = -x2cout - dx2;
      pixin = prof->pix;
      for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0)
        {
        x1t = cd12*x2 + cd11*x1;
        x2t = cd22*x2 + cd21*x1;
        for (ix1=profit->modnaxisn[0]; ix1--;)
          {
          r2 = x1t*x1t+x2t*x2t;
          if (r2<r2minxout)
            {
            x1in = x1t*ca - x2t*sa;
            x2in = invwidth2*(x1t*sa + x2t*ca);
            *(pixin++) = (r2>r2minxin) ?
				(r2minxout - r2)*invr2xdif*expf(-x2in*x2in)
				: expf(-x2in*x2in);
            }
          else
            *(pixin++) = 0.0;
          x1t += cd11;
          x2t += cd21;
          }
        }
      prof->lostfluxfrac = 0.0;
      threshflag = 1;
      break;
    case PROF_INRING:
      rmin = *prof->featstart;
      r2minxin = *prof->featstart-4.0**prof->featwidth;
      if (r2minxin < 0.0)
        r2minxin = 0.0;
      r2minxin *= r2minxin;
      r2minxout = *prof->featstart+4.0**prof->featwidth;
      r2minxout *= r2minxout;
      invwidth2 = 0.5 / (*prof->featwidth**prof->featwidth);
      cd22 /= *prof->feataspect;
      cd21 /= *prof->feataspect;
      x1 = -x1cout - dx1;
      x2 = -x2cout - dx2;
      pixin = prof->pix;
      for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0)
        {
        x1t = cd12*x2 + cd11*x1;
        x2t = cd22*x2 + cd21*x1;
        for (ix1=profit->modnaxisn[0]; ix1--;)
          {
          r2 = x1t*x1t+x2t*x2t;
          if (r2>r2minxin && r2<r2minxout)
            {
            r = sqrt(r2) - rmin;
            *(pixin++) = expf(-invwidth2*r*r);
            }
          else
            *(pixin++) = 0.0;
          x1t += cd11;
          x2t += cd21;
          }
        }
      prof->lostfluxfrac = 0.0;
      threshflag = 1;
      break;
    case PROF_OUTRING:
      rmin = *prof->featstart;
      r2minxin = *prof->featstart-4.0**prof->featwidth;
      if (r2minxin < 0.0)
        r2minxin = 0.0;
      r2minxin *= r2minxin;
      r2minxout = *prof->featstart+4.0**prof->featwidth;
      r2minxout *= r2minxout;
      invwidth2 = 0.5 / (*prof->featwidth**prof->featwidth);
      x1 = -x1cout - dx1;
      x2 = -x2cout - dx2;
      pixin = prof->pix;
      for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0)
        {
        x1t = cd12*x2 + cd11*x1;
        x2t = cd22*x2 + cd21*x1;
        for (ix1=profit->modnaxisn[0]; ix1--;)
          {
          r2 = x1t*x1t+x2t*x2t;
          if (r2>r2minxin && r2<r2minxout)
            {
            r = sqrt(r2) - rmin;
            *(pixin++) = expf(-invwidth2*r*r);
            }
          else
            *(pixin++) = 0.0;
          x1t += cd11;
          x2t += cd21;
          }
        }
      prof->lostfluxfrac = 0.0;
      threshflag = 1;
      break;
    default:
/*---- Tabulated profile: remap each pixel */
/*---- Initialize multi-dimensional counters */
     for (d=0; d<2; d++)
        {
        posout[d] = 1.0;	/* FITS convention */
        dnaxisn[d] = profit->modnaxisn[d] + 0.99999;
        }

      for (e=0; e<prof->naxis - 2; e++)
        {
        d = 2+e;
/*------ Compute position along axis */
        posin[d] = (*prof->extra[e]-prof->extrazero[e])/prof->extrascale[e]+1.0;
/*------ Keep position within boundaries and let interpolation do the rest */
        if (posin[d] < 0.99999)
          {
          if (prof->extracycleflag[e])
            posin[d] += (float)prof->naxisn[d];
          else
            posin[d] = 1.0;
          }
        else if (posin[d] > (float)prof->naxisn[d])
          {
          if (prof->extracycleflag[e])
          posin[d] = (prof->extracycleflag[e])?
		  fmod(posin[d], (float)prof->naxisn[d])
		: (float)prof->naxisn[d];
          }
        }
      x1cin = (float)(prof->naxisn[0]/2);
      x2cin = (float)(prof->naxisn[1]/2);
      pixin = prof->pix;
      for (i=npix; i--;)
        {
        x1 = posout[0] - x1cout - 1.0 - dx1;
        x2 = posout[1] - x2cout - 1.0 - dx2;
        posin[0] = cd11*x1 + cd12*x2 + x1cin + 1.0;
        posin[1] = cd21*x1 + cd22*x2 + x2cin + 1.0;
        *(pixin++) = prof_interpolate(prof, posin);
        for (d=0; d<2; d++)
          if ((posout[d]+=1.0) < dnaxisn[d])
            break;
          else
            posout[d] = 1.0;
        }
      prof->lostfluxfrac = 0.0;
      threshflag = 1;
    break;
    }

/* For complex profiles, threshold to the brightest pixel value at border */
  if (threshflag)
    {
/*-- Now find truncation threshold */
/*-- Find the shortest distance to a vignet border */
    rmax = x1cout;
    if (rmax > (r = x2cout))
      rmax = r;
    rmax += 0.01;
    if (rmax<1.0)
      rmax = 1.0;
    r2max = rmax*rmax;
    rmin = rmax - 1.0;
    r2min = rmin*rmin;

/*-- Find best threshold (the max around the circle with radius rmax */
    dx2 = -x2cout;
    pixin = prof->pix;
    thresh = -BIG;
    for (ix2=profit->modnaxisn[1]; ix2--; dx2 += 1.0)
      {
      dx1 = -x1cout;
      for (ix1=profit->modnaxisn[0]; ix1--; dx1 += 1.0)
        if ((val=*(pixin++))>thresh && (r2=dx1*dx1+dx2*dx2)>r2min && r2<r2max)
          thresh = val;
      }

/*-- Threshold and measure the flux */
    flux = 0.0;
    pixin = prof->pix;
    for (i=npix; i--; pixin++)
      if (*pixin >= thresh)
        flux += (double)*pixin;
      else
        *pixin = 0.0;
    }
  else
    {
    flux = 0.0;
    pixin = prof->pix;
    for (i=npix; i--;)
      flux += (double)*(pixin++);
    }

  if (extfluxfac_flag)
    fluxfac = prof->fluxfac;
  else
    {
    if (prof->lostfluxfrac < 1.0)
      flux /= (1.0 - prof->lostfluxfrac);

    prof->fluxfac = fluxfac = fabs(flux)>0.0? profit->fluxfac/fabs(flux) : 0.0;
    }

  pixin = prof->pix;
  for (i=npix; i--;)
    *(pixin++) *= fluxfac;

/* Correct final flux */
  fluxfac = *prof->flux;
  pixin = prof->pix;
  pixout = profit->modpix;
  for (i=npix; i--;)
    *(pixout++) += fluxfac**(pixin++);

  return *prof->flux;
  }


/****** prof_moments **********************************************************
PROTO	int	prof_moments(profitstruct *profit, profstruct *prof)
PURPOSE	Computes (analytically or numerically) the 2nd moments of a profile.
INPUT	Profile-fitting structure,
	profile structure,
	optional pointer to 3xnparam Jacobian matrix.
OUTPUT	Index to the profile flux for further processing.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	20/08/2010
 ***/
int	prof_moments(profitstruct *profit, profstruct *prof, double *jac)
  {
   double	*dmx2,*dmy2,*dmxy,
		m20, a2, ct,st, mx2fac, my2fac,mxyfac, dc2,ds2,dcs,
		bn,bn2, n,n2, nfac,nfac2, hscale2, dmdn;
   int		nparam, index;

  if (jac)
/*-- Clear output Jacobian */
    {
    nparam = profit->nparam;
    memset(jac, 0, nparam*3*sizeof(double));
    dmx2 = jac;
    dmy2 = jac + nparam;
    dmxy = jac + 2*nparam;
    }
  else
    dmx2 = dmy2 = dmxy = NULL;		/* To avoid gcc -Wall warnings */

  m20 = 0.0;				/* to avoid gcc -Wall warnings */
  index = 0;				/* to avoid gcc -Wall warnings */


  if (prof->posangle)
    {
    a2 = *prof->aspect**prof->aspect;
    ct = cos(*prof->posangle*DEG);
    st = sin(*prof->posangle*DEG);
    mx2fac = ct*ct + st*st*a2;
    my2fac = st*st + ct*ct*a2;
    mxyfac = ct*st * (1.0 - a2);
    if (jac)
      {
      dc2 = -2.0*ct*st*DEG;
      ds2 =  2.0*ct*st*DEG;
      dcs = (ct*ct - st*st)*DEG;
      }
    else
      dc2 = ds2 = dcs = 0.0;		/* To avoid gcc -Wall warnings */
    switch(prof->code)
      {
      case PROF_SERSIC:
        n = fabs(*prof->extra[0]);
        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 */
        nfac  = prof_gamma(4.0*n) / (prof_gamma(2.0*n)*pow(bn, 2.0*n));
        hscale2 = 0.5 * *prof->scale**prof->scale;
        m20 = hscale2 * nfac;
        if (jac)
          {
          dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
			= *prof->scale * nfac * mx2fac;
          dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
			= *prof->scale * nfac * my2fac;
          dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
			= *prof->scale * nfac * mxyfac;
          n2 = n+0.01;
          bn2 = 2.0*n2 - 1.0/3.0 + 4.0/(405.0*n2) + 46.0/(25515.0*n2*n2)
		+ 131.0/(1148175*n2*n2*n2);	/* Ciotti & Bertin 1999 */
          nfac2 = prof_gamma(4.0*n2) / (prof_gamma(2.0*n2)*pow(bn2, 2.0*n2));
          dmdn = 100.0 * hscale2 * (nfac2-nfac);
          dmx2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mx2fac;
          dmy2[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * my2fac;
          dmxy[profit->paramindex[PARAM_SPHEROID_SERSICN]] = dmdn * mxyfac;
          dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= 2.0 * m20 * st*st * *prof->aspect;
          dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= 2.0 * m20 * ct*ct * *prof->aspect;
          dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= -2.0 * m20 * ct*st * *prof->aspect;
          dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
          dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
          dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
          }
        index = profit->paramindex[PARAM_SPHEROID_FLUX];
        break; 
      case PROF_DEVAUCOULEURS:
        m20 = 10.83995 * *prof->scale**prof->scale;
        if (jac)
          {
          dmx2[profit->paramindex[PARAM_SPHEROID_REFF]]
			= 21.680 * *prof->scale * mx2fac;
          dmy2[profit->paramindex[PARAM_SPHEROID_REFF]]
			= 21.680 * *prof->scale * my2fac;
          dmxy[profit->paramindex[PARAM_SPHEROID_REFF]]
			= 21.680 * *prof->scale * mxyfac;
          dmx2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= 2.0 * m20 * st*st * *prof->aspect;
          dmy2[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= 2.0 * m20 * ct*ct * *prof->aspect;
          dmxy[profit->paramindex[PARAM_SPHEROID_ASPECT]]
			= -2.0 * m20 * ct*st * *prof->aspect;
          dmx2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (dc2+ds2*a2);
          dmy2[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (ds2+dc2*a2);
          dmxy[profit->paramindex[PARAM_SPHEROID_POSANG]] = m20 * (1.0-a2)*dcs;
          }
        index = profit->paramindex[PARAM_SPHEROID_FLUX];
        break;
      case PROF_EXPONENTIAL:
        m20 = 3.0 * *prof->scale**prof->scale;
        if (jac)
          {
          dmx2[profit->paramindex[PARAM_DISK_SCALE]]
			= 6.0 * *prof->scale * mx2fac;
          dmy2[profit->paramindex[PARAM_DISK_SCALE]]
			= 6.0 * *prof->scale * my2fac;
          dmxy[profit->paramindex[PARAM_DISK_SCALE]]
			= 6.0 * *prof->scale * mxyfac;
          dmx2[profit->paramindex[PARAM_DISK_ASPECT]]
			= 2.0 * m20 * st*st * *prof->aspect;
          dmy2[profit->paramindex[PARAM_DISK_ASPECT]]
			= 2.0 * m20 * ct*ct * *prof->aspect;
          dmxy[profit->paramindex[PARAM_DISK_ASPECT]]
			= -2.0 * m20 * ct*st * *prof->aspect;
          dmx2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (dc2 + ds2*a2);
          dmy2[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (ds2 + dc2*a2);
          dmxy[profit->paramindex[PARAM_DISK_POSANG]] = m20 * (1.0 - a2) * dcs;
          }
        index = profit->paramindex[PARAM_DISK_FLUX];
        break;
      case PROF_ARMS:
        m20 = 1.0;
        index = profit->paramindex[PARAM_ARMS_FLUX];
        break;
      case PROF_BAR:
        m20 = 1.0;
        index = profit->paramindex[PARAM_BAR_FLUX];
        break;
      case PROF_INRING:
        m20 = 1.0;
        index = profit->paramindex[PARAM_INRING_FLUX];
        break;
      case PROF_OUTRING:
        m20 = 1.0;
        index = profit->paramindex[PARAM_OUTRING_FLUX];
        break;
      default:
        error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ",
		"prof_moments()");
      break;
      }

    prof->mx2 = m20*mx2fac;
    prof->my2 = m20*my2fac;
    prof->mxy = m20*mxyfac;
    
    }
  else
    prof->mx2 = prof->my2 = prof->mxy = 0.0;

  return index;
  }


/****** prof_interpolate ******************************************************
PROTO	float	prof_interpolate(profstruct *prof, float *posin)
PURPOSE	Interpolate a multidimensional model profile at a given position.
INPUT	Profile structure,
	input position vector.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	10/12/2006
 ***/
static float	prof_interpolate(profstruct *prof, float *posin)
  {
   float		dpos[2+PROFIT_MAXEXTRA],
			kernel_vector[INTERP_MAXKERNELWIDTH],
			*kvector, *pixin,*pixout,
			val;
   long			step[2+PROFIT_MAXEXTRA],
			start, fac;
   int			linecount[2+PROFIT_MAXEXTRA],