/*** File libwcs/imgetwcs.c
 *** June 24, 2016
 *** By Jessica Mink, jmink@cfa.harvard.edu (remotely based on UIowa code)
 *** Harvard-Smithsonian Center for Astrophysics
 *** Copyright (C) 1996-2016
 *** 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "wcs.h"
#include "lwcs.h"

/* Get the C* WCS fields in  a FITS header based on a reference catalog
 * do it by finding stars in the image and in the reference catalog and
 * finding the rotation and offsets which result in a best-fit.
 * verbose generates extra info on stderr.
 * try using deeper reference star catalog searches if there is trouble.
 * return 1 if all ok, else 0
 */

/* These parameters can be set on the command line */
static double secpix0 = PSCALE;		/* Set image scale--override header */
static double secpix2 = PSCALE;		/* Set image scale 2--override header */
static double *cd0 = NULL;		/* Set CD matrix--override header */
static double rot0 = 361.0;		/* Initial image rotation */
static int comsys = WCS_J2000;		/* Command line center coordinte system */
static int wp0 = 0;			/* Initial width of image */
static int hp0 = 0;			/* Initial height of image */
static double ra0 = -99.0;		/* Initial center RA in degrees */
static double dec0 = -99.0;		/* Initial center Dec in degrees */
static double xref0 = -99999.0;		/* Reference pixel X coordinate */
static double yref0 = -99999.0;		/* Reference pixel Y coordinate */
static int ptype0 = -1;			/* Projection type to fit */
static int  nctype = NWCSTYPE;		/* Number of possible projections */
static char ctypes[NWCSTYPE][4];	/* 3-letter codes for projections */
static int usecdelt = 0;		/* Use CDELT if 1, else CD matrix */
static char *dateobs0 = NULL;		/* Initial DATE-OBS value in FITS date format */

struct WorldCoor *ChangeFITSWCS();

/* Set a nominal world coordinate system from image header info.
 * If the image center is not FK5 (J2000) equinox, convert it
 * Return a WCS structure if OK, else return NULL
 */

struct WorldCoor *
GetFITSWCS (filename, header, verbose, cra, cdec, dra, ddec, secpix, wp, hp,
	    sysout, eqout)

