/* @(#)scfaux.c	17.1.1.1 (ESO-DMD) 01/25/02 17:36:27 */
/*===========================================================================
  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
===========================================================================*/

/*+++++++++++++++++++++++ extension module SCFAUX +++++++++++++++++++++++
.LANGUAGE C
.IDENTIFICATION Module SCFAUX
.AUTHOR         K. Banse 		ESO - Garching
.KEYWORDS       standard interfaces, bulk data files
.ENVIRONMENT    VMS and UNIX
.COMMENTS 
holds auxiliary routines Fextcor, Fcnvsng, Fextfr, Finsfr
.VERSION  [1.00] 890929: creation
 
 011121		last modif

----------------------------------------------------------------------------*/

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

#include <fileexts.h>

#define DISK_REC    512


static int   naxis, npix[3];
static int   newnax, newpix[3];

static double     newstart[3];
static double     start[3], step[3];

static struct FCT_STRUCT  *fctpntr;

/*

*/
 
int Fextcor(fno,string,sta,end,size)
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE
from string [xa,ya,za:xe,ye,ze] or @no  extract the start + end frame pixels
.ALGORITHM
use Fcnvcor + CGN_CNVT to do the job
.RETURNS
return status ( 0 = o.k )
-----------------------------------------------------------------------*/

int fno;     	/* IN: input file id  */
char	*string;     	/* IN: string with coords: [xa,ya,za:xe,ye,ze]  */
int   *sta;	/* OUT: start frame pixels  */
int   *end;	/* OUT: end frame pixels  */
int   *size;	/* OUT: size of subframe  */


{
int    nout, off, n1, k, iav;
int    knul = 0, kuni = 0;
int    ibuf;

float  rbuf;

double dbuf;

char   work[80], test;


for (k=0; k<3; k++)
   {
   npix[k] = 1;
   sta[k] = 1;
   end[k] = 1;
   }

/*  obtain necessary info from input frame */
 
(void) SCDRDI(fno,"NAXIS",1,1,&iav,&naxis,&kuni,&knul);
(void) SCDRDI(fno,"NPIX",1,naxis,&iav,npix,&kuni,&knul);
(void) SCDRDD(fno,"START",1,naxis,&iav,start,&kuni,&knul);
(void) SCDRDD(fno,"STEP",1,naxis,&iav,step,&kuni,&knul);
   

nout = 1;				/*  determine true NAXIS  */
for (n1=naxis-1; n1>=0; n1--)
   {
   if (npix[n1] > 1)
      {
      nout = n1 + 1;
      break;
      }
   }
newnax = nout;


if (string[0] == '@')		/* handle [...:...] and @... differently  */
   {
   if (newnax == 3)
      {
      end[0] = npix[0];
      end[1] = npix[1];
      }
   else if (newnax == 2)
      end[0] = npix[0];
   else
      goto inv_syntax;

   k = CGN_CNVT(&string[1],1,1,&ibuf,&rbuf,&dbuf);
   if (k != 1) goto inv_syntax;

   sta[newnax-1] = ibuf;
   end[newnax-1] = ibuf;
   }

else
   {
   (void) strcpy(work,string);		/* copy string, could be a constant */
   nout = 0;  off = 1;  n1 = naxis - 1;

st_loop:		/*  loop over start coords  */
   if (nout == n1)
      test = ':';
   else
      test= ',';

   k = CGN_INDEXC(&work[off],test);
   if (k <= 0) goto inv_syntax;

   work[off+k] = '\0';
   sta[nout] = Fcnvsng(&work[off],nout);
   if (sta[nout] < 1) goto inv_syntax;

   off += (k+1);
   if (++nout < naxis) goto st_loop;

   nout = 0;

end_loop:			/* loop over end coords  */
   if (nout == n1)
      test = ']';
   else
      test= ',';

   k = CGN_INDEXC(&work[off],test);
   if (k <= 0) goto inv_syntax;

   work[off+k] = '\0';
   end[nout] = Fcnvsng(&work[off],nout);
   if (end[nout] < 1) goto inv_syntax;

   if (++nout < naxis) 
      {
      off += (k+1);
      goto end_loop;
      }
   }

*size = 1;
newnax = 1;
for (k=0; k<3; k++)
   {
   newpix[k] = end[k] - sta[k] + 1;
   if (newpix[k] > 1) newnax = k + 1;
   *size *= newpix[k];
   }

return ERR_NORMAL;

inv_syntax:
return ERR_INPINV;
}

/*

*/

int Fcnvsng(stri,no)
char    *stri;
int     no;

