profit.c

/* Don't go further if out of frame!! */
  ix = profit->ix;
  iy = profit->iy;
  if (ix<0 || ix>=field->width || iy<field->ymin || iy>=field->ymax)
    return 0;

  backnoise2 = field->backsig*field->backsig;
  sn = (int)profit->subsamp;
  sflag = (sn>1);
  w = profit->objnaxisn[0]*sn;
  h = profit->objnaxisn[1]*sn;
  if (sflag)
    backnoise2 *= (PIXTYPE)sn;
  invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0;
  satlevel = field->satur_level - profit->obj->bkg;
  rad2 = h/2.0;
  if (rad2 > w/2.0)
    rad2 = w/2.0;
  rad2 *= rad2;

/* Set the image boundaries */
  pixout = profit->objpix;
  wpixout = profit->objweight;
  ymin = iy-h/2;
  ymax = ymin + h;
  if (ymin<field->ymin)
    {
    off = (field->ymin-ymin-1)/sn + 1;
    pixout += off*profit->objnaxisn[0];
    wpixout += off*profit->objnaxisn[0];
    ymin += off*sn;
    }
  if (ymax>field->ymax)
    ymax -= ((ymax-field->ymax-1)/sn + 1)*sn;

  xmin = ix-w/2;
  xmax = xmin + w;
  dw = 0;
  if (xmax>field->width)
    {
    off = (xmax-field->width-1)/sn + 1;
    dw += off;
    xmax -= off*sn;
    }
  if (xmin<0)
    {
    off = (-xmin-1)/sn + 1;
    pixout += off;
    wpixout += off;
    dw += off;
    xmin += off*sn;
    }
/* Make sure the input frame size is a multiple of the subsampling step */
  if (sflag)
    {
/*
    if (((rem=ymax-ymin)%sn))
      {
      ymin += rem/2;
      ymax -= (rem-rem/2);
      }
    if (((rem=xmax-xmin)%sn))
      {
      xmin += rem/2;
      pixout += rem/2;
      wpixout += rem/2;
      dw += rem;
      xmax -= (rem-rem/2);
      }
*/
    sw = field->width;
    }

/* Copy the right pixels to the destination */
  npix = 0;
  if (wfield)
    {
    wthresh = wfield->weight_thresh;
    gainflag = prefs.weightgain_flag;
    if (sflag)
      {
/*---- Sub-sampling case */
      for (y=ymin; y<ymax; y+=sn, pixout+=dw,wpixout+=dw)
        {
        for (x=xmin; x<xmax; x+=sn)
          {
          pix = wpix = 0.0;
          badflag = 0;
          for (sy=0; sy<sn; sy++)
            {
            dy2 = (y+sy-iy);
            dy2 *= dy2;
            dx = (x-ix);
            spixin = &PIX(field, x, y+sy);
            swpixin = &PIX(wfield, x, y+sy);
            for (sx=sn; sx--;)
              {
              dr2 = dy2 + dx*dx;
              dx++;
              spix = *(spixin++);
              swpix = *(swpixin++);
              if (dr2<rad2 && spix>-BIG && spix<satlevel && swpix<wthresh)
                {
                pix += spix;
                wpix += swpix;
                }
              else
                badflag=1;
              }
            }
          *(pixout++) = pix;
          if (!badflag)	/* A single bad pixel ruins is all (saturation, etc.)*/
            {
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
		(gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
            npix++;
            }
          else
            *(wpixout++) = 0.0;
          }
        }
      }
    else
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
        {
        dy2 = y-iy;
        dy2 *= dy2;
        pixin = &PIX(field, xmin, y);
        wpixin = &PIX(wfield, xmin, y);
        for (x=xmin; x<xmax; x++)
          {
          dx = x-ix;
          dr2 = dy2 + dx*dx;
          pix = *(pixin++);
          wpix = *(wpixin++);
          if (dr2<rad2 && pix>-BIG && pix<satlevel && wpix<wthresh)
            {
            *(pixout++) = pix;
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
		(gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
            npix++;
            }
          else
            *(pixout++) = *(wpixout++) = 0.0;
          }
        }
    }
  else
    {
    if (sflag)
      {
/*---- Sub-sampling case */
      for (y=ymin; y<ymax; y+=sn, pixout+=dw, wpixout+=dw)
        {
        for (x=xmin; x<xmax; x+=sn)
          {
          pix = 0.0;
          badflag = 0;
          for (sy=0; sy<sn; sy++)
            {
            dy2 = y+sy-iy;
            dy2 *= dy2;
            dx = x-ix;
            spixin = &PIX(field, x, y+sy);
            for (sx=sn; sx--;)
              {
              dr2 = dy2 + dx*dx;
              dx++;
              spix = *(spixin++);
              if (dr2<rad2 && spix>-BIG && spix<satlevel)
                pix += spix;
              else
                badflag=1;
              }
            }
          *(pixout++) = pix;
          if (!badflag)	/* A single bad pixel ruins is all (saturation, etc.)*/
            {
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain:0.0));
            npix++;
            }
          else
            *(wpixout++) = 0.0;
          }
        }
      }
    else
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
        {
        dy2 = y-iy;
        dy2 *= dy2;
        pixin = &PIX(field, xmin, y);
        for (x=xmin; x<xmax; x++)
          {
          dx = x-ix;
          dr2 = dy2 + dx*dx;
          pix = *(pixin++);
          if (dr2<rad2 && pix>-BIG && pix<satlevel)
            {
            *(pixout++) = pix;
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain : 0.0));
            npix++;
            }
          else
            *(pixout++) = *(wpixout++) = 0.0;
          }
        }
    }
 
  return npix;
  }


