text: /*** File libwcs/gscread.c
   *** September 22, 2009
   *** By Jessica Mink, jmink@cfa.harvard.edu
   *** Harvard-Smithsonian Center for Astrophysics
   *** Copyright (C) 1996-2009
   *** 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: jmink@cfa.harvard.edu
             Postal address: Jessica Mink
                             Smithsonian Astrophysical Observatory
                             60 Garden St.
                             Cambridge, MA 02138 USA
   */
  
  #include <stdlib.h>
  #include <unistd.h>
  #include <stdio.h>
  #include <math.h>
  #include <string.h>
  #include "fitsfile.h"
  #include "wcs.h"
  #include "wcscat.h"
  
  /* Pathname of northern hemisphere GSC CDROM  or search engine URL */
  static char cdn[64]="/data/astrocat/gsc1";
  
  /* Uncomment following line to use ESO GSC server for GSC
  static char cdn[64]="http://archive.eso.org/skycat/servers/gsc-server";
   */
  
  /* Pathname of southern hemisphere GSC CDROM */
  static char cds[64]="/data/astrocat/gsc2";
  
  /* Pathname of northern hemisphere GSC-ACT CDROM  or search engine URL */
  static char cdna[64]="/data/astrocat/gscact1";
  
  /* Pathname of southern hemisphere GSC-ACT CDROM */
  static char cdsa[64]="/data/astrocat/gscact2";
  
  static void gscpath();
  static int gscreg();
  
  static char *table = NULL;	/* FITS table buffer */
  static int ltab = 0;		/* Length of FITS table buffer */
  static double *gdist = NULL;	/* Array of distances to stars */
  static int ndist = 0;
  
  void
  gscfree()
  { free ((void *)table); ltab = 0;
    free ((void *)gdist); ndist = 0;
    return; }
  
  static int classd = -1; /* Desired object class
  			   (-2=all objects, -1=all, 0=stars, 3=nonstars) */
  void
  setgsclass (class)
  int class;
  { classd = class; return; }
  
  /* GSCREAD -- Read HST Guide Star Catalog stars from CDROM */
  
  int
  gscread (refcat,cra,cdec,dra,ddec,drad,dradi,distsort,sysout,eqout,epout,
  	 mag1,mag2,nstarmax,gnum,gra,gdec,gmag,gtype,nlog)
  
  int	refcat;		/* Catalog code (GSC or GSCACT) */
  double	cra;		/* Search center J2000 right ascension in degrees */
  double	cdec;		/* Search center J2000 declination in degrees */
  double	dra;		/* Search half width in right ascension in degrees */
  double	ddec;		/* Search half-width in declination in degrees */
  double	drad;		/* Limiting separation in degrees (ignore if 0) */
  double	dradi;		/* Inner separation in degrees for annulus (ignore if 0) */
  int	distsort;	/* 1 to sort stars by distance from center */
  int	sysout;		/* Search coordinate system */
  double	eqout;		/* Search coordinate equinox */
  double	epout;		/* Proper motion epoch (0.0 for no proper motion) */
  double	mag1,mag2;	/* Limiting magnitudes (none if equal) */
  int	nstarmax;	/* Maximum number of stars to be returned */
  double	*gnum;		/* Array of Guide Star numbers (returned) */
  double	*gra;		/* Array of right ascensions (returned) */
  double	*gdec;		/* Array of declinations (returned) */
  double	**gmag;		/* Array of magnitudes (returned) */
  int	*gtype;		/* Array of object types (returned) */
  int	nlog;		/* 1 for diagnostics */
  {
      double ra1,ra2;	/* Limiting right ascensions of region in degrees */
      double dec1,dec2;	/* Limiting declinations of region in degrees */
      double dist = 0.0;  /* Distance from search center in degrees */
      double faintmag=0.0; /* Faintest magnitude */
      double maxdist=0.0; /* Largest distance */
      int	faintstar=0;	/* Faintest star */
      int	farstar=0;	/* Most distant star */
      int magsort=0;
      int nreg;		/* Number of input FITS tables files */
      double xnum;		/* Guide Star number */
      int rlist[100];	/* List of input FITS tables files */
      char inpath[64];	/* Pathname for input FITS table file */
      char entry[100];	/* Buffer for FITS table row */
      int class;		/* Object class (0>star, 3>other) */
      int sysref=WCS_J2000;	/* Catalog coordinate system */
      double eqref=2000.0;	/* Catalog equinox */
      double epref=2000.0;	/* Catalog epoch */
      struct Keyword kw[8];	/* Keyword structure */
      struct Keyword *kwn;
  
      int verbose;
      int pass;
      int rnum, num0, num, itot,ireg;
      int ik,nk,itable,ntable,jstar;
      int nbline,npos,nbhead;
      int nbr,nrmax,nstar,i;
      int ift;
      int wrap;
      double ra,rasum,dec,decsum,perr,perr2,perrsum,msum;
      double mag,merr,merr2,merrsum;
      int band0 = 0;
      int band = 0;
      int class0 = 0;
      double merr0 = 0.0;
      double mag0 = 0.0;
      double perr0 = 0.0;
      double ra0 = 0.0;
      double dec0 = 0.0;
      double rra1, rra2, rdec1, rdec2;
      double rdist, ddist;
      char *str;
      char *url;
      char *title;
      char cstr[32], numstr[32], rastr[32], decstr[32], catid[16];
  
      itot = 0;
      if (nlog == 1)
  	verbose = 1;
      else
  	verbose = 0;
      verbose = nlog;
  
      /* If root pathname is a URL, search and return */
      url = cdn;
      if (refcat == GSC) {
  	url = cdn;
  	if ((str = getenv("GSC_NORTH")) == NULL)
  	    str = getenv ("GSC_PATH");
  	if (str != NULL) url = str;
  	}
      if (refcat == GSCACT) {
  	url = cdna;
  	if ((str = getenv("GSCACT_NORTH")) == NULL)
  	    str = getenv ("GSCACT_PATH");
  	if (str != NULL) url = str;
  	}
      if (!strncmp (url, "http:",5))
  	return (webread (url,"gsc",distsort,cra,cdec,dra,ddec,drad,dradi,
  			 sysout,eqout,epout,mag1,mag2,magsort,nstarmax,
  			 gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog));
  
      /* Allocate FITS table buffer which is saved between calls */
      if (ltab < 1) {
  	ltab = 10000;
  	table = (char *)calloc (ltab, sizeof (char));
  	if (table == NULL) {
  	    fprintf (stderr, "GSCREAD: cannot allocate FITS table buffer\n");
  	    return (0);
  	    }
  	}
  
      for (i = 0; i < 100; i++)
  	entry[i] = 0;
  
      /* Set path to Guide Star Catalog */
      if (refcat == GSCACT) {
  	if ((str = getenv("GSCACT_NORTH")) != NULL )
  	    strcpy (cdna,str);
  	if ((str = getenv("GSCACT_SOUTH")) != NULL )
  	    strcpy (cdsa,str);
  	}
      else {
  	if ((str = getenv("GSC_NORTH")) != NULL )
  	    strcpy (cdn,str);
  	if ((str = getenv("GSC_SOUTH")) != NULL )
  	    strcpy (cds,str);
  	}
  
      wcscstr (cstr, sysout, eqout, epout);
  
      SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose);
  
      /* Make mag1 always the smallest magnitude */
      if (mag2 < mag1) {
  	mag = mag2;
  	mag2 = mag1;
  	mag1 = mag;
  	}
  
      /* Find Guide Star Catalog regions in which to search */
      nrmax = 100;
      rra1 = ra1;
      rra2 = ra2;
      rdec1 = dec1;
      rdec2 = dec2;
      RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,0.0,
  	    &rra1, &rra2, &rdec1, &rdec2,  &wrap, verbose);
      nreg = gscreg (refcat,rra1,rra2,rdec1,rdec2,table,nrmax,rlist,verbose);
      if (nreg <= 0) {
  	fprintf (stderr,"GSCREAD:  no Guide Star regions found\n");
  	return (0);
  	}
  
      /* Allocate array for distances from search center */
      if (nstarmax > ndist) {
  	if (ndist > 0)
  	    free ((void *)gdist);
  	gdist = (double *) malloc (nstarmax * sizeof (double));
  	if (gdist == NULL) {
  	    fprintf (stderr,"GSCREAD:  cannot allocate separation array\n");
  	    return (0);
  	    }
  	ndist = nstarmax;
  	}
  
      /* Set keyword list */
      nk = 8;
      strcpy (kw[0].kname,"GSC_ID");
      strcpy (kw[1].kname,"RA_DEG");
      strcpy (kw[2].kname,"DEC_DEG");
      strcpy (kw[3].kname,"POS_ERR");
      strcpy (kw[4].kname,"MAG");
      strcpy (kw[5].kname,"MAG_ERR");
      strcpy (kw[6].kname,"MAG_BAND");
      strcpy (kw[7].kname,"CLASS");
      for (ik = 0; ik < nk; ik++) {
  	kw[ik].lname = (int) strlen (kw[ik].kname);
  	kw[ik].kn = 0;
  	kw[ik].kf = 0;
  	kw[ik].kl = 0;
  	}
      nstar = 0;
  
      /* Write header if printing star entries as found */
      if (nstarmax < 1) {
  	char *revmessage;
  	revmessage = getrevmsg();
  	title = CatName (refcat, NULL);
  	printf ("catalog	%s\n", title);
  	free (title);
  	ra2str (rastr, 31, cra, 3);
  	printf ("ra	%s\n", rastr);
  	dec2str (decstr, 31, cdec, 2);
  	printf ("dec	%s\n", decstr);
  	if (drad != 0.0) {
  	    printf ("radmin	%.1f\n", drad*60.0);
  	    if (dradi > 0)
  		printf ("radimin	%.1f\n", dradi*60.0);
  	    }
  	else {
  	    printf ("dramin	%.1f\n", dra*60.0* cosdeg (cdec));
  	    printf ("ddecmin	%.1f\n", ddec*60.0);
  	    }
  	printf ("radecsys	%s\n", cstr);
  	printf ("equinox	%.3f\n", eqout);
  	printf ("epoch	%.3f\n", epout);
  	printf ("program	scat %s\n", revmessage);
  	CatID (catid, refcat);
  	printf ("%s	ra          	dec         	", catid);
  	printf ("magv 	class	band	n	arcmin\n");
  	printf ("---------	------------	------------	");
  	printf ("-----	-----	----	-	------\n");
  	}
  
      /* Loop through region list */
      for (ireg = 0; ireg < nreg; ireg++) {
  	gscpath (refcat, rlist[ireg], inpath);
  
      /* Read size and keyword info from FITS table header */
  	kwn = kw;
  	ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead);
  
  	rnum = rlist[ireg];
  	num0 = 0;
  	rasum = 0.0;
  	decsum = 0.0;
  	msum = 0.0;
  	perrsum = 0.0;
  	merrsum = 0.0;
  	npos = 0;
  	num = 0;
  	fitsrtlset();
  	jstar = 0;
  	class = 0;
  
  	/* Loop through FITS table for this region */
  	for (itable = 0; itable <= ntable; itable++) {
  
  	    if (itable < ntable) {
  		nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry);
  		if (nbr < nbline) {
  		    fprintf (stderr,"GSCREAD: %d / %d bytes read, line %d / %d, region %d\n",
  			      nbr,nbline,itable,ntable,rnum);
  		    break;
  		    }
  
  		/* Extract selected fields */
  
  		/* Star number within region */
  		num0 = ftgeti4 (entry, &kw[0]);
  
  		/* Right ascension in degrees */
  		ra0 = ftgetr8 (entry, &kw[1]);
  
  		/* Declination in degrees */
  		dec0 = ftgetr8 (entry, &kw[2]);
  
  		/* Position error */
  		perr0 = ftgetr8 (entry, &kw[3]);
  
  		/* Magnitude */
  		mag0 = ftgetr8 (entry, &kw[4]);
  
  		/* Magnitude error */
  		merr0 = ftgetr8 (entry, &kw[5]);
  
  		/* Bandpass code */
  		band0 = ftgeti4 (entry, &kw[6]);
  
  		/* Object class code */
  		class0 = ftgeti4 (entry, &kw[7]);
  		}
  	    else
  		num0 = 0;
  
  	/* Compute mean position and magnitude for object */
  	    if (itable > 0 && npos > 0 &&
  		((classd < -1 && band != band0) ||
  		(classd < -1 && class != class0) || 
  		num != num0)) {
  
  		pass = 1;
  		if (perrsum == 0.0 || merrsum == 0)
  		    pass = 0;
  		else {
  		    ra = rasum / perrsum;
  		    dec = decsum / perrsum;
  		    mag = msum / merrsum;
  		    }
  
  		if (pass > 0 && classd > -1 && class != classd)
  		    pass = 0;
  
  		/* Check magnitude and position limits */
  		if (pass > 0 && mag1 != mag2 && (mag < mag1 || mag > mag2))
  		    pass = 0;
  
  		if (pass) {
  		    wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout);
  		
  		    /* Compute distance from search center */
  		    if (drad > 0 || distsort)
  			dist = wcsdist (cra,cdec,ra,dec);
  		    else
  			dist = 0.0;
  
  		    /* Check position limits */
  		    if (drad > 0) {
  			if (dist > drad)
  			    pass = 0;
  			if (dradi > 0.0 && dist < dradi)
  			    pass = 0;
  			}
  		    else {
  			ddist = wcsdist (cra,cdec,cra,dec);
  			if (ddist > ddec)
  			    pass = 0;
  			rdist = wcsdist (cra,dec,ra,dec);
  			if (rdist > dra)
  			    pass = 0;
  			}
  		    }
  
  		if (pass) {
  		    class = class + (band * 100) + (npos * 10000);
  		    xnum = (double)rnum + (0.0001 * (double) num);
  
  		    /* Write star position and magnitudes to stdout */
  		    if (nstarmax < 1) {
  			CatNum (refcat, -9, 0, xnum, numstr);
  			ra2str (rastr, 31, ra, 3);
  			dec2str (decstr, 31, dec, 2);
  			dist = wcsdist (cra,cdec,ra,dec) * 60.0;
  			printf ("%s	%s	%s", numstr,rastr,decstr);
  			printf ("	%.2f	%d	%d	%d	%.2f\n",
  				mag, class-(band*100)-(npos*10000), band, npos, dist);
  			}
  
  		    /* Save star position in table */
  		    else if (nstar < nstarmax) {
  			gnum[nstar] = xnum;
  			gra[nstar] = ra;
  			gdec[nstar] = dec;
  			gmag[0][nstar] = mag;
  			gtype[nstar] = class;
  			gdist[nstar] = dist;
  			if (dist > maxdist) {
  			    maxdist = dist;
  			    farstar = nstar;
  			    }
  			if (mag > faintmag) {
  			    faintmag = mag;
  			    faintstar = nstar;
  			    }
  			}
  
  		    /* If too many stars and distance sorting,
  			replace furthest star */
  		    else if (distsort) {
  			if (dist < maxdist) {
  			    gnum[farstar] = xnum;
  			    gra[farstar] = ra;
  			    gdec[farstar] = dec;
  			    gmag[0][farstar] = mag;
  			    gtype[farstar] = class;
  			    gdist[farstar] = dist;
  			    maxdist = 0.0;
  
  			    /* Find new farthest star */
  			    for (i = 0; i < nstarmax; i++) {
  				if (gdist[i] > maxdist) {
  				    maxdist = gdist[i];
  				    farstar = i;
  				    }
  				}
  			    }
  			}
  
  		    /* If too many stars, replace faintest star */
  		    else if (mag < faintmag) {
  			gnum[faintstar] = xnum;
  			gra[faintstar] = ra;
  			gdec[faintstar] = dec;
  			gmag[0][faintstar] = mag;
  			gtype[faintstar] = class;
  			gdist[faintstar] = dist;
  			faintmag = 0.0;
  
  			/* Find new faintest star */
  			for (i = 0; i < nstarmax; i++) {
  			    if (gmag[0][i] > faintmag) {
  				faintmag = gmag[0][i];
  				faintstar = i;
  				}
  			    }
  			}
  		    nstar++;
  		    jstar++;
  		    if (nlog == 1)
  			fprintf (stderr,"GSCREAD: %04d.%04d: %9.5f %9.5f %s %5.2f %d %d\n",
  				rnum,num,ra,dec,cstr,mag,class,npos);
  		    }
  
  	/* Reset star position for averaging */
  		rasum = 0.0;
  		decsum = 0.0;
  		msum = 0.0;
  		perrsum = 0.0;
  		merrsum = 0.0;
  		npos = 0;
  		}
  
  	/* Add information from current line to current object */
  
  	    /* Check object class */
       	    if ((classd > -1 && class0 == classd) ||
  		classd < -2 || (classd < 0 && class0 != 5)) {
  		perr = perr0;
  		perr2 = perr * perr;
  		if (perr2 <= 0.0) perr2 = 0.01;
  		rasum = rasum + (ra0 / perr2);
  		decsum = decsum + (dec0 / perr2);
  		perrsum = perrsum + (1.0 / perr2);
  		if (merr0 <= 0.0) merr0 = 0.01;
  		merr = merr0;
  		merr2 = merr * merr;
  		msum = msum + (mag0 / merr2);
  		merrsum = merrsum + (1.0 / merr2);
  		num = num0;
  		class = class0;
  		band = band0;
  		npos++;
  		}
  
  	    /* Log operation */
  	    if (nlog > 0 && itable%nlog == 0)
  		fprintf (stderr,"GSCREAD: %4d / %4d: %5d / %5d  / %5d sources, region %4d.%04d\r",
  			 ireg,nreg,jstar,itable,ntable,rlist[ireg],num0);
  
  	/* End of region */
  	    }
  
  	/* Close region input file */
  	(void) close (ift);
  	itot = itot + itable;
  	if (nlog > 0)
  	    fprintf (stderr,"GSCREAD: %4d / %4d: %5d / %5d  / %5d sources from region %4d    \n",
  		     ireg+1,nreg,jstar,itable,ntable,rlist[ireg]);
  	}
  
      /* Close output file and summarize transfer */
      if (nlog > 0) {
  	if (nreg > 1)
  	    fprintf (stderr,"GSCREAD: %d regions: %d / %d found\n",nreg,nstar,itot);
  	else
  	    fprintf (stderr,"GSCREAD: 1 region: %d / %d found\n",nstar,itable);
  	if (nstar > nstarmax )
  	    fprintf (stderr,"GSCREAD: %d stars found; only %d returned\n",
  		     nstar,nstarmax);
  	}
      return (nstar);
  }
  
  /* GSCRNUM -- Read HST Guide Star Catalog stars from CDROM */
  
  int
  gscrnum (refcat, nstars, sysout, eqout, epout, gnum,gra,gdec,gmag,gtype,nlog)
  
  int	refcat;		/* Catalog code (GSC or GSCACT) */
  int	nstars;		/* Number of stars to find */
  int	sysout;		/* Search coordinate system */
  double	eqout;		/* Search coordinate equinox */
  double	epout;		/* Proper motion epoch (0.0 for no proper motion) */
  double	*gnum;		/* Array of Guide Star numbers (returned) */
  double	*gra;		/* Array of right ascensions (returned) */
  double	*gdec;		/* Array of declinations (returned) */
  double	**gmag;		/* Array of magnitudes (returned) */
  int	*gtype;		/* Array of object types (returned) */
  int	nlog;		/* 1 for diagnostics */
  {
      char *table;		/* FITS table */
      char inpath[64];		/* Pathname for input FITS table file */
      char entry[100];		/* Buffer for FITS table row */
      int class;			/* Object class (0>star, 3>other) */
      int sysref=WCS_J2000;	/* Catalog coordinate system */
      double eqref=2000.0;	/* Catalog equinox */
      struct Keyword kw[8];	/* Keyword structure */
      struct Keyword *kwn;
  
      int rnum, num0, num, itot;
      int ik,nk,itable,ntable,jstar;
      int nbline,npos,nbhead;
      int nbr,nstar,i, snum;
      int ift;
      double dnum;
      double ra,rasum,dec,decsum,perr,perr2,perrsum,msum;
      double mag,merr,merr2,merrsum;
      int band0 = 0;
      int band = 0;
      int class0 = 0;
      double merr0 = 0.0;
      double mag0 = 0.0;
      double perr0 = 0.0;
      double ra0 = 0.0;
      double dec0 = 0.0;
      char *str;
      char *url;
  
      itot = 0;
  
      /* If root pathname is a URL, search and return */
      url = cdn;
      if (refcat == GSC) {
  	url = cdn;
  	if ((str = getenv("GSC_NORTH")) == NULL)
  	    str = getenv ("GSC_PATH");
  	if (str != NULL) url = str;
  	if (!strncmp (url, "http:",5))
  	    return (webrnum (str,"gsc",nstars,sysout,eqout,epout,1,
  			     gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog));
  	}
      if (refcat == GSCACT) {
  	url = cdna;
  	if ((str = getenv("GSCACT_NORTH")) == NULL)
  	    str = getenv ("GSCACT_PATH");
  	if (str != NULL) url = str;
  	}
      if (!strncmp (url, "http:",5))
  	return (webrnum (str,"gsc",nstars,sysout,eqout,epout,1,
  			 gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog));
  
      if (ltab < 1) {
  	ltab = 10000;
  	table = (char *)calloc (ltab, sizeof (char));
  	}
  
      for (i = 0; i < 100; i++)
  	entry[i] = 0;
  
      /* Set path to Guide Star Catalog */
      if ((str = getenv("GSC_NORTH")) != NULL )
  	strcpy (cdn,str);
      if ((str = getenv("GSC_SOUTH")) != NULL )
  	strcpy (cds,str);
  
      /* Set keyword list */
      nk = 8;
      strcpy (kw[0].kname,"GSC_ID");
      strcpy (kw[1].kname,"RA_DEG");
      strcpy (kw[2].kname,"DEC_DEG");
      strcpy (kw[3].kname,"POS_ERR");
      strcpy (kw[4].kname,"MAG");
      strcpy (kw[5].kname,"MAG_ERR");
      strcpy (kw[6].kname,"MAG_BAND");
      strcpy (kw[7].kname,"CLASS");
      for (ik = 0; ik < nk; ik++) {
  	kw[ik].lname = (int) strlen (kw[ik].kname);
  	kw[ik].kn = 0;
  	kw[ik].kf = 0;
  	kw[ik].kl = 0;
  	}
      nstar = 0;
  
      /* Loop through star list */
      for (jstar = 0; jstar < nstars; jstar++) {
  	rnum = (int) gnum[jstar];
  	gscpath (refcat, rnum, inpath);
  	dnum = (gnum[jstar] - (double)rnum) * 10000.0;
  	snum = (int) (dnum + 0.5);
  
      /* Read size and keyword info from FITS table header */
  	kwn = kw;
  	ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead);
  	if (ift < 0) {
  	    fprintf (stderr,"GSCRNUM: File %s not found\n",inpath);
  	    return (0);
  	    }
  	num0 = 0;
  	rasum = 0.0;
  	decsum = 0.0;
  	msum = 0.0;
  	perrsum = 0.0;
  	merrsum = 0.0;
  	npos = 0;
  	num = 0;
  	fitsrtlset();
  	class = 0;
  
  	/* Loop through FITS table for this region */
  	for (itable = 0; itable <= ntable; itable++) {
  
  	    if (itable < ntable) {
  		nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry);
  		if (nbr < nbline) {
  		    fprintf (stderr,"GSCRNUM: %d / %d bytes read, line %d / %d, region %d\n",
  			      nbr,nbline,itable,ntable,rnum);
  		    break;
  		    }
  
  	 /* Extract selected fields */
  
  		/* Star number within region */
  		num0 = ftgeti4 (entry, &kw[0]);
  		if (num0 < snum)
  		    continue;
  		else if (num == 0)
  		    num = num0;
  
  		/* Right ascension in degrees */
  		ra0 = ftgetr8 (entry, &kw[1]);
  
  		/* Declination in degrees */
  		dec0 = ftgetr8 (entry, &kw[2]);
  
  		/* Position error */
  		perr0 = ftgetr8 (entry, &kw[3]);
  
  		/* Magnitude */
  		mag0 = ftgetr8 (entry, &kw[4]);
  
  		/* Magnitude error */
  		merr0 = ftgetr8 (entry, &kw[5]);
  
  		/* Bandpass code */
  		band0 = ftgeti4 (entry, &kw[6]);
  
  		/* Object class code */
  		class0 = ftgeti4 (entry, &kw[7]);
  		}
  	    else
  		num0 = 0;
  
  	    /* Compute mean position and magnitude for object */
  	    if (num != num0 && itable > 0 && npos > 0) {
  		ra = rasum / perrsum;
  		dec = decsum / perrsum;
  		wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout);
  		mag = msum / merrsum;
  		class = class + (band * 100) + (npos * 10000);
  
  		/* Save star position in table */
  		gra[nstar] = ra;
  		gdec[nstar] = dec;
  		gmag[0][nstar] = mag;
  		gtype[nstar] = class;
  
  		nstar++;
  		if (nlog == 1)
  		    fprintf (stderr,"GSCRNUM: %04d.%04d: %9.5f %9.5f %5.2f %d %d\n",
  			     rnum,num,ra,dec,mag,class,npos);
  
  		/* Reset star position for averaging */
  		rasum = 0.0;
  		decsum = 0.0;
  		msum = 0.0;
  		perrsum = 0.0;
  		merrsum = 0.0;
  		npos = 0;
  		break;
  		}
  
  	    /* Add information from current line to current object */
  	    perr = perr0;
  	    perr2 = perr * perr;
  	    if (perr2 <= 0.0) perr2 = 0.01;
  	    rasum = rasum + (ra0 / perr2);
  	    decsum = decsum + (dec0 / perr2);
  	    perrsum = perrsum + (1.0 / perr2);
  	    if (merr0 <= 0.0) merr0 = 0.01;
  	    merr = merr0;
  	    merr2 = merr * merr;
  	    msum = msum + (mag0 / merr2);
  	    merrsum = merrsum + (1.0 / merr2);
  	    num = num0;
  	    class = class0;
  	    band = band0;
  	    npos = npos + 1;
  
  	    /* Log operation */
  	    if (nlog > 0 && itable%nlog == 0)
  		fprintf (stderr,"GSCRNUM: %4d / %4d: %5d / %5d sources, region %4d.%04d\r",
  			 jstar+1,nstars,itable,ntable,rnum,snum);
  
  	/* End of region */
  	    }
  
  	/* Close region input file */
  	(void) close (ift);
  	itot = itot + itable;
  	if (nlog > 0)
  	    fprintf (stderr,"GSCRNUM: %4d / %4d: %5d / %5d sources, region %4d.%04d\n",
  		     jstar+1,nstars,itable,ntable,rnum,snum);
  	}
  
      return (nstars);
  }
  
  
  /* GSCBIN -- Fill FITS WCS image with HST Guide Star Catalog objects */
  
  int
  gscbin (refcat, wcs, header, image, mag1, mag2, magscale, nlog)
  
  int	refcat;		/* Catalog code (GSC or GSCACT) */
  struct WorldCoor *wcs;	/* World coordinate system for image */
  char	*header;	/* FITS header for output image */
  char	*image;		/* Output FITS image */
  double	mag1,mag2;	/* Limiting magnitudes (none if equal) */
  double	magscale;	/* Scaling factor for magnitude to pixel flux
  			 * (number of catalog objects per bin if 0) */
  int	nlog;		/* 1 for diagnostics */
  {
      double cra;		/* Search center J2000 right ascension in degrees */
      double cdec;	/* Search center J2000 declination in degrees */
      double dra;		/* Search half width in right ascension in degrees */
      double ddec;	/* Search half-width in declination in degrees */
      int	sysout = wcs->syswcs;	/* Image coordinate system */
      double eqout = wcs->equinox; /* Image coordinate equinox */
      double epout = wcs->epoch;	/* Image epoch */
      double ra1,ra2;	/* Limiting right ascensions of image in degrees */
      double dec1,dec2;	/* Limiting declinations of image in degrees */
      int nreg;		/* Number of input FITS tables files */
      double xnum;		/* Guide Star number */
      int rlist[100];	/* List of input FITS tables files */
      char inpath[64];	/* Pathname for input FITS table file */
      char entry[100];	/* Buffer for FITS table row */
      int class, class0;	/* Object class (0>star, 3>other) */
      int sysref=WCS_J2000;	/* Catalog coordinate system */
      double eqref=2000.0;	/* Catalog equinox */
      struct Keyword kw[8];	/* Keyword structure */
      struct Keyword *kwn;
  
      int verbose;
      int pass;
      int wrap;
      int rnum, num0, num, itot,ireg;
      int ik,nk,itable,ntable,jstar;
      int nbline,npos,nbhead;
      int nbr,nrmax,nstar,i;
      int ift, band0, band;
      double ra,ra0,rasum,dec,dec0,decsum,perr,perr0,perr2,perrsum,msum;
      double mag,mag0,merr,merr0,merr2,merrsum;
      double rra1, rra2, rdec1, rdec2;
      double rdist, ddist;
      char *str;
      char cstr[32];
      int bitpix, w, h;	/* Image bits/pixel and pixel width and height */
      double logt = log(10.0);
      double xpix, ypix, flux;
      int ix, iy;
      int offscl;
  
      itot = 0;
      if (nlog == 1)
  	verbose = 1;
      else
  	verbose = 0;
      verbose = nlog;
  
      /* Set image parameters */
      bitpix = 0;
      (void)hgeti4 (header, "BITPIX", &bitpix);
      w = 0;
      (void)hgeti4 (header, "NAXIS1", &w);
      h = 0;
      (void)hgeti4 (header, "NAXIS2", &h);
  
      if (refcat == GSCACT) {
  	if ((str = getenv("GSCACT_NORTH")) == NULL)
  	    str = getenv ("GSCACT_PATH");
  	}
  
      /* Allocate FITS table buffer which is saved between calls */
      if (ltab < 1) {
  	ltab = 10000;
  	table = (char *)calloc (ltab, sizeof (char));
  	if (table == NULL) {
  	    fprintf (stderr, "GSCBIN: cannot allocate FITS table buffer\n");
  	    return (0);
  	    }
  	}
  
      for (i = 0; i < 100; i++)
  	entry[i] = 0;
  
      /* Set path to Guide Star Catalog */
      if (refcat == GSCACT) {
  	if ((str = getenv("GSCACT_NORTH")) != NULL )
  	    strcpy (cdna,str);
  	if ((str = getenv("GSCACT_SOUTH")) != NULL )
  	    strcpy (cdsa,str);
  	}
      else {
  	if ((str = getenv("GSC_NORTH")) != NULL )
  	    strcpy (cdn,str);
  	if ((str = getenv("GSC_SOUTH")) != NULL )
  	    strcpy (cds,str);
  	}
  
      wcscstr (cstr, sysout, eqout, epout);
  
      wcssize (wcs, &cra, &cdec, &dra, &ddec);
      SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose);
  
      /* Make mag1 always the smallest magnitude */
      if (mag2 < mag1) {
  	mag = mag2;
  	mag2 = mag1;
  	mag1 = mag;
  	}
  
      /* Find Guide Star Catalog regions in which to search */
      nrmax = 100;
      rra1 = ra1;
      rra2 = ra2;
      rdec1 = dec1;
      rdec2 = dec2;
      RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epout,0.0,
  	    &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose);
      nreg = gscreg (refcat,rra1,rra2,rdec1,rdec2,table,nrmax,rlist,verbose);
      if (nreg <= 0) {
  	fprintf (stderr,"GSCBIN:  no Guide Star regions found\n");
  	return (0);
  	}
  
      /* Set keyword list */
      nk = 8;
      strcpy (kw[0].kname,"GSC_ID");
      strcpy (kw[1].kname,"RA_DEG");
      strcpy (kw[2].kname,"DEC_DEG");
      strcpy (kw[3].kname,"POS_ERR");
      strcpy (kw[4].kname,"MAG");
      strcpy (kw[5].kname,"MAG_ERR");
      strcpy (kw[6].kname,"MAG_BAND");
      strcpy (kw[7].kname,"CLASS");
      for (ik = 0; ik < nk; ik++) {
  	kw[ik].lname = (int) strlen (kw[ik].kname);
  	kw[ik].kn = 0;
  	kw[ik].kf = 0;
  	kw[ik].kl = 0;
  	}
      nstar = 0;
  
      /* Loop through region list */
      for (ireg = 0; ireg < nreg; ireg++) {
  	gscpath (refcat, rlist[ireg], inpath);
  
      /* Read size and keyword info from FITS table header */
  	kwn = kw;
  	ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead);
  
  	rnum = rlist[ireg];
  	num0 = 0;
  	rasum = 0.0;
  	decsum = 0.0;
  	msum = 0.0;
  	perrsum = 0.0;
  	merrsum = 0.0;
  	npos = 0;
  	num = 0;
  	fitsrtlset();
  	jstar = 0;
  	class = 0;
  
  	/* Loop through FITS table for this region */
  	for (itable = 0; itable <= ntable; itable++) {
  
  	    if (itable < ntable) {
  		nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry);
  		if (nbr < nbline) {
  		    fprintf (stderr,"GSCBIN: %d / %d bytes read, line %d / %d, region %d\n",
  			      nbr,nbline,itable,ntable,rnum);
  		    break;
  		    }
  
  		/* Extract selected fields */
  
  		/* Star number within region */
  		num0 = ftgeti4 (entry, &kw[0]);
  
  		/* Right ascension in degrees */
  		ra0 = ftgetr8 (entry, &kw[1]);
  
  		/* Declination in degrees */
  		dec0 = ftgetr8 (entry, &kw[2]);
  
  		/* Position error */
  		perr0 = ftgetr8 (entry, &kw[3]);
  
  		/* Magnitude */
  		mag0 = ftgetr8 (entry, &kw[4]);
  
  		/* Magnitude error */
  		merr0 = ftgetr8 (entry, &kw[5]);
  
  		/* Bandpass code */
  		band0 = ftgeti4 (entry, &kw[6]);
  
  		/* Object class code */
  		class0 = ftgeti4 (entry, &kw[7]);
  		}
  	    else
  		num0 = 0;
  
  	/* Compute mean position and magnitude for object */
  	    if (itable > 0 && npos > 0 &&
  		((classd < -1 && band != band0) ||
  		(classd < -1 && class != class0) || 
  		num != num0)) {
  
  		pass = 1;
  		if (perrsum == 0.0 || merrsum == 0)
  		    pass = 0;
  		else {
  		    ra = rasum / perrsum;
  		    dec = decsum / perrsum;
  		    mag = msum / merrsum;
  		    }
  
  		if (pass > 0 && classd > -1 && class != classd)
  		    pass = 0;
  
  		/* Check magnitude and position limits */
  		if (pass > 0 && mag1 != mag2 && (mag < mag1 || mag > mag2))
  		    pass = 0;
  
  		if (pass) {
  		    wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout);
  		
  		    /* Check position limits */
  		    ddist = wcsdist (cra,cdec,cra,dec);
  		    if (ddist > ddec)
  			pass = 0;
  		    rdist = wcsdist (cra,dec,ra,dec);
  		    if (rdist > dra)
  			pass = 0;
  		    }
  
  		/* Save star in FITS image */
  		if (pass) {
  		    class = class + (band * 100) + (npos * 10000);
  		    xnum = (double)rnum + (0.0001 * (double) num);
  		    wcs2pix (wcs, ra, dec, &xpix, &ypix, &offscl);
  		    if (!offscl) {
  			if (magscale > 0.0)
  			    flux = magscale * exp (logt * (-mag / 2.5));
  			else
  			    flux = 1.0;
  			ix = (int) (xpix + 0.5);
  			iy = (int) (ypix + 0.5);
  			addpix1 (image, bitpix, w,h, 0.0,1.0, xpix,ypix, flux);
  			nstar++;
  			jstar++;
  			}
  		    else {
  			ix = 0;
  			iy = 0;
  			}
  		    if (nlog == 1) {
  			fprintf (stderr,"GSCBIN: %04d.%04d: %9.5f %9.5f %s %d %d",
  				 rnum, num, ra, dec, cstr, class, npos);
  			if (magscale > 0.0)
  			    fprintf (stderr, " %5.2f", mag);
  			if (!offscl)
  			    flux = getpix1 (image, bitpix, w, h, 0.0, 1.0, ix, iy);
  			else
  			    flux = 0.0;
  			fprintf (stderr," (%d,%d): %f\n", ix, iy, flux);
  			}
  		    }
  
  	/* Reset star position for averaging */
  		rasum = 0.0;
  		decsum = 0.0;
  		msum = 0.0;
  		perrsum = 0.0;
  		merrsum = 0.0;
  		npos = 0;
  		}
  
  	/* Add information from current line to current object */
  
  	    /* Check object class */
       	    if ((classd > -1 && class0 == classd) ||
  		classd < -2 || (classd < 0 && class0 != 5)) {
  		perr = perr0;
  		perr2 = perr * perr;
  		if (perr2 <= 0.0) perr2 = 0.01;
  		rasum = rasum + (ra0 / perr2);
  		decsum = decsum + (dec0 / perr2);
  		perrsum = perrsum + (1.0 / perr2);
  		if (merr0 <= 0.0) merr0 = 0.01;
  		merr = merr0;
  		merr2 = merr * merr;
  		msum = msum + (mag0 / merr2);
  		merrsum = merrsum + (1.0 / merr2);
  		num = num0;
  		class = class0;
  		band = band0;
  		npos++;
  		}
  
  	    /* Log operation */
  	    if (nlog > 0 && itable%nlog == 0)
  		fprintf (stderr,"GSCBIN: %4d / %4d: %5d / %5d  / %5d sources, region %4d.%04d\r",
  			 ireg,nreg,jstar,itable,ntable,rlist[ireg],num0);
  
  	/* End of region */
  	    }
  
  	/* Close region input file */
  	(void) close (ift);
  	itot = itot + itable;
  	if (nlog > 0)
  	    fprintf (stderr,"GSCBIN: %4d / %4d: %5d / %5d  / %5d sources from region %4d    \n",
  		     ireg+1,nreg,jstar,itable,ntable,rlist[ireg]);
  	}
  
      /* Close output file and summarize transfer */
      if (nlog > 0) {
  	if (nreg > 1)
  	    fprintf (stderr,"GSCBIN: %d regions: %d / %d found\n",nreg,nstar,itot);
  	else
  	    fprintf (stderr,"GSCBIN: 1 region: %d / %d found\n",nstar,itable);
  	}
      return (nstar);
  }
  
  
  /* First region in each declination zone */
  int zreg1[24]={1,594,1178,1729,2259,2781,3246,3652,4014,4294, 4492,4615,
  	      4663,5260,5838,6412,6989,7523,8022,8464,8840,9134,9346,9490};
  
  /* Last region in each declination zone */
  int zreg2[24]={593,1177,1728,2258,2780,3245,3651,4013,4293,4491,4614,4662,
  	       5259,5837,6411,6988,7522,8021,8463,8839,9133,9345,9489,9537};
  
  /* Directory for each declination zone */
  char zdir[24][8]={"n0000","n0730","n1500","n2230","n3000","n3730","n4500",
  		"n5230","n6000","n6730","n7500","n8230","s0000","s0730",
  		"s1500","s2230","s3000","s3730","s4500","s5230","s6000",
  		"s6730","s7500","s8230"};
  
  static struct Keyword rkw[15];
  static int nrkw = 13;
  
  /* GSCREG -- search the HST Guide Star Catalog index table for fields
   * in the specified range of coordinates and magnitudes.  Build a
   * list containing the pathnames of the files on the cd-rom.
   */
  
  static int
  gscreg (refcat, ra1, ra2, dec1, dec2, table, nrmax, rgns, verbose)
  
  int	refcat;		/* Catalog code (GSC or GSCACT) */
  double	ra1, ra2;	/* Right ascension limits in degrees */
  double	dec1, dec2; 	/* Declination limits in degrees */
  char	*table;		/* Table data buffer */
  int	nrmax;		/* Maximum number of regions to find */
  int	*rgns;		/* Region numbers (returned)*/
  int	verbose;	/* 1 for diagnostics */
  
  {
      int nrgn;		/* Number of regions found (returned) */
      char tabpath[64];	/* Pathname for regions table */
      int nrows;		/* Number of entries in table */
      int nchar;		/* Number of characters per line in table */
      int nwrap;		/* 1 if 0h included in RA span*/
      int iwrap;
      struct Keyword *kwn;
      double gscra(), gscdec();
      int gsczone();
  
      char fitsline[120];
      int irow,iz1,iz2,ir1,ir2,jr1,jr2,i;
      int nsrch,nsrch1,nbhead,nbr;
      int ift;
      double ralow, rahi;
      double declow, dechi, decmin, decmax;
      int regnum;
  
      /* Set up keyword list for table entries to extract */
      strcpy (rkw[0].kname,"REG_NO");
      strcpy (rkw[1].kname,"RA_H_LOW");
      strcpy (rkw[2].kname,"RA_M_LOW");
      strcpy (rkw[3].kname,"RA_S_LOW");
      strcpy (rkw[4].kname,"RA_H_HI");
      strcpy (rkw[5].kname,"RA_M_HI");
      strcpy (rkw[6].kname,"RA_S_HI");
      strcpy (rkw[7].kname,"DECSI_LO");
      strcpy (rkw[8].kname,"DEC_D_LO");
      strcpy (rkw[9].kname,"DEC_M_LO");
      strcpy (rkw[10].kname,"DECSI_HI");
      strcpy (rkw[11].kname,"DEC_D_HI");
      strcpy (rkw[12].kname,"DEC_M_HI");
      rkw[13].kname[0] = 0;
  
      /* Add lengths of keywords to keyword structures */
      for (i = 0; i < 13; i++)
  	rkw[i].lname = strlen (rkw[i].kname);
  
      /* Initialize region list to zeroes */
      for (i = 0; i < nrmax; i++)
  	rgns[i] = 0;
  
      /* Initialize FITS table line buffer to zeroes */
      for (i = 0; i < 120; i++)
  	fitsline[i] = (char) 0;
  
      nrgn = 0;
  
      /* Set pathnames to Guide Star Catalog CDROM */
      if (refcat == GSCACT)
  	strcpy (tabpath,cdna);
      else
  	strcpy (tabpath,cdn);
  
      /* Set pathname for region table file */
      strcat (tabpath,"/tables/regions.tbl");
  
      /* Open the index table */
      kwn = rkw;
      ift = fitsrtopen (tabpath,&nrkw,&kwn, &nrows, &nchar, &nbhead);
  
      /* If the northern hemisphere CDROM cannot be read, try the southern */
      if (ift < 0) {
  	if (refcat == GSCACT)
  	    strcpy (tabpath,cdsa);
  	else
  	    strcpy (tabpath,cds);
  	strcat (tabpath,"/tables/regions.tbl");
  	ift = fitsrtopen (tabpath,&nrkw,&kwn,&nchar,&nrows,&nbhead);
  	if (ift < 0) {
  	    fprintf (stderr,"GSCREG:  error reading region table %s\n",tabpath);
  	    return (0);
  	    }
  	}
  
      /* Find region range to search based on declination */
      iz1 = gsczone (dec1);
      iz2 = gsczone (dec2);
      jr1 = 0;
      jr2 = 0;
      nwrap = 1;
  
      /* Search region northern hemisphere or only one region */
      if (iz2 >= iz1) {
  	ir1 = zreg1[iz1];
  	ir2 = zreg2[iz2];
  	}
  
      /* Search region in southern hemisphere with multiple regions */
      else if (dec1 < 0 && dec2 < 0) {
  	ir1 = zreg1[iz2];
  	ir2 = zreg2[iz1];
  	}
  
      /* Search region spans equator */
      else if (dec1 < 0 && dec2 >= 0) {
  	ir1 = zreg1[12];
  	ir2 = zreg2[iz1];
  	jr1 = zreg1[0];
  	jr2 = zreg2[iz2];
  	nwrap = 2;
  	}
      else {
  	ir1 = 1;
  	ir2 = 0;
  	}
  
      nsrch = ir2 - ir1 + 1;
      if (verbose)
  	fprintf (stderr,"GSCREG: searching %d regions: %d - %d\n",nsrch,ir1,ir2);
      if (jr1 > 0) {
  	nsrch1 = jr2 - jr1 + 1;
  	if (verbose)
  	    fprintf (stderr,"GSCREG: searching %d regions: %d - %d\n",nsrch1,jr1,jr2);
  	}
      if (verbose)
  	fprintf(stderr,"GSCREG: RA: %.5f - %.5f, Dec: %.5f - %.5f\n",ra1,ra2,dec1,dec2);
  
      if (nsrch < 1) {
  	return (0);
  	}
  
      nrgn = 0;
  
      for (iwrap = 0; iwrap < nwrap; iwrap++) {
  
  	for (irow = ir1 - 1; irow < ir2; irow++) {
  
  	/* Read next line of region table */
  	    nbr = fitsrtline (ift,nbhead,ltab,table,irow,nchar,fitsline);
  	    if (nbr < nchar) {
  		fprintf (stderr,"GSREG: %d / %d bytes read for row %d\n",nbr,nchar,irow);
  		break;
  		}
  
  	/* Declination range of the gs region */
  	/* note:  southern dechi and declow are reversed */
  	    dechi = gscdec (fitsline, 10, 11, 12);
  	    declow = gscdec (fitsline, 7, 8, 9);
  	    if (dechi > declow) {
  		decmin = declow;
  		decmax = dechi;
  		}
  	    else {
  		decmax = declow;
  		decmin = dechi;
  		}
  
  	    if (decmax >= dec1 && decmin <= dec2) {
  
  	    /* Right ascension range of the Guide Star Catalog region */
  		ralow = gscra (fitsline, 1, 2, 3);
  		rahi = gscra (fitsline, 4, 5, 6);
  		if (rahi <= 0.0) rahi = 360.0;
  
  	    /* Check RA if 0h RA not between region RA limits */
  		if (ra1 < ra2) {
  		    if (ralow <= ra2 && rahi >= ra1) {
  
  			/* Get region number from FITS table */
  			regnum = ftgeti4 (fitsline, &rkw[0]);
  			if (verbose)
  			    fprintf (stderr,"GSCREG: Region %d added to search\n",regnum);
  
  			/* Add this region to list, if there is space */
  			if (nrgn < nrmax) {
  			    rgns[nrgn] = regnum;
  			    nrgn = nrgn + 1;
  			    }
  			}
  		    }
  
  	    /* Check RA if 0h RA is between region RA limits */
  		else {
  		    if (ralow > rahi) rahi = rahi + 360.0;
  		    if (ralow <= ra2 || rahi >= ra1) {
  
  		    /* Get region number from FITS table */
  			regnum = ftgeti4 (fitsline, &rkw[0]);
  			if (verbose)
  			    fprintf (stderr,"GSCREG: Region %d added to search\n",regnum);
  
  		    /* Add this region to list, if there is space */
  			if (nrgn < nrmax) {
  			    rgns[nrgn] = regnum;
  			    nrgn = nrgn + 1;
  			    }
  			}
  		    }
  		}
  	    }
  
  	/* Handle wrap-around through the equator */
  	ir1 = jr1;
  	ir2 = jr2;
  	jr1 = 0;
  	jr2 = 0;
  	}
  
      (void) close (ift);
      return (nrgn);
  }
  
  
  /* GSCRA -- returns right ascension in degrees from the GSC index table
   *  This is converted from the hours, minutes, and seconds columns.
   */
  
  double
  gscra (fitsline, hcol, mcol, scol)
  
  char *fitsline;		/* Index table line */
  int hcol;		/* Column index for hours */
  int mcol;		/* Column index for minutes */
  int scol;		/* Column index for seconds */
  
  {
      double ra;		/* Right ascension in fractional degrees */
      double hrs;		/* Hours of right ascension */
      double min;		/* Minutes of right ascension */
      double sec;		/* Seconds of right ascension */
  
  /*  hours of right ascension */
      hrs = ftgetr8 (fitsline, &rkw[hcol]);
  
  /* minutes of right ascension */
      min = ftgetr8 (fitsline, &rkw[mcol]);
  
  /* seconds of right ascension */
      sec = ftgetr8 (fitsline, &rkw[scol]);
  
  /* right ascension in degrees */
      ra = hrs + (min / 60.0) + (sec / 3600.0);
      ra = ra * 15.0;
  
      return (ra);
  }
  
  
  /*  GSCDEC -- returns the declination in degrees from the GSC index table.
   *  This is converted from the sign, degrees, minutes, and seconds columns.
   */
  
  double
  gscdec (fitsline, sgncol, dcol, mcol)
  
  char *fitsline;		/* Index table line */
  int sgncol;		/* Column index for sign */
  int dcol;		/* Column index for degrees */
  int mcol;		/* Column index for minutes */
  
  {
  double	dec;		/* Declination in fractional degrees */
  
  char sgn[4];		/* Sign of declination */
  double deg;		/* Degrees of declination*/
  double min;		/* Minutes of declination */
  
      /* Get declination sign from table */
      (void) ftgetc (fitsline, &rkw[sgncol], sgn, 3);
  
      /* Get degrees of declination from table */
      deg = ftgetr8 (fitsline, &rkw[dcol]);
  
      /* Get minutes of declination from table */
      min = ftgetr8 (fitsline, &rkw[mcol]);
  
      dec = deg + (min / 60.0);
  
      /* Negative declination */
      if (strchr (sgn, '-') != NULL)
  	dec = -dec;
  
      return (dec);
  }
  
  
  /*  GSCZONE -- find the zone number where a declination can be found */
  
  int
  gsczone (dec)
  
  double dec;	/* declination in degrees */
  
  {
  int zone;		/* gsc zone (returned) */
  double	zonesize;
  int ndeczones = 12;	/* number of declination zones per hemisphere */
  
      /* Width of declination zones */
      zonesize = 90.0 / ndeczones;
  
      zone = ((int) (dec / zonesize));
      if (dec < 0)
  	zone = ndeczones - zone;
  	
      return (zone);
  }
  
  
  /* GSCPATH -- Get HST Guide Star Catalog region FITS file pathname */
  
  static void
  gscpath (refcat, regnum, path)
  
  int	refcat;		/* Catalog code (GSC or GSCACT) */
  int	regnum;		/* Guide Star Catalog region number */
  char	*path;		/* Pathname of GSC region FITS file */
  
  {
      int zone = 0;	/* Name of Guide Star Catalog zone directory */
      int i;
  
      /* Get zone directory name given region number */
      for (i = 0; i < 24; i++) {
  	if (regnum >= zreg1[i] && regnum <= zreg2[i]) {
  	    zone = i;
  	    break;
  	    }
  	}
  
      /* Set the pathname using the appropriate GSC CDROM directory */
  
      /* Northern hemisphere disk (volume 1) */
      if (regnum < zreg1[13]) {
  	if (refcat == GSCACT)
  	    sprintf (path,"%s/%s/%04d.gsc", cdna, zdir[zone], regnum);
  	else
  	    sprintf (path,"%s/gsc/%s/%04d.gsc", cdn, zdir[zone], regnum);
  	}
  
      /* Southern hemisphere disk (volume 2) */
      else {
  	if (refcat == GSCACT)
  	    sprintf (path,"%s/%s/%04d.gsc", cdsa, zdir[zone], regnum);
  	else
  	    sprintf (path,"%s/gsc/%s/%04d.gsc", cds, zdir[zone], regnum);
  	}
  
      return;
  }
  
  /* Feb 16 1996	New program
   * May 14 1996	Change arguments to FITSRTOPEN
   * May 17 1996	Fix bug so equal magnitude limits accepts anything
   * May 17 1996	Use new FITSRTOPEN which internalizes FITS header
   * May 24 1996	Fix string decoding bug and region search
   * May 31 1996	Use stream I/O
   * Jun 10 1996	Remove unused variables after using lint
   * Jul  1 1996	Fix GSC pathname
   * Aug  6 1996	Minor changes after lint
   * Sep 17 1996	Fix bug causing incomplete region access; improve logging
   * Oct 17 1996	Change FITS table entry calls
   * Nov 13 1996	Return no more than maximum star number
   * Nov 13 1996	Write all error messages to stderr with subroutine names
   * Nov 15 1996	Implement search radius; change input arguments
   * Dec 12 1996	Make logging run on single line per region
   * Dec 17 1996  Add code to keep brightest or closest stars if too many found
   * Dec 18 1996	Add code to get star information by number
   *
   * Mar 20 1997	Remove unused variables and fixed logging in GSCRNUM after lint
   * Oct 10 1997	Use standard I/O instead of stream I/O when reading FITS files
   * Nov  6 1997	Do not print overrun unless logging
   *
   * May 13 1998	Do not set classd in gscread()
   * May 13 1998	Print all stars if classd is < -1
   * May 27 1998	Include fitsio.h instead of fitshead.h
   * Jun 24 1998	Add string lengths to ra2str() and dec2str() calls
   * Aug  6 1998	Change fitsio.h to fitsfile.h
   * Sep 16 1998	Use limiting radius, if present
   * Sep 22 1998	Convert to desired output coordinate system
   * Oct 26 1998	Fix bug in region selection
   * Oct 29 1998	Correctly assign numbers when too many stars are found
   *
   * May 25 1999	Allocate table buffer only once
   * Jun 16 1999	Use SearchLim()
   * Aug 16 1999	Add RefLim() to get converted search coordinates right
   * Aug 16 1999	Fix bug to fix failure to search across 0:00 RA
   * Aug 25 1999	Return real number of stars from gscread()
   * Sep 10 1999	Set class selection with subroutine, not argument
   * Sep 16 1999	Fix bug which didn't always return closest stars
   * Sep 16 1999	Add distsort argument so brightest stars in circle works, too
   * Sep 22 1999	Rewrite table allocation so it works; make ltab static
   * Oct 21 1999	Include wcscat.h, unistd.h
   *
   * Mar 27 2000	Drop unused variables after lint
   * Mar 28 2000	Make default to read all classes
   * Jun 26 2000	Add coordinate system to SearchLim() arguments
   * Nov 29 2000	Add option to read cataog across the web
   * Dec 11 2000	Allow search engine URL in cdn[]
   * Dec 12 2000	Fix wrong web subroutine in gscrnum()
   *
   * Apr 24 2001	Add bandpass code and number of entries to object class
   * Apr 24 2001	Return individual entries if class is < -1
   * May 23 2001	Add support for GSC-ACT
   * May 30 2001	Use GSC_NORTH and GSCACT_NORTH instead of *_PATH for consistency
   * Jun 27 2001  Print stars as found in gscread() if nstarmax < 1
   * Jun 27 2001	Allocate distance array only if larger one is needed
   * Jun 27 2001	Add gscfree() to free table buffer and distance array
   * Sep 11 2001	Use single magnitude argument to gscread() and webread()
   * Sep 21 2001	Clean up web interface
   * Nov 20 2001	Change cos(degrad)) to cosdeg()
   *
   * Mar 28 2002	Change pathnames to /data/astrocat
   * Jul 31 2002	Always check for GSCACT_PATH, not GSCA_PATH
   * Oct  2 2002	Fix pass-through tab table header
   *
   * Mar 10 2003	Improve position limits
   * Apr  3 2003	Drop wrap in gscread(); add test for npos in gscread()
   * Apr 14 2003	Explicitly get revision date if nstarmax < 1
   * Aug 22 2003	Add dradi for inner edge of annulus search
   * Sep 25 2003	Add gscbin() to fill an image with sources
   * Oct  6 2003	Update gscread() and gscbin() for improved RefLim()
   * Nov 18 2003	Initialize image size and bits/pixel from header in gscbin()
   * Dec  1 2003	Add missing tab to n=-1 header
   * Dec 12 2003	Fix bug in wcs2pix() call in gscbin()
   *
   * Aug 27 2004	Include math.h
   *
   * Jun 20 2006	Initialize uninitialized variables
   * Nov 16 2006	Fix binning
   *
   * Jan 10 2007	Add match=1 argument to webrnum()
   * Jan 10 2007	Rewrite web access in gscread() and gscrnum() to reduce code
   *
   * Sep 22 2009	Initialize lengths of FITS table columns
   * Sep 22 2009	Change region table keywords from DEC*LOW to DEC*LO
   */
  