/*
		     Copyright (c) 1991
		   Space Research Groningen
		    All Rights Reserved.
	P.O. Box 800, 9700 AV GRONINGEN. The Netherlands


@ subroutine irds_merge( character, character, integer )

#>            irds_merge.dc2

Function:     IRDS_MERGE

Purpose:      merge IRDSX in IRDS

Category:     IRAS

File:         irds_merge.c

Author:       Do Kester

Use:          call irds_merge( 
                        irdsx,    I     character*(*)
                        irds,     I/O   character*(*)
                        status )  I/O   integer
        irdsx   name of irds to merge with IRDS
        irds    name of irds to merge in
        status   0 : OK. 
		-1 : both nonexistent
                -2 : the two IRDS's dont pertain to the same plate
		-3 : not the same instrument
		-4 : partly or non existing IRDS

description:    
	IRDSX and IRDS should pertain to the same plate i.e. all the
	output parameters obtained with irds_enquire should be the
	same. If not the subroutine returns with status = -2.
	The snips in the two IRDS's which have the same sop/att
	combination are glued into one; blank values are filled 
	in between if the two are not adjacent in time. 
	The final IRDS will have a minimal snip duration, necessary 
	to fit all the snips in IRDS.
	The snips remain in increasing order of sop and att.

externals:      

updates:        23/02/90 DK, document creation
		13/11/90 DK, code
		06/02/91 DK, make global variables and functions `static'
		08/03/91 DK, copy NOISE descriptors too
		12/02/92 FL, also copy BPHF data
		26/02/92 DK, copy SCANTYPE
		03/03/92 DRB,copy OBS
		18/03/92 FL, new enquire_snip and extend, so copying OBS,
		             SCANTYPE and PLATE obsolete.
		07/05/92 DK, flag error for cripple IRDS's
		11/03/93 DK, delete partial irds and use userfio
#< 
*/

#include "string.h"
#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "gipsyc.h"
#include "irds_basic.h"
#include "ircc_rate.h"
#include "irds_rd_samples.h"
#include "irds_wr_samples.h"
#include "irds_rd_bphf.h"
#include "irds_wr_bphf.h"
#include "gds_rename.h"
#include "gds_delete.h"
#include "userfio.h"
#include "nelc.h"
#include "gdsd_rchar.h"
#include "gdsd_wchar.h"
#include "gdsd_rreal.h"
#include "gdsd_wreal.h"
#include "gdsd_rint.h"
#include "gdsd_wint.h"
#include "gdsc_word.h"
#include "gdsc_size.h"
#include "gdsc_range.h"
#include "gdsi_read.h"
#include "gdsi_write.h"

static void irds_fetch_snips( fchar, fchar, fint*, fint*, fint*, fint*, 
fint*, fint* );
static void irds_copy_snip( fchar, fchar, fint*, fint*, fint*, fint*, 
fint*, fint*, fint* );
static fint irdsd_noise( fchar, fchar, fchar, fint*, fint*, fint*, fint*,
fint*, fint* );

#define		FOREVER		for ( ; ; )
/* allocate n floats to pointer x, either new or overwriting old */
#define		fmalloc( x, n ) \
			( float * )(x) ? realloc( (x), (n) ) : malloc( n )
#define		fsize( n )	( (n) * sizeof( float ) )
#define		fallocmem( x, n )	( x = fmalloc( x, fsize( n ) ) )
/* allocate n doubles to pointer x, either new or overwriting old */
#define		dmalloc( x, n ) \
			( double * )(x) ? realloc( (x), (n) ) : malloc( n )
#define		dsize( n )	( (n) * sizeof( double ) )
#define		dallocmem( x, n )	( x = dmalloc( x, dsize( n ) ) )