char	*filename;	/* FITS or IRAF file name */
char	*header;	/* Image FITS header */
int	verbose;	/* Extra printing if =1 */
double	*cra;		/* Center right ascension in degrees (returned) */
double	*cdec;		/* Center declination in degrees (returned) */
double	*dra;		/* Right ascension half-width in degrees (returned) */
double	*ddec;		/* Declination half-width in degrees (returned) */
double	*secpix;	/* Arcseconds per pixel (returned) */
int	*wp;		/* Image width in pixels (returned) */
int	*hp;		/* Image height in pixels (returned) */
int	*sysout;	/* Coordinate system to return (0=image, returned) */
double	*eqout;		/* Equinox to return (0=image, returned) */
{
    int naxes;
    double eq1, x, y;
    double ra1, dec1, dx, dy;
    double xmin, xmax, ymin, ymax, ra2, dec2, ra3, dec3, ra4, dec4;
    double dra0, dra1, dra2, dra3, dra4;
    struct WorldCoor *wcs;
    char rstr[64], dstr[64], cstr[16];

    /* Initialize WCS structure from possibly revised FITS header */
    wcs = ChangeFITSWCS (filename, header, verbose);
    if (wcs == NULL) {
	return (NULL);
	}
    *hp = (int) wcs->nypix;
    *wp = (int) wcs->nxpix;

    /* If incomplete WCS in header, drop out */
    if (nowcs (wcs)) {
	setwcsfile (filename);
	/* wcserr(); */
	if (verbose)
	    fprintf (stderr,"Insufficient information for initial WCS\n");
	return (NULL);
	}

    /* If in linear coordinates, do not print as sexigesimal */
    if (wcs->sysout < 1 || wcs->sysout == 6 || wcs->sysout == 10)
	wcs->degout = 1;

    /* Set flag to get appropriate equinox for catalog search */
    if (!*sysout)
	*sysout = wcs->syswcs;
    if (*eqout == 0.0)
	*eqout = wcs->equinox;
    eq1 = wcs->equinox;
    if (wcs->coorflip) {
	ra1 = wcs->crval[1];
	dec1 = wcs->crval[0];
	}
    else {
	ra1 = wcs->crval[0];
	dec1 = wcs->crval[1];
	}

    /* Print reference pixel position and value */
    if (verbose && (eq1 != *eqout || wcs->syswcs != *sysout)) {
	if (wcs->degout) {
	    deg2str (rstr, 32, ra1, 6);
	    deg2str (dstr, 32, dec1, 6);
	    }
	else {
	    ra2str (rstr, 32, ra1, 3);
	    dec2str (dstr, 32, dec1, 2);
	    }
	wcscstr (cstr, wcs->syswcs, wcs->equinox, wcs->epoch);
	fprintf (stderr,"Reference pixel (%.2f,%.2f) %s %s %s\n",
		 wcs->xrefpix, wcs->yrefpix, rstr, dstr, cstr);
	}

    /* Get coordinates of corners for size for catalog searching */
    dx = wcs->nxpix;
    dy = wcs->nypix;
    xmin = 0.5;
    ymin = 0.5;
    xmax = 0.5 + dx;
    ymax = 0.5 + dy;
    pix2wcs (wcs, xmin, ymin, &ra1, &dec1);
    pix2wcs (wcs, xmin, ymax, &ra2, &dec2);
    pix2wcs (wcs, xmax, ymin, &ra3, &dec3);
    pix2wcs (wcs, xmax, ymax, &ra4, &dec4);

    /* Convert search corners to output coordinate system and equinox */
    if (wcs->syswcs > 0 && wcs->syswcs != 6 && wcs->syswcs != 10) {
	wcscon (wcs->syswcs,*sysout,wcs->equinox,*eqout,&ra1,&dec1,wcs->epoch);
	wcscon (wcs->syswcs,*sysout,wcs->equinox,*eqout,&ra2,&dec2,wcs->epoch);
	wcscon (wcs->syswcs,*sysout,wcs->equinox,*eqout,&ra3,&dec3,wcs->epoch);
	wcscon (wcs->syswcs,*sysout,wcs->equinox,*eqout,&ra4,&dec4,wcs->epoch);
	}

    /* Find center and convert to output coordinate system and equinox */
    x = 0.5 + (dx * 0.5);
    y = 0.5 + (dy * 0.5);
    pix2wcs (wcs, x, y, cra, cdec);
    if (wcs->syswcs > 0 && wcs->syswcs != 6 && wcs->syswcs != 10)
	wcscon (wcs->syswcs,*sysout,wcs->equinox,*eqout,cra,cdec,wcs->epoch);

    /* Find maximum half-width in declination */
    *ddec = fabs (dec1 - *cdec);
    if (fabs (dec2 - *cdec) > *ddec)
	*ddec = fabs (dec2 - *cdec);
    if (fabs (dec3 - *cdec) > *ddec)
	*ddec = fabs (dec3 - *cdec);
    if (fabs (dec4 - *cdec) > *ddec)
	*ddec = fabs (dec4 - *cdec);

    /* Find maximum half-width in right ascension */
    dra0 = (dx / dy) * (*ddec / cos (*cdec));
    dra1 = ra1 - *cra;
    dra2 = ra2 - *cra;
    if (*cra < 0 && *cra + dra0 > 0.0) {
	dra1 = -(dra1 - 360.0);
	dra2 = -(dra2 - 360.0);
	}
    if (dra1 > 180.0)
	dra1 = dra1 - 360.0;
    else if (dra1 < -180.0)
	dra1 = dra1 + 360.0;
    else if (dra1 < 0.0)
	dra1 = -dra1;
    if (dra2 > 180.0)
	dra2 = dra2 - 360.0;
    else if (dra2 < -180.0)
	dra2 = dra2 + 360.0;
    else if (dra2 < 0.0)
	dra2 = -dra2;
    dra3 = *cra - ra3;
    dra4 = *cra - ra4;
    if (*cra > 0 && *cra - dra0 < 0.0) {
	dra3 = dra3 + 360.0;
	dra4 = dra4 + 360.0;
	}
    if (dra3 > 180.0)
	dra3 = dra3 - 360.0;
    else if (dra3 < -180.0)
	dra3 = dra3 + 360.0;
    else if (dra3 < 0.0)
	dra3 = -dra3;
    if (dra4 > 180.0)
	dra4 = dra4 - 360.0;
    else if (dra4 < -180.0)
	dra4 = dra4 + 360.0;
    else if (dra4 < 0.0)
	dra4 = -dra4;
    *dra = dra1;
    if (dra2 > *dra)
	*dra = dra2;
    if (dra3 > *dra)
	*dra = dra3;
    if (dra4 > *dra)
	*dra = dra4;

    /* wcssize (wcs, cra, cdec, dra, ddec); */

    /* Set reference pixel to center of image if it has not been set */
    if (wcs->xref == -999.0 && wcs->yref == -999.0) {
	wcs->xref = *cra;
	wcs->cel.ref[0] = *cra;
	wcs->crval[0] = *cra;
	wcs->yref = *cdec;
	wcs->cel.ref[1] = *cdec;
	wcs->crval[1] = *cdec;
	ra1 = *cra;
	dec1 = *cdec;
	if (wcs->xrefpix == 0.0 && wcs->yrefpix == 0.0) {
	    wcs->xrefpix = 0.5 + (double) wcs->nxpix * 0.5;
	    wcs->yrefpix = 0.5 + (double) wcs->nypix * 0.5;
	    }
	wcs->xinc = *dra * 2.0 / (double) wcs->nxpix;
	wcs->yinc = *ddec * 2.0 / (double) wcs->nypix;
	/* hchange (header,"PLTRAH","PLT0RAH");
	wcs->plate_fit = 0; */
	}

    /* Convert center to desired coordinate system */
    else if (wcs->syswcs != *sysout && wcs->equinox != *eqout) {
	wcscon (wcs->syswcs, *sysout, wcs->equinox, *eqout, &ra1, &dec1, wcs->epoch);
	if (wcs->coorflip) {
	    wcs->yref = ra1;
	    wcs->xref = dec1;
	    }
	else {
	    wcs->xref = ra1;
	    wcs->yref = dec1;
	    }
	}

    /* Compute plate scale to return if it was not set on the command line */
    if (secpix0 <= 0.0) {
	pix2wcs (wcs, wcs->xrefpix-0.5, wcs->yrefpix, &ra1, &dec1);
	pix2wcs (wcs, wcs->xrefpix+0.5, wcs->yrefpix, &ra2, &dec2);
	*secpix = 3600.0 * wcsdist (ra1, dec1, ra2, dec2);
	}

    wcs->crval[0] = wcs->xref;
    wcs->crval[1] = wcs->yref;
    if (wcs->coorflip) {
	wcs->cel.ref[0] = wcs->crval[1];
	wcs->cel.ref[1] = wcs->crval[0];
	}
    else {
	wcs->cel.ref[0] = wcs->crval[0];
	wcs->cel.ref[1] = wcs->crval[1];
	}

    if (wcs->syswcs > 0 && wcs->syswcs != 6 && wcs->syswcs != 10) {
	wcs->cel.flag = 0;
	wcs->wcsl.flag = 0;
	}
    else {
	wcs->lin.flag = LINSET;
	wcs->wcsl.flag = WCSSET;
	}

    wcs->equinox = *eqout;
    wcs->syswcs = *sysout;
    wcs->sysout = *sysout;
    wcs->eqout = *eqout;
    wcs->sysin = *sysout;
    wcs->eqin = *eqout;
    wcscstr (cstr,*sysout,*eqout,wcs->epoch);
    strcpy (wcs->radecsys, cstr);
    strcpy (wcs->radecout, cstr);
    strcpy (wcs->radecin, cstr);
    wcsininit (wcs, wcs->radecsys);
    wcsoutinit (wcs, wcs->radecsys);

    naxes = wcs->naxis;
    if (naxes < 1 || naxes > 9) {
	naxes = wcs->naxes;
	wcs->naxis = naxes;
	}

    if (usecdelt) {
	hputnr8 (header, "CDELT1", 9, wcs->xinc);
	if (naxes > 1) {
	    hputnr8 (header, "CDELT2", 9, wcs->yinc);
	    hputnr8 (header, "CROTA2", 9, wcs->rot);
	    }
	hdel (header, "CD1_1");
	hdel (header, "CD1_2");
	hdel (header, "CD2_1");
	hdel (header, "CD2_2");
	}
    else {
	hputnr8 (header, "CD1_1", 9, wcs->cd[0]);
	if (naxes > 1) {
	    hputnr8 (header, "CD1_2", 9, wcs->cd[1]);
	    hputnr8 (header, "CD2_1", 9, wcs->cd[2]);
	    hputnr8 (header, "CD2_2", 9, wcs->cd[3]);
	    }
	}

    /* Print reference pixel position and value */
    if (verbose) {
	if (wcs->degout) {
	    deg2str (rstr, 32, ra1, 6);
            deg2str (dstr, 32, dec1, 6);
	    }
	else {
	    ra2str (rstr, 32, ra1, 3);
            dec2str (dstr, 32, dec1, 2);
	    }
	wcscstr (cstr,*sysout,*eqout,wcs->epoch);
	fprintf (stderr,"Reference pixel (%.2f,%.2f) %s %s %s\n",
		wcs->xrefpix, wcs->yrefpix, rstr, dstr, cstr);
	}

    /* Image size for catalog search */
    if (verbose) {
	if (wcs->degout) {
	    deg2str (rstr, 32, *cra, 6);
            deg2str (dstr, 32, *cdec, 6);
	    }
	else {
	    ra2str (rstr, 32, *cra, 3);
	    dec2str (dstr, 32, *cdec, 2);
	    }
	wcscstr (cstr, *sysout, *eqout, wcs->epoch);
	fprintf (stderr,"Search at %s %s %s", rstr, dstr, cstr);
	if (wcs->degout) {
	    deg2str (rstr, 32, *dra, 6);
            deg2str (dstr, 32, *ddec, 6);
	    }
	else {
	    ra2str (rstr, 32, *dra, 3);
	    dec2str (dstr, 32, *ddec, 2);
	    }
	fprintf (stderr," +- %s %s\n", rstr, dstr);
	fprintf (stderr,"Image width=%d height=%d, %g arcsec/pixel\n",
				*wp, *hp, *secpix);
	}

    return (wcs);
}


