text: /*** File libwcs/matchstar.c
   *** January 11, 2007
   *** By Doug Mink, dmink@cfa.harvard.edu
   *** Harvard-Smithsonian Center for Astrophysics
   *** Copyright (C) 1996-2007
   *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
  
      This library is free software; you can redistribute it and/or
      modify it under the terms of the GNU Lesser General Public
      License as published by the Free Software Foundation; either
      version 2 of the License, or (at your option) any later version.
  
      This library is distributed in the hope that it will be useful,
      but WITHOUT ANY WARRANTY; without even the implied warranty of
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      Lesser General Public License for more details.
      
      You should have received a copy of the GNU Lesser General Public
      License along with this library; if not, write to the Free Software
      Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  
      Correspondence concerning WCSTools should be addressed as follows:
             Internet email: dmink@cfa.harvard.edu
             Postal address: Doug Mink
                             Smithsonian Astrophysical Observatory
                             60 Garden St.
                             Cambridge, MA 02138 USA
   */
  
  /* StarMatch (ns, sx, sy, ng, gra, gdec, goff, gx, gy, tol, wcs, nfit, debug)
   *  Find shift, scale, and rotation of image stars to best-match reference stars
   *
   * ReadMatch (filename, sx, sy, gra, gdec, debug)
   *  Read in x, y, RA, and Dec of pre-match stars in image
   *
   * WCSMatch (nmatch, sbx, sby, gbra, gbdec, debug)
   *  Find shift, scale, and rotation of image stars to best-match reference stars
   *
   * FitMatch (ns, sx, sy, ng, gra, gdec, gx, gy, tol, wcs, nfit, debug)
   *  Fit shift, scale, and rotation of image stars to RA/Dec/X/Y matches
   *
   * wcs_amoeba (wcs0) Set up temp arrays and call multivariate solver
   * chisqr (v) Compute the chisqr of the vector v
   * amoeba (p, y, ndim, ftol, itmax, funk, nfunk)
   *    Multivariate solver from Numerical Recipes
   */
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <string.h>
  #include "wcs.h"
  #include "lwcs.h"
  #include "fitsfile.h"
  #include "wcscat.h"
  
  #define NPAR 8
  #define NPAR1 9
  
  #define ABS(a) ((a) < 0 ? (-(a)) : (a))
  
  static void wcs_amoeba ();
  extern void setnofit();
  extern int getfilelines();
  
  /* Statics used by the chisqr evaluator */
  static double	*sx_p;
  static double	*sy_p;
  static double	*gra_p;
  static double	*gdec_p;
  static double	xref_p, yref_p;
  static double	xrefpix, yrefpix;
  static int	nbin_p;
  static int	nfit;	/* Number of parameters to fit */
  static int	pfit0 = 0;	/* List of parameters to fit, 1 per digit */
  static int	cdfit = 0;	/* 1 if CD matrix has been fit */
  static int	resid_refine = 0;
  static int	minbin=2;	/* Minimum number of coincidence hits needed */
  static int	minmatch0 = MINMATCH;	/* matches to drop out of loop */
  static int	nitmax0 = NMAX;		/* max iterations to stop fit */
  static int	binarray;	/* If =1, bin matched stars */
  static int	vfit[NPAR1]; /* Parameters being fit: index to value vector
  				1= RA,		  2= Dec,
  				3= X plate scale, 4= Y plate scale
  				5= rotation,	  6= second rotation (skew),
  				7= optical axis X,8= optical axis Y */
  
  /* Find shift, scale, and rotation of image stars to best-match reference stars
   * Get best match by finding which offsets between pairs of s's and g's
   * work for the most other pairs of s's and g's
   * N.B. we assume rotation will be "small enough" so that initial guesses can
   *   be done using just shifts.
   * Return count of total coincidences found, else 0 if none or -1 if trouble.
   */
  
  int
  StarMatch (ns,sx,sy,refcat,ng,gnum,gra,gdec,goff,gx,gy,tol,wcs,debug)
  
  int	ns;		/* Number of image stars */
  double	*sx;		/* Image star X coordinates in pixels */
  double	*sy;		/* Image star Y coordinates in pixels */
  int	refcat;		/* Reference Catalog code */
  int	ng;		/* Number of reference stars */
  double	*gnum;		/* Reference star catalog numbers */
  double	*gra;		/* Reference star right ascensions in degrees */
  double	*gdec;		/* Reference star right ascensions in degrees */
  int	*goff;		/* Reference star offscale flags */
  double	*gx;		/* Reference star X coordinates in pixels */
  double	*gy;		/* Reference star Y coordinates in pixels */
  double	tol;		/* +/- this many pixels is a hit */
  struct WorldCoor *wcs;	/* World coordinate structure (fit returned) */
  int	debug;
  
  {
      double dx, bestdx, dxi;
      double dy, bestdy, dyi;
      double dx2, dy2, dxy, dxys, dxs, dys, dxsum, dysum;
      double *mx, *my, *mxy;
      int nmatch;
      int s, g, si, gi, igs;
      int nbin;
      double *sbx, *sby;	/* malloced array of s stars in best bin */
      double *gbra, *gbdec;	/* malloced array of g stars in best bin */
      int peaks[NPEAKS+1];	/* history of bin counts */
      int dxpeaks[NPEAKS+1], dypeaks[NPEAKS+1]; /* history of dx/dy at peaks */
      int npeaks;		/* entries in use in peaks[] */
      int maxnbin, i, nmatchd;
      int minmatch;
      int *is, *ig, *ibs, *ibg;
      char rastr[32], decstr[32], numstr[32];
      double xref0, yref0, xinc0, yinc0, rot0, xrefpix0, yrefpix0, cd0[4];
      int bestbin;	/* Number of coincidences for refit */
      int pfit;		/* List of parameters to fit, 1 per digit */
      char vpar[16];	/* List of parameters to fit */
      char *vi;
      char vc;
      int ParamFit();
      double tol2 = tol * tol;
      double maxnum;
      int nnfld = 0;
  
      /* Set minimum number of matches between image and reference stars to fit */
      if (ns > ng) {
  	minmatch = 0.5 * ng;
  	if (minmatch > minmatch0)
  	    minmatch = 0.25 * ng;
  	if (minmatch > minmatch0)
  	    minmatch = minmatch0;
  	}
      else {
  	minmatch = 0.5 * ns;
  	if (minmatch > minmatch0)
  	    minmatch = 0.25 * ns;
  	if (minmatch > minmatch0)
  	    minmatch = minmatch0;
  	}
  
      /* Set format for numbers, if listed */
      if (debug) {
  	maxnum = gnum[0];
  	for (gi = 1; gi < ng; gi++) {
  	    if (gnum[gi] > maxnum)
  		maxnum = gnum[gi];
  	    }
  	nnfld = CatNumLen (refcat, maxnum, 0);
          }
  
      /* Set maximum number of matches and allocate match indices */
      if (ng > ns)
  	maxnbin = (int) ((double) ng * 1.25);
      else
  	maxnbin = (int) ((double) ns * 1.25);
      if (debug)
  	fprintf (stderr,"Match history: nim=%d nref=%d tol=%3.0f minbin=%d minmatch=%d):\n",
  		 ns, ng, tol, minbin, minmatch);
  
      /* Allocate arrays in which to save match information */
      is = (int *) calloc (maxnbin, sizeof(int));
      ig = (int *) calloc (maxnbin, sizeof(int));
      ibs = (int *) calloc (maxnbin, sizeof(int));
      ibg = (int *) calloc (maxnbin, sizeof(int));
  
      /* Try matching stars using the current WCS first */
      nmatch = 0;
      bestdx = 0.0;
      bestdy = 0.0;
      dxsum = 0.0;
      dysum = 0.0;
      dxs = 0.0;
      dys = 0.0;
  
      /* Build 2-d histogram of offset distribution */
      if (binarray) {
  	int nbx, nnbx, nbx2, nby, nnby, nby2;
  	int idxmax, idymax, idx, idy, tol2;
  	int *obin, *obini, nset, npos;
  
  	/* Loop through image stars */
  	tol2 = 2 * (int) tol;
  	nbx = wcs->nxpix / tol2;
  	nnbx = -nbx;
  	nbx2 = 2 * nbx;
  	nby = wcs->nypix / tol2;
  	nby2 = 2 * nby;
  	nnby = -nby;
  	obin = calloc (nbx2 * nby2, sizeof(int));
  	nset = 0;
  	npos = 0;
  	for (s = 0; s < ns; s++) {
  
  	    /* Loop through reference catalog stars */
  	    for (g = 0; g < ng; g++) {
  		dx  = gx[g] - sx[s];
  		dy = gy[g] - sy[s];
  
  		/* Add to number in this offset bin */
  		idx = (int) (dx / tol2);
  		idy = (int) (dy / tol2);
  		if (idx > nnbx && idx < nbx && idy > nnby && idy < nby) {
  		    obin[((idy+nby) * nbx2) + idx + nbx]++;
  		    nset++;
  		    }
  		npos++;
  		}
  	    }
  
  	/* Find offset bin with maximim number of entries */
  	idxmax = 0;
  	idymax = 0;
  	nmatch = 0;
  	obini = obin;
  	for (idy = nnby; idy < nby; idy++) {
  	    for (idx = nnbx; idx < nbx; idx++) {
  		if (*obini > nmatch) {
  		    idxmax = idx;
  		    idymax = idy;
  		    nmatch = *obini;
  		    if (debug)
  			fprintf (stderr, "%5d at offset %5d, %5d\n", *obini,
  				 idx*tol2, idy*tol2);
  		    }
  		obini++;
  		}
  	    }
  	/* If we found enough matches, we can proceed with this offset */
  	if (nmatch >= minmatch) {
  	    bestdx = tol2 * (double) (idxmax - nbx);
  	    bestdy = tol2 * (double) (idymax - nby);
  	    if (debug)
  		fprintf (stderr, "%d matches found at mean offset %6.3f %6.3f\n",
  			nmatch, bestdx, bestdy);
  	    }
  	free (obin);
  	}
  
      /* Vote for closest match */
      else {
  	mx = (double *) calloc (maxnbin, sizeof(double));
  	my = (double *) calloc (maxnbin, sizeof(double));
  	mxy = (double *) calloc (maxnbin, sizeof(double));
  
  	/* Loop through image stars */
  	for (s = 0; s < ns; s++) {
  	    dxys = tol2;
  	    igs = -1;
  
  	    /* Loop through reference catalog stars */
  	    for (g = 0; g < ng; g++) {
  
  		/* Try reference catalog star only if it is on the image */
  		dx = gx[g] - sx[s];
  		dy = gy[g] - sy[s];
  		dx2 = dx * dx;
  		dy2 = dy * dy;
  		dxy = dx2 + dy2;
  
  		/* Check offset less than tolerance or this star's closest match */
  		if (dxy < dxys) {
  		    dxys = dxy;
  		    dxs = dx;
  		    dys = dy;
  		    igs = g;
  		    ibs[nmatch] = s;
  		    ibg[nmatch] = g;
  		    }
  		}
  
  	    /* If a match was found */
  	    if (igs > -1) {
  
  		/* if new match is closer than old match, replace it */
  		if (mxy[igs] > 0.0) {
  		    if (dxy < mxy[igs]) {
  			dxsum = dxsum - mx[igs];
  			dysum = dysum - my[igs];
  			dxsum = dxsum + dxs;
  			dysum = dysum + dys;
  			if (debug) {
  			    CatNum (refcat, nnfld, 0, gnum[ibg[nmatch]], numstr);
  			    ra2str (rastr, 31, gra[ibg[nmatch]], 3);
  			    dec2str (decstr, 31, gdec[ibg[nmatch]], 2);
  			    fprintf (stderr, "*%3d %s %s %s %7.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f\n",
  				 nmatch, numstr, rastr, decstr, 
  				 gx[ibg[nmatch]], gy[ibg[nmatch]], 
  				 sx[ibs[nmatch]], sy[ibs[nmatch]], 
  				 dxs, dys, sqrt (dxys));
  			    }
  			}
  		    }
  	
  		/* If not matched before, use new match */
  		else {
  		    dxsum = dxsum + dxs;
  		    dysum = dysum + dys;
  		    if (debug) {
  			CatNum (refcat, nnfld, 0, gnum[ibg[nmatch]], numstr);
  			ra2str (rastr, 31, gra[ibg[nmatch]], 3);
  			dec2str (decstr, 31, gdec[ibg[nmatch]], 2);
  			fprintf (stderr, " %3d %s %s %s %7.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f\n",
  			     nmatch, numstr, rastr, decstr, 
  			     gx[ibg[nmatch]], gy[ibg[nmatch]], 
  			     sx[ibs[nmatch]], sy[ibs[nmatch]], 
  			     dxs, dys, sqrt (dxys));
  			}
  		    nmatch++;
  		    mx[igs] = dxs;
  		    my[igs] = dys;
  		    mxy[igs] = dxy;
  		    }
  		}
  	    }
  	free (mxy);
  	free (mx);
  	free (my);
  
  	/* If we found enough matches, we can proceed with this offset */
  	if (nmatch >= minmatch) {
  	    bestdx = dxsum / (double) nmatch;
  	    bestdy = dysum / (double) nmatch;
  	    if (debug)
  		fprintf (stderr, "%d matches found at mean offset %6.3f %6.3f\n",
  			nmatch, bestdx, bestdy);
  	    }
  	}
  
      /* Otherwise, we will look for a coarse alignment assuming no additional rotation.
       * This will allow us to collect a set of stars that correspond and
       * establish an initial guess of the solution.
       */
      if (nmatch < minmatch) {
  	if (debug)
  	    fprintf (stderr, "%d matches found  less than %d minimum\n",
  			nmatch, minmatch);
  	npeaks = 0;
  	nmatch = 0;
  	for (i = 0; i < NPEAKS; i++) {
  	    peaks[i] = 0;
  	    dxpeaks[i] = 0;
  	    dypeaks[i] = 0;
  	    }
  	bestdx = 0.0;
  	bestdy = 0.0;
  	for (s = 0; s < ns; s++) {
  	    for (g = 0; g < ng; g++) {
  		dx = gx[g] - sx[s];
  		dy = gy[g] - sy[s];
  		nbin = 0;
  		for (gi = 0; gi < ng; gi++) {
  		    for (si = 0; si < ns; si++) {
  			dxi = gx[gi] - sx[si] - dx;
  			if (dxi < 0)
  			    dxi = -dxi;
  			dyi = gy[gi] - sy[si] - dy;
  			if (dyi < 0)
  			    dyi = -dyi;
  			if (dxi <= tol && dyi <= tol) {
  			/* if (debug)
  			    fprintf (stderr,"%d %d %d %d %5.1f %5.1f %5.1f %5.1f\n",
  				     g,s,gi,si,dx,dy,dxi,dyi); */
  			    is[nbin] = si;
  			    ig[nbin] = gi;
  			    nbin++;
  			    }
  			}
  		    }
  		/* if (debug)
  		    fprintf (stderr,"%d %d %d %d %d\n", g,s,gi,si,nbin); */
  		if (nbin > 1 && nbin >= nmatch) {
  		    int i;
  		    nmatch = nbin;
  		    bestdx = (double) dx;
  		    bestdy = (double) dy;
  		    for (i = 0; i < nbin; i++) {
  			ibs[i] = is[i];
  			ibg[i] = ig[i];
  			}
  	
  		    /* keep last NPEAKS nmatchs, dx and dy;
  		     * put newest first in arrays */
  		    if (npeaks > 0) {
  			for (i = npeaks; i > 0; i--) {
  			    peaks[i] = peaks[i-1];
  			    dxpeaks[i] = dxpeaks[i-1];
  			    dypeaks[i] = dypeaks[i-1];
  			    }
  			}
  		    peaks[0] = nmatch;
  		    if (bestdx > 0.0)
  			dxpeaks[0] = (int) (bestdx + 0.5);
  		    else
  			dxpeaks[0] = (int) (bestdx - 0.5);
  		    if (bestdy > 0)
  			dypeaks[0] = (int) (bestdy + 0.5);
  		    else
  			dypeaks[0] = (int) (bestdy - 0.5);
  		    if (npeaks < NPEAKS)
  			npeaks++;
  		    if (debug)
  			fprintf (stderr,"%d: %d/%d matches at image %d cat %d: dx= %d dy= %d\n",
  				npeaks, nmatch, minmatch, s, g, dxpeaks[0], dypeaks[0]);
  		    }
  		if (nmatch > minmatch)
  		    break;
  		}
  	    if (nmatch > minmatch)
  		break;
  	    }
  
  	/* if (debug) {
  	    int i;
  	    for (i = 0; i < npeaks; i++)
  		fprintf (stderr," %d bins at dx=%d dy=%d\n",
  			 peaks[i], dxpeaks[i], dypeaks[i]);
  	    } */
  
  	/* peak is broad */
  	if (npeaks < 2 || peaks[1] == peaks[0]) {
  	    if (debug)
  		fprintf (stderr,"  Broad peak of %d bins at dx=%.0f dy=%.0f\n",
  		         peaks[0], bestdx, bestdy);
  	    }
  	}
  
      /* too few hits */
      if (nmatch < minbin)
  	return (nmatch);
  
      /* Get X and Y coordinates of matches from best binning */
      nmatchd = nmatch * sizeof (double);
      if (!(sbx = (double *) malloc (nmatchd)))
  	fprintf (stderr," Could not allocate %d bytes for SBX\n", nmatchd);
      if (!(sby = (double *) malloc (nmatchd)))
  	fprintf (stderr," Could not allocate %d bytes for SBY\n", nmatchd);
      if (!(gbra = (double *) malloc (nmatchd)))
  	fprintf (stderr," Could not allocate %d bytes for GBRA\n", nmatchd);
      if (!(gbdec = (double *) malloc (nmatchd)))
  	fprintf (stderr," Could not allocate %d bytes for GBDEC\n", nmatchd);
      for (i = 0; i < nmatch; i++) {
  	sbx[i] = sx[ibs[i]];
  	sby[i] = sy[ibs[i]];
  	gbra[i] = gra[ibg[i]];
  	gbdec[i] = gdec[ibg[i]];
  	}
  
      /* Reset image center based on star matching */
      wcs->xref = wcs->xref + (bestdx * wcs->xinc);
      if (wcs->xref < 0.0) wcs->xref = 360.0 + wcs->xref;
      wcs->yref = wcs->yref + (bestdy * wcs->yinc);
  
      /* Fit WCS to matched stars */
  
      /* Provide non-parametric access to the star lists */
      sx_p = sbx;
      sy_p = sby;
      gra_p = gbra;
      gdec_p = gbdec;
      xref_p = wcs->xref;
      yref_p = wcs->yref;
      xrefpix = wcs->xrefpix;
      yrefpix = wcs->yrefpix;
      nbin_p = nmatch;
  
      /* Number of parameters to fit from command line or number of matches */
      pfit = ParamFit (nmatch);
  
      /* Get parameters to fit from digits of pfit */
      sprintf (vpar, "%d", pfit);
      nfit = 0;
      vfit[0] = -1;
      for (i = 1; i < NPAR1; i++) {
  	vc = i + 48;
  	vi = strchr (vpar, vc);
  	if (vi != NULL) {
  	    vfit[i] = vi - vpar;
  	    nfit++;
  	    }
  	else
  	    vfit[i] = -1;
  	}
  
      /* Set initial guesses for parameters which are being fit */
      xref0 = wcs->xref;
      yref0 = wcs->yref;
      xinc0 = wcs->xinc;
      yinc0 = wcs->yinc;
      rot0 = wcs->rot;
      xrefpix0 = wcs->xrefpix;
      yrefpix0 = wcs->yrefpix;
      cd0[0] = wcs->cd[0];
      cd0[1] = wcs->cd[1];
      cd0[2] = wcs->cd[2];
      cd0[3] = wcs->cd[3];
      if (vfit[6] > -1)
  	cdfit = 1;
      else
  	cdfit = 0;
  
      /* Fit image star coordinates to reference star positions */
      wcs_amoeba (wcs);
  
      if (debug) {
  	fprintf (stderr,"\nAmoeba fit:\n");
  	ra2str (rastr, 31, xref0, 3);
  	dec2str (decstr, 31, yref0, 2);
  	fprintf (stderr,"   initial guess:\n");
  	if (vfit[6] > -1)
  	    fprintf (stderr," cra= %s cdec= %s cd = %9.7f,%9.7f,%9.7f,%9.7f ",
  		     rastr, decstr, cd0[0], cd0[1], cd0[2], cd0[3]);
  	else
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		     rastr, decstr, xinc0*3600.0, yinc0*3600.0, rot0);
  	fprintf (stderr,"(%8.2f,%8.2f\n", xrefpix0, yrefpix0);
  
  	ra2str (rastr, 31, wcs->xref, 3);
  	dec2str (decstr, 31, wcs->yref, 2);
  	fprintf (stderr,"\nfirst solution:\n");
  	if (vfit[6] > -1)
  	    fprintf (stderr," cra= %s cdec= %s cd = %9.7f,%9.7f,%9.7f,%9.7f ",
  		     rastr,decstr,wcs->cd[0],wcs->cd[1],wcs->cd[2],wcs->cd[3]);
  	else
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		     rastr,decstr,3600.0*wcs->xinc,3600.0*wcs->yinc,wcs->rot);
  	fprintf (stderr,"(%8.2f,%8.2f)\n", wcs->xrefpix, wcs->yrefpix);
  	}
  
      /* If we have extra bins, repeat with the best ones */
      bestbin = nfit + 1;
      if (resid_refine && nmatch > bestbin) {
  	double *resid = (double *) malloc (nmatch * sizeof(double));
  	double *xe = (double *) malloc (nmatch * sizeof(double));
  	double *ye = (double *) malloc (nmatch * sizeof(double));
  	int i, j;
  	double xmean, ymean, rmean, xsumsq, ysumsq, diff;
  	double mx, my, xsig, ysig, rsig, siglim;
  	char wcstring[64];
  	double xsum = 0.0;
  	double ysum = 0.0;
  	double rsum = 0.0;
  	double dmatch = (double)nmatch;
  	double dmatch1 = (double)(nmatch - 1);
  
  	/* Compute residuals at each star location */
  	for (i = 0; i < nmatch; i++) {
  	    pix2wcs (wcs, sbx[i], sby[i], &mx, &my);
  	    xe[i] = (mx - gbra[i]) * 3600.0;
  	    ye[i] = (my - gbdec[i]) * 3600.0;
  	    resid[i] = sqrt (xe[i]*xe[i] + ye[i]*ye[i]);
  	    if (debug) {
  		pix2wcst (wcs, sbx[i], sby[i], wcstring, 64);
  		fprintf (stderr,"%3d (%8.3f,%8.3f) -> %s %6.3f %6.3f %6.3f\n",
  		    i, sbx[i], sby[i], wcstring, xe[i], ye[i], resid[i]);
  		}
  	    xsum = xsum + xe[i];
  	    ysum = ysum + ye[i];
  	    rsum = rsum + resid[i];
  	    }
  
  	/* Compute means and standard deviations */
  	xmean = xsum / dmatch;
  	ymean = ysum / dmatch;
  	rmean = rsum / dmatch;
  	xsumsq = 0.0;
  	ysumsq = 0.0;
  	for (i = 0; i < nmatch; i++) {
  	    diff = xe[i] - xmean;
  	    xsumsq = xsumsq + (diff * diff);
  	    diff = ye[i] - ymean;
  	    ysumsq = ysumsq + (diff * diff);
  	    }
  	xsig = sqrt (xsumsq / dmatch1);
  	ysig = sqrt (ysumsq / dmatch1);
  	rsig = sqrt ((xsumsq + ysumsq)/ dmatch1);
  	siglim = 2.0 * rsig;
  	if (debug) {
  	    fprintf (stderr,"Mean x: %6.3f/%6.3f y: %6.3f/%6.3f r: %6.3f/%6.3f\n",
  		    xmean, xsig, ymean, ysig, rmean, rsig);
  	    }
  
  	/* sort by increasing total residual */
  	for (i = 0; i < nmatch-1; i++) {
  	    for (j = i+1; j < nmatch; j++) {
  		if (resid[j] < resid[i]) {
  		    double tmp;
  
  		    tmp = sbx[i]; sbx[i] = sbx[j]; sbx[j] = tmp;
  		    tmp = sby[i]; sby[i] = sby[j]; sby[j] = tmp;
  		    tmp = gbra[i]; gbra[i] = gbra[j]; gbra[j] = tmp;
  		    tmp = gbdec[i]; gbdec[i] = gbdec[j]; gbdec[j] = tmp;
  		    tmp = resid[i]; resid[i] = resid[j]; resid[j] = tmp;
  		    }
  		}
  	    }
  
  	/* Cut off points at residual of two sigma */
  	for (i = 0; i < nmatch; i++) {
  	    if (resid[i] > siglim) {
  		if (i > bestbin) bestbin = i - 1;
  		break;
  		}
  	    }
  
  	xref_p = wcs->xref;
  	if (xref_p < 0.0) xref_p = 360.0 + xref_p;
  	yref_p = wcs->yref;
  	xrefpix = wcs->xrefpix;
  	yrefpix = wcs->yrefpix;
  	nbin_p = bestbin;
  	wcs_amoeba (wcs);
  
  	if (debug) {
  	    ra2str (rastr, 31, wcs->xref, 3);
  	    dec2str (decstr, 31, wcs->yref, 2);
  	    fprintf (stderr,"\nresid solution:\n");
  	    fprintf (stderr,"\n%d points < %.3f arcsec residuals refit\n",
  			    bestbin, siglim);
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		 rastr, decstr, 3600.0*wcs->xinc, 3600.0*wcs->yinc, wcs->rot);
  	    fprintf (stderr,"(%8.2f,%8.2f)\n", wcs->xrefpix, wcs->yrefpix);
  	    }
  	free (resid);
  	free (xe);
  	free (ye);
  	}
  
      free (sbx);
      free (sby);
      free (gbra);
      free (gbdec);
      free (is);
      free (ig);
      free (ibs);
      free (ibg);
  
      return (nmatch);
  }
  
  int
  ParamFit (nbin)
  
  int	nbin;	/* Number of point to be fit */
  {
      int pfit;
  
      if (pfit0 != 0) {
  	if (pfit0 < 3)
  	    pfit = 12;
  	else if (pfit0 == 3)	/* Fit center and plate scale */
  	    pfit = 123;
  	else if (pfit0 == 4)	/* Fit center, plate scale, rotation */
  	    pfit = 1235;
  	else if (pfit0 == 5)	/* Fit center, x&y plate scales, rotation */
  	    pfit = 12345;
  	else if (pfit0 == 6)	/* Fit center, x&y plate scales, x&y rotations */
  	    pfit = 123456;
  	else if (pfit0 == 7)	/* Fit center, x&y plate scales, rotation, refpix */
  	    pfit = 1234578;
  	else if (pfit0 == 8)	/* Fit center, x&y plate scales, x&y rotation, refpix */
  	    pfit = 12345678;
  	else
  	    pfit = pfit0;
  	}
      else if (nbin < 4)
  	pfit = 12;
      else if (nbin < 6)
  	pfit = 123;
      else
  	pfit = 12345;
      return (pfit);
  }
  
  
  int
  NParamFit (nbin)
  
  int	nbin;	/* Number of point to be fit */
  {
      int pfit;
  
      pfit = ParamFit (nbin);
      if (pfit < 1)
  	return (0);
      else if (pfit < 10)
  	return (1);
      else if (pfit < 100)
  	return (2);
      else if (pfit < 1000)
  	return (3);
      else if (pfit < 10000)
  	return (4);
      else if (pfit < 100000)
  	return (5);
      else if (pfit < 1000000)
  	return (6);
      else if (pfit < 10000000)
  	return (7);
      else
  	return (8);
  }
  
  
  int
  ReadMatch (filename, sx, sy, sra, sdec, debug)
  
  char	*filename;	/* Name of file containing matches */
  double	**sx;		/* Image star X coordinates in pixels */
  double	**sy;		/* Image star Y coordinates in pixels */
  double	**sra;		/* Probable image star right ascensions in degrees */
  double	**sdec;		/* Probable image star declinations in degrees */
  int	debug;		/* Printed debugging information if not zero */
  
  {
      int nbytes, nread, ir, ntok, itok, iytok;
      double *tx, *ty, *tra, *tdec, ra, dec, x, y;
      int ndec;
      int nmatch = 0;	/* Number of matches read from file */
      char rastr[32], decstr[32];
  
      /* If tab file, read from ra, dec, x, y  columns */
      if (istab (filename)) {
  	}
  
      /* Otherwise, assume first 4 columns are x, y, ra, dec */
      else {
  	char line[1025];
  	char *nextline, *lastchar;
  	FILE *fd;
  	struct Tokens tokens;	/* Token structure */
  	char *cwhite;		/* additional whitespace characters */
  	char token[256];
  
  	cwhite = NULL;
  
  	/* Open input file */
  	if (!strcmp (filename, "stdin")) {
  	    fd = stdin;
  	    nread = 1000;
  	    }
  	else {
  	    nread = getfilelines (filename);
  	    if (!(fd = fopen (filename, "r"))) {
  	    fprintf (stderr, "SetWCSFITS: Match file %s could not be opened\n",
  		     filename);
  		return (0);
  		}
  	    }
  	nbytes = nread * sizeof (double);
  	if (!(tra = (double *) calloc (nread, sizeof(double))))
  	    fprintf (stderr, "Could not calloc %d bytes for gra\n", nbytes);
  	if (!(tdec = (double *) calloc (nread, sizeof(double))))
  	    fprintf (stderr, "Could not calloc %d bytes for gdec\n", nbytes);
  	if (!(tx = (double *) calloc (nread, sizeof(double))))
  	    fprintf (stderr, "Could not calloc %d bytes for sx\n", nbytes);
  	if (!(ty = (double *) calloc (nread, sizeof(double))))
  	    fprintf (stderr, "Could not calloc %d bytes for sy\n", nbytes);
  	*sra = tra;
  	*sdec = tdec;
  	*sx = tx;
  	*sy = ty;
  
  	nmatch = 0;
  	nextline = line;
          for (ir = 0; ir < nread; ir++) {
  	    if (fgets (line, 1024, fd) == NULL)
  		break;
  
  	    /* Skip lines with comments */
  	    if (line[0] == '#')
  		continue;
  
  	    /* Drop linefeeds */
  	    lastchar = nextline + strlen(nextline) - 1;
  	    if (*lastchar < 32)
  		*lastchar = (char) 0;
  
  	    /* Read X, Y, RA, and Dec from each line,
  		skipping line if all four are not present and numbers */
  	    ntok = setoken (&tokens, line, cwhite);
  	    if (ntok < 1)
  		break;
  	    if (ntok < 4)
  		continue;
  
  	    /* if (debug)
  		fprintf (stderr, "%d: %s\n", nmatch, line); */
  
  	    /* Image X coordinate or RA */
  	    itok = 1;
  	    if (getoken(&tokens, itok, token, 256)) {
  
  		/* Read RA, Dec, X, Y if first token has : in it */
  		if (strchr (token, ':') != NULL) {
  		    ra = str2ra (token);
  		    iytok = 4;
  		    if (getoken(&tokens, 2, token, 256))
  			dec = str2dec (token);
  		    if (getoken(&tokens, 3, token, 256)) {
  			if (isnum (token))
  			    x = atof (token);
  			else {
  			    iytok = 5;
  			    if (getoken(&tokens, 4, token, 256)) {
  				if (isnum (token))
  				    x = atof (token);
  				else
  				    continue;
  				}
  			    }
  			}
  		    if (getoken(&tokens, iytok, token, 256)) {
  			if (isnum (token))
  			    y = atof (token);
  			else
  			    continue;
  			}
  		    tx[nmatch] = x;
  		    ty[nmatch] = y;
  		    tra[nmatch] = ra;
  		    tdec[nmatch] = dec;
  		    nmatch++;
  		    continue;
  		    }
  		if (isnum (token))
  		    x = atof (token);
  		else
  		    continue;
  		}
  	    else
  		continue;
  
  	    /* Image Y coordinate */
  	    itok++;
  	    if (getoken(&tokens, itok, token, 256)) {
  		if (isnum (token))
  		    y = atof (token);
  		else
  		    continue;
  		}
  	    else
  		continue;
  
  	    /* Right ascension */
  	    itok++;
  	    if (getoken(&tokens, itok, token, 256)) {
  		if (isnum (token) == 1) {
  		    ra = atof (token);
  		    itok++;
  		    if (getoken(&tokens, itok, token, 256)) {
  			if (isnum (token) == 2)
  			    ra = ra + (atof (token) / 60.0);
  			else if (isnum (token) == 1) {
  			    ra = ra + (atof (token) / 60.0);
  			    itok++;
  			    if (getoken(&tokens, itok, token, 256)) {
  				if (isnum (token))
  				    ra = ra + (atof (token) / 3600.0);
  				}
  			    }
  			}
  		    ra = ra * 15.0;
  		    }
  		else
  		    ra = str2ra (token);
  		}
  	    else
  		continue;
  
  	    /* Declination */
  	    itok++;
  	    if (getoken(&tokens, itok, token, 256)) {
  		if (isnum (token) == 1) {
  		    dec = atof (token);
  		    itok++;
  		    if (strchr (token, '-') != NULL)
  			ndec = 1;
  		    else
  			ndec = 0;
  		    if (getoken(&tokens, itok, token, 256)) {
  			if (isnum (token) == 2) {
  			    if (ndec)
  				dec = dec - (atof (token) / 60.0);
  			    else
  				dec = dec + (atof (token) / 60.0);
  			    }
  			else if (isnum (token) == 1) {
  			    if (ndec)
  				dec = dec - (atof (token) / 60.0);
  			    else
  				dec = dec + (atof (token) / 60.0);
  			    itok++;
  			    if (getoken(&tokens, itok, token, 256)) {
  				if (isnum (token)) {
  				    if (ndec)
  					dec = dec - (atof (token) / 3600.0);
  				    else
  					dec = dec + (atof (token) / 3600.0);
  				    }
  				}
  			    }
  			}
  		    }
  		else
  		    dec = str2dec (token);
  		}
  	    else
  		continue;
  	    tx[nmatch] = x;
  	    ty[nmatch] = y;
  	    tra[nmatch] = ra;
  	    tdec[nmatch] = dec;
  	    if (debug) {
  		ra2str (rastr, 31, tra[nmatch], 3);
  		dec2str (decstr, 31, tdec[nmatch], 2);
  		fprintf (stderr, "%d: %8.3f %8.3f %s %s\n", nmatch,
  			 tx[nmatch], ty[nmatch], rastr, decstr);
  		}
  	    nmatch++;
  	    }
  	}
  
      return (nmatch);
  }
  
  
  /* Find shift, scale, and rotation of image stars to best-match reference stars
   */
  
  void
  WCSMatch (nmatch, sbx, sby, gbra, gbdec, debug)
  
  int	nmatch;		/* Number of matched stars */
  double	*sbx;		/* Image star X coordinates in pixels */
  double	*sby;		/* Image star Y coordinates in pixels */
  double	*gbra;		/* Reference star right ascensions in degrees */
  double	*gbdec;		/* Reference star right ascensions in degrees */
  int	debug;		/* Printed debugging information if not zero */
  
  {
      int i;
      double xdiff, ydiff;
      int nsc, j;
      double tx = 0.0;
      double ty = 0.0;
      double tra = 0.0;
      double tdec = 0.0;
      double tdiff = 0.0;
      double cra, cdec, cx, cy, scale;
      double dmatch;
      double skydiff, imdiff;
      extern double getsecpix();
      extern void getcenter(),getrefpix(),setdcenter(),setrefpix(),setsecpix();
  
      dmatch = (double) nmatch;
  
      if (debug) {
  	fprintf (stderr,"%d matched stars:\n", nmatch);
  	}
  
      /* too few hits */
      if (nmatch < 2)
  	return;
  
      /* Compute plate scale and center of stars */
      nsc = 0;
      for (i = 0; i < nmatch; i++) {
  	tx = tx + sbx[i];
  	ty = ty + sby[i];
  	tra = tra + gbra[i];
  	tdec = tdec + gbdec[i];
  	for (j = i+1; j < nmatch; j++) {
  	    skydiff = wcsdist (gbra[i], gbdec[i], gbra[j], gbdec[j]);
  	    xdiff = sbx[j] - sbx[i];
  	    ydiff = sby[j] - sby[i];
  	    imdiff = sqrt ((xdiff * xdiff) + (ydiff * ydiff));
  	    scale = skydiff / imdiff;
  	    tdiff = tdiff + scale;
  	    nsc++;
  	    if (debug) {
  		fprintf (stderr,"%d %d: sky: %8g, image: %8g, %8g deg/pix", 
  			i, j, skydiff, imdiff, scale);
  		fprintf (stderr," = %8g arcsec/pix\n", scale * 3600.0);
  		}
  	    }
  	}
      
      /* Reset image center based on star matching */
      getcenter (&cra, &cdec);
      if (cra == -99.0 && cdec == -99.0) {
  	cra = tra / dmatch;
  	cdec = tdec / dmatch;
  	setdcenter (cra, cdec);
  	}
      getrefpix (&cx, &cy);
      if (cx == -99999.0) {
  	cx = tx / dmatch;
  	cy = ty / dmatch;
  	setrefpix (cx, cy);
  	}
      scale = getsecpix();
      if (scale == 0.0) {
  	scale = tdiff / (double) nsc;
  	setsecpix (3600.0 * scale);
  	}
      if (debug)
  	fprintf (stderr,"scale = %8g deg/pix = %8g arcsec/pix\n", scale, scale*3600.0);
      return;
  }
  
  
  /* Find shift, scale, and rotation of image stars to best-match reference stars
   * Return count of total coincidences found, else 0 if none or -1 if trouble.
   */
  
  int
  FitMatch (nmatch, sbx, sby, gbra, gbdec, wcs, debug)
  
  int	nmatch;		/* Number of matched stars */
  double	*sbx;		/* Image star X coordinates in pixels */
  double	*sby;		/* Image star Y coordinates in pixels */
  double	*gbra;		/* Reference star right ascensions in degrees */
  double	*gbdec;		/* Reference star right ascensions in degrees */
  struct WorldCoor *wcs;	/* World coordinate structure (fit returned) */
  int	debug;		/* Printed debugging information if not zero */
  
  {
      int i;
      char rastr[32], decstr[32];
      double xref0, yref0, xinc0, yinc0, rot0, xrefpix0, yrefpix0, cd0[4];
      int bestbin;	/* Number of coincidences for refit */
      int pfit;		/* List of parameters to fit, 1 per digit */
      char vpar[16];	/* List of parameters to fit */
      double xdiff, ydiff;
      char *vi;
      char vc;
      int nsc, j;
  /*    double equinox = wcs->equinox; */
      double tx = 0.0;
      double ty = 0.0;
      double tra = 0.0;
      double tdec = 0.0;
      double tdiff = 0.0;
      double scale;
  /*    double dmatch; */
      double skydiff, imdiff;
  
  /*    dmatch = (double) nmatch; */
  
      if (debug) {
  	fprintf (stderr,"%d matched stars:\n", nmatch);
  	}
  
      /* too few hits */
      if (nmatch < minbin)
  	return (nmatch);
  
      /* Compute plate scale and center of stars */
      nsc = 0;
      for (i = 0; i < nmatch; i++) {
  	tx = tx + sbx[i];
  	ty = ty + sby[i];
  	tra = tra + gbra[i];
  	tdec = tdec + gbdec[i];
  	for (j = i+1; j < nmatch; j++) {
  	    skydiff = wcsdist (gbra[i], gbdec[i], gbra[j], gbdec[j]);
  	    xdiff = sbx[j] - sbx[i];
  	    ydiff = sby[j] - sby[i];
  	    imdiff = sqrt ((xdiff * xdiff) + (ydiff * ydiff));
  	    scale = skydiff / imdiff;
  	    tdiff = tdiff + scale;
  	    nsc++;
  	    if (debug) {
  		fprintf (stderr,"%d %d: sky: %8g, image: %8g, %8g deg/pix", 
  			i, j, skydiff, imdiff, scale);
  		fprintf (stderr," = %8g arcsec/pix\n", scale * 3600.0);
  		}
  	    }
  	}
      
      /* Reset image center in WCS data structure based on star matching */
      /* cra = tra / dmatch;
      cdec = tdec / dmatch;
      cx = tx / dmatch;
      cy = ty / dmatch;
      scale = tdiff / (double) nsc;
      if (debug)
  	fprintf (stderr,"scale = %8g deg/pix = %8g arcsec/pix\n", scale, scale*3600.0);
      wcsreset (wcs, cx, cy, cra, cdec, scale, 0.0, 0.0, NULL, equinox); */
  
      /* Provide non-parametric access to the star lists */
      sx_p = sbx;
      sy_p = sby;
      gra_p = gbra;
      gdec_p = gbdec;
      xref_p = wcs->xref;
      if (xref_p < 0.0) xref_p = 360.0 + xref_p;
      yref_p = wcs->yref;
      xrefpix = wcs->xrefpix;
      yrefpix = wcs->yrefpix;
      nbin_p = nmatch;
  
      /* Number of parameters to fit from command line or number of matches */
      if (pfit0 != 0) {
  	if (pfit0 < 3)
  	    pfit = 12;
  	else if (pfit0 == 3)	/* Fit center and plate scale */
  	    pfit = 123;
  	else if (pfit0 == 4)	/* Fit center, plate scale, rotation */
  	    pfit = 1235;
  	else if (pfit0 == 5)	/* Fit center, x&y plate scales, rotation */
  	    pfit = 12345;
  	else if (pfit0 == 6)	/* Fit center, x&y plate scales, x&y rotations */
  	    pfit = 123456;
  	else if (pfit0 == 7)	/* Fit center, x&y plate scales, rotation, refpix */
  	    pfit = 1234578;
  	else if (pfit0 == 8)	/* Fit center, x&y plate scales, x&y rotation, refpix */
  	    pfit = 12345678;
  	else
  	    pfit = pfit0;
  	}
      else if (nmatch < 4)
  	pfit = 12;
      else if (nmatch < 6)
  	pfit = 123;
      else
  	pfit = 12345;
  
      /* Get parameters to fit from digits of pfit */
      sprintf (vpar, "%d", pfit);
      nfit = 0;
      vfit[0] = -1;
      for (i = 1; i < NPAR1; i++) {
  	vc = i + 48;
  	vi = strchr (vpar, vc);
  	if (vi != NULL) {
  	    vfit[i] = vi - vpar;
  	    nfit++;
  	    }
  	else
  	    vfit[i] = -1;
  	}
  
      /* Set initial guesses for parameters which are being fit */
      xref0 = wcs->xref;
      yref0 = wcs->yref;
      xinc0 = wcs->xinc;
      yinc0 = wcs->yinc;
      rot0 = wcs->rot;
      xrefpix0 = wcs->xrefpix;
      yrefpix0 = wcs->yrefpix;
      cd0[0] = wcs->cd[0];
      cd0[1] = wcs->cd[1];
      cd0[2] = wcs->cd[2];
      cd0[3] = wcs->cd[3];
      if (vfit[6] > -1)
  	cdfit = 1;
      else
  	cdfit = 0;
  
      /* Fit image star coordinates to reference star positions */
      wcs_amoeba (wcs);
  
      if (debug) {
  	fprintf (stderr,"\nAmoeba fit:\n");
  	ra2str (rastr, 31, xref0, 3);
  	dec2str (decstr, 31, yref0, 2);
  	fprintf (stderr,"   initial guess:\n");
  	if (vfit[6] > -1)
  	    fprintf (stderr," cra= %s cdec= %s cd = %9.7f,%9.7f,%9.7f,%9.7f ",
  		     rastr, decstr, cd0[0], cd0[1], cd0[2], cd0[3]);
  	else
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		     rastr, decstr, xinc0*3600.0, yinc0*3600.0, rot0);
  	fprintf (stderr,"(%8.2f,%8.2f\n", xrefpix0, yrefpix0);
  
  	ra2str (rastr, 31, wcs->xref, 3);
  	dec2str (decstr, 31, wcs->yref, 2);
  	fprintf (stderr,"\nfirst solution:\n");
  	if (vfit[6] > -1)
  	    fprintf (stderr," cra= %s cdec= %s cd = %9.7f,%9.7f,%9.7f,%9.7f ",
  		     rastr,decstr,wcs->cd[0],wcs->cd[1],wcs->cd[2],wcs->cd[3]);
  	else
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		     rastr,decstr,3600.0*wcs->xinc,3600.0*wcs->yinc,wcs->rot);
  	fprintf (stderr,"(%8.2f,%8.2f)\n", wcs->xrefpix, wcs->yrefpix);
  	}
  
      /* If we have extra bins, repeat with the best ones */
      bestbin = nfit + 1;
      if (resid_refine && nmatch > bestbin) {
  	double *resid = (double *) malloc (nmatch * sizeof(double));
  	double *xe = (double *) malloc (nmatch * sizeof(double));
  	double *ye = (double *) malloc (nmatch * sizeof(double));
  	int i, j;
  	double xmean, ymean, rmean, xsumsq, ysumsq, diff;
  	double mra, mdec, xsig, ysig, rsig, siglim;
  	char wcstring[64];
  	double xsum = 0.0;
  	double ysum = 0.0;
  	double rsum = 0.0;
  	double dmatch = (double)nmatch;
  	double dmatch1 = (double)(nmatch - 1);
  
  	/* Compute residuals at each star location */
  	for (i = 0; i < nmatch; i++) {
  	    pix2wcs (wcs, sbx[i], sby[i], &mra, &mdec);
  	    xe[i] = (mra - gbra[i]) * 3600.0;
  	    ye[i] = (mdec - gbdec[i]) * 3600.0;
  	    resid[i] = sqrt (xe[i]*xe[i] + ye[i]*ye[i]);
  	    if (debug) {
  		pix2wcst (wcs, sbx[i], sby[i], wcstring, 64);
  		fprintf (stderr,"%3d (%8.3f,%8.3f) -> %s %6.3f %6.3f %6.3f\n",
  		    i, sbx[i], sby[i], wcstring, xe[i], ye[i], resid[i]);
  		}
  	    xsum = xsum + xe[i];
  	    ysum = ysum + ye[i];
  	    rsum = rsum + resid[i];
  	    }
  
  	/* Compute means and standard deviations */
  	xmean = xsum / dmatch;
  	ymean = ysum / dmatch;
  	rmean = rsum / dmatch;
  	xsumsq = 0.0;
  	ysumsq = 0.0;
  	for (i = 0; i < nmatch; i++) {
  	    diff = xe[i] - xmean;
  	    xsumsq = xsumsq + (diff * diff);
  	    diff = ye[i] - ymean;
  	    ysumsq = ysumsq + (diff * diff);
  	    }
  	xsig = sqrt (xsumsq / dmatch1);
  	ysig = sqrt (ysumsq / dmatch1);
  	rsig = sqrt ((xsumsq + ysumsq)/ dmatch1);
  	siglim = 2.0 * rsig;
  	if (debug) {
  	    fprintf (stderr,"Mean x: %6.3f/%6.3f y: %6.3f/%6.3f r: %6.3f/%6.3f\n",
  		    xmean, xsig, ymean, ysig, rmean, rsig);
  	    }
  
  	/* sort by increasing total residual */
  	for (i = 0; i < nmatch-1; i++) {
  	    for (j = i+1; j < nmatch; j++) {
  		if (resid[j] < resid[i]) {
  		    double tmp;
  
  		    tmp = sbx[i]; sbx[i] = sbx[j]; sbx[j] = tmp;
  		    tmp = sby[i]; sby[i] = sby[j]; sby[j] = tmp;
  		    tmp = gbra[i]; gbra[i] = gbra[j]; gbra[j] = tmp;
  		    tmp = gbdec[i]; gbdec[i] = gbdec[j]; gbdec[j] = tmp;
  		    tmp = resid[i]; resid[i] = resid[j]; resid[j] = tmp;
  		    }
  		}
  	    }
  
  	/* Cut off points at residual of two sigma */
  	for (i = 0; i < nmatch; i++) {
  	    if (resid[i] > siglim) {
  		if (i > bestbin) bestbin = i - 1;
  		break;
  		}
  	    }
  
  	xref_p = wcs->xref;
  	if (xref_p < 0.0) xref_p = 360.0 + xref_p;
  	yref_p = wcs->yref;
  	xrefpix = wcs->xrefpix;
  	yrefpix = wcs->yrefpix;
  	nbin_p = bestbin;
  	wcs_amoeba (wcs);
  
  	if (debug) {
  	    ra2str (rastr, 31, wcs->xref, 3);
  	    dec2str (decstr, 31, wcs->yref, 2);
  	    fprintf (stderr,"\nresid solution:\n");
  	    fprintf (stderr,"\n%d points < %.3f arcsec residuals refit\n",
  			    bestbin, siglim);
  	    fprintf (stderr," cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ",
  		 rastr, decstr, 3600.0*wcs->xinc, 3600.0*wcs->yinc, wcs->rot);
  	    fprintf (stderr,"(%8.2f,%8.2f)\n", wcs->xrefpix, wcs->yrefpix);
  	    }
  	free (resid);
  	free (xe);
  	free (ye);
  	}
  
      return (nmatch);
  }
  
  struct WorldCoor *wcsf;
  
  static double wcs_chisqr ();
  
  /* From Numerical Recipes */
  void amoeba();
  static double amotry();
  
  
  /* Set up the necessary temp arrays and call the amoeba() multivariate solver */
  
  static void
  wcs_amoeba (wcs0)
  
  struct WorldCoor *wcs0;
  
  {
      double *p[NPAR1];				  /* used as p[NPAR1][NPAR] */
      double vguess[NPAR], vp[NPAR], vdiff[NPAR];
      double p0[NPAR], p1[NPAR], p2[NPAR], p3[NPAR], p4[NPAR],
  	   p5[NPAR], p6[NPAR], p7[NPAR], p8[NPAR]; /* used as px[0..NPAR-1] */
      double y[NPAR1];				  /* used as y[1..NPAR] */
      double xinc1, yinc1, xrefpix1, yrefpix1, rot, cd[4];
      double sumx, sumy, sumr;
      int iter;
      int i, j;
      int nfit1;
      char rastr[32],decstr[32];
      int nitmax;
  
      nitmax = nitmax0;
      if (nfit > NPAR)
  	nfit = NPAR;
      nfit1 = nfit + 1;
      wcsf = wcs0;
  
  /* Initialize guess and difference vectors to zero */
      for (i = 0; i < NPAR; i++) {
  	vguess[i] = 0.0;
  	vdiff[i] = 0.0;
  	}
  
      /* Optical axis center (RA and Dec degrees) */
      if (vfit[1] > -1) {
  	vguess[vfit[1]] = 0.0;
  	vdiff[vfit[1]] = 5.0 * wcsf->xinc;
  	}
      if (vfit[2] > -1) {
  	vguess[vfit[2]] = 0.0;
  	vdiff[vfit[2]] = 5.0 * wcsf->yinc;
  	}
      /* Second rotation about optical axis (degrees) -> CD matrix */
      if (vfit[6] > -1) {
  	wcsf->rotmat = 1;
  	vguess[vfit[3]] = wcsf->cd[0];
  	vdiff[vfit[3]] = wcsf->xinc * 0.03;
  	vguess[vfit[4]] = wcsf->cd[1];
  	vdiff[vfit[4]] = wcsf->yinc * 0.03;
  	vguess[vfit[5]] = wcsf->cd[2];
  	vdiff[vfit[5]] = wcsf->xinc * 0.03;
  	vguess[vfit[6]] = wcsf->cd[3];
  	vdiff[vfit[6]] = wcsf->yinc * 0.03;
  	}
  
      else {
  
      /* Plate scale at optical axis right ascension or both (degrees/pixel) */
  	if (vfit[3] > -1) {
  	    vguess[vfit[3]] = wcsf->xinc;
  	    vdiff[vfit[3]] = wcsf->xinc * 0.03;
  	    }
  
      /* Plate scale in declination at optical axis (degrees/pixel) */
  	if (vfit[4] > -1) {
  	    vguess[vfit[4]] = wcsf->yinc;
  	    vdiff[vfit[4]] = wcsf->yinc * 0.03;
  	    }
  
      /* Rotation about optical axis in degrees */
  	if (vfit[5] > -1) {
  	    vguess[vfit[5]] = wcsf->rot;
  	    vdiff[vfit[5]] = 0.5;
  	    }
  	}
  
  /* Reference pixel (optical axis) */
      if (vfit[7] > -1) {
  	vguess[vfit[7]] = 0.0;
  	vdiff[vfit[7]] = 10.0;
  	}
      if (vfit[8] > -1) {
  	vguess[vfit[8]] = 0.0;
  	vdiff[vfit[8]] = 10.0;
  	}
  
  /* Set up matrix of nfit+1 initial guesses.
   * The supplied guess, plus one for each parameter altered by a small amount
   */
      p[0] = p0;
      if (nfit > 0) p[1] = p1;
      if (nfit > 1) p[2] = p2;
      if (nfit > 2) p[3] = p3;
      if (nfit > 3) p[4] = p4;
      if (nfit > 4) p[5] = p5;
      if (nfit > 5) p[6] = p6;
      if (nfit > 6) p[7] = p7;
      if (nfit > 7) p[8] = p8;
      for (i = 0; i <= nfit; i++) {
  	for (j = 0; j < nfit; j++)
  	    p[i][j] = vguess[j];
  	if (i > 0 && i <= nfit)
  	    p[i][i-1] = vguess[i-1] + vdiff[i-1];
  	 y[i] = wcs_chisqr (p[i], -i);
  	}
  
  #define	PDUMP
  #ifdef	PDUMP
      fprintf (stderr,"Before:\n");
      for (i = 0; i < nfit1; i++) {
  	if (vfit[1] > -1)
  	    ra2str (rastr, 31, p[i][vfit[1]] + xref_p, 3);
  	else
  	    ra2str (rastr, 31, wcsf->xref, 3);
  	if (vfit[2] > -1)
  	    dec2str (decstr, 16, p[i][vfit[2]]+yref_p, 2);
  	else
  	    dec2str (decstr, 16, wcsf->yref, 2);
  	if (vfit[6] > -1) {
  	    cd[0] = p[i][vfit[3]];
  	    cd[1] = p[i][vfit[4]];
  	    cd[2] = p[i][vfit[5]];
  	    cd[3] = p[i][vfit[6]];
  	    fprintf (stderr,"%d: %s %s CD: %7.5f,%7.5f,%7.5f,%7.5f ",
  		    i, rastr, decstr, cd[0],cd[1],cd[2],cd[3]);
  	    }
  	else {
  	    if (vfit[3] > -1)
  		xinc1 = p[i][vfit[3]];
  	    else
  		xinc1 = wcsf->xinc;
  	    if (vfit[4] > -1)
  		yinc1 = p[i][vfit[4]];
  	    else if (vfit[3] > -1) {
  		if (xinc1 < 0)
  		    yinc1 = -xinc1;
  		else
  		    yinc1 = xinc1;
  		}
  	    else
  		yinc1 = wcsf->yinc;
  	    if (vfit[5] > -1)
  		rot = p[i][vfit[5]];
  	    else
  		rot = wcsf->rot;
  	    fprintf (stderr,"%d: %s %s del=%6.4f,%6.4f rot=%5.3f ",
  		    i, rastr, decstr, 3600.0*xinc1, 3600.0*yinc1, rot);
  	    }
  
  	if (vfit[7] > -1)
  	    xrefpix1 = xrefpix + p[i][vfit[7]];
  	else
  	    xrefpix1 = wcsf->xrefpix;
  	if (vfit[8] > -1)
  	    yrefpix1 = yrefpix + p[i][vfit[8]];
  	else
  	    yrefpix1 = wcsf->yrefpix;
  	fprintf (stderr,"(%8.2f,%8.2f) y=%g\n", xrefpix1, yrefpix1, y[i]);
  	}
  #endif
  
      amoeba (p, y, nfit, FTOL, nitmax, wcs_chisqr, &iter);
  
  #define	PDUMP
  #ifdef	PDUMP
      fprintf (stderr,"\nAfter:\n");
      for (i = 0; i < nfit1; i++) {
  	if (vfit[1] > -1)
  	    ra2str (rastr, 31, p[i][vfit[1]] + xref_p, 3);
  	else
  	    ra2str (rastr, 31, wcsf->xref, 3);
  	if (vfit[2] > -1)
  	    dec2str (decstr, 31, p[i][vfit[2]]+yref_p, 2);
  	else
  	    dec2str (decstr, 31, wcsf->yref, 2);
  	if (vfit[6] > -1) {
  	    cd[0] = p[i][vfit[3]];
  	    cd[1] = p[i][vfit[4]];
  	    cd[2] = p[i][vfit[5]];
  	    cd[3] = p[i][vfit[6]];
  	    fprintf (stderr,"%d: %s %s CD: %7.5f,%7.5f,%7.5f,%7.5f ",
  		    i, rastr, decstr, cd[0],cd[1],cd[2],cd[3]);
  	    }
  	else {
  	    if (vfit[3] > -1)
  		xinc1 = p[i][vfit[3]];
  	    else
  		xinc1 = wcsf->xinc;
  	    if (vfit[4] > -1)
  		yinc1 = p[i][vfit[4]];
  	    else if (vfit[3] > -1) {
  		if (xinc1 < 0)
  		    yinc1 = -xinc1;
  		else
  		    yinc1 = xinc1;
  		}
  	    else
  		yinc1 = wcsf->yinc;
  	    if (vfit[5] > -1)
  		rot = p[i][vfit[5]];
  	    else
  		rot = wcsf->rot;
  	    fprintf (stderr,"%d: %s %s del=%6.4f,%6.4f rot=%5.3f ",
  		    i,rastr,decstr, 3600.0*xinc1, 3600.0*yinc1, rot);
  	    }
  	if (vfit[7] > -1)
  	    xrefpix1 = xrefpix + p[i][vfit[7]];
  	else
  	    xrefpix1 = wcsf->xrefpix;
  	if (vfit[8] > -1)
  	    yrefpix1 = yrefpix + p[i][vfit[8]];
  	else
  	    yrefpix1 = wcsf->yrefpix;
  	fprintf (stderr,"(%8.2f,%8.2f) y=%g\n", xrefpix1, yrefpix1, y[i]);
  	}
  #endif
  
      /* On return, all entries in p[1..NPAR] are within FTOL;
       * Return the average, though you could just pick the first one
       */
      for (j = 0; j < nfit; j++) {
  	double sum = 0.0;
  	for (i = 0; i < nfit1; i++)
  	    sum += p[i][j];
  	vp[j] = sum / (double)nfit1;
  	}
      if (vfit[1] > -1) {
  	wcsf->xref = xref_p + vp[vfit[1]];
  	if (wcsf->xref < 0.0) wcsf->xref = 360.0 + wcsf->xref;
  	}
      if (vfit[2] > -1)
  	wcsf->yref = yref_p + vp[vfit[2]];
      if (vfit[6] > -1) {
  	wcsf->cd[0] = vp[vfit[3]];
  	wcsf->cd[1] = vp[vfit[4]];
  	wcsf->cd[2] = vp[vfit[5]];
  	wcsf->cd[3] = vp[vfit[6]];
  	}
      else {
  	if (vfit[3] > -1)
  	    wcsf->xinc = vp[vfit[3]];
  	if (vfit[4] > -1)
  	    wcsf->yinc = vp[vfit[4]];
  	else if (vfit[3] > -1) {
  	    if (wcsf->xinc < 0)
  		wcsf->yinc = -wcsf->xinc;
  	    else
  		wcsf->yinc = wcsf->xinc;
  	    }
  	if (vfit[5] > -1)
  	    wcsf->rot = vp[vfit[5]];
  	}
      if (vfit[7] > -1)
  	wcsf->xrefpix = xrefpix + vp[vfit[7]];
      if (vfit[8] > -1)
  	wcsf->yrefpix = yrefpix + vp[vfit[8]];
  
  #define RESIDDUMP
  #ifdef RESIDDUMP
      ra2str (rastr, 31, wcsf->xref, 3);
      dec2str (decstr, 31, wcsf->yref, 2);
  
      if (vfit[6] > -1)
  	fprintf (stderr,"iter=%d\n cra= %s cdec= %s CD=%9.7f,%9.7f,%9.7f,%9.7f ", iter,
  		rastr, decstr, wcsf->cd[0], wcsf->cd[1], wcsf->cd[2],
  		wcsf->cd[3]);
      else
  	fprintf (stderr,"iter=%d\n cra= %s cdec= %s del=%7.4f,%7.4f rot=%7.4f ", iter,
  		rastr, decstr, wcsf->xinc*3600.0, wcsf->yinc*3600.0, wcsf->rot);
      fprintf (stderr,"(%8.2f,%8.2f)\n", wcsf->xrefpix, wcsf->yrefpix);
      sumx = 0.0;
      sumy = 0.0;
      sumr = 0.0;
      for (i = 0; i < nbin_p; i++) {
  	double mra, mdec, ex, ey, er;
  	char rastr[32], decstr[32];
  
  	pix2wcs (wcsf, sx_p[i], sy_p[i], &mra, &mdec);
  	ex = 3600.0 * (mra - gra_p[i]);
  	ey = 3600.0 * (mdec - gdec_p[i]);
  	er = sqrt (ex * ex + ey * ey);
  	sumx = sumx + ex;
  	sumy = sumy + ey;
  	sumr = sumr + er;
  
  	ra2str (rastr, 31, gra_p[i], 3);
  	dec2str (decstr, 31, gdec_p[i], 2);
  	fprintf (stderr,"%2d: c: %s %s ", i+1, rastr, decstr);
  	ra2str (rastr, 31, mra, 3);
  	dec2str (decstr, 31, mdec, 2);
  	fprintf (stderr, "i: %s %s %6.3f %6.3f %6.3f\n",
  		rastr, decstr, 3600.0*ex, 3600.0*ey,
  		3600.0*sqrt(ex*ex + ey*ey));
  	}
      sumx = sumx / (double)nbin_p;
      sumy = sumy / (double)nbin_p;
      sumr = sumr / (double)nbin_p;
      fprintf (stderr,"mean dra: %6.3f, ddec: %6.3f, dr = %6.3f\n", sumx, sumy, sumr);
  #endif
  }
  
  
  /* Compute the chisqr of the vector v, where
   * v[0]=cra, v[1]=cdec, v[2]=ra deg/pix, v[3]=dec deg/pix,
   * v[4]=rotation, v[5]=2nd rotation->CD matrix, v[6]=ref x, and v[7] = ref y
   * chisqr is in arcsec^2
   */
  
  static double
  wcs_chisqr (v, iter)
  
  double	*v;	/* Vector of parameter values */
  int	iter;	/* Number of iterations */
  
  {
      double chsq;
      char rastr[32],decstr[32];
      double xmp, ymp, dx, dy, cd[4], *cdx;
      double crval1, crval2, cdelt1, cdelt2, crota, crpix1, crpix2;
      int i, offscale;
  
      /* Set WCS parameters from fit parameter vector */
  
      /* Sky coordinates at optical axis (degrees) */
      if (vfit[1] > -1)
  	crval1 = xref_p + v[vfit[1]];
      else
  	crval1 = wcsf->xref;
      if (vfit[2] > -1)
  	crval2 = yref_p + v[vfit[2]];
      else
  	crval2 = wcsf->yref;
  
      /* CD matrix */
      if (vfit[6] > -1) {
  	cdelt1 = 0.0;
  	cdelt2 = 0.0;
  	crota = 0.0;
  	cd[0] = v[vfit[3]];
  	cd[1] = v[vfit[4]];
  	cd[2] = v[vfit[5]];
  	cd[3] = v[vfit[6]];
  	cdx = cd;
  	}
  
      else {
  	/* Plate scale (degrees/pixel) */
  	if (vfit[3] > -1)
  	    cdelt1 = v[vfit[3]];
  	else
  	    cdelt1 = wcsf->xinc;
  	if (vfit[4] > -1)
  	    cdelt2 = v[vfit[4]];
  	else if (vfit[3] > -1) {
  	    if (cdelt1 < 0)
  		cdelt2 = -cdelt1;
  	    else
  		cdelt2 = cdelt1;
  	    }
  	else
  	    cdelt2 = wcsf->yinc;
  
  	/* Rotation angle (degrees) */
  	if (vfit[5] > -1)
  	    crota = v[vfit[5]];
  	else
  	    crota = wcsf->rot;
  	cdx = NULL;
  	}
  
      /* Optical axis pixel coordinates */
      if (vfit[7] > -1)
  	crpix1 = xrefpix + v[vfit[7]];
      else
  	crpix1 = wcsf->xrefpix;
      if (vfit[8] > -1)
  	crpix2 = yrefpix + v[vfit[8]];
      else
  	crpix2 = wcsf->yrefpix;
      if (wcsreset (wcsf,crpix1,crpix2,crval1,crval2,cdelt1,cdelt2,crota,cdx)) {
  	fprintf (stderr,"CHISQR: Cannot reset WCS!\n");
  	return (0.0);
  	}
  
      /* Compute sum of squared residuals for these parameters */
      chsq = 0.0;
      for (i = 0; i < nbin_p; i++) {
  	wcs2pix (wcsf, gra_p[i], gdec_p[i], &xmp, &ymp, &offscale);
  	/* if (!offscale) { */
  	    dx = xmp - sx_p[i];
  	    dy = ymp - sy_p[i];
  	    chsq += dx*dx + dy*dy;
  	    /* } */
  	}
  
  #define TRACE_CHSQR
  #ifdef TRACE_CHSQR
      ra2str (rastr, 31, wcsf->xref, 3);
      dec2str (decstr, 31, wcsf->yref, 2);
      if (vfit[6] > -1)
  	fprintf (stderr,"%4d: %s %s CD: %9.7f,%9.7f,%9.7f,%9.7f ",
  		iter, rastr, decstr, wcsf->cd[0],wcsf->cd[1],wcsf->cd[2],
  		wcsf->cd[3]);
      else
  	fprintf (stderr,"%4d: %s %s %9.7f,%9.7f %8.5f ",
  		iter, rastr, decstr, wcsf->xinc*3600.0, wcsf->yinc*3600.0,
  		wcsf->rot);
      fprintf (stderr,"(%8.2f,%8.2f) -> %f\r",
  	     wcsf->xrefpix, wcsf->yrefpix, chsq);
  #endif
      return (chsq);
  }
  
  /* The following subroutines are based on those in Numerical Recipes in C */
  
  /* amoeba.c */
  
  #define ALPHA 1.0
  #define BETA 0.5
  #define GAMMA 2.0
  
  void
  amoeba (p, y, ndim, ftol, itmax, funk, nfunk)
  
  double	**p;
  double	y[];
  double	ftol;
  int	itmax;
  double	(*funk)();
  int	ndim;
  int	*nfunk;
  
  {
  int i,j,ilo,ihi,inhi,ndim1=ndim+1;
  double ytry,ysave,sum,rtol,*psum;
  
      psum = (double *) malloc ((unsigned)ndim * sizeof(double));
      *nfunk = 0;
      for (j=0; j<ndim; j++) {
  	for (i=0,sum=0.0; i<ndim1; i++)
  	    sum += p[i][j]; psum[j]=sum;
  	}
      for (;;) {
  	ilo=1;
  	if (y[0] > y[1]) {
  	    inhi = 1;
  	    ihi = 0;
  	    }
  	else {
  	    inhi = 0;
  	    ihi = 1;
  	    }
  	for (i = 0; i < ndim1; i++) {
  	    if (y[i] < y[ilo])
  		ilo=i;
  	    if (y[i] > y[ihi]) {
  		inhi=ihi;
  		ihi=i;
  		}
  	    else if (y[i] > y[inhi])
  		if (i != ihi)
  		    inhi=i;
  	    }
  	rtol = 2.0 * fabs(y[ihi]-y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]));
  	if (rtol < ftol)
  	    break;
  	if (*nfunk >= itmax) {
  	    fprintf (stderr,"Too many iterations in amoeba fit %d > %d",*nfunk,itmax);
  	    return;
  	    }
  	ytry = amotry (p, y, psum, ndim, funk, ihi, nfunk, -ALPHA);
  	if (ytry <= y[ilo])
  	    ytry = amotry (p, y, psum, ndim, funk, ihi, nfunk, GAMMA);
  	else if (ytry >= y[inhi]) {
  	    ysave = y[ihi];
  	    ytry = amotry (p,y,psum,ndim,funk,ihi,nfunk,BETA);
  	    if (ytry >= ysave) {
  		for (i = 0; i < ndim1; i++) {
  		    if (i != ilo) {
  			for (j = 0; j < ndim; j++) {
  			    psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
  			    p[i][j] = psum[j];
  			    }
  			y[i]=(*funk)(psum, *nfunk);
  			}
  		    }
  		*nfunk += ndim;
  		for (j=0; j<ndim; j++) {
  		    for (i=0,sum=0.0; i<ndim1; i++)
  			sum += p[i][j]; psum[j]=sum;
  		    }
  		}
  	    }
  	}
      free (psum);
      return;
  }
  
  
  static double
  amotry (p, y, psum, ndim, funk, ihi, nfunk, fac)
  
  double	**p;
  double	*y;
  double	*psum;
  double	(*funk)();
  double	fac;
  int	ndim;
  int	ihi;
  int	*nfunk;
  
  {
      int j;
      double fac1,fac2,ytry,*ptry;
  
      ptry = (double *) malloc ((unsigned) ndim * sizeof(double));
      fac1 = (1.0 - fac) / ndim;
      fac2 = fac1 - fac;
      for (j = 0; j < ndim; j++)
  	ptry[j] = psum[j] * fac1 - p[ihi][j] * fac2;
      ytry = (*funk)(ptry, *nfunk);
      ++(*nfunk);
      if (ytry < y[ihi]) {
  	y[ihi] = ytry;
  	for (j = 0; j < ndim; j++) {
      	    psum[j] +=  ptry[j] - p[ihi][j];
      	    p[ihi][j] = ptry[j];
  	    }
  	}
      free (ptry);
      return ytry;
  }
  
  void
  setbin (binflag)
  int binflag;
  { binarray = binflag; return;}
  
  void
  setresid_refine (refine)
  int refine;
  { resid_refine = refine; return; }
  
  int
  getresid_refine ()
  { return (resid_refine); }
  
  void
  setnfit (nfit)
  int nfit;
  {
      if (nfit == 0)
  	setnofit();
      else if (nfit < 0) {
  	pfit0 = -nfit;
  	resid_refine = 1;
  	}
      else {
  	pfit0 = nfit;
  	resid_refine = 0;
  	}
      return;
  }
  
  int
  getnfit ()
  { return (pfit0); }
  
  int
  iscdfit ()
  { return (cdfit); }
  
  void
  setminmatch (minmatch)
  int minmatch;
  { minmatch0 = minmatch; return; }
  
  void
  setminbin (minbin1)
  int minbin1;
  { minbin = minbin1; return; }
  
  void
  setnitmax (nitmax)
  int nitmax;
  { nitmax0 = nitmax; return; }
  
  /* Aug  6 1996	New subroutine
   * Sep  1 1996	Move constants to lwcs.h
   * Sep  3 1996	Use offscale pixels for chi^2 computation
   * Sep  3 1996	Overprint chi^2 in verbose mode
   * Oct 15 1996	Fix am* subroutine declarations
   * Nov 19 1996	Fix bug regarding rotation
   *
   * Jul 21 1997	Add reference pixel position fitting
   * Aug  4 1997	Increase maximum iterations from 750 to 1000 in lwcs.h
   * Aug 28 1997	Fix VGUESS dimension bug
   * Sep  9 1997	Print RA and Dec offsets in residual listing
   * Sep  9 1997	Turn on resid_refinement if number of parameters to fit negated
   * Sep  9 1997	Fit separate horizontal and vertical plate scales if nfit=5
   * Sep  9 1997	Fix bugs associated with fitting optical axis
   * Sep 12 1997	Add chip rotation instead of second plate scale
   * Oct  2 1997	Keep second plate scale AND chip rotation
   * Oct 16 1997	Try to deal with reference pixel position correctly
   * Nov  5 1997	Select parameters one at a time, in any order
   * Nov 12 1997	Add PFIT=3 to fit center and plate scale only
   * Dec 15 1997	Fix minor bugs after lint
   *
   * Jan 26 1998	Remove chip rotation code
   * Jan 29 1998	Streamline initialization code
   * Feb 19 1998	Fix bug in initialization code
   * Mar  3 1998	Fix residual-refining code
   * Mar 20 1998	Add option to fit CD matrix
   * Mar 25 1998	Make amoeba() externally callable
   * Mar 26 1998	Return instead of crashing when too many iterations
   * Apr 21 1998	Drop out of loop if more than half of stars are matched
   * Apr 27 1998	Fix bug handling nfit=8
   * Jun 24 1998	Fix bug summing unitialized values for mean after fit
   * Jun 24 1998	Add string lengths to ra2str() and dec2str() calls
   * Oct  8 1998	Initialize bestdx and bestdy to zero
   * Dec  8 1998	Fix declaration of amotry()
   *
   * Apr 21 1999	Add subroutines to set and retrieve resid_refine independently
   * Jul 21 1999	Add FitMatch() to fit WCS to already-matched stars
   * Sep  8 1999	Fix bug found by Jean-Baptiste Marquette
   * Oct  1 1999	Add ReadMatch() to read a set of matches from a file
   * Oct 20 1999	Include wcscat.h
   *
   * Feb 15 2000	Add iscdfit() to return whether CD matrix is being fit
   * Mar 10 2000	Add debug statement to list max matches as they are found
   * Mar 10 2000	Change loop order to image stars first
   * Dec 18 2000	Write half of ReadMatch() to deal with ASCII files
   *
   * Jan  2 2001	Modify ReadMatch() to read hh mm ss dd mm ss, too
   * Jan  9 2001	Work on FitMatch()
   * Jan 11 2001	All diagnostic printing goes to stderr
   * Feb 28 2001	Ignore coordinate system if present after match file coordinates
   * Jun 18 2001	Add maximum length of returned string to getoken()
   * Aug  2 2001	Separate parameter listing and counting into subroutines
   * Sep 19 2001	Drop fitshead.h; it is in wcs.h
   * Sep 24 2001	Ease match numeric criterium if half num is > 40
   * Oct 15 2001	Simplify error message
   * Oct 16 2001	Read minimum match to drop out of loop from lwcs.h
   * Oct 31 2001	Simplify innermost loop to try for more speed
   * Nov  1 2001	Add goff to StarMatch() arguments
   * Nov  5 2001	Use current WCS with no offset before trying offset matching
   * Nov  6 2001	Add setnitmax() to set maximum number of amoeba iterations
   * Nov  7 2001	Add setminbin to set minimum number of matches for fit
   * Nov 16 2001	Allocate slightly more than maxbin to handle dense fields
   *
   * Jul 31 2002	Add getnfit() to return current number of parameters being fit
   * Aug 30 2002	Fix WCSMatch() to set scale in arcsec, not degrees
   *
   * Jan 30 2003	Remove uninitialized variable in WCSMatch()
   * Mar 13 2003	Do not include malloc.h on Apples and Convexes
   * Apr  3 2003	Clean up code with lint
   * Nov 18 2003	Drop include of malloc.h; it is in stdlib.h
   *
   * Aug 30 2004	Declare void various external set*() calls
   *
   * Jun 19 2006	Initialize unitialized variables dxs and dys
   * Sep 26 2006	Increase length of rastr and destr from 16 to 32
   * Oct 23 2006	Add bin array option to see if it works better
   *
   * Jan  8 2007	Drop unused variables; fix two format disagreements
   * Jan 11 2007	Include fitsfile.h
   */ 
  