#ifndef	NULL				/* stupid */
#define		NULL		0	/* null pointer */
#endif
#define	 	MAXAXIS		4
#define	 	MAXTAPES	10	/* number of tapes */
#define		LL		2	/* longitude and latitude */
#define		LP		21	/* length of character item in GDS */
#define		LI		20	/* length of character item in GDS */
#define		SNIPS		3	/* axis of the snips */
#define		DETS		2	/* axis of the detectors */
#define		TICKS		1	/* axis of the satcals */
#define		SAMPLES		0	/* axis of the smaples */
#define		OVERMAXSOPATT	999999  /* larger than any sop*1000+att */
#define		ERR_NONEXIST	0	/* irds's are non-existent */
/* error numbers 2 - 8 as in irds_enquire */
#define		ERR_INSTR	8	/* different instruments */
#define		ERR_CENTER	9	/* different centers */
#define		ERR_SIZE	10	/* different sizes */
#define		ERR_COOR	11	/* different coordinate systems*/
#define		ERR_IMPROPER	12	/* improper irds */
/* define error messages */
static char *errmes[] = {	
	"Both IRDS's are non-existent.",
	"it contains not enough (>=4) axes.",
	"it misses INSTRUME.",
	"it misses PRCLON and/or PRCLAT.",
	"it misses SIZELON and/or SIZELAT.",
	"it misses SKYSYS.",
	"it misses OBJECT.",
	"it misses EPOCHE.",    
	"The IRDS's have different instruments (INSTRUME).",
	"The IRDS's have different centers (PRCLON, PRCLAT).",
	"The IRDS's have different sizes (SIZELON, SIZELAT).",
	"The IRDS's have different coor. systems (SKYSYS,EPOCHE).",
	"is not a proper IRDS:", 
};
			

/* define global variables in this file */
static char	fc[21] = "\0" ;
static fint	sop1, obs1, att1, scancal1, scandur1, snipcal1, snipdur1;
static fchar	scantype1 = { fc, 20 } ;
static float	psi1, psirate1, theta1;
static fint	sop2, obs2, att2, scancal2, scandur2, snipcal2, snipdur2;
static fchar	scantype2 = { fc, 20 } ;
static float	psi2, psirate2, theta2;

/************************************************************************/

void irds_merge_c(
	fchar	irds1,
	fchar	irds2,
	fint 	*status )

