/*
 				field.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	Part of:	SWarp
*
*	Author:		E.BERTIN (IAP)
*
*	Contents:	Handling of field structures.
*
*	Last modify:	11/03/2001
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

#include <math.h>
#include <pwd.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <sys/types.h>

#include "define.h"
#include "globals.h"
#include "fitscat.h"
#include "fitswcs.h"
#include "back.h"
#include "coadd.h"
#include "data.h"
#include "field.h"
#include "key.h"
#include "prefs.h"
#include "wcs/wcs.h"

/****** load_field ************************************************************
PROTO	fieldstruct *load_field(char *filename, int frameno, int fieldno)
PURPOSE	Initialize a field structure (in read mode)
INPUT	File name,
	FITS extension number in file (0=primary)
	Field number in coaddition (for various config settings)
	Field flags,
OUTPUT	The new field pointer if OK, NULL otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	21/05/2000
 ***/
fieldstruct	*load_field(char *filename, int frameno, int fieldno)

  {
   fieldstruct	*field;
   tabstruct	*tab;
   char		*pstr;
   int		i;

/* First allocate memory for the new field (and nullify pointers) */
  QCALLOC(field,fieldstruct, 1);
  field->flags = FIELD_READ;
/* Set conversion flags */
  field->cflags = (prefs.interp_flag[fieldno]?CONVERT_INTERP:0)
                |(prefs.subback_flag[fieldno]?CONVERT_BACKSUB:0)
                |(prefs.dyncomp_flag[fieldno]?CONVERT_DYNCOMP:0);
  field->frameno = frameno;
  strcpy (field->filename, filename);
/* A short, "relative" version of the filename */
  if (!(field->rfilename = strrchr(field->filename, '/')))
    field->rfilename = field->filename;
  else
    field->rfilename++;

/* Create a file name with a "header" extension */
  strcpy(field->hfilename, filename);
  if (!(pstr = strrchr(field->hfilename, '.')))
    pstr = field->hfilename+strlen(field->hfilename);
  sprintf(pstr, "%s", prefs.head_suffix);

  sprintf(gstr, "Looking for %s", field->rfilename);
  NFPRINTF(OUTPUT, gstr);

/* Check the image exists and read important info (image size, etc...) */
  field->cat = read_cat(filename);
  if (!field->cat)
    error(EXIT_FAILURE, "*Internal Error*: Can't read ", filename);
  tab=field->cat->tab;

/* Select the desired FITS extension (if frameno != 0) */
  if (frameno >= field->cat->ntab)
    error(EXIT_FAILURE, "*Internal Error*: FITS extension unavailable in ",
		filename);
  for (i=frameno; i--;)
    tab = tab->nexttab;
  field->tab = tab;

/* Insert additional header informations from the "header" file */
  if (!read_aschead(field->hfilename, frameno, tab))
    warning("FITS header data read in ", field->hfilename);

  if (tab->naxis<1)
    error(EXIT_FAILURE, "*Error*: Zero-dimensional table in ", filename);

/* Force data to be at least 2D */
  if (tab->naxis<2)
    {
    tab->naxis = 2;
    QREALLOC(tab->naxisn, int, 2);
    tab->naxisn[1] = 1;
    }

/* Read WCS information in FITS header */
  field->wcs = read_wcs(tab);

/* Set field width and field height (the latter can be "virtual") */
  field->width = tab->naxisn[0];
  field->height = 1;
  for (i=1; i<tab->naxis; i++)
    field->height *= tab->naxisn[i];
  field->npix = field->width*field->height;

/*-- Background */
  field->backw = prefs.back_size[fieldno]<field->width ?
					prefs.back_size[fieldno]
					: field->width;
  field->backh = prefs.back_size[fieldno]<field->height ?
					prefs.back_size[fieldno]
					: field->height;
  field->nbackp = field->backw * field->backh;
  if ((field->nbackx = (field->width-1)/field->backw + 1) < 1)
    field->nbackx = 1;
  if ((field->nbacky = (field->height-1)/field->backh + 1) < 1)
    field->nbacky = 1;
  field->nback = field->nbackx * field->nbacky;
  field->nbackfx = field->nbackx>1 ? prefs.back_fsize[fieldno] : 1;
  field->nbackfy = field->nbacky>1 ? prefs.back_fsize[fieldno] : 1;
/* Set the back_type flag if absolute background is selected */
  field->back_type = prefs.back_type[fieldno];
  field->backdefault = prefs.back_default[fieldno];

/* Set the flux scale */
  field->fscale = prefs.fscale_default[fieldno];
  fitsread(tab->headbuf,prefs.fscale_keyword,&field->fscale,H_FLOAT,T_DOUBLE);
  if (field->fscale == 0.0)
    {
    warning(field->filename, " has flux scale = 0: I will take 1 instead");
    field->fscale = 1.0;
    }
/* This scaling will now be applied to all data from this field */
  tab->bscale *= field->fscale;
/* Set the dynamic compression factor */
  field->gain = prefs.gain_default[fieldno];
  fitsread(tab->headbuf, prefs.gain_keyword, &field->gain, H_FLOAT, T_DOUBLE);
/* gain = 0 means gain = infinite and thus no dynamic compression */
  if (field->gain == 0.0)
    field->cflags &= ~CONVERT_DYNCOMP;
  else
    field->gain /= field->fscale;

  return field;
  }