/* ChangeFITSWCS: Modify FITS WCS header from command line values */

struct WorldCoor *
ChangeFITSWCS (filename, header, verbose)

char	*filename;	/* FITS or IRAF file name */
char	*header;	/* Image FITS header */
int	verbose;	/* Extra printing if =1 */
{
    int nax, i, hp, wp;
    double xref, yref, degpix, secpix;
    struct WorldCoor *wcs;
    char temp[16];
    char *cwcs;

    /* Set the world coordinate system from the image header */
    if (strlen (filename) > 0) {
	cwcs = strchr (filename, '%');
	if (cwcs != NULL)
	    cwcs++;
	}
    if (!strncmp (header, "END", 3)) {
	cwcs = NULL;
	for (i = 0; i < 2880; i++)
	    header[i] = (char) 32;
	hputl (header, "SIMPLE", 1);
	hputi4 (header, "BITPIX", 0);
	hputi4 (header, "NAXIS", 2);
	hputi4 (header, "NAXIS1", 1);
	hputi4 (header, "NAXIS2", 1);
	}

    /* Set image dimensions */
    nax = 0;
    if (hp0 > 0 || wp0 > 0) {
	hp = hp0;
	wp = wp0;
	if (hp > 0 && wp > 0)
	    nax = 2;
	else
	    nax = 1;
	hputi4 (header, "NAXIS", nax);
	hputi4 (header, "NAXIS1", wp);
	hputi4 (header, "NAXIS2", hp);
	}
    else if (hgeti4 (header,"NAXIS",&nax) < 1 || nax < 1) {
	if (hgeti4 (header, "WCSAXES", &nax) < 1)
	    return (NULL);
	else {
	    if (hgeti4 (header, "IMAGEW", &wp) < 1)
		return (NULL);
	    if (hgeti4 (header, "IMAGEH", &hp) < 1)
		return (NULL);
	    }
	}
    else {
	if (hgeti4 (header,"NAXIS1",&wp) < 1)
	    return (NULL);
	if (hgeti4 (header,"NAXIS2",&hp) < 1)
	    return (NULL);
	}

    /* Set plate center from command line, if it is there */
    if (ra0 > -99.0 && dec0 > -99.0) {
	hputnr8 (header, "CRVAL1" ,8,ra0);
	hputnr8 (header, "CRVAL2" ,8,dec0);
	hputra (header, "RA", ra0);
	hputdec (header, "DEC", dec0);
	if (comsys == WCS_B1950) {
	    hputi4 (header, "EPOCH", 1950);
	    hputi4 (header, "EQUINOX", 1950);
	    hputs (header, "RADECSYS", "FK4");
	    }
	else {
	    hputi4 (header, "EPOCH", 2000);
	    hputi4 (header, "EQUINOX", 2000);
	    if (comsys == WCS_GALACTIC)
		hputs (header, "RADECSYS", "GALACTIC");
	    else if (comsys == WCS_ECLIPTIC)
		hputs (header, "RADECSYS", "ECLIPTIC");
	    else if (comsys == WCS_ICRS)
		hputs (header, "RADECSYS", "ICRS");
	    else
		hputs (header, "RADECSYS", "FK5");
	    }
	if (hgetr8 (header, "SECPIX", &secpix)) {
	    degpix = secpix / 3600.0;
	    hputnr8 (header, "CDELT1", 8, -degpix);
	    hputnr8 (header, "CDELT2", 8, degpix);
	    hdel (header, "CD1_1");
	    hdel (header, "CD1_2");
	    hdel (header, "CD2_1");
	    hdel (header, "CD2_2");
	    } 
	}
    if (ptype0 > -1 && ptype0 < nctype) {
	strcpy (temp,"RA---");
	strcat (temp, ctypes[ptype0]);
	hputs (header, "CTYPE1", temp);
	strcpy (temp,"DEC--");
	strcat (temp, ctypes[ptype0]);
	hputs (header, "CTYPE2", temp);
	}

    /* Set reference pixel from command line, if it is there */
    if (xref0 > -99999.0 && yref0 > -99999.0) {
	hputr8 (header, "CRPIX1", xref0);
	hputr8 (header, "CRPIX2", yref0);
	}
    else if (hgetr8 (header, "CRPIX1", &xref) < 1) {
	xref = 0.5 + (double) wp / 2.0;
	yref = 0.5 + (double) hp / 2.0;
	hputnr8 (header, "CRPIX1", 3, xref);
	hputnr8 (header, "CRPIX2", 3, yref);
	}

    /* Set plate scale from command line, if it is there */
    if (secpix0 != 0.0 || cd0 != NULL) {
	if (secpix2 != 0.0) {
	    secpix = 0.5 * (secpix0 + secpix2);
	    hputnr8 (header, "SECPIX1", 5, secpix0);
	    hputnr8 (header, "SECPIX2", 5, secpix2);
	    degpix = -secpix0 / 3600.0;
	    hputnr8 (header, "CDELT1", 8, degpix);
	    degpix = secpix2 / 3600.0;
	    hputnr8 (header, "CDELT2", 8, degpix);
	    hdel (header, "CD1_1");
	    hdel (header, "CD1_2");
	    hdel (header, "CD2_1");
	    hdel (header, "CD2_2");
	    }
	else if (secpix0 != 0.0) {
	    secpix = secpix0;
	    hputnr8 (header, "SECPIX", 5, secpix);
	    degpix = secpix / 3600.0;
	    hputnr8 (header, "CDELT1", 8, -degpix);
	    hputnr8 (header, "CDELT2", 8, degpix);
	    hdel (header, "CD1_1");
	    hdel (header, "CD1_2");
	    hdel (header, "CD2_1");
	    hdel (header, "CD2_2");
	    }
	else {
	    hputr8 (header, "CD1_1", cd0[0]);
	    hputr8 (header, "CD1_2", cd0[1]);
	    hputr8 (header, "CD2_1", cd0[2]);
	    hputr8 (header, "CD2_2", cd0[3]);
	    hdel (header, "CDELT1");
	    hdel (header, "CDELT2");
	    hdel (header, "CROTA1");
	    hdel (header, "CROTA2");
	    }
	if (!ksearch (header,"CRVAL1")) {
	    hgetra (header, "RA", &ra0);
	    hgetdec (header, "DEC", &dec0);
	    hputnr8 (header, "CRVAL1", 8, ra0);
	    hputnr8 (header, "CRVAL2", 8, dec0);
	    }
	if (!ksearch (header,"CRPIX1")) {
	    xref = (double) wp / 2.0;
	    yref = (double) hp / 2.0;
	    hputnr8 (header, "CRPIX1", 3, xref);
	    hputnr8 (header, "CRPIX2", 3, yref);
	    }
	if (!ksearch (header,"CTYPE1")) {
	    if (comsys == WCS_GALACTIC) {
		hputs (header, "CTYPE1", "GLON-TAN");
		hputs (header, "CTYPE2", "GLAT-TAN");
		}
	    else {
		hputs (header, "CTYPE1", "RA---TAN");
		hputs (header, "CTYPE2", "DEC--TAN");
		}
	    }
	}

    /* Set rotation angle from command line, if it is there */
    if (rot0 < 361.0) {
	hputnr8 (header, "CROTA1", 5, rot0);
	hputnr8 (header, "CROTA2", 5, rot0);
	}

    /* Set observation date for epoch, if it is there */
    if (dateobs0 != NULL)
	hputs (header, "DATE-OBS", dateobs0);

    /* Initialize WCS structure from FITS header */
    wcs = wcsinitn (header, cwcs);

    /* If incomplete WCS in header, drop out */
    if (nowcs (wcs)) {
	setwcsfile (filename);
	/* wcserr(); */
	if (verbose)
	    fprintf (stderr,"Insufficient information for initial WCS\n");
	return (NULL);
	}
    return (wcs);
}