{
	fint		naxis1, naxis2, axes1[MAXAXIS], axes2[MAXAXIS];
	fint		status1, status2;
	double		center1[LL], center2[LL], size1[LL], size2[LL];
	fchar		object1, object2, instru1, instru2, coor1, coor2;
	char		o1[LP], o2[LP], i1[LP], i2[LP], c1[LP], c2[LP];
	fint		sop, obs, att, scancal, scandur, snipcal, snipdur;
	fchar		scantype = { fc, 20 } ;
	float		psi, psirate, theta, epoch1, epoch2;
	fint		snip1, snip2, done1, done2, todo1, todo2;
	fint		snipdone1, snipdone2 ;
	fint		longsnip, axes[MAXAXIS], level, error;
	char		irm[15] = "MERGESCRATCH  " ;
	fchar		irmerge, observer;
	fint		startsnip, snip, rate;
	fint		noiser = 0;
	fint		no ;

	status1 = status2 = 0;
	object2.a = o2; object2.l = LI;
	instru2.a = i2; instru2.l = LI;
	coor2.a = c2; coor2.l = LI;

	irds_enquire_c( irds2, object2, instru2, &naxis2, axes2, center2, 
			size2, coor2, &epoch2, &status2 );

	object1.a = o1; object1.l = LI;
	instru1.a = i1; instru1.l = LI;
	coor1.a = c1; coor1.l = LI;
	irds_enquire_c( irds1, object1, instru1, &naxis1, axes1, center1, 
			size1, coor1, &epoch1, &status1 );

/* check for existence of at least one IRDS */
	if ( status1 == -1 && status2 == -1 ) {
		errorf( SERIOUS, errmes[ERR_NONEXIST] ) ;
		*status = -1;
		return;			/* both non-existent */
	}

/* check proper structure of IRDS's */
	if ( status1 < -2 || status2 < -2 ) {
		if ( status1 ) {
			errorf( SERIOUS, "%.*s %s", nelc_c( irds1 ),
					irds1.a, errmes[ERR_IMPROPER] ) ;
			errorf( SERIOUS, errmes[ 1 - status1 ] ) ;
		}
		if ( status2 )  {
			errorf( SERIOUS, "%.*s %s", nelc_c( irds1 ),
					irds2.a, errmes[ERR_IMPROPER] ) ;
			errorf( SERIOUS, errmes[ abs( status2 ) ] ) ;
		}
		*status = -2;
		return;			/* invalid IRDS('s) */
	}

/* check whether the IRDS's pertain to the same instrument and area */
	if ( !status2 && !status2 ) {
		if ( strncmp( instru1.a, instru2.a, LI ) ) {
			errorf( SERIOUS, errmes[ERR_INSTR] ) ;
			*status = -3;
			return;		/* different instruments */
		}

/* I do not think the following errors invalidate an irds completely; 
   if the user wants to mix data from different areas that is his/her choice.
   Give a warning message only. */
		if ( center1[1] != center2[1] || center1[0] != center2[0] )
			errorf( WARNING, errmes[ERR_CENTER] ) ;
		if ( size1[1] != size2[1] || size1[0] != size2[0] )
			errorf( WARNING, errmes[ERR_SIZE] ) ;
		if ( strncmp( coor1.a, coor2.a, LI ) || epoch1 != epoch2 )
			errorf( WARNING, errmes[ERR_COOR] ) ;
	}

/***
    use snipdone1/2 in stead of snip1/2 as the variables for the loop over
    the snip-axes. This to ensure that you do not lose the last snip of one
    of the two sets, which otherwise may happen once in a while.

    FL.
          If you do not understand, ask me.
***/
	snip1 = snip2 = snipdone1 = snipdone2 = 0;
	done1 = done2 = 1;
	todo1 = ( !status1 && naxis1 == 4 ) ? axes1[SNIPS] : 0;
	todo2 = ( !status2 && naxis2 == 4 ) ? axes2[SNIPS] : 0;
	longsnip = 0;
	while ( snipdone1 < todo1 || snipdone2 < todo2 ) {
		irds_fetch_snips( irds1, irds2, &snip1, &snip2, 
			&snipcal, &snipdur, &done1, &done2 );
		if (done1) snipdone1++;
		if (done2) snipdone2++;
		if ( snipdur > longsnip ) longsnip = snipdur;
	}

/* create a scratch IRDS of the proper length */
	axes[SAMPLES] = axes2[SAMPLES]; 
	axes[TICKS] = longsnip;
	axes[DETS] = axes2[DETS];
	irmerge.a = irm ; irmerge.l = 14 ;
	*status = 0; error = 0;
	level = 0;
	observer.a = o1; observer.l = LI;
	gdsd_rchar_c( irds2, tofchar( "OBSERVER" ), &level, observer, &error );
/* delete if an irds with this name already exists */
	irds_delete_c( irmerge, status );
	irds_create_c( irmerge, instru2, axes, center2, size2, coor2, &epoch2,
		object2, observer, status );
/* set BUNIT to the value in irds2 */
	error = 0;
	gdsd_rchar_c( irds2, tofchar( "BUNIT" ), &level, observer, &error );
	gdsd_wchar_c( irmerge, tofchar( "BUNIT" ), &level, observer, &error );

/***
    use snipdone1/2 in stead of snip1/2 as the variables for the loop over
    the snip-axes. This to ensure that you do not lose the last snip of one
    of the two sets, which otherwise may happen once in a while.

    FL.
          If you do not understand, ask me.
***/
	snip = snip1 = snip2 = snipdone1 = snipdone2 = 0;
	done1 = done2 = 1;
	rate = ircc_rate_c( instru2 );
	while ( snipdone1 < todo1 || snipdone2 < todo2 ) {
		snip++;
		irds_fetch_snips( irds1, irds2, &snip1, &snip2,
			&snipcal, &snipdur, &done1, &done2 );
		if ( done2 ) {		/* take remaining info from irds2 */
			sop = sop2; obs = obs2; att = att2;
			scancal = scancal2; scandur = scandur2;
			for( no = 0 ; no < scantype.l ; 
				scantype.a[no] = scantype2.a[no++] ) ;
			psi = psi2; psirate = psirate2; theta = theta2;
		} else {		/* take remaining info from irds1 */
			sop = sop1; obs = obs1 ; att = att1;
			scancal = scancal1; scandur = scandur1;
			for( no = 0 ; no < scantype.l ; 
				scantype.a[no] = scantype1.a[no++] ) ;
			psi = psi1; psirate = psirate1; theta = theta1;
		}

		irds_extend_c( irmerge, &sop, &obs, &att, scantype, &scancal, 
			&scandur, &snipcal, &snipdur, &psi, &psirate, &theta, 
			status );
		if ( done1 ) {
			snipdone1++ ;
			startsnip = snipcal1 - snipcal + 1;
			statusf( "%.*s snip %4d; "
				"copied to snip %4d at tick %d ", 
				nelc_c( irds1 ), irds1.a, 
				snip1, snip, startsnip );
			irds_copy_snip( irmerge, irds1, &snip, &snip1, 
				&startsnip, &snipdur1, &rate, &axes[DETS], 
				&error );
		}
		if ( done2 ) {
			snipdone2++ ;
			startsnip = snipcal2 - snipcal + 1;
			statusf( "%.*s snip %4d; "
				"copied to snip %4d at tick %d ", 
				nelc_c( irds2 ), irds2.a, 
				snip2, snip, startsnip );
			irds_copy_snip( irmerge, irds2, &snip, &snip2, 
				&startsnip, &snipdur2, &rate, &axes[DETS], 
				&error );
		}
		if ( ! noiser ) {
			noiser = irdsd_noise( irds1, irds2, irmerge,
				&snip1, &snip2, &snip, &done1, &done2, 
				&axes[DETS] );
		}
	}

	gds_delete_c( irds2, &error ) ;

	if ( gds_rename_c( irmerge, irds2 ) != 0 ) {
	  errorf( SERIOUS, "IRDS exists only partly" ) ;
	  *status = -4 ;
	}

	anyoutf( ANYOUT_DEF, "The new irds `%.*s' contains %d snips", 
		nelc_c( irds2 ), irds2.a, snip );
	anyoutf( ANYOUT_DEF, "and it has a sniplength of %d ticks.", longsnip );
	if ( snip == 0 ) {
	  anyoutf( ANYOUT_DEF, "No data found; the (partial) irds is deleted" );
	  error = 0 ;
	  gds_delete_c( irds2, &error ) ;
	}

	return;
}