{
int    m, nn;

float  rbuf, rm;

double dbuf;


if (*stri == '<')
   m = 1;

else if (*stri == '>')
   m = npix[no];

else if (*stri == '@')
   {
   nn = CGN_CNVT(&stri[1],1,1,&m,&rbuf,&dbuf);
   if (nn != 1) m = -1;
   }

else
   {
   nn = CGN_CNVT(stri,2,1,&m,&rbuf,&dbuf);
   if (nn != 1)
      m = -1;
   else
      {
      rm = (float) ((rbuf - start[no]) / step[no]);
      m = (int) rm;
      m++;
      }
   }

return (m);
}

/*

*/
 
int Fextfr(fno,sta,end,tno)
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE
extract subframe (tno) from father (fno)
.ALGORITHM
use conv_dat or rddisk to get stuff in
then write compact subframe out via wrdisk
.RETURNS
return status ( 0 = o.k )
-----------------------------------------------------------------------*/
int fno;        /* IN: father file id  */
int *sta;	/* IN: start coords */
int *end;	/* IN: end coords */
int tno;	/* IN: son file id */

{
register int nr;
int   ysize, felemz;
int   xsize, felemi, felemo, status;
int   savlink, dattyp, forma, nbyt, kpix, totsize, actsz;
int   kuni = 0;

char   *wpntra, *wp[2];


/* update descriptors for subframe  */

totsize = 1;
for (nr=0; nr<naxis; nr++)
   {
   totsize *= newpix[nr];
   newstart[nr] = start[nr] + ((sta[nr]-1) * step[nr]);
   }

/* and write descriptors  */

fctpntr = FCT.ENTRIES +  tno;		/* son: clear LINK first, */
savlink = fctpntr->LINK[0];		/* so we can */
fctpntr->LINK[0] = 0;			/* write into the extracted frame */
(void) SCDWRI(tno,"NAXIS",&newnax,1,1,&kuni);
(void) SCDWRI(tno,"NPIX",newpix,1,naxis,&kuni);
(void) SCDWRD(tno,"START",newstart,1,naxis,&kuni);
(void) SCDWRD(tno,"STEP",step,1,naxis,&kuni);
(void) SCDWRI(tno,"__NAXIS",&naxis,1,1,&kuni);
(void) SCDWRI(tno,"__NPIX",npix,1,3,&kuni);
(void) SCDWRI(tno,"__STPIX",sta,1,3,&kuni);
(void) SCDWRI(tno,"__ENPIX",end,1,3,&kuni);
fctpntr->LINK[0] = savlink;			/* reset LINK */
if (savlink == 3) return ERR_NORMAL;		/* no data extraction */


felemz = sta[0] + (sta[1]-1)*npix[0] + (sta[2]-1)*npix[0]*npix[1];
felemi = felemz;
felemo = 1;
xsize = newpix[0];
ysize = newpix[1];

fctpntr = FCT.ENTRIES +  fno;		/* use DATTYP of father */
dattyp = fctpntr->DATTYP;
nbyt = get_byte(dattyp);
forma = fctpntr->FORMAT;

fctpntr = FCT.ENTRIES +  tno;
fctpntr->DATTYP = fctpntr->FORMAT = dattyp;	/* that's equal for the son */
fctpntr->NOBYTE = nbyt;
fctpntr->PIXPBL =  DISK_REC / fctpntr->NOBYTE;
fctpntr = FCT.ENTRIES +  fno;


kpix = xsize * ysize;		/* allocate virtual memory for file copy  */
wpntra = malloc((unsigned int)(kpix*nbyt));
if (wpntra == (char *) 0) return ERR_INPINV;


if (fctpntr->FILTYP > 0)                /* FITS file */
   {
  zf_loop:					/* reading loop  for planes */
   wp[0] = wpntra;
   for(nr=0; nr<ysize; nr++)
      {				/* must pass via SCFGET for FITS files ... */
      status = SCFGET(fno,felemi,xsize,&actsz,wp[0]);
      if (status != ERR_NORMAL) return status;

      wp[0] += xsize*nbyt;
      felemi += npix[0];				/* offset in father */
      }

   status = wrdisk(tno,felemo,kpix,wpntra);     /* write plane to disk */
   if (status != ERR_NORMAL) return status;	/* because it's bdf */
   felemo += kpix;

   if (felemo <= totsize)
      {
      felemz += (npix[0] * npix[1]);
      felemi = felemz;
      goto zf_loop;
      }
   }
else
   {
  z_loop:					/* reading loop  for planes */
   wp[0] = wpntra;
   for(nr=0; nr<ysize; nr++)
      {
      if (dattyp != forma)
         status = conv_dat(3,fno,felemi,xsize,&actsz,wp);
      else
         status = rddisk(fno,felemi,xsize,&actsz,wp[0]);
      if (status != ERR_NORMAL) return status;

      wp[0] += xsize*nbyt;
      felemi += npix[0];				/* offset in father */
      }

   status = wrdisk(tno,felemo,kpix,wpntra);	/* write plane to disk */
   if (status != ERR_NORMAL) return status;
   felemo += kpix;

   if (felemo <= totsize)
      {
      felemz += (npix[0] * npix[1]);
      felemi = felemz;
      goto z_loop;
      }
   }


free(wpntra);
return ERR_NORMAL;
}