/****** read_aschead ********************************************************
PROTO	int	read_aschead(char *filename, int frameno, tabstruct *tab)
PURPOSE	Read a ASCII header file and update the current field's tab
INPUT	Name of the ASCII file,
	Frame number (if extensions),
	Tab structure.
OUTPUT	RETURN_OK if the file was found, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	04/03/2001
 ***/
int	read_aschead(char *filename, int frameno, tabstruct *tab)
  {
   static char	keyword[MAXCHAR],data[MAXCHAR],comment[MAXCHAR];
   FILE		*file;
   h_type	htype;
   t_type	ttype;
   int		i, pos;

  if (file=fopen(filename, "r"))
    {
/*- Skip previous ENDs in multi-FITS extension headers */
    for (i=(frameno?(frameno-1):0); i--;)
      while (fgets(gstr, MAXCHAR, file) && strncmp(gstr,"END     ",8));
    pos = fitsfind(tab->headbuf, "END     ");
    memset(gstr, ' ', 80);
    while (fgets(gstr, MAXCHAR, file) && strncmp(gstr,"END     ",8))
      {
      fitspick(gstr, keyword, data, &htype, &ttype, comment);
/*---- Block critical keywords */
      if (!wstrncmp(keyword, "SIMPLE  ", 8)
	||!wstrncmp(keyword, "BITPIX  ", 8)
	||!wstrncmp(keyword, "NAXIS?? ", 8)
	||!wstrncmp(keyword, "BSCALE  ", 8)
	||!wstrncmp(keyword, "BZERO   ", 8))
        continue;
/*---- Always keep a one-line margin */
      if ((pos+1)*80>=tab->headnblock*FBSIZE)
        {
        tab->headnblock++;
        QREALLOC(tab->headbuf, char, tab->headnblock*FBSIZE);
        }
      if (fitsadd(tab->headbuf, keyword, comment)>=pos)
        pos++;
      fitswrite(tab->headbuf, keyword, data, htype, ttype);
      memset(gstr, ' ', 80);
      }
    fclose(file);
/*-- Update the tab data */
    readbasic_head(tab);
    return RETURN_OK;
    }
  else
    return RETURN_ERROR;
  }


/****** inherit_field ********************************************************
PROTO	fieldstruct *inherit_field(char *filename, fieldstruct *reffield,
				int fflags)
PURPOSE	Make a copy of a field structure.
INPUT	Reference field pointer.
	flags.
OUTPUT	The new field pointer if OK, NULL otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	07/07/2000
 ***/
fieldstruct	*inherit_field(char *filename, fieldstruct *reffield,
				int flags)
  {
   fieldstruct	*field;

/* First allocate memory for the new field (and nullify pointers) */
  QCALLOC(field, fieldstruct, 1);

/* Copy what is important and reset the remaining */
  *field = *reffield;
  field->cat = new_cat(1);
  inherit_cat(reffield->cat, field->cat);
  field->tab = field->cat->tab;
  field->tab->cat = field->cat;
  field->flags = flags;
  field->reffield = reffield;
/* We don't need this */
  field->back = NULL;
  field->dback = NULL;
  field->sigma = NULL;
  field->dsigma = NULL;
  field->pix = NULL;
  field->backline = NULL;
  field->wcs = NULL;
  field->rawmin = NULL;
  field->rawmax = NULL;
  field->reffield = NULL;

  strcpy(field->filename, filename);
  strcpy(field->cat->filename, filename);

  return field;
  }