/*********************************************************************/
static void irds_fetch_snips(
	fchar	irds1,
	fchar	irds2,
	fint	*snip1,
	fint	*snip2,
	fint	*snipcal,
	fint	*snipdur,
	fint	*done1,
	fint 	*done2 )

{
	fint		stat1, stat2, last, last1, last2;
	static int	sopatt1, sopatt2;

	stat1 = stat2 = 0;
	if ( *done1 ) {		/* fetch a new snip of irds1 */
		(*snip1)++;
		irds_enquire_snip_c( irds1, snip1, &sop1, &obs1, &att1, 
			scantype1, &scancal1,
 			&scandur1, &snipcal1, &snipdur1, &psi1, &psirate1,
			&theta1, &stat1 );
		if ( !stat1 ) sopatt1 = 1000 * sop1 + att1;
		else{
			anyoutf( ANYOUT_TST,
			"irds_fetch_snips: error in enquire_snip; set %.*s %d",
			nelc_c(irds2), irds2.a, stat1 ) ;
			sopatt1 = OVERMAXSOPATT;
		}
	}
	if ( *done2 ) {		/* fetch a new snip of irds2 */
		(*snip2)++;
		irds_enquire_snip_c( irds2, snip2, &sop2, &obs2, &att2, 
			scantype2, &scancal2,
			&scandur2, &snipcal2, &snipdur2, &psi2, &psirate2,
			&theta2, &stat2 );
		if ( !stat2 ) sopatt2 = 1000 * sop2 + att2;
		else{
			anyoutf( ANYOUT_TST,
			"irds_fetch_snips: error in enquire_snip; set %.*s %d",
			nelc_c(irds2), irds2.a, stat2 ) ;
			sopatt2 = OVERMAXSOPATT;
		}
	}
	anyoutf( ANYOUT_TST, "sopoatt1/soppatt2 %d %d", sopatt1, sopatt2 ) ;
	if ( sopatt1 == sopatt2 ){
		*snipcal = ( snipcal1 < snipcal2 ) ? snipcal1 : snipcal2;
		last1 = snipcal1 + snipdur1;
		last2 = snipcal2 + snipdur2;
		last = ( last1 > last2 ) ? last1 : last2 ;
		*snipdur = last - *snipcal;
		*done1 = 1; *done2 = 1;
	}
	else{
		if ( sopatt1 < sopatt2 ) {
			*snipdur = snipdur1;
			*snipcal = snipcal1;
			*done1 = 1; *done2 = 0;
		}
		else{
			*snipdur = snipdur2;
			*snipcal = snipcal2;
			*done1 = 0; *done2 = 1;
		}
	}
	return;
}

/**********************************************************************/

static void irds_copy_snip( 
	fchar	irmerge, 		/* to write in */
	fchar	irds, 			/* to read from */
	fint	*wsnip, 		/* snip number to write */
	fint	*rsnip, 		/* snip number to read */
	fint	*wstart, 		/* satcal to start writing */
	fint	*rlength,		/* ticks to read */
	fint	*rate,			/* samples per tick */
	fint	*maxdet,		/* number of detectors */		
	fint	*error )

