psf.c

              ppix = *(pd++)/pflux[j];
              ppix -=  psfmasks[j][p]**pw;
              chi2 += ppix*ppix;
              if (chi2>1E29) chi2=1E28;
            }
          obj2->chi2_psf = pflux[j]>0?
		chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;

        }
      
    }
  thepsfit->niter = thedpsfit->niter;
  thepsfit->npsf = npsf;

  for (j=0; j<npsf; j++)
    {
      thedpsfit->x[j] = ix+pdeltax[j]+1.0;
      thedpsfit->y[j] = iy+pdeltay[j]+1.0;
      thedpsfit->flux[j] = pflux[j];
      thedpsfit->fluxerr[j] = pfluxerr[j];
      thepsfit->x[j] = ix+pdeltax[j]+1.0;
      thepsfit->y[j] = iy+pdeltay[j]+1.0;
      thepsfit->flux[j] = pflux[j];
      thepsfit->fluxerr[j] = pfluxerr[j];
    }
    
  
  for (i=0; i<npsfmax; i++)
    {
      QFREE(psfmasks[i]);
      QFREE(psfmaskx[i]);
      QFREE(psfmasky[i]);
    }

  
  QFREE(psfmasks);
  QFREE(psfmaskx);
  QFREE(psfmasky);
  QFREE(pdatah);
  QFREE(pdata);
  QFREE(pdata2);
  QFREE(pdata3);
  QFREE(pweighth);
  QFREE(pweight);
  QFREE(pdata);
  QFREE(pmat);
   
  if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
      || prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
      || prefs.check[CHECK_PCOPROTOS])
    {
      QFREE(checkmask);
    }
  return;
}

/******************************* psf_build **********************************/
/*
Build the local PSF (function of "context").
*/
void	psf_build(psfstruct *psf)
  {
   static double	pos[POLY_MAXDIM];
   double	*basis, fac;
   float	*ppc, *pl;
   int		i,n,p, ndim, npix;

  if (psf->build_flag)
    return;

  npix = psf->masksize[0]*psf->masksize[1];

/* Reset the Local PSF mask */
  memset(psf->maskloc, 0, npix*sizeof(float));

/* Grab the context vector */
  ndim = psf->poly->ndim;
  for (i=0; i<ndim; i++)
    {
    ttypeconv(psf->contextindex[i]<0? (char *)psf->context[i]
		: *((char **)psf->context[i])
			+ psf->contextindex[i]*t_size[psf->contexttyp[i]],
		&pos[i], psf->contexttyp[i],T_DOUBLE);
    pos[i] = (pos[i] - psf->contextoffset[i]) / psf->contextscale[i];
    }
  poly_func(psf->poly, pos);

  basis = psf->poly->basis;

  ppc = psf->maskcomp;
/* Sum each component */
  for (n = (psf->maskdim>2?psf->masksize[2]:1); n--;)
    {
    pl = psf->maskloc;
    fac = *(basis++);
    for (p=npix; p--;)
      *(pl++) +=  fac**(ppc++);
    }

  psf->build_flag = 1;

  return;
  }


/******************************** psf_fwhm **********************************/
/*
Return the local PSF FWHM.
*/
double	psf_fwhm(psfstruct *psf)
  {
   float	*pl,
		val, max;
   int		n,p, npix;

  if (!psf->build_flag)
    psf_build(psf);

  npix = psf->masksize[0]*psf->masksize[1];
  max = -BIG;
  pl = psf->maskloc;
  for (p=npix; p--;)
    if ((val=*(pl++)) > max)
      max = val;
  pl = psf->maskloc;
  max /= 2.0;
  n = 0;
  for (p=npix; p--;)
    if (*(pl++) >= max)
      n++;

  return 2.0*sqrt(n/PI)*psf->pixstep;
  }


/*****************************compute_gradient*********************************/

double *compute_gradient(float *weight,int width, int height,
                         float *masks,float *maskx,float *masky
                        ,double *m)
{
  int x,y;
  float	*w, *ps,*px,*py;
    
  /*------ copy of the (weighted) PSF, with outer ring set to zero */
      ps = masks;
      w = weight;
      for (y=0; y<height; y++)
        for (x=0; x<width; x++, ps++, w++)
          *(m++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*ps**w):0)):0;
      /*------ (weighted) PSF gradient in x (kernel for first moment in x) */
      ps = masks;
      px = maskx;
      w = weight;
      for (y=0; y<height; y++)
        for (x=0; x<width; x++, ps++, w++)
          *(m++) = ((*px++) = (x?(x>=(width-1)?0:*(ps+1)-*(ps-1)):0))**w/2;
      /*------ (weighted) PSF gradient in y (kernel for first moment in y) */
      ps = masks; 
      py = masky;
      w = weight;
      for (y=0; y<height; y++)
        for (x=0; x<width; x++, ps++, w++)
          *(m++) = (*(py++)=(y?(y>=(height-1)?0:*(ps+width)-*(ps-width)):0))
            **w/2;
  return m;
}