/****** end_field ************************************************************
PROTO	void end_field(fieldstruct *field)
PURPOSE	Free a field structure.
INPUT	field structure pointer.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	16/04/2000
 ***/
void	end_field(fieldstruct *field)

  {
/* Check first that a tab structure is present */
  if (field->tab)
    {
/*-- Terminate astrometry */
    end_wcs(field->wcs);
/*-- End memory mapping (if allocated) */
    free_body(field->tab);

/*-- Close cat if not already done */
    if (field->cat)
      free_cat(field->cat, 1);
    field->tab = NULL;
    }

  end_back(field);
  field->pix = NULL;
  free(field);

  return;
  }


/****** printinfo_field ******************************************************
PROTO	void printinfo_field(fieldstruct *field)
PURPOSE	Print info about a field
INPUT	Pointer to the field.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	11/03/2001
 ***/
void	printinfo_field(fieldstruct *field, fieldstruct *wfield)

  {
   wcsstruct		*wcs;
   static double	pixpos[NAXIS],
			wcspos[NAXIS];
   double		pixscale;
   int			i;
  
/* Information about the file */
  if (field->frameno)
      sprintf(gstr, "Extension #%d:", field->frameno);
    else
      sprintf(gstr, "");
  QPRINTF(OUTPUT, "  %s   \"%.20s\"   %s   %dx%d   %d bits (%s)\n",
	gstr, *field->ident? field->ident: "no ident",
	wfield? "weighted" : "unweighted",	
	field->width, field->height, field->tab->bytepix*8,
	field->tab->bitpix>0?
		(field->tab->compress_type!=COMPRESS_NONE ?
			"compressed":"integers") : "floats");

/* Astrometry */
  wcs = field->wcs;

/* Find field center */
  for (i=0; i<wcs->naxis; i++)
    pixpos[i] = (wcs->naxisn[i]+1.0)/2.0;
  raw_to_wcs(wcs, pixpos, wcspos);
  if (wcs->lat != wcs->lng)
    {
    pixscale = wcs->wcsscale[wcs->lng]*DEG;
    QPRINTF(OUTPUT, "    Center: %s %s   %.3g'x%.3g'  Scale: %.4g ''/pixel\n",
	degtosexal(wcspos[wcs->lng]),
	degtosexde(wcspos[wcs->lat]),
	wcs->naxisn[wcs->lng]*pixscale/ARCMIN,
	wcs->naxisn[wcs->lat]*pixscale/ARCMIN,
	pixscale/ARCSEC);
    }
  else if (wcs->naxis >= 2)
    {
    QPRINTF(OUTPUT,
	"    Center: %.3g,%.3g   %.3gx%.3g  Scale: %.4gx%.4g /pixel\n",
	wcspos[0],
	wcspos[1],
	wcs->naxisn[0]*wcs->wcsscale[0],
	wcs->naxisn[1]*wcs->wcsscale[1],
	wcs->wcsscale[0],
	wcs->wcsscale[1]);
    }
  else
    QPRINTF(OUTPUT, "    Center: %.3g   %.3g  Scale: %.3g /pixel\n",
	wcspos[0],
	wcs->naxisn[0]*wcs->wcsscale[0],
	wcs->wcsscale[0]);

    for (i=0; i<wcs->naxis; i++)
    if (i==wcs->lat || i==wcs->lng)
      {
      pixscale = wcs->wcsscale[i]*DEG;
      break;
      }

/* Photometry */
  QPRINTF(OUTPUT, "    Gain: %.4g e-/ADU   Flux scaling: %.4g X\n",
	field->gain,
	field->fscale);

  return;
  }


/******* init_field ********************************************************
PROTO	fieldstruct *init_field(fieldstruct **infield, fieldstruct **inwfield,
			int ninput, char *filename)

PURPOSE	Automatically set appropriate output field parameters according to the
	prefs and a set of input fields.
INPUT	Input field ptr array,
	Input weight field ptr array,
	number of input fields,
	Filename.
OUTPUT	Pointer to the new output field.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 11/03/2001
 ***/