/****** profit_spiralindex ****************************************************
PROTO	float profit_spiralindex(profitstruct *profit)
PURPOSE	Compute the spiral index of a galaxy image (positive for arms
	extending counter-clockwise and negative for arms extending CW, 0 for
	no spiral pattern).
INPUT	Profile-fitting structure.
OUTPUT	Vector of residuals.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	18/05/2011
 ***/
float profit_spiralindex(profitstruct *profit)
  {
   objstruct	*obj;
   obj2struct	*obj2;
   float	*dx,*dy, *fdx,*fdy, *gdx,*gdy, *gdxt,*gdyt, *pix,
		fwhm, invtwosigma2, hw,hh, ohw,ohh, x,y,xstart, tx,ty,txstart,
		gx,gy, r2, spirindex, invsig, val, sep;
   PIXTYPE	*fpix;
   int		i,j, npix;

  npix = profit->objnaxisn[0]*profit->objnaxisn[1];

  obj = profit->obj;
  obj2 = profit->obj2;
/* Compute simple derivative vectors at a fraction of the object scale */
  fwhm = profit->guessradius * 2.0 / 4.0;
  if (fwhm < 2.0)
    fwhm = 2.0;
  sep = 2.0;

  invtwosigma2 = -(2.35*2.35/(2.0*fwhm*fwhm));
  hw = (float)(profit->objnaxisn[0]/2);
  ohw = profit->objnaxisn[0] - hw;
  hh = (float)(profit->objnaxisn[1]/2);
  ohh = profit->objnaxisn[1] - hh;
  txstart = -hw;
  ty = -hh;
  QMALLOC(dx, float, npix);
  pix = dx;
  for (j=profit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=profit->objnaxisn[0]; i--; tx+=1.0)
      {
      x = tx < -0.5? tx + hw : tx - ohw;
      *(pix++) = exp(invtwosigma2*((x+sep)*(x+sep)+y*y))
		- exp(invtwosigma2*((x-sep)*(x-sep)+y*y));
      }
    }
  QMALLOC(dy, float, npix);
  pix = dy;
  ty = -hh;
  for (j=profit->objnaxisn[1]; j--; ty+=1.0)
    {
    tx = txstart;
    y = ty < -0.5? ty + hh : ty - ohh;
    for (i=profit->objnaxisn[0]; i--; tx+=1.0)
      {
      x = tx < -0.5? tx + hw : tx - ohw;
      *(pix++) = exp(invtwosigma2*(x*x+(y+sep)*(y+sep)))
		- exp(invtwosigma2*(x*x+(y-sep)*(y-sep)));
      }
    }

  QMALLOC(gdx, float, npix);
  gdxt = gdx;
  fpix = profit->objpix;
  invsig = npix/profit->sigma;
  for (i=npix; i--; fpix++)
    {
    val = *fpix > -1e29? *fpix*invsig : 0.0;
    *(gdxt++) = (val>0.0? log(1.0+val) : -log(1.0-val));
    }
  gdy = NULL;			/* to avoid gcc -Wall warnings */
  QMEMCPY(gdx, gdy, float, npix);
  fdx = fft_rtf(dx, profit->objnaxisn);
  fft_conv(gdx, fdx, profit->objnaxisn);
  fdy = fft_rtf(dy, profit->objnaxisn);
  fft_conv(gdy, fdy, profit->objnaxisn);