void
setcdelt()			/* Set flag to use CDELTn, not CD matrix */
{usecdelt = 1; return;}

void
setnpix (nx, ny)		/* Set image size */
int nx, ny;
{ wp0 = nx; hp0 = ny; return; }

void
setrot (rot)
double rot;
{ rot0 = rot; return; }

void
setsecpix (secpix)		/* Set first axis arcseconds per pixel */
double secpix;
{ secpix0 = secpix; return; }

double
getsecpix ()		/* Return first axis arcseconds per pixel */
{ return (secpix0); }

void
setsecpix2 (secpix)		/* Set second axis arcseconds per pixel */
double secpix;
{ secpix2 = secpix; return; }

void
setcd (cd)		/* Set initial CD matrix */
double *cd;
{ int i;
  if (cd0 != NULL) free (cd0);
  cd0 = (double *) calloc (4, sizeof (double));
  for (i = 0; i < 4; i++) cd0[i] = cd[i];
  return; }

void
setsys (comsys0)		/* Set WCS coordinates as FK4 */
int comsys0;
{ comsys = comsys0; return; }

void
setcenter (rastr, decstr)	/* Set center sky coordinates in strings */
char *rastr, *decstr;
{
    ra0 = str2ra (rastr);
    dec0 = str2dec (decstr);
    return;
}

