/* @(#)wrldco.c	17.1.1.1 (ESO-DMD) 01/25/02 17:40:37 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program 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 General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
  MA 02139, USA.
 
  Correspondence concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT:  Copyright (c) 1987 European Southern Observatory,
                                         all rights reserved
.IDENTIFIER  module wrldco.c
.LANGUAGE    C
.AUTHOR      K. Banse			IPG-ESO Garching
.KEYWORDS    
.COMMENTS    holds fp2wc 

.ENVIRONment MIDAS
             #include <midas_def.h>      Prototypes for MIDAS interfaces
 
.VERSIONS    1.00       951215  build on `fp2wc' from cutil.c  KB

010213			last modif
 
------------------------------------------------------------*/

/* Define _POSIX_SOURCE to indicate that this is a POSIX program */

#define  _POSIX_SOURCE 1

#include <stdio.h>
#include <stdlib.h>

#include <midas_def.h>
#include <wcs.h>

#define MAXDIM      4		/* synchronize with error message below */
#define MAXDIMSQ    MAXDIM*MAXDIM

/*

*/
 
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE     convert from Pixel to World Coordinate and vice versa
.INPUT/OUTPUT
  call as    fp2wc(flag,imno,coordin,coordout)

  input:
          int flag         = 0 for initialization,
          		   = 1 for conversion fp -> wc,
          		   = -1 for conversion wc -> fp
          int imno         image no. for flag = 0
                           unused if flag != 0
  in/output:
          double *coordin     frame pixel (flag = 1)
			      world coord (flag = -1)
          double *coordout    world coord (flag = 1)
                              frame pixel (flag = -1)

.RETURNS    status  0 = o.k. for flag != 0
                    0 or -1 for flag=0 to tell, if non- or linear coord system
                    > 0 if error for all flags
------------------------------------------------------------*/

#ifdef __STDC__
int  fp2wc(int flag,int imno,double *coordin,double *coordout)

#else
int  fp2wc(flag,imno,coordin,coordout)
int flag, imno;
double *coordin;
double *coordout;
#endif