/*****************************compute_gradient_phot*****************************
****/

double *compute_gradient_phot(float  *pweight,int width, int height,
                         float *pmasks,double *pm)

{
  int x,y;
  float  *pw, *pps;
    
  /*------ copy of the (weighted) PSF, with outer ring set to zero */
      pps = pmasks;
      pw = pweight;
      for (y=0; y<height; y++)
        for (x=0; x<width; x++, pps++, pw++)
          *(pm++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*pps**pw):0)):0;

  return pm;
}

/**************************compute_pos********************************/

void compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,double radmin2,
                         double radmax2,double r2,double *sol,double *flux 
                        ,double *deltax,double *deltay,double *pdx,double *pdy)
{
  int j,k,convflag,npsfflag,npsf; 
  double dx,dy;

  dx=*pdx;
  dy=*pdy;
  convflag=*pconvflag;
  npsfflag=*pnpsfflag;
  npsf=*pnpsf;
  for (j=0; j<npsf; j++)
    {
      flux[j] = sol[j*PSF_NA];
      /*------ Update the PSF shifts */
      if (fabs(flux[j])>0.0)
        {
          dx = -sol[j*PSF_NA+1]/((npsf>1?2:1)*flux[j]);
          dy = -sol[j*PSF_NA+2]/((npsf>1?2:1)*flux[j]);
        }
      
      deltax[j] += dx;
      deltay[j] += dy;
      /*------ Continue until all PSFs have come to a complete stop */
      if ((dx*dx+dy*dy) > radmin2)
        convflag = 1;
    }
  for (j=0; j<npsf; j++)
    {
      /*------ Exit if too much decentering or negative flux */
      for (k=j+1; k<npsf; k++)
        {
          dx = deltax[j]-deltax[k];
          dy = deltay[j]-deltay[k];
          if (dx*dx+dy*dy<r2/4.0)
            {
              flux[j] = -BIG;
              break;
            }
        }
      if (flux[j]<0.0
          || (deltax[j]*deltax[j] + deltay[j]*deltay[j]) > radmax2)
        {
          npsfflag = 0;
          convflag = 0;
          npsf--;
          break;
        }
    }
  *pdx=dx;
  *pdy=dy;
  *pconvflag=convflag;
  *pnpsfflag= npsfflag;
  *pnpsf=npsf;
  return;
}

/**************************compute_pos_phot********************************/

void compute_pos_phot(int *pnpsf,double *sol,double *flux)
{
  int j,npsf;   
  npsf=*pnpsf;
  for (j=0; j<npsf; j++)
    {
      flux[j] = sol[j];     
    }
  *pnpsf=npsf;
  return;
}


/************************************compute_poserr*****************************
*********/

void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
                    double *y2,double *xy, int npsf)
{
  double vara,covab,varb, f2;

  /*------ Variances and covariance along x and y */
  vara = *(var += (PSF_NA*npsf+1)*(j*PSF_NA+1));
  covab = *(++var);
  varb = *(var += PSF_NA*npsf);
  f2 = sol[PSF_NA*j];
  f2 *= f2;
  obj2->poserrmx2_psf = vara/f2;
  obj2->poserrmy2_psf = varb/f2;
  obj2->poserrmxy_psf = covab/f2;

  /*------ If requested, translate variances to major and minor error axes... */
  if (FLAG(obj2.poserra_psf))
    {
      double    pmx2,pmy2,temp,theta;
      
      if (fabs(temp=obj2->poserrmx2_psf-obj2->poserrmy2_psf) > 0.0)
        theta = atan2(2.0 * obj2->poserrmxy_psf,temp) / 2.0;
      else
        theta = PI/4.0;
      
      temp = sqrt(0.25*temp*temp+obj2->poserrmxy_psf*obj2->poserrmxy_psf);
      pmy2 = pmx2 = 0.5*(obj2->poserrmx2_psf+obj2->poserrmy2_psf);
      pmx2+=temp;
      pmy2-=temp;

      obj2->poserra_psf = (float)sqrt(pmx2);
      obj2->poserrb_psf = (float)sqrt(pmy2);
      obj2->poserrtheta_psf = theta*180.0/PI;
    }
  
  /*------ ...Or ellipse parameters */
  if (FLAG(obj2.poserr_cxx))
    {
      double    xm2,ym2, xym, temp;
      
      xm2 = obj2->poserrmx2_psf;
      ym2 = obj2->poserrmy2_psf;
      xym = obj2->poserrmxy_psf;
      obj2->poserrcxx_psf = (float)(ym2/(temp=xm2*ym2-xym*xym));
      obj2->poserrcyy_psf = (float)(xm2/temp);
      obj2->poserrcxy_psf = (float)(-2*xym/temp);
    }
  return;
}