{
	static	fint	prevticks = 0 ;
	fint		status = 0 ;
	fint		rtick, sdet, samps;
	static float	*reals = 0;
	static double	*srlon = NULL, *esrlon = NULL;
	static double	*srlat = NULL, *esrlat = NULL;
	static double	*twist = NULL;
	double		lngsun, sunrate ;

	if( *rlength > prevticks ){
		free( srlon  ) ; free( esrlon ) ; free( srlat  ) ;
		free( esrlat ) ; free( twist  ) ;
		if( (srlon = malloc(*rlength*sizeof(double))) == NULL ||
		    (esrlon = malloc(*rlength*sizeof(double))) == NULL ||
		    (srlat = malloc(*rlength*sizeof(double))) == NULL ||
		    (esrlat = malloc(*rlength*sizeof(double))) == NULL ||
		    (twist = malloc(*rlength*sizeof(double))) == NULL
		){
			*error = SERIOUS ;
			errorf( SERIOUS, "irds_merge: not enough memory" ) ;
			return ;
		}
		prevticks = *rlength ;
	}
	samps = *rate * *rlength ;
	if ( ! fallocmem( reals, samps ) ) {
		*error = SERIOUS ;
		errorf( SERIOUS, "irds_merge: not enough memory" ) ;
		return ;
	}

	rtick = 1;
	for ( sdet = 1; sdet <= *maxdet; sdet++ ) {
		irds_rd_samples_c( irds, rsnip, &sdet, &rtick, 
				reals, &samps, error );
		irds_wr_samples_c( irmerge, wsnip, &sdet, wstart,
				reals, &samps, error );
	}
	irds_rd_bphf_c( irds, rsnip, &rtick, 
				srlon, esrlon, srlat, esrlat, twist,
				&lngsun, &sunrate, rlength, &status ) ;
	if( status == 0 ) irds_wr_bphf_c( irmerge, wsnip, wstart,
				srlon, esrlon, srlat, esrlat, twist,
				&lngsun, &sunrate, rlength, &status ) ;
	return;
}

/************************************************************************/

/* combine the noise at snip-detector level in new irmerge */

static fint irdsd_noise( 
	fchar	irds1, 		/* In	irds to merge */
	fchar	irds2, 		/* In	irds to merge */
	fchar	irmerge, 	/* Out	irds to write */
	fint	*snip1, 	/* In	snip number in irds1 */
	fint	*snip2, 	/* In	snip number in irds2 */
	fint	*snip, 		/* In   snip number in irmerge */
	fint	*done1, 	/* In	if true do snip1 */
	fint	*done2,		/* In   if true do snip2 */
	fint	*maxdet )	/* In	number of detectors */

{
	fint		level1, level2, levelm, level, detaxis, snipaxis ;
	fint		err, sdet ;
	float		noise, noise1, noise2 ;
	
	err = level = level1 = level2 = levelm = 0 ;
	detaxis = DETS + 1 ;
	snipaxis = SNIPS + 1 ;
	if ( *done1 ) {
		level1 = gdsc_word_c( irds1, &snipaxis, snip1, &level1, &err ) ;
	}
	if ( *done2 ) {
		level2 = gdsc_word_c( irds2, &snipaxis, snip2, &level2, &err ) ;
	}
	levelm = gdsc_word_c( irmerge, &snipaxis, snip, &levelm, &err ) ;
	for ( sdet = 1; sdet <= *maxdet; sdet++ ) {
		noise = 0 ;
		if ( *done1 ) {
			err = 0 ;
			level = gdsc_word_c( irds1, &detaxis, &sdet,
						&level1, &err ) ;
			gdsd_rreal_c( irds1, tofchar( "NOISE" ), &level, 
						&noise1, &err ) ;
			if ( err != level ) return( err ) ;
			noise += noise1 * noise1 ;
		}
		if ( *done2 ) {
			err = 0 ;
			level = gdsc_word_c( irds2, &detaxis, &sdet,
						&level2, &err ) ;
			gdsd_rreal_c( irds2, tofchar( "NOISE" ), &level, 
						&noise2, &err ) ;
			if ( err != level ) return( err ) ;
			noise += noise2 * noise2 ;
		}
		noise = sqrt( noise / ( *done1 + *done2 ) ) ;
		err = 0 ;
		level = gdsc_word_c( irmerge, &detaxis, &sdet, &levelm, &err ) ;
		gdsd_wreal_c( irmerge, tofchar( "NOISE" ), &level, 
				&noise, &err ) ;
	}
	return( 0 ) ;
}