/*

*/
 
int Finsfr(fno,tno)
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE
reinsert subframe (tno) into father (fno)
.ALGORITHM
get compact subframe in via rddisk
then write it out to father via conv_dat or wrdisk
.RETURNS
return status ( 0 = o.k )
-----------------------------------------------------------------------*/
int fno;        /* IN: father file id  */
int tno;        /* IN: son file id */

{
register int  nr;
int   xsize, felemi, felemo, status;
int   sta[3], end[3];
int   savlink, iav, kpix, nbyt, totsize, actsz;
int   felemz, ysize, dattyp, forma;
int   knul = 0, kuni = 0;

char  *wpntra, *wp[2];


/* get descriptors with values of father  */

fctpntr = FCT.ENTRIES +  tno;		/* son: clear LINK first, */
savlink = fctpntr->LINK[0];		/* so we can read */
fctpntr->LINK[0] = 0;			/* descr.s from extracted frame */
(void) SCDRDI(tno,"__NPIX",1,3,&iav,npix,&kuni,&knul);
(void) SCDRDI(tno,"__STPIX",1,3,&iav,sta,&kuni,&knul);
(void) SCDRDI(tno,"__ENPIX",1,3,&iav,end,&kuni,&knul);
fctpntr->LINK[0] = savlink;			/* reset LINK */

felemz = sta[0] + (sta[1]-1)*npix[0] + (sta[2]-1)*npix[0]*npix[1];
felemo = felemz;
felemi = 1;
xsize = end[0] - sta[0] + 1;
ysize = end[1] - sta[1] + 1;
totsize = xsize * ysize * (end[2] - sta[2] + 1);


fctpntr = FCT.ENTRIES + fno;	/* allocate virtual memory for file copy  */
dattyp = fctpntr->DATTYP;		/* same for father + son */
nbyt = get_byte(dattyp);
forma = fctpntr->FORMAT;

kpix = xsize * ysize;			/* pull in a full squeezed plane */
wpntra = malloc((unsigned int)(kpix*nbyt));
if (wpntra == (char *) 0) return ERR_INPINV;


if (fctpntr->FILTYP > 0)                /* FITS file */
   {
  zf_loop:                                      /* reading loop  for planes */
   wp[0] = wpntra;
   status = rddisk(tno,felemi,kpix,&actsz,wpntra);
   if (status != ERR_NORMAL) return status;

   for(nr=0; nr<ysize; nr++)
      {                         /* must pass via SCFGET for FITS files ... */
      status = SCFPUT(fno,felemo,xsize,wp[0]);
      if (status != ERR_NORMAL) return status;

      wp[0] += xsize*nbyt;
      felemo += npix[0];                                /* offset in father */
      }

   felemi += kpix;
   if (felemi <= totsize)
      {
      felemz += (npix[0] * npix[1]);
      felemo = felemz;
      goto zf_loop;
      }
   }
else
   {
  z_loop:                                      /* reading loop for planes */
   wp[0] = wpntra;
   status = rddisk(tno,felemi,kpix,&actsz,wpntra);
   if (status != ERR_NORMAL) return status;

   for (nr=0; nr<ysize; nr++)
      {
      if (dattyp != forma)
         status = conv_dat(2,fno,felemo,xsize,&actsz,wp);
      else
         status = wrdisk(fno,felemo,xsize,wp[0]);
      if (status != ERR_NORMAL) return status;

      wp[0] += xsize*nbyt;
      felemo += npix[0];
      }

   felemi += kpix;
   if (felemi <= totsize) 
      {
      felemz += (npix[0] * npix[1]);
      felemo = felemz;
      goto z_loop;		/* get next plane */
      }
   }


free(wpntra);
return ERR_NORMAL;
}