fieldstruct *init_field(fieldstruct **infield, fieldstruct **inwfield,
			int ninput, char *filename)
  {
   extern char		celsysname[][2][8];
   static float		scale[MAXINFIELD];

   fieldstruct		*field;
   tabstruct		*tab;
   wcsstruct		*wcs;
   double		wcsmin,wcsmax,wcsmin1,wcsmax1,wcsmin2,wcsmax2,
			pixscale;
   char			*pstr;
   int			i,j,n, naxis, lat,lng, countmin0,countmax0,
			countmin,countmax;

/* First allocate memory for the new field */
  QCALLOC(field,fieldstruct, 1);
  field->flags = FIELD_WRITE;

  strcpy (field->filename, filename);
/* A short, "relative" version of the filename */
  if (!(field->rfilename = strrchr(field->filename, '/')))
    field->rfilename = field->filename;
  else
    field->rfilename++;
/* Create a file name with a "header" extension */
  strcpy(field->hfilename, filename);
  if (!(pstr = strrchr(field->hfilename, '.')))
    pstr = field->hfilename+strlen(field->hfilename);
  sprintf(pstr, "%s", prefs.head_suffix);

  field->cat = new_cat(1);
  init_cat(field->cat);
  strcpy(field->cat->filename, field->filename);
  field->tab = tab = field->cat->tab;
  tab->cat = field->cat;

/*-------------------------- Set the astrometry -----------------------------*/
/* Check that angular axes are the same and adopt them for output */
  lat = lng = -1;
  for (i=0; i<ninput; i++)
    {
    if (lng == -1)
      lng = infield[i]->wcs->lng;
    if (lat == -1)
      lat = infield[i]->wcs->lat;
    if (infield[i]->wcs->lng != lng || infield[i]->wcs->lat != lat)
      error(EXIT_FAILURE, "*Error*: Mismatched axis type in ",
				infield[i]->filename); 
    if (!i)
      naxis = infield[i]->wcs->naxis;
    else if (infield[i]->wcs->naxis != naxis)
      error(EXIT_FAILURE, "*Error*: Mismatched number of axes in ",
				infield[i]->filename); 
    }

  QMALLOC(tab->naxisn, int, naxis);

/* Fill a new WCS structure */
  QCALLOC(wcs, wcsstruct, 1);
  field->wcs = wcs;
  wcs->naxis = tab->naxis = naxis;
  QCALLOC(wcs->projp, double, naxis*100);

  wcs->lng = lng;
  wcs->lat = lat;
/* Copy the types and units */
  for (i=0; i<naxis; i++)
    {
    strncpy(wcs->cunit[i], infield[0]->wcs->cunit[i], 8);
    strncpy(wcs->ctype[i], infield[0]->wcs->ctype[i], 8);
    }

/* Change the Celestial system if needed */
  if (lng!=lat && prefs.celsys_type != CELSYS_NATIVE)
    {
    n = prefs.celsys_type - 1;	/* Drop the ``NATIVE'' option */
    strncpy(wcs->ctype[lng], celsysname[n][0], 4);
    strncpy(wcs->ctype[lat], celsysname[n][1], 4);
    }

  for (n=0; n<naxis; n++)
    { 
    switch (prefs.center_type[n])
      {
      case CENTER_MOST:
/*------ Find (in a dumb way) the densest part of the coadd */
        countmin0 = countmax0 = 0;
        for (j=0; j<ninput; j++)
          {
          wcsmin1 = infield[j]->wcs->wcsmin[n];
          wcsmax1 = infield[j]->wcs->wcsmax[n];
          countmin = countmax = 0;
          for (i=0; i<ninput; i++)
            {
            wcsmin2 = infield[i]->wcs->wcsmin[n];
            wcsmax2 = infield[i]->wcs->wcsmax[n];
/*---------- Test for the lower limit */
            if (wcsmin2<=wcsmin1)
              countmin++;
            if (wcsmax2>wcsmin)
              countmin++;
/*---------- Test for the upper limit */
            if (wcsmax2>=wcsmax1)
              countmax++;
            if (wcsmin2<wcsmax1)
              countmax++;
            }

/*-------- Update the "most intersecting" limits */
          if (countmin>countmin0)
            {
            wcsmin = wcsmin1;
            countmin0 = countmin;
            }
          if (countmax>countmax0)
            {
            wcsmax = wcsmax1;
            countmax0 = countmax;
            }
          }
        wcs->crval[n] = (wcsmin+wcsmax)/2.0;
        break;

      case CENTER_ALL:
        wcsmin = BIG;
        wcsmax = -BIG;
        for (i=0; i<ninput; i++)
          {
          wcsmin2 = infield[i]->wcs->wcsmin[n];
          wcsmax2 = infield[i]->wcs->wcsmax[n];
/*------- Test for the lower limit */
          if (wcsmin2<wcsmin)
            wcsmin = wcsmin2;
/*------- Test for the upper limit */
          if (wcsmax2>wcsmax)
            wcsmax = wcsmax2;
          }
        wcs->crval[n] = (wcsmin+wcsmax)/2.0;
        break;

      case CENTER_MANUAL:
        pstr = prefs.image_center[n];
	wcs->crval[n] = strchr(pstr, ':') ?
			  (n==lng?sextodegal(pstr):sextodegde(pstr))
			: atof(pstr);
        break;

      default:
        error(EXIT_FAILURE,
	      "*Internal Error*: Unknown area type in ", "init_field()");
        break;
      }

    switch (prefs.pixscale_type[n])
      {
      case PIXSCALE_MIN:
        pixscale = BIG;
        for (i=0; i<ninput; i++)
          if (infield[i]->wcs->wcsscale[n] < pixscale)
            pixscale = infield[i]->wcs->wcsscale[n];
        break;
      case PIXSCALE_MAX:
        pixscale = -BIG;
        for (i=0; i<ninput; i++)
          if (infield[i]->wcs->wcsscale[n] > pixscale)
            pixscale = infield[i]->wcs->wcsscale[n];
        break;
      case PIXSCALE_MEDIAN:
        for (i=0; i<ninput; i++)
          scale[i] = (float)infield[i]->wcs->wcsscale[n];
        pixscale = (double)hmedian(scale, ninput);
        break;
      case PIXSCALE_FIT:
        if (prefs.center_type[n] == CENTER_MANUAL)
          error(EXIT_FAILURE,
	      "*Error*: cannot choose an output pixel scale if ",
		"CENTER_TYPE is in MANUAL mode");
        if (!prefs.image_size[n])
          error(EXIT_FAILURE,
	      "*Error*: cannot choose an output pixel scale if ",
		"IMAGE_SIZE is not provided");

        pixscale = (wcsmax - wcsmin) / (double)prefs.image_size[n];
        break;

      case PIXSCALE_MANUAL:
        pixscale = prefs.pixscale[n];
        if ((n==lng || n==lat) && lat!=lng)
          pixscale *= (ARCSEC/DEG);
        break;

      default:
        error(EXIT_FAILURE,
	      "*Internal Error*: Unknown pixel scale type in ","init_field()");
        break;
      }

    if (!pixscale)
      error(EXIT_FAILURE, "*Error*: null output pixel size ", "");

    wcs->cd[n*(naxis+1)] = wcs->cdelt[n] = pixscale;

/*-- Compute image size. It is necessarily > 0 */
    if (prefs.image_size[n])
      tab->naxisn[n] = wcs->naxisn[n] = prefs.image_size[n];
    else
      {
/*---- Add a 5% margin in field size */
      if ((n==lng || n==lat) && lat!=lng)
        tab->naxisn[n] = wcs->naxisn[n]
		= (int)(fmod(wcsmax-wcsmin+360.0, 360.0)*1.05/pixscale) + 1;
      else
        tab->naxisn[n] = wcs->naxisn[n]
		= (int)((wcsmax-wcsmin)*1.05/pixscale) + 1;
      }

    wcs->crpix[n] = (wcs->naxisn[n]+1.0)/2.0;
    }

/* The special case of longitude */

  if (lat!=lng)
    {
    if (strcmp(prefs.projection_name, "NONE"))
      {
      strcpy(wcs->ctype[lng]+5, prefs.projection_name);
      strcpy(wcs->ctype[lat]+5, prefs.projection_name);
      }
    else
      {
      strcpy(wcs->ctype[lng], "");
      strcpy(wcs->ctype[lat], "");
      }
    if (!prefs.image_size[lng])
      {
      tab->naxisn[lng] = wcs->naxisn[lng]
		= (int)((wcs->naxisn[lng]-1)*cos(wcs->crval[lat]*DEG))+1;
      wcs->crpix[lng] = (wcs->naxisn[lng]+1.0)/2.0;
      }
    else if (prefs.pixscale_type[lng] == PIXSCALE_FIT)
      wcs->cdelt[lng] *= cos(wcs->crval[lat]*DEG);
/*-- Make pixel scales equal in alpha and delta */
    if (wcs->cdelt[lng] < wcs->cdelt[lat])
      wcs->cdelt[lng] = wcs->cdelt[lat];
    else
      wcs->cdelt[lat] = wcs->cdelt[lng];
/*-- No negative longitude!! */
    wcs->crval[lng] = fmod(wcs->crval[lng]+360.0, 360.0);
    }

/* Default equinox, RA-DEC-sys, longpole and latpole */
  wcs->equinox = 2000.0;
  wcs->radecsys = RDSYS_FK5;
  wcs->longpole = 999.0;
  wcs->latpole = 999.0;

  write_wcs(tab, wcs);

/* Insert additional header informations from the "header" file */
  if (read_aschead(field->hfilename, 0, tab))
    {
    QCALLOC(wcs->wcsprm, struct wcsprm, 1);

/*-- Test if the WCS is recognized and a celestial pair is found */
    wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);

/*-- Initialize other WCS structures */
    init_wcs(wcs);
/*-- Find the range of coordinates */
    range_wcs(wcs);
    }
  else
    {
    warning("FITS header data read in ", field->hfilename);
/*-- Drop the current WCS structure and update from what's in the header */
    end_wcs(wcs);
    field->wcs = wcs = read_wcs(tab);
    }

  field->width = tab->naxisn[0];
  field->height = 1;
  for (i=1; i<tab->naxis; i++)
    field->height *= tab->naxisn[i];
  field->npix = field->width*field->height;