/* Compute estimator */
  invtwosigma2 = -1.18*1.18 / (2.0*profit->guessradius*profit->guessradius);
  xstart = -hw - obj->mx + (int)(obj->mx+0.49999);
  y = -hh -  obj->my + (int)(obj->my+0.49999);;
  spirindex = 0.0;
  gdxt = gdx;
  gdyt = gdy;
  for (j=profit->objnaxisn[1]; j--; y+=1.0)
    {
    x = xstart;
    for (i=profit->objnaxisn[0]; i--; x+=1.0)
      {
      gx = *(gdxt++);
      gy = *(gdyt++);
      if ((r2=x*x+y*y)>0.0)
        spirindex += (x*y*(gx*gx-gy*gy)+gx*gy*(y*y-x*x))/r2
			* exp(invtwosigma2*r2);
      }
    }

  free(dx);
  free(dy);
  free(fdx);
  free(fdy);
  free(gdx);
  free(gdy);

  return spirindex;
  }


/****** profit_moments ****************************************************
PROTO	void profit_moments(profitstruct *profit, obj2struct *obj2)
PURPOSE	Compute the 2nd order moments from the unconvolved object model.
INPUT	Profile-fitting structure,
	Pointer to obj2 structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	22/04/2011
 ***/
void	 profit_moments(profitstruct *profit, obj2struct *obj2)
  {
   profstruct	*prof;
   double	dpdmx2[6], cov[4],
		*jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
		*dmx2,*dmy2,*dmxy,
		m0,invm0, mx2,my2,mxy, den,invden,
		temp, temp2,invtemp2,invstemp2,
		pmx2,theta, flux, dval;
   float	 *covart;
   int		findex[MODEL_NMAX],
		i,j,p, nparam;

/*  hw = (float)(profit->modnaxisn[0]/2);*/
/*  hh = (float)(profit->modnaxisn[1]/2);*/
/*  r2max = hw<hh? hw*hw : hh*hh;*/
/*  xstart = -hw;*/
/*  y = -hh;*/
/*  pix = profit->modpix;*/
/*  mx2 = my2 = mxy = mx = my = sum = 0.0;*/
/*  for (iy=profit->modnaxisn[1]; iy--; y+=1.0)*/
/*    {*/
/*    x = xstart;*/
/*    for (ix=profit->modnaxisn[0]; ix--; x+=1.0)*/
/*      if (y*y+x*x <= r2max)*/
/*        {*/
/*        val = *(pix++);*/
/*        sum += val;*/
/*        mx  += val*x;*/
/*        my  += val*y;*/
/*        mx2 += val*x*x;*/
/*        mxy += val*x*y;*/
/*        my2 += val*y*y;*/
/*        }*/
/*      else*/
/*        pix++;*/
/*    }*/

/*  if (sum <= 1.0/BIG)*/
/*    sum = 1.0;*/
/*  mx /= sum;*/
/*  my /= sum;*/
/*  obj2->prof_mx2 = mx2 = mx2/sum - mx*mx;*/
/*  obj2->prof_my2 = my2 = my2/sum - my*my;*/
/*  obj2->prof_mxy = mxy = mxy/sum - mx*my;*/

  nparam = profit->nparam;
  if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
    {
/*-- Set up Jacobian matrices */
    QCALLOC(jac, double, nparam*3);
    QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
    QMALLOC(dcovar, double, nparam*nparam);
    dcovart = dcovar;
    covart = profit->covar;
    for (i=nparam*nparam; i--;)
      *(dcovart++) = (double)(*(covart++));
    dmx2 = jac;
    dmy2 = jac+nparam;
    dmxy = jac+2*nparam;
    }
  else
    jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;

  m0 = mx2 = my2 = mxy = 0.0;
  for (p=0; p<profit->nprof; p++)
    {
    prof = profit->prof[p];
    findex[p] = prof_moments(profit, prof, pjac);
    flux = *prof->flux;
    m0 += flux;
    mx2 += prof->mx2*flux;
    my2 += prof->my2*flux;
    mxy += prof->mxy*flux;
    if (jac)
      {
      jact = jac;
      pjact = pjac;
      for (j=nparam*3; j--;)
        *(jact++) += flux * *(pjact++);
      }
    }
  invm0 = 1.0 / m0;
  obj2->prof_mx2 = (mx2 *= invm0);
  obj2->prof_my2 = (my2 *= invm0);
  obj2->prof_mxy = (mxy *= invm0);
/* Complete the flux derivative of moments */
  if (jac)
    {
    for (p=0; p<profit->nprof; p++)
      {
      prof = profit->prof[p];
      dmx2[findex[p]] = prof->mx2 - mx2;
      dmy2[findex[p]] = prof->my2 - my2;
      dmxy[findex[p]] = prof->mxy - mxy;
      }
    jact = jac;
    for (j=nparam*3; j--;)
      *(jact++) *= invm0;
    }

/* Handle fully correlated profiles (which cause a singularity...) */
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
    {
    mx2 += 0.0833333;
    my2 += 0.0833333;
    temp2 = mx2*my2-mxy*mxy;
    }

/* Use the Jacobians to compute the moment covariance matrix */
  if (jac)
    propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
						pjac);	/* We re-use pjac */

  if (FLAG(obj2.prof_pol1))
    {
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
    if (mx2+my2 > 1.0/BIG)
      {
      obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
      obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
      if (FLAG(obj2.prof_pol1err))
        {
/*------ Compute the Jacobian of polarisation */
        invden = 1.0/(mx2+my2);
        dpdmx2[0] =  2.0*my2*invden*invden;
        dpdmx2[1] = -2.0*mx2*invden*invden;
        dpdmx2[2] =  0.0;
        dpdmx2[3] = -2.0*mxy*invden*invden;
        dpdmx2[4] = -2.0*mxy*invden*invden;
        dpdmx2[5] =  2.0*invden;

/*------ Use the Jacobian to compute the polarisation covariance matrix */
        propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
						pjac);	/* We re-use pjac */
        obj2->prof_pol1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
        obj2->prof_pol2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
        obj2->prof_pol12corr = (dval=cov[0]*cov[3]) > 0.0?
					(float)(cov[1]/sqrt(dval)) : 0.0;
        }
      }
    else
      obj2->prof_pol1 = obj2->prof_pol2
	= obj2->prof_pol1err = obj2->prof_pol2err = obj2->prof_pol12corr = 0.0;
    }

  if (FLAG(obj2.prof_e1))
    {
/*--- "Ellipticity", i.e. module = (a-b)/(a+b) */
    if (mx2+my2 > 1.0/BIG)
      {
      den = (temp2>=0.0) ? mx2+my2+2.0*sqrt(temp2) : mx2+my2;
      invden = 1.0/den;
      obj2->prof_e1 = (float)(invden * (mx2 - my2));
      obj2->prof_e2 = (float)(2.0 * invden * mxy);
      if (FLAG(obj2.prof_e1err))
        {
/*------ Compute the Jacobian of ellipticity */
        invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0;
        dpdmx2[0] = ( den - (1.0+my2*invstemp2)*(mx2-my2))*invden*invden;
        dpdmx2[1] = (-den - (1.0+mx2*invstemp2)*(mx2-my2))*invden*invden;
        dpdmx2[2] = 2.0*mxy*invstemp2*(mx2-my2)*invden*invden;
        dpdmx2[3] = -2.0*mxy*(1.0+my2*invstemp2)*invden*invden;
        dpdmx2[4] = -2.0*mxy*(1.0+mx2*invstemp2)*invden*invden;
        dpdmx2[5] =  (2.0*den+4.0*mxy*mxy*invstemp2)*invden*invden;

/*------ Use the Jacobian to compute the ellipticity covariance matrix */
        propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
					pjac);	/* We re-use pjac */
        obj2->prof_e1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
        obj2->prof_e2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
        obj2->prof_e12corr = (dval=cov[0]*cov[3]) > 0.0?
					(float)(cov[1]/sqrt(dval)) : 0.0;
        }
      }
    else
      obj2->prof_e1 = obj2->prof_e2
	= obj2->prof_e1err = obj2->prof_e2err = obj2->prof_e12corr = 0.0;
    }

  if (FLAG(obj2.prof_cxx))
    {
    invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
    obj2->prof_cxx = (float)(my2*invtemp2);
    obj2->prof_cyy = (float)(mx2*invtemp2);
    obj2->prof_cxy = (float)(-2*mxy*invtemp2);
    }

  if (FLAG(obj2.prof_a))
    {
    if ((fabs(temp=mx2-my2)) > 0.0)
      theta = atan2(2.0 * mxy,temp) / 2.0;
    else
      theta = PI/4.0;

    temp = sqrt(0.25*temp*temp+mxy*mxy);
    pmx2 = 0.5*(mx2+my2);
    obj2->prof_a = (float)sqrt(pmx2 + temp);
    obj2->prof_b = (float)sqrt(pmx2 - temp);
    obj2->prof_theta = theta*180.0/PI;
    }

