/* @(#)spmerge.c	17.1.1.1 (ES0-DMD) 01/25/02 17:56:15 */
/*===========================================================================
  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 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding 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
===========================================================================*/

/* Program : spmerge.c
/* Author  : P. Ballester   -     ESO Garching                    */
/* Date    : 08.09.94      Version 1.0                            */
/*           17.09.99      do not divide by 0.0 otherwise you may */
/*                         get problems with STAT/IMA (SW).       */
/*                                                                */
/* Purpose :                                                      */
/*                                                                */
/* Merging two spectra (ramp, optimal, constant weights)          */ 
/*                                                                */
/* Input:                                                         */  
/*             - name of input spectrum 1       : IN_A            */
/*             - name of input spectrum 2       : IN_B            */
/*             - name of input std deviation 1  : OUT_A           */
/*             - name of input std deviation 2  : OUT_B           */
/*             - Width of zeroed interval       : INPUTI(1)       */
/*             - name of output                 : INPUTC          */
/*                                                                */
/* Notes: Purified ...                                            */

#include <math.h>
#include <midas_def.h>
#define  ipos2D(col,row,npix) row*npix[0]+col

main ()

{
    char    inspec1[100], inspec2[100], weight1[100];
    char    weight2[100], outframe[100];
    char    cunit[64], ident[72], mode[12];
    int     actvals, x, naxis, row, col, unit, null;
    int     imns1, imns2, imnw1, imnw2, imnout;
    float   *pntrs1, *pntrs2, *pntrw1, *pntrw2, *pntrout, v1, v2, v12;
    double  end[2], zeroed, lambda;
    int     pixel[3], interv, limit;
    int     npixs1[2], npixs2[2], npixw1[2], npixw2[2], npix[2];
    double  starts1[2], starts2[2], startw1[2], startw2[2], start[2];
    double  steps1[2], steps2[2], stepw1[2], stepw2[2], step[2];
    int     xs1, xs2, xw1, xw2;
    int     cnt;

    SCSPRO ("splocext");

    strcpy(cunit,""), strcpy(ident,"");

    SCKGETC ("IN_A",   1, 60, &actvals, inspec1);
    SCKGETC ("IN_B",   1, 60, &actvals, inspec2);
    SCKGETC ("OUT_A",  1, 60, &actvals, weight1);
    SCKGETC ("OUT_B",  1, 60, &actvals, weight2);
    SCKGETC ("INPUTC", 1, 60, &actvals, outframe);
    SCKRDI  ("INPUTI", 1, 1,  &actvals, &limit, &unit, &null);

    SCIGET(inspec1, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 1, 
           &naxis, npixs1,starts1, steps1, ident, cunit, (char **)&pntrs1, &imns1);

    SCIGET(inspec2, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 1, 
           &naxis, npixs2, starts2, steps2, ident, cunit, (char **)&pntrs2, &imns2);

    SCIGET(weight1, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 1, 
           &naxis, npixw1, startw1, stepw1, ident, cunit, (char **)&pntrw1, &imnw1);

    SCIGET(weight2, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 1, 
           &naxis, npixw2, startw2, stepw2, ident, cunit, (char **)&pntrw2, &imnw2);

    strcpy(ident,"Optimal merging of spectra  ");
    strcat(ident,inspec1);
    strcat (ident," and ");
    strcat(ident,inspec2);

    if (starts1[0] > starts2[0]) 
    SCETER(9,
    "MERGE/SPEC: Spectrum 1 must precede spectrum 2 in wavelength ...\n");
    if (steps1[0] != steps2[0] || 
        steps1[0] != stepw1[0] || 
        steps1[0] != stepw2[0])
    SCETER(9,
    "MERGE/SPEC: Steps must be equal for all input images...Exit.\n");

    if (naxis == 1) npix[1]=1;
    end[0]   = starts1[0] + (npixs1[0] -1.)*steps1[0];
    end[1]   = starts2[0] + (npixs2[0] -1.)*steps2[0];

    start[0] = starts1[0], step[0] = steps1[0];
    npix[0] = (int)((end[1]-starts1[0])/step[0]) + 1;
    if (naxis == 1) npix[1]=1;
    else            start[1]=starts1[1], step[1]=steps1[1], npix[1]=npixs1[1];

    zeroed   = limit*step[0];

    SCIPUT(outframe, D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
           npix, start, step, ident, cunit, (char **)&pntrout, &imnout);

    cnt = 0;
    for (row=0; row<npix[1]; row++) {
    for (col=0; col<npix[0]; col++) {
        lambda = start[0] + col*step[0];  /* Current wavelength */
        xs1    = (lambda - starts1[0])/steps1[0]; /* Position in spectrum 1 */
        xs2    = (lambda - starts2[0])/steps2[0]; /* Position in spectrum 2 */
        xw2    = (lambda - startw2[0])/stepw2[0]; /* Position in weight 2 */
        xw1    = (lambda - startw1[0])/stepw1[0]; /* Position in weight 1 */

        if (lambda < (start[0]+zeroed))
            pntrout[ipos2D(col,row,npix)] = 0.;
        if (lambda >= (start[0]+zeroed) && lambda < (starts2[0]+zeroed))
            pntrout[ipos2D(col,row,npix)] = pntrs1[ipos2D(xs1,row,npixs1)];
        if (lambda > (end[0]-zeroed) && lambda <= (end[1]-zeroed))
            pntrout[ipos2D(col,row,npix)] = pntrs2[ipos2D(xs2,row,npixs2)];
        if (lambda > (end[1]-zeroed)) 
            pntrout[ipos2D(col,row,npix)] = 0.; 

        if (lambda >= (starts2[0]+zeroed) && lambda <= (end[0]-zeroed))
	  {
	    v1 = pow(pntrw1[ipos2D(xw1,row,npixw1)], 2.);
	    v2 = pow(pntrw2[ipos2D(xw2,row,npixw2)], 2.);
	    v12 = v1 + v2;
	    if (v12 != 0.0)
	      pntrout[ipos2D(col,row,npix)] = 
		(v1 * pntrs1[ipos2D(xs1,row,npixs1)] + 
		 v2 * pntrs2[ipos2D(xs2,row,npixs2)]) / v12;
	    else
	      {
		pntrout[ipos2D(col,row,npix)] = 0.0;
		cnt++;
	      }
	  }
	}
     } 
    if (cnt > 0) {
      sprintf(outframe,"%d undefined pixels ... set to 0.0\n",cnt);
      SCTPUT(outframe);
    }

    SCSEPI();

}