/* Default flux scale and gain*/
  field->fscale = 1.0;
  field->gain = 0.0;

/* Add relevant information to FITS header */
  writefitsinfo_field(field);

  return field;
  }


/******* scale_field *********************************************************
PROTO	void scale_field(fieldstruct *field, fieldstruct *reffield)
PURPOSE	Compute the flux-scaling factor for each input field.
INPUT	Field ptr,
	Reference field ptr.
OUTPUT	-.
NOTES   -.
AUTHOR  E. Bertin (IAP)
VERSION 07/03/2001
 ***/
void	scale_field(fieldstruct *field, fieldstruct *reffield)
  {
   wcsstruct		*wcs;
   static double	raw[NAXIS], wcspos2[NAXIS];
   double		*wcspos,*wcsscale,
			inscale,outscale;
   int			i, naxis, lng,lat;

  wcs = reffield->wcs;
  naxis = wcs->naxis;
  lng = wcs->lng;
  lat = wcs->lat;
/* ``Go'' to the position where scale has been computed on input image */
  wcspos = field->wcs->wcsscalepos;
  wcs_to_raw(wcs, wcspos, raw);
  wcsscale = field->wcs->wcsscale;
/* Compute scaling factors for input and output images */
  inscale = outscale = 1.0;
  for (i=0; i<naxis; i++)
    {
    if ((i==lng || i==lat) && lng!=lat)
      outscale *= sqrt(wcs_scale(wcs, raw));
    else
      {
      raw[i] += 1.0;
      raw_to_wcs(wcs, raw, wcspos2);
      outscale *= fabs(wcspos2[i] - wcspos[i]);
      raw[i] -= 1.0;
      }
    inscale *= wcsscale[i];
    }

  if (inscale!=0.0 && outscale!=0.0)
    {
    field->fscale *= outscale/inscale;
    field->tab->bscale *= outscale/inscale;
    }

  return;
  }