/* Free memory used by Jacobians */
  free(jac);
  free(pjac);
  free(dcovar);

  return;
  }


/****** profit_convmoments ****************************************************
PROTO	void profit_convmoments(profitstruct *profit, obj2struct *obj2)
PURPOSE	Compute the 2nd order moments of the convolved object model.
INPUT	Profile-fitting structure,
	Pointer to obj2 structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	12/04/2011
 ***/
void	 profit_convmoments(profitstruct *profit, obj2struct *obj2)
  {
   double	hw,hh, r2max, x,xstart,y, mx2,my2,mxy,mx,my,sum, dval,
		temp,temp2,invtemp2, pmx2, theta;
   PIXTYPE	*pix;
   int		ix,iy, w,h;

  w = profit->modnaxisn[0];
  h = profit->modnaxisn[1];
  hw = (double)(w/2);
  hh = (double)(h/2);

  r2max = hw<hh? hw*hw : hh*hh;
  xstart = -hw;
  y = -hh;
  pix = profit->cmodpix;
  mx2 = my2 = mxy = mx = my = sum = 0.0;
  for (iy=h; iy--; y+=1.0)
    {
    x = xstart;
    for (ix=w; ix--; x+=1.0)
      if (y*y+x*x <= r2max)
        {
        dval = *(pix++);
        sum += dval;
        mx  += dval*x;
        my  += dval*y;
        mx2 += dval*x*x;
        mxy += dval*x*y;
        my2 += dval*y*y;
        }
      else
        pix++;
    }

  if (sum <= 1.0/BIG)
    sum = 1.0;
  mx /= sum;
  my /= sum;
  obj2->prof_convmx2 = (mx2 = mx2/sum - mx*mx)*profit->pixstep*profit->pixstep;
  obj2->prof_convmy2 = (my2 = my2/sum - my*my)*profit->pixstep*profit->pixstep;
  obj2->prof_convmxy = (mxy = mxy/sum - mx*my)*profit->pixstep*profit->pixstep;

/* Handle fully correlated profiles (which cause a singularity...) */
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
    {
    mx2 += 0.0833333;
    my2 += 0.0833333;
    temp2 = mx2*my2-mxy*mxy;
    }

  temp2 *= profit->pixstep*profit->pixstep;

  if (FLAG(obj2.prof_convcxx))
    {
    invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
    obj2->prof_convcxx = (float)(my2*invtemp2);
    obj2->prof_convcyy = (float)(mx2*invtemp2);
    obj2->prof_convcxy = (float)(-2*mxy*invtemp2);
    }

  if (1 /*FLAG(obj2.prof_conva)*/)
    {
    if ((fabs(temp=mx2-my2)) > 0.0)
      theta = atan2(2.0 * mxy,temp) / 2.0;
    else
      theta = PI/4.0;

    temp = sqrt(0.25*temp*temp+mxy*mxy);
    pmx2 = 0.5*(mx2+my2);
    obj2->prof_conva = (float)sqrt(pmx2 + temp)*profit->pixstep;
    obj2->prof_convb = (float)sqrt(pmx2 - temp)*profit->pixstep;
    obj2->prof_convtheta = theta/DEG;
    }

  return;
  }