{
int  naxis, stat, retstat;
int  ioff, mr, null, nval, unit;
int  e_c=1, e_log=0, e_disp=0;
register int  nr;
static int npix[MAXDIM];

static char  ctype[6][9];

double   tmpr[MAXDIM], phi=0, theta=0;
static double  start[MAXDIM], step[MAXDIM], refpix[MAXDIM], crval[MAXDIM];
static double  end[MAXDIM], pcval[MAXDIMSQ];

static struct wcsprm  wcs;
static struct linprm  lin;
static struct prjprm  prj;
static struct celprm  cel;



retstat = 0;

if (flag == 0)
   {
   char     cunit[100], pcnam[12];
   int      rstat, ln, ln1, ln2;
   double   longpole, latpole;


   (void) SCDRDI(imno,"NAXIS",1,1,&nval,&naxis,&unit,&null);
   if (naxis > MAXDIM) 
      {
      SCTPUT("WCO: max. 4 dimensions supported...");
      return (1);
      }

   (void) SCDRDI(imno,"NPIX",1,naxis,&nval,npix,&unit,&null);
   (void) SCDRDD(imno,"START",1,naxis,&nval,start,&unit,&null);
   (void) SCDRDD(imno,"STEP",1,naxis,&nval,step,&unit,&null);

   nval = (naxis+1)*16;
   (void) SCDRDC(imno,"CUNIT",1,1,nval,&ln,cunit,&unit,&null);
   mr = 16;
   for (nr=0; nr<naxis; nr++)
      {
      (void) strncpy(&ctype[nr][0],&cunit[mr],8);
      ctype[nr][8] = '\0';
      mr += 16;
      }

   SCECNT("GET",&ln,&ln1,&ln2);                    /* disable error abort */
   SCECNT("PUT",&e_c,&e_log,&e_disp);

   ioff = 0;
   for (nr=1; nr<(naxis+1); nr++)
      {
      for (mr=1; mr<(naxis+1); mr++)
         {
         (void) sprintf(pcnam,"PC00%d00%d",nr,mr);
         stat = SCDRDD(imno,pcnam,1,1,&nval,&pcval[ioff],&unit,&null);
         if (stat != ERR_NORMAL)
            {
            if (nr == mr)
               pcval[ioff] = 1.0;                  /* defaults of PC matrix */
            else
               pcval[ioff] = 0.0;
            }
         ioff++;
         }
      }

   stat = SCDRDD(imno,"LONGPOLE",1,1,&nval,&longpole,&unit,&null);
   if (stat != ERR_NORMAL) longpole = 999.0;
   stat = SCDRDD(imno,"LATPOLE",1,1,&nval,&latpole,&unit,&null);
   if (stat != ERR_NORMAL) latpole = 999.0;
   rstat = SCDRDD(imno,"REFPIX",1,naxis,&nval,refpix,&unit,&null);

   SCECNT("PUT",&ln,&ln1,&ln2);


   if (naxis == 1) 
      wcs.flag = 999;

   else
      {
      wcs.flag = 0;                           /* setup wcsprm structure */
      stat =  wcsset(naxis,ctype,&wcs);
      if (stat != 0)
         {
         SCTPUT("WCO: inconsistent projection types...");
         return (2);	/* inconsistent projection types */
         }
      }

   if (rstat != ERR_NORMAL)
      {					/* we didn't have descr. REFPIX */
      for (nr=0; nr<naxis; nr++)
         refpix[nr] = npix[nr]/2;
      if (wcs.flag != 999)
         SCTPUT
("Warning: No descr. REFPIX (FITS keyword CRVALn) found -> use center pixels");
      }

   for (nr=0; nr<naxis; nr++)	/* recompute CRVALn FITS keyword */
      {
      crval[nr] = start[nr] + (refpix[nr]-1)*step[nr];
      if (npix[nr] == 1)
         {
         if (fabs(step[nr]) < 1.e-35) 
            {
            step[nr] = 1.0;
            crval[nr] = start[nr];
            }
         else
            crval[nr] = start[nr] + (refpix[nr]-1)*step[nr];
         }
      }


   lin.flag = 0;                           /* setup linprm structure */
   lin.naxis = naxis;
   lin.crpix = refpix;
   lin.pc = pcval;
   lin.cdelt = step;

   cel.flag = 0;                           /* setup celprm structure */
   cel.ref[0] = crval[0];
   cel.ref[1] = crval[1];
   cel.ref[2] = longpole;
   cel.ref[3] = latpole;

   prj.flag = 0;                           /* setup prjprm structure */
   for (nr=0; nr<10; nr++)
      {
      prj.p[nr] = 0.0;
      prj.w[nr] = 0.0;
      }
   prj.r0 = 0.0;
 
   if (wcs.flag == 999) 
      {					/* get endpoints for later check */
      for (nr=0; nr<naxis; nr++)
         end[nr] = start[nr] + ((npix[nr]-1) * step[nr]);
      return (-1);		/* linear coord_sys */
      }
   }

else
   {
   if (flag > 0)		/* pixel coords => world coords */
      {				/* 1, 2, ...  counting */
      if (wcs.flag == 999)
         {				/* no projection type specified */
         for (nr=0; nr<lin.naxis; nr++)
            {
            coordout[nr] = start[nr] + ((coordin[nr]-1.0)*step[nr]);
            if ((coordin[nr] <= 0.0) || 
                (coordin[nr] > (double) npix[nr])) retstat = 55;
            }
         }

      else
         {
         stat = wcsrev(ctype,&wcs,coordin,&lin,tmpr,&prj,&phi,&theta,crval,
                       &cel,coordout);
         if (stat != 0) 
            {
            SCTPUT("WCO: problems with wcsrev routine...");
            retstat = 5;
            }
         }
      }

   else				/* world coords => pixel coords */
      {
      if (wcs.flag == 999)
         {				/* no projection type specified */
         for (nr=0; nr<lin.naxis; nr++) 
            {
            coordout[nr] = (coordin[nr] - start[nr])/step[nr] + 1.0;
            mr = CGN_DNINT(coordout[nr]);
            if ((mr < 1) || (mr > npix[nr])) retstat = 55;
            }
         }

      else
         {
         stat = wcsfwd(ctype,&wcs,coordin,crval,&cel,&phi,&theta,&prj,tmpr,
                       &lin,coordout);
         if (stat != 0) 
            {
            SCTPUT("WCO: problems with wcsfwd routine...");
            retstat = 5;
            }
         }
      }
   }

return retstat;
}