/******************************** svdfit ************************************/
/*
General least-square fit A.x = b, based on Singular Value Decomposition (SVD).
Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
Note: the a and v matrices are transposed with respect to the N.R. convention.
*/
void svdfit(double *a, float *b, int m, int n, double *sol,
	double *vmat, double *wmat)
  {
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))
#define	PYTHAG(a,b)	((at=fabs(a)) > (bt=fabs(b)) ? \
				  (ct=bt/at,at*sqrt(1.0+ct*ct)) \
				: (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define	TOL		1.0e-11

   int			flag,i,its,j,jj,k,l,nm,mmi,nml;
   double		c,f,h,s,x,y,z,
			anorm, g, scale,
			at,bt,ct,maxarg1,maxarg2,
			thresh, wmax,
			*w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
			*tmpp, *rv1,*tmp;
   float		*bp;

  anorm = g = scale = 0.0;
  if (m < n)
    error(EXIT_FAILURE, "*Error*: Not enough rows for solving the system ",
	"in svdfit()");
  
  QMALLOC(rv1, double, n);
  QMALLOC(tmp, double, n);
  l = nm = nml = 0;			/* To avoid gcc -Wall warnings */
  for (i=0;i<n;i++)
    {
    l = i+1;
    nml = n-l;
    rv1[i] = scale*g;
    g = s = scale = 0.0;
    if ((mmi = m - i) > 0)
      {
      ap = ap0 = a+i*(m+1);
      for (k=mmi;k--;)
        scale += fabs(*(ap++));
      if (scale)
        {
        for (ap=ap0,k=mmi; k--; ap++)
          {
          *ap /= scale;
          s += *ap**ap;
          }
        f = *ap0;
        g = -SIGN(sqrt(s),f);
        h = f*g-s;
        *ap0 = f-g;
        ap10 = a+l*m+i;
        for (j=nml;j--; ap10+=m)
          {
          for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
            s += *(ap1++)**(ap++);
          f = s/h;
          for (ap=ap0,ap1=ap10,k=mmi; k--;)
            *(ap1++) += f**(ap++);
          }
        for (ap=ap0,k=mmi; k--;)
          *(ap++) *= scale;
        }
      }
    wmat[i] = scale*g;
    g = s = scale = 0.0;
    if (i < m && i+1 != n)
      {
      ap = ap0 = a+i+m*l;
      for (k=nml;k--; ap+=m)
        scale += fabs(*ap);
      if (scale)
        {
        for (ap=ap0,k=nml;k--; ap+=m)
          {
          *ap /= scale;
          s += *ap**ap;
          }
        f=*ap0;
        g = -SIGN(sqrt(s),f);
        h=f*g-s;
        *ap0=f-g;
        rv1p = rv1+l;
        for (ap=ap0,k=nml;k--; ap+=m)
          *(rv1p++) = *ap/h;
        ap10 = a+l+m*l;
        for (j=m-l; j--; ap10++)
          {
          for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
            s += *ap1**ap;
          rv1p = rv1+l;
          for (ap1=ap10,k=nml;k--; ap1+=m)
            *ap1 += s**(rv1p++);
          }
        for (ap=ap0,k=nml;k--; ap+=m)
          *ap *= scale;
        }
      }
    anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
    }

  for (i=n-1;i>=0;i--)
    {
    if (i < n-1)
      {
      if (g)
        {
        ap0 = a+l*m+i;
        vp0 = vmat+i*n+l;
        vp10 = vmat+l*n+l;
        g *= *ap0;
        for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
          *(vp++) = *ap/g;
        for (j=nml; j--; vp10+=n)
          {
          for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
            s += *ap**(vp1++);
          for (vp=vp0,vp1=vp10,k=nml; k--;)
            *(vp1++) += s**(vp++);
          }
        }
      vp = vmat+l*n+i;
      vp1 = vmat+i*n+l;
      for (j=nml; j--; vp+=n)
        *vp = *(vp1++) = 0.0;
      }
    vmat[i*n+i]=1.0;
    g=rv1[i];
    l=i;
    nml = n-l;
    }

  for (i=(m<n?m:n); --i>=0;)
    {
    l=i+1;
    nml = n-l;
    mmi=m-i;
    g=wmat[i];
    ap0 = a+i*m+i;
    ap10 = ap0 + m;
    for (ap=ap10,j=nml;j--;ap+=m)
      *ap=0.0;
    if (g)
      {
      g=1.0/g;
      for (j=nml;j--; ap10+=m)
        {
        for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
              s += *(++ap)**(++ap1);
        f = (s/(*ap0))*g;
        for (ap=ap0,ap1=ap10,k=mmi;k--;)
          *(ap1++) += f**(ap++);
        }
      for (ap=ap0,j=mmi;j--;)
        *(ap++) *= g;
      }
    else
      for (ap=ap0,j=mmi;j--;)
        *(ap++)=0.0;
    ++(*ap0);
    }

  for (k=n; --k>=0;)
      {
      for (its=0;its<100;its++)
        {
        flag=1;
        for (l=k;l>=0;l--)
          {
          nm=l-1;
          if (fabs(rv1[l])+anorm == anorm)
            {
            flag=0;
            break;
            }
          if (fabs(wmat[nm])+anorm == anorm)
            break;
          }
        if (flag)
          {
          c=0.0;
          s=1.0;
          ap0 = a+nm*m;
          ap10 = a+l*m;
          for (i=l; i<=k; i++,ap10+=m)
            {
            f=s*rv1[i];
            if (fabs(f)+anorm == anorm)
              break;
            g=wmat[i];
            h=PYTHAG(f,g);
            wmat[i]=h;
            h=1.0/h;
            c=g*h;
            s=(-f*h);
            for (ap=ap0,ap1=ap10,j=m; j--;)
              {
              z = *ap1;
              y = *ap;
              *(ap++) = y*c+z*s;
              *(ap1++) = z*c-y*s;
              }
            }
          }
        z=wmat[k];
        if (l == k)
          {
          if (z < 0.0)
            {
            wmat[k] = -z;
            vp = vmat+k*n;
            for (j=n; j--; vp++)
              *vp = (-*vp);
            }
          break;
          }
        if (its == 99)
          error(EXIT_FAILURE, "*Error*: No convergence in 100 SVD iterations ",
		"in svdfit()");
        x=wmat[l];
        nm=k-1;
        y=wmat[nm];
        g=rv1[nm];
        h=rv1[k];
        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
        g=PYTHAG(f,1.0);
        f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
        c=s=1.0;
        ap10 = a+l*m;
        vp10 = vmat+l*n;
        for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
          {
          i=j+1;
          g=rv1[i];
          y=wmat[i];
          h=s*g;
          g=c*g;
          z=PYTHAG(f,h);
          rv1[j]=z;
          c=f/z;
          s=h/z;
          f=x*c+g*s;
          g=g*c-x*s;
          h=y*s;
          y=y*c;
          for (vp=(vp1=vp10)+n,jj=n; jj--;)
            {
            z = *vp;
            x = *vp1;
            *(vp1++) = x*c+z*s;
            *(vp++) = z*c-x*s;
            }
          z=PYTHAG(f,h);
          wmat[j]=z;
          if (z)
            {
            z=1.0/z;
            c=f*z;
            s=h*z;
            }
          f=c*g+s*y;
          x=c*y-s*g;
          for (ap=(ap1=ap10)+m,jj=m; jj--;)
            {
            z = *ap;
            y = *ap1;
            *(ap1++) = y*c+z*s;
            *(ap++) = z*c-y*s;
            }
          }
        rv1[l]=0.0;
        rv1[k]=f;
        wmat[k]=x;
        }
      }

  wmax=0.0;
  w = wmat;
  for (j=n;j--; w++)
    if (*w > wmax)
      wmax=*w;
  thresh=TOL*wmax;
  w = wmat;
  for (j=n;j--; w++)
    if (*w < thresh)
      *w = 0.0;

  w = wmat;
  ap = a;
  tmpp = tmp;
  for (j=n; j--; w++)
    {
    s=0.0;
    if (*w)
      {
      bp = b;
      for (i=m; i--;)
        s += *(ap++)**(bp++);
      s /= *w;
      }
    else
      ap += m;
    *(tmpp++) = s;
    }

  vp0 = vmat;
  for (j=0; j<n; j++,vp0++)
    {
    s=0.0;
    tmpp = tmp;
    for (vp=vp0,jj=n; jj--; vp+=n)
      s += *vp**(tmpp++);
    sol[j]=s;
    }

/* Free temporary arrays */
  free(tmp);
  free(rv1);

  return;
  }

#undef SIGN
#undef MAX
#undef PYTHAG
#undef TOL

/******************************** svdvar ************************************/
/*
Computation of the covariance matrix from the SVD vmat and wmat matrices.A
dapted from Numerical Recipes in C, 2nd Ed. (p. 679).
*/
void svdvar(double *v, double *w, int n, double *cov)
  {
   static double	wti[PSF_NTOT];
   double		sum;
   int			i,j,k;

  for (i=0; i<n; i++)
    wti[i] = w[i]? 1.0/(w[i]*w[i]) : 0.0;

  for (i=0; i<n; i++)
    for (j=0; j<=i; j++)
      {
      for (sum=0.0,k=0; k<n; k++)
        sum += v[k*n+i]*v[k*n+j]*wti[k];
      cov[j*n+i] = cov[i*n+j] = sum;
      }

  return;
  }