/*
                                  makeit.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*       Part of:        Swarp
*
*       Author:         E.BERTIN (IAP)
*
*       Contents:       Main loop
*
*       Last modify:    04/03/2001
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

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

#include "define.h"
#include "globals.h"
#include "fitscat.h"
#include "fitswcs.h"
#include "back.h"
#include "data.h"
#include "prefs.h"
#include "field.h"

/********************************** makeit ***********************************/
void	makeit(void)
  {
   fieldstruct		**infield,**inwfield, *outfield,*outwfield,
			*field, *wfield;
   catstruct		*cat;
   tabstruct		*tab;
   double		w,w1,w2, invscale;
   static int	       	list[MAXINFIELD],next[MAXINFIELD];
   int			i,j,k, ninfield, ntinfield;

/* Install the signal-catching routines for temporary file cleanup */
  install_cleanup();
/* Load input images */
  ninfield = prefs.ninfield;
/* First check input files and count FITS extensions when available */
  ntinfield = k = 0;
  for (i=0; i<ninfield; i++)
    {
    if (!(cat=read_cat(prefs.infield_name[i])))
      {
      sprintf(gstr, "*Error*: %s not found", prefs.infield_name[i]);
      error(EXIT_FAILURE, gstr,"");
      }
/*-- Examine all extensions */
    tab=cat->tab;
    for (j=0; j<cat->ntab; j++)
      {
      if (tab->naxis
	&& (!tab->tfields))
        {
        list[k++] = j;
        next[i]++;
        }
      tab = tab->nexttab;
      }
    ntinfield += next[i];
    if (!next[i])
      warning("No suitable data found in ", cat->filename);
    free_cat(cat, 1);
    }

/* Now load the ``fields'' */
  QMALLOC(infield, fieldstruct *, ntinfield);
  QMALLOC(inwfield, fieldstruct *, ntinfield);
  NFPRINTF(OUTPUT, "Examining input data...")
  k = 0;
  for (i=0; i<ninfield; i++)
    {
    for (j=0; j<next[i]; j++, k++)
      {
      infield[k] = load_field(prefs.infield_name[i], list[k], i);
/*---- Set flags */
      inwfield[k] = (prefs.weight_type[i])?
	load_weight(prefs.inwfield_name[i],infield[k], list[k], i,
	prefs.weight_type[i]) : NULL;
/*---- Now we can close input files */
      close_cat(infield[k]->cat);
      if (inwfield[k])
        close_cat(inwfield[k]->cat);
      }
    }

/* Load or create output image */
  outwfield = NULL;
  NFPRINTF(OUTPUT, "Creating NEW output image...")
  outfield = init_field(infield, inwfield, ntinfield, prefs.outfield_name);
  NFPRINTF(OUTPUT, "")
  QPRINTF(OUTPUT, "------- Output File %s:\n", outfield->rfilename);
  printinfo_field(outfield, outwfield);
  NFPRINTF(OUTPUT, "Creating NEW weight-map...")
  outwfield = init_weight(prefs.outwfield_name, outfield);

/* Compute projected limits and scaling in output frame (avoid unnecessary */
/* computations) */
  NFPRINTF(OUTPUT, "Framing and scaling images...")
  for (i=0; i<ntinfield; i++)
    {
    frame_wcs(infield[i]->wcs, outfield->wcs);
    if (prefs.fscalastro_type==FSCALASTRO_FIXED)
      scale_field(infield[i], outfield);
    }

/* Read and transform the data */
  NFPRINTF(OUTPUT, "Loading input data...")
  k = 0;
  w1 = w2 = invscale = 0.0;
  for (i=0; i<ninfield; i++)
    {
    for (j=0; j<next[i]; j++, k++)
      {
/*---- Display some info */
      if (!j)
        {
        NFPRINTF(OUTPUT, "")
        QPRINTF(OUTPUT, "-------------- File %s:\n", infield[k]->rfilename);
        }
      printinfo_field(infield[k], inwfield[k]);
/*---- Open input files */
      if (open_cat(infield[k]->cat, READ_ONLY) != RETURN_OK)
        error(EXIT_FAILURE, "*Error*: Cannot re-open ", infield[k]->filename);
      if (inwfield[k])
        {
        if (open_cat(inwfield[k]->cat, READ_ONLY) != RETURN_OK)
          error(EXIT_FAILURE, "*Error*: Cannot re-open ",
		inwfield[k]->filename);
        }
/*---- Pre-compute the background map */
      make_back(infield[k], inwfield[k]);
      if (inwfield[k])
        sprintf(gstr, "   Weight scale: %.7g", inwfield[k]->sigfac);
      else
        sprintf(gstr, "");
      NPRINTF(OUTPUT, "    Background: %.7g   RMS: %.7g%s\n",
        infield[k]->backmean, infield[k]->backsig, gstr);
/*---- Update gain information */
      w = 1/(infield[k]->backsig*infield[k]->backsig);
      w1 += w;
      if (infield[k]->gain)
        w2 += w/infield[k]->gain;
      invscale += 1/infield[k]->fscale;
/*---- Read (and convert) the weight data */
      if (inwfield[k])
        {
        sprintf(gstr, "Reading %s", inwfield[k]->filename);
        NFPRINTF(OUTPUT, gstr)
        read_weight(inwfield[k]);
        }
/*---- Read (and convert) the data */
      sprintf(gstr, "Reading %s", infield[k]->filename);
      NFPRINTF(OUTPUT, gstr)
      read_data(infield[k], inwfield[k]);
/*---- Resample the data (no need to close catalogs) */
      sprintf(gstr, "Resampling %s", infield[k]->filename);
      NFPRINTF(OUTPUT, gstr)
      resample_field(&infield[k], &inwfield[k], outfield, outwfield,
	prefs.resamp_type);
      }
    }

/* Approximation to the final equivalent gain */
  if (w2)
    outfield->gain = (w1/w2)*invscale;
  else
    outfield->gain = 0.0;

/* Go! */
  coadd_fields(infield, inwfield, ntinfield, outfield, outwfield,
		prefs.coadd_type, BIG);

/* Close files and free memory */
  NFPRINTF(OUTPUT, "Closing files...")
  for (k=0; k<ntinfield; k++)
    {
    end_field(infield[k]);
    if (inwfield[k])
      end_field(inwfield[k]);
    }
  free(infield);
  free(inwfield);

  end_field(outfield);
  end_field(outwfield);
  cleanup_files();

  return;
  }