/****** profit_surface ****************************************************
PROTO	void profit_surface(profitstruct *profit, obj2struct *obj2)
PURPOSE	Compute surface brightnesses from the unconvolved object model.
INPUT	Pointer to the profile-fitting structure,
	Pointer to obj2 structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	24/01/2011
 ***/
void	 profit_surface(profitstruct *profit, obj2struct *obj2)
  {
   profitstruct	hdprofit;
   double	dsum,dhsum,dsumoff, dhval, frac, seff;
   float	*spix, *spixt,
		val,vmax,
		scalefac, imsizefac, flux, lost, sum, lostfluxfrac;
   int		i,p, imax, npix, neff;

/* Allocate "high-definition" raster only to make measurements */
  hdprofit.modnaxisn[0] = hdprofit.modnaxisn[1] = PROFIT_HIDEFRES;
  npix = hdprofit.nmodpix = hdprofit.modnaxisn[0]*hdprofit.modnaxisn[1];
/* 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 = scalefac*scalefac;
  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];
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));


int r,t;
double ratio,ratio0,ang,ang0, x,x0,y,y0;
list = profit->paramlist;
index = profit->paramindex;
for (p=0; p<nparam; p++)
param[p] = profit->paraminit[p];

ratio0 = profit->paraminit[index[PARAM_SPHEROID_ASPECT]];
ang0 = profit->paraminit[index[PARAM_SPHEROID_POSANG]];
x0 = profit->paraminit[index[PARAM_X]];
y0 = profit->paraminit[index[PARAM_Y]];
for (r=0;r<check->height;r++)
for (t=0; t<check->width;t++)
{
//x = (r-10.0)/100.0 + x0;
//y = (t-10.0)/100.0 + y0;
ratio = ratio0*exp((r-10.0)/400.0);
ang = ang0+(t-10.0)/3.0;

for (i=0; i<PARAM_NPARAM; i++)
{
//if (list[i] && i==PARAM_X)
//param[index[i]] = x;
//if (list[i] && i==PARAM_Y)
//param[index[i]] = y;
if (list[i] && i==PARAM_SPHEROID_ASPECT)
param[index[i]] = ratio;
if (list[i] && i==PARAM_SPHEROID_POSANG)
param[index[i]] = ang;
//if (list[i] && i==PARAM_SPHEROID_REFF)
//param[index[i]] = profit->paraminit[index[i]]*sqrt(ratio0/ratio);
}
profit_residuals(profit,field,wfield, PROFIT_DYNPARAM, param,profit->resi);
*((float *)check->pix + t + r*check->width) = profit->chi2;
}
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 (overwrites oversampled model!!) */
  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;
/*-- Recompute hi-def model raster without oversampling */
/*-- and with the same flux correction factor */
    memset(hdprofit.modpix,0, npix*sizeof(float));
    for (p=0; p<profit->nprof; p++)
      prof_add(&hdprofit, profit->prof[p], 1);
    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	18/05/2011
 ***/
void	profit_resetparam(profitstruct *profit, paramenum paramtype)
  {
   objstruct	*obj;
   obj2struct	*obj2;
   float	param, parammin,parammax, range;
   parfitenum	fittype;

  obj = profit->obj;
  obj2 = profit->obj2;
  param = parammin = parammax = 0.0;	/* Avoid gcc -Wall warnings*/

  switch(paramtype)
    {
    case PARAM_BACK:
      fittype = PARFIT_LINBOUND;
      param = 0.0;
      parammin = -6.0*obj->sigbkg;
      parammax =  6.0*obj->sigbkg;
      break;
    case PARAM_X:
      fittype = PARFIT_LINBOUND;
      param = obj->mx - (int)(obj->mx+0.49999);
      range = profit->guessradius*4.0;
      if (range>profit->objnaxisn[0]*2.0)
        range = profit->objnaxisn[0]*2.0;
      parammin = -range;
      parammax =  range;
      break;
    case PARAM_Y:
      fittype = PARFIT_LINBOUND;
      param = obj->my - (int)(obj->my+0.49999);
      range = profit->guessradius*4.0;
      if (range>profit->objnaxisn[1]*2)
        range = profit->objnaxisn[1]*2;
      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(obj->a/obj->b);
      parammin = 0.01;
      parammax = param * 10.0;
      break;
    case PARAM_SPHEROID_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = FLAG(obj2.prof_disk_flux)? 1.0 : obj->b/obj->a;
      parammin = FLAG(obj2.prof_disk_flux)? 0.5 : 0.01;
      parammax = FLAG(obj2.prof_disk_flux)? 2.0 : 100.0;
      break;
    case PARAM_SPHEROID_POSANG:
      fittype = PARFIT_UNBOUND;
      param = obj->theta;
      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;
      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(obj->a/obj->b);
      parammin = 0.01/1.67835;
      parammax = param * 10.0;
      break;
    case PARAM_DISK_ASPECT:
      fittype = PARFIT_LOGBOUND;
      param = obj->b/obj->a;
      parammin = 0.01;
      parammax = 100.0;
      break;
    case PARAM_DISK_POSANG:
      fittype = PARFIT_UNBOUND;
      param = obj->theta;
      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;