void
setdcenter (ra, dec)		/* Set center sky coordinates in degrees */
double ra, dec;
{ ra0 = ra; dec0 = dec; return; }

void
getcenter (ra, dec)		/* Return initial reference sky coordinates */
double *ra, *dec;
{ *ra = ra0; *dec = dec0; return; }

void
setrefpix (x, y)		/* Set reference pixel image coordinates */
double x, y;
{ xref0 = x; yref0 = y; return; }

void
getrefpix (x, y)		/* Return initial ref pixel image coordinates */
double *x, *y;
{ *x = xref0; *y = yref0; return; }

void
setproj (ptype)
char*	ptype;
{
    int i;

    /* Set up array of projection types */
    strcpy (ctypes[0], "LIN");
    strcpy (ctypes[1], "AZP");
    strcpy (ctypes[2], "SZP");
    strcpy (ctypes[3], "TAN");
    strcpy (ctypes[4], "SIN");
    strcpy (ctypes[5], "STG");
    strcpy (ctypes[6], "ARC");
    strcpy (ctypes[7], "ZPN");
    strcpy (ctypes[8], "ZEA");
    strcpy (ctypes[9], "AIR");
    strcpy (ctypes[10], "CYP");
    strcpy (ctypes[11], "CAR");
    strcpy (ctypes[12], "MER");
    strcpy (ctypes[13], "CEA");
    strcpy (ctypes[14], "COP");
    strcpy (ctypes[15], "COD");
    strcpy (ctypes[16], "COE");
    strcpy (ctypes[17], "COO");
    strcpy (ctypes[18], "BON");
    strcpy (ctypes[19], "PCO");
    strcpy (ctypes[20], "SFL");
    strcpy (ctypes[21], "PAR");
    strcpy (ctypes[22], "AIT");
    strcpy (ctypes[23], "MOL");
    strcpy (ctypes[24], "CSC");
    strcpy (ctypes[25], "QSC");
    strcpy (ctypes[26], "TSC");
    strcpy (ctypes[27], "NCP");
    strcpy (ctypes[28], "GLS");
    strcpy (ctypes[29], "DSS");
    strcpy (ctypes[30], "PLT");
    strcpy (ctypes[31], "TNX");
    strcpy (ctypes[32], "ZPX");
    strcpy (ctypes[33], "TPV");

    ptype0 = -1;
    ptype[3] = (char) 0;
    for (i = 0; i < nctype; i++) {
	if (!strcasecmp (ptype, ctypes[i]))
	    ptype0 = i;
	}
    return;
}