/****** writefitsinfo_field **************************************************
PROTO	void writefitsinfo_field(fieldstruct *field)
PURPOSE	Add relevant information to a field FITS header.
INPUT	Pointer to the field.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	04/09/2000
 ***/
void	writefitsinfo_field(fieldstruct *field)

  {
   struct passwd	*pwd;
   struct tm		*tm;
   time_t		thetime;
   tabstruct		*tab;
   extern pkeystruct	key[];
   extern char		keylist[][16];
   char			*pstr;
   int			d, n;

  tab = field->tab;
/* Keep some margin */
  QREALLOC(tab->headbuf, char,
	80*(tab->headnblock*36+15+tab->naxis*3+prefs.ninfield*9));

  fitsadd(tab->headbuf, "COMMENT ", "");
  fitsadd(tab->headbuf, "SOFTNAME", "The software that processed those data");
  fitswrite(tab->headbuf, "SOFTNAME", BANNER, H_STRING, T_STRING);
  fitsadd(tab->headbuf, "SOFTVERS", "Version of the software");
  fitswrite(tab->headbuf, "SOFTVERS", VERSION, H_STRING, T_STRING);
  fitsadd(tab->headbuf, "SOFTDATE", "Release date of the software");
  fitswrite(tab->headbuf, "SOFTDATE", DATE, H_STRING, T_STRING);
  fitsadd(tab->headbuf, "SOFTAUTH", "Maintainer of the software");
  fitswrite(tab->headbuf, "SOFTAUTH", COPYRIGHT, H_STRING, T_STRING);
  fitsadd(tab->headbuf, "SOFTINST", "Institute");
  fitswrite(tab->headbuf, "SOFTINST", INSTITUTE, H_STRING, T_STRING);
  fitsadd(tab->headbuf, "COMMENT  ", "");

/* User name */
  if (pwd=getpwuid(geteuid()))
    {
    fitsadd(tab->headbuf, "AUTHOR", "Who ran the software");
    fitswrite(tab->headbuf, "AUTHOR", strtok(pwd->pw_gecos,","),
		H_STRING, T_STRING);
    }

/* Host name */
  if (!gethostname(gstr, 80))
    {
    fitsadd(tab->headbuf, "ORIGIN", "Where it was done");
    fitswrite(tab->headbuf, "ORIGIN", gstr, H_STRING, T_STRING);
    }

/* Obs dates */
  thetime = time(NULL);
  tm = gmtime(&thetime);
  sprintf(gstr,"%04d-%02d-%02dT%02d:%02d:%02d",
                        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday,
                        tm->tm_hour, tm->tm_min, tm->tm_sec);
  fitsadd(tab->headbuf, "DATE    ", "When it was started (GMT)");
  fitswrite(tab->headbuf, "DATE    ", gstr, H_STRING, T_STRING);

/* Config parameters */
/* COMBINE_TYPE */
  fitsadd(tab->headbuf,"COMBINET","COMBINE_TYPE config parameter for " BANNER);
  fitswrite(tab->headbuf, "COMBINET", key[findkeys("COMBINE_TYPE", keylist,
	FIND_STRICT)].keylist[prefs.coadd_type], H_STRING, T_STRING);

/* Axis-dependent config parameters */
  fitsadd(tab->headbuf, "COMMENT  ", "");
  fitsadd(tab->headbuf, "COMMENT  ", "Axis-dependent config parameters");
  for (d=0; d<tab->naxis; d++)
    {
/*-- RESAMPLING_TYPE */
    sprintf(gstr, "RESAMPT%1d", d+1);
    fitsadd(tab->headbuf, gstr,
		"RESAMPLING_TYPE config parameter");
    fitswrite(tab->headbuf, gstr, key[findkeys("RESAMPLING_TYPE", keylist,
	FIND_STRICT)].keylist[prefs.resamp_type[d]], H_STRING, T_STRING);
/*-- CENTER_TYPE */
    sprintf(gstr, "CENTERT%1d", d+1);
    fitsadd(tab->headbuf, gstr, "CENTER_TYPE config parameter");
    fitswrite(tab->headbuf, gstr, key[findkeys("CENTER_TYPE", keylist,
	FIND_STRICT)].keylist[prefs.center_type[d]], H_STRING, T_STRING);
/*-- PIXELSCALE_TYPE */
    sprintf(gstr, "PSCALET%1d", d+1);
    fitsadd(tab->headbuf, gstr, "PIXELSCALE_TYPE config parameter");
    fitswrite(tab->headbuf, gstr, key[findkeys("PIXELSCALE_TYPE", keylist,
	FIND_STRICT)].keylist[prefs.pixscale_type[d]], H_STRING, T_STRING);
    }

/* Image-dependent config parameters */
  if (prefs.writefileinfo_flag)
    {
    fitsadd(tab->headbuf, "COMMENT  ", "");
    fitsadd(tab->headbuf, "COMMENT  ", "Image-dependent config parameters");
    for (n=0; n<prefs.ninfield; n++)
      {
/*---- Input images */
      sprintf(gstr, "FILE%04d", n+1);
      fitsadd(tab->headbuf, gstr, "Input filename");
/*---- A short, "relative" version of the filename */
      if (!(pstr = strrchr(prefs.infield_name[n], '/')))
        pstr = prefs.infield_name[n];
      else
        pstr++;
      fitswrite(tab->headbuf, gstr, pstr, H_STRING, T_STRING);
/*---- Input weights */
      if (n<prefs.ninwfield)
        {
        sprintf(gstr, "WGHT%04d", n+1);
        fitsadd(tab->headbuf, gstr, "Input weight-map");
/*------ A short, "relative" version of the filename */
        if (!(pstr = strrchr(prefs.inwfield_name[n], '/')))
          pstr = prefs.inwfield_name[n];
        else
          pstr++;
        fitswrite(tab->headbuf, gstr, pstr, H_STRING, T_STRING);
/*------ WEIGHT_TYPE */
        sprintf(gstr, "WGHT%04d", n+1);
        fitsadd(tab->headbuf, gstr,
		"WEIGHT_TYPE config parameter");
        fitswrite(tab->headbuf, gstr, key[findkeys("WEIGHT_TYPE", keylist,
	  FIND_STRICT)].keylist[prefs.weight_type[n]], H_STRING, T_STRING);
        }
/*---- Interpolation flag */
      sprintf(gstr, "INTF%04d", n+1);
      fitsadd(tab->headbuf, gstr, "INTERPOLATE config flag");
      fitswrite(tab->headbuf, gstr, &prefs.interp_flag[n], H_BOOL, T_STRING);
/*---- Background-subtraction flag */
      sprintf(gstr, "SUBF%04d", n+1);
      fitsadd(tab->headbuf, gstr, "SUBTRACT_BACK config flag");
      fitswrite(tab->headbuf, gstr, &prefs.subback_flag[n], H_BOOL, T_STRING);
/*---- BACK_TYPE */
      sprintf(gstr, "BCKT%04d", n+1);
      fitsadd(tab->headbuf, gstr, "BACK_TYPE config parameter");
      fitswrite(tab->headbuf, gstr, key[findkeys("BACK_TYPE", keylist,
        FIND_STRICT)].keylist[prefs.back_type[n]], H_STRING, T_STRING);
/*---- BACK_DEFAULT */
      if (prefs.back_type[n]==BACK_ABSOLUTE)
        {
        sprintf(gstr, "BCKD%04d", n+1);
        fitsadd(tab->headbuf, gstr,"BACK_DEFAULT config parameter");
        fitswrite(tab->headbuf, gstr, &prefs.back_default[n], H_EXPO,T_DOUBLE);
        }
/*---- BACK_SIZE */
      sprintf(gstr, "BCKS%04d", n+1);
      fitsadd(tab->headbuf, gstr,"BACK_SIZE config parameter");
      fitswrite(tab->headbuf, gstr, &prefs.back_size[n], H_INT, T_LONG);
/*---- BACK_FILTERSIZE */
      sprintf(gstr, "BKFS%04d", n+1);
      fitsadd(tab->headbuf, gstr,"BACK_FILTERSIZE config parameter");
      fitswrite(tab->headbuf, gstr, &prefs.back_fsize[n], H_INT, T_LONG);
      }
    }
/* Keep only what's necessary */
  tab->headnblock = ((fitsfind(tab->headbuf, "END     ")+36)*80)/FBSIZE;
  QREALLOC(tab->headbuf, char, tab->headnblock*FBSIZE);

  return;
  }