void
setdateobs (dateobs)
char *dateobs;
{ dateobs0 = calloc (strlen (dateobs), sizeof (char));
  strcpy (dateobs0, dateobs);
  return; }


/* Feb 29 1996	New program
 * May 23 1996	Use pre-existing WCS for center, if it is present
 * May 29 1996	Simplify program by always using WCS structure
 * Jun 12 1996	Be more careful with nominal WCS setting
 * Jul  3 1996	Set epoch from old equinox if not already set
 * Jul 19 1996	Set image center in WCS if DSS WCS
 * Aug  5 1996	Check for SECPIX1 as well as SECPIX
 * Aug  7 1996	Save specified number of decimal places in header parameters
 * Aug  7 1996	Rename old center parameters
 * Aug 26 1996	Decrease default pixel tolerance from 20 to 10
 * Sep  1 1996	Set plate scale default in lwcs.h
 * Sep  3 1996	Fix bug to set plate scale from command line
 * Oct 15 1996	Break off from imsetwcs.c
 * Oct 16 1996	Clean up center setting so eqref is used
 * Oct 17 1996	Do not print error messages unless verbose is set
 * Oct 30 1996	Keep equinox from image if EQREF is zero
 * Nov  1 1996	Declare undeclared subroutines; remove unused variables
 * Nov  4 1996	Add reference pixel and projection to wcsset() call
 * Nov 14 1996	Add GetLimits() to deal with search limits around the poles
 * Nov 15 1996	Drop GetLimits(); code moved to individual catalog routines
 * Dec 10 1996	Fix precession and make equinox double

 * Feb 19 1997	If eqout is 0, use equinox of image
 * Feb 24 1997	Always convert center to output equinox (bug fix)
 * Mar 20 1997	Declare EQ2 double instead of int, fixing a bug
 * Jul 11 1997	Allow external (command line) setting of reference pixel coords
 * Sep 26 1997	Set both equinox and epoch if input center coordinates
 * Nov  3 1997	Separate WCS reference pixel from search center
 * Dec  8 1997	Set CDELTn using SECPIX if it is in the header
 *
 * Jan  6 1998	Do not print anything unless verbose is set
 * Jan 29 1998	Use flag to allow AIPS classic WCS subroutines
 * Mar  1 1998	Set x and y plate scales from command line if there are two
 * Mar  2 1998	Do not reset plate solution switch
 * Mar  6 1998	Add option to reset center sky coordinates in degrees
 * Mar  6 1998	Add option to set projection type
 * Mar 18 1998	Initialize reference pixel if CDELTn's being set
 * Mar 20 1998	Only set CTYPEn if CDELTn or CRVALn are set
 * Apr 10 1998	Add option to set search area to circle or box
 * Apr 20 1998	Move GetArea() to scat.c
 * Apr 24 1998	Always convert image reference coordinate to catalog equinox
 * Apr 28 1998	Change coordinate system flags to WCS_*
 * Jun  1 1998	Print error message if WCS cannot be initialized
 * Jun 24 1998	Add string lengths to ra2str() and dec2str() calls
 * Jun 25 1998	Leave use of AIPS wcs to wcs.c file
 * Sep 17 1998	Add sysout to argument list and use scscon() for conversion
 * Sep 25 1998	Make sysout==0 indicate output in image coordinate system
 * Oct 28 1998	Set coordinate system properly to sysout/eqout in GetFITSWCS()
 *
 * Apr  7 1999	Add filename argument to GetFITSWCS
 * Apr 29 1999	Add option to set image size
 * Jun  2 1999	Fix sign of CDELT1 if secpix2 and secpix0 are set
 * Jul  7 1999	Fix conversion of center coordinates to refsys
 * Jul  9 1999	Fix bug which reset command-line-set reference pixel coordinate
 * Oct 21 1999	Fix declarations after lint
 * Nov  1 1999	Add option to write CD matrix
 * Nov  1 1999	If CDELTn set from command line delete previous header CD matrix
 * Nov 12 1999	Add galactic coordinates as command line option
 * Nov 16 1999	Set radecsys correctly for command line galactic
 *
 * Feb 15 2000	Add option to override the header CD matrix (like CDELTs)
 * Feb 29 2000	Fix bug, converting reference pixel WCS coordinates everywhere
 * Mar 27 2000	Drop unused subroutine setradius()
 * May 24 2000	Print degrees in debugging messages if output format
 * Sep 14 2000	Set xinc and yinc correctly if center pixel in header
 *
 * Jan 11 2001	All printing to stderr
 * Sep 19 2001	Drop fitshead.h; it is in wcs.h
 * Oct 19 2001	Allow DATE-OBS to be set
 *
 * Apr  3 2002	Update projection types to match list in wcs.h and wcs.c
 * Aug  2 2002	Add getsecpix(), getrefpix(), getcneter() to return presets
 *
 * Mar 25 2003	Fix GetFITSWCS() to return correct search area for rotated image
 * Mar 25 2003	Write out CTYPEn with quotes
 * Mar 27 2003	Fix half-pixel bug in computation of image center
 * Apr  3 2003	Drop unused variables after lint
 * Jun  6 2003	Set xref and yref to center if -999, not 0
 * Jul 21 2003	Fix bug setting secpix if it was not set on the command line
 *		(found by Takehiko Wada, ISAS)
 * Sep 23 2003	In setproj(), use strcasecmp() instead of strcmp()
 * Sep 26 2003	If reference pixel not set, set center correctly
 * Oct  6 2003	Change wcs->naxes to wcs->naxis to match WCSLIB 3.2
 * Dec  3 2003	Add wcs->naxes back as an alternative
 *
 * Jul 26 2004	Fix image size when wrapping around RA=0:00
 * Sep 16 2004	Fix verbose mode search size in GetFITSWCS()
 * Oct 29 2004	Fix problem setting RA size from limits
 *
 * Jan 20 2005	Fix cel.ref assignment if axis are switched
 * Jul 20 2005	Fix bug which reversed dimensions when setting image size
 * Jul 21 2005	Fix bug which caused bad results at RA ~= 0.0
 * Aug 30 2005	Implement multiple WCS's, though not modification thereof, yet
 * Nov  1 2005	Set RADECSYS to ICRS if appropriate
 *
 * Oct 30 2006	Do not precess LINEAR or XY coordinates
 *
 * Jun  1 2007	In GetFITSWCS, deal with no input file
 * Jun  5 2007	Add ChangeFITSWCS to set header WCS arguments and WCS
 * Jul  3 2007	Fix bug by setting hp and wp
 * Jul 26 2007	If first line of header is END, initialize other needed values
 * Oct 19 2007	Return NULL from GetFITSWCS() immediately if no WCS in header
 *
 * Mar 24 2009	Set dimensions from IMAGEW and IMAGEH if WCSAXES > 0
 *
 * Apr 06 2010	Set hp from IMAGEH in ChangeFITSWCS() (from Paul Liptack)
 * Apr  7 2010	In ChangeFITSWCS() set number of WCS projections from NWCSTYPE
 *
 * Sep  1 2011	Add ZPX and TPV projections to setproj()
 *
 * Jun 24 2016	ptype contains only first 3 letters of projection codes
 */