/****************************************************************/
/*      splitcomb						*/
/****************************************************************/
/*	Short Description :					*/
/*	Splits a ser file with interleaved data and combines	*/
/*	them with Hadamard coefficients.			*/
/****************************************************************/
/*	Syntax:							*/
/*	SPLITCOMB subspectra startexpno rowsperinc		*/
/*								*/
/*	subspectra: 	The number of subspectra (2)		*/
/*	startexpno:	Expno to put the first subspectra into,	*/
/*			defaults to expno*100+1			*/
/*	rowsperinc:	The number of rows (with TD-datapoints) */
/*			per increment. Depending on type of	*/
/*			FnMODE and order of the interleaved	*/
/*			data this can be 1 or 2. Defaults to 1.	*/
/*								*/
/*	Usage:							*/
/*	1	splitcomb 2					*/
/*			Splits to two expnos:			*/
/* 			Row 1+2,3+4,.. goes to the first 	*/
/* 			Row 1-2,3-4,.. goes to the second 	*/
/*								*/
/*	2	splitcomb 2 11 2				*/
/*			Splits to expnos 11 and 12:		*/
/* 			Row 1+3,2+4,.. goes to the first 	*/
/* 			Row 1-3,2-4,.. goes to the second 	*/
/*								*/
/*	3	splitcomb (without parameters)			*/
/*			Opens up a dialog window		*/
/*								*/
/* 	A Hadamard matrix is a square matrix whose entries are	*/
/*	either +1 or -1 and whose rows are mutually orthogonal.	*/
/*								*/
/*	The Hadamard coefficients (n=8) are:			*/
/*								*/
/*	+1 +1 +1 +1 +1 +1 +1 +1					*/
/*	+1 -1 +1 -1 +1 -1 +1 -1					*/
/*	+1 +1 -1 -1 +1 +1 -1 -1					*/
/*	+1 -1 -1 +1 +1 -1 -1 +1					*/
/*	+1 +1 +1 +1 -1 -1 -1 -1					*/
/*	+1 -1 +1 -1 -1 +1 -1 +1					*/
/*	+1 +1 -1 -1 -1 -1 +1 +1					*/
/*	+1 -1 -1 +1 -1 +1 +1 -1					*/
/*								*/
/*	This program supports up to eight subspectra. 		*/
/*								*/
/****************************************************************/
/*	Nils Nyberg, CRC, 2005-07-18				*/
/****************************************************************/
/*	Changelog						*/
/*	2005-09-16	Changed default rowsperinc to 1		*/ 
/****************************************************************/

#include <fcntl.h>
#include <limits.h>

#define MAXSIZE 16384
#define MAXOUTFILES 8

char	infile[PATH_MAX],
	outfile[MAXOUTFILES][PATH_MAX]; 
	
char	*outfilepointer[MAXOUTFILES];

int td1, td2, td1s, parmode, mc2, byteorder;
int i, j, k, l;
int subspectra, startexpno, rowsperinc;
int nargs;


int serindata[MAXOUTFILES][MAXSIZE];
int seroutdata[MAXOUTFILES][MAXSIZE];

int H[8][8] = {
	{+1, +1, +1, +1, +1, +1, +1, +1},
	{+1, -1, +1, -1, +1, -1, +1, -1},
	{+1, +1, -1, -1, +1, +1, -1, -1},
	{+1, -1, -1, +1, +1, -1, -1, +1},
	{+1, +1, +1, +1, -1, -1, -1, -1},
	{+1, -1, +1, -1, -1, +1, -1, +1},
	{+1, +1, -1, -1, -1, -1, +1, +1},
	{+1, -1, -1, +1, -1, +1, +1, -1}
	};


FILE *fpin, *fpout[MAXOUTFILES];


/* The program only works on 2D datasets */
GETCURDATA
FETCHPARS("PARMODE",&parmode);
if (parmode != 1) {STOPMSG( "Not a 2D dataset" );}

/* Get the processing parameter MC2 */
FETCHPAR1("MC2",&mc2);

/* Get the acq byte order */
FETCHPARS("BYTORDA",&byteorder);


/* Get arguments */

nargs = sscanf(cmd,"%d %d %d",&subspectra,&startexpno,&rowsperinc);

switch (nargs)
	{
	default:
	case -1:
	case 0:
		subspectra = 2;
		GETINT("Enter number of subspectra: ", subspectra);
		
		startexpno = 1 + expno * 100;
		GETINT("Enter start-expno", startexpno);
		
		rowsperinc = 1;
		GETINT("Enter number of rows per increment: ", rowsperinc);
		break;

	case 1:
		startexpno = 1 + expno * 100;
		rowsperinc = 1;
		break;
		
	case 2:
		rowsperinc = 1;
		break;
		
	case 3:
		break;
	}
(void)sprintf(text,"%d %d %d",subspectra, startexpno, rowsperinc);
Show_status(text);



/* Get and set dataset dimension */
FETCHPARS("TD",&td2)
td2 = ( (td2 + 255) / 256 ) * 256;

FETCHPAR1("TD",&td1)
FETCHPAR1S("TD",&td1s)

/* Set ser-file name (infile) */
(void)sprintf(infile,"%s/data/%s/nmr/%s/%d/ser",disk,user,name,expno);

/* Set output ser files (outfile) */
for (i=0; i < subspectra; i++)
	{
	(void)sprintf(outfile[i],"%s/data/%s/nmr/%s/%d/ser",disk,user,name,startexpno+i);
	outfilepointer[i] = outfile[i];
	}
	
/* Copy raw data to new datasets */
for (i=0; i < subspectra; i++)
	{
	WRA(startexpno+i);
	}

/* Open files for read and write */
fpin=fopen(infile,"rb");
for (i=0; i < subspectra; i++)
	{
	fpout[i]=fopen(outfilepointer[i],"wb");
	}
	
/* Do the actual splitting and combinations */
for (i=0; i < td1s/(rowsperinc * subspectra); i++)
	{
	/* Read in data */
	for (j=0; j < subspectra; j++)
		{
		fread(serindata[j],sizeof(int),td2*rowsperinc,fpin);
		local_swap4(serindata[j],sizeof(int)*td2*rowsperinc,byteorder);
		}
	/* Do the combinations */
	for (j=0; j < subspectra; j++)
		{
		for (k=0; k < subspectra; k++)
			{
			for (l=0; l < td2*rowsperinc; l++)
				{
				if (k==0)
					{
					seroutdata[j][l]  = serindata[k][l] * H[j][k];
					}
				else
					{
					seroutdata[j][l] += serindata[k][l] * H[j][k];
					}
				}
			}
		}
	/* Write the data back */
	for (j=0; j < subspectra; j++)
		{
		local_swap4(seroutdata[j],sizeof(int)*td2*rowsperinc,byteorder);
		fwrite(seroutdata[j],sizeof(int),td2*rowsperinc,fpout[j % subspectra]);
		}
	}
	
/* Close the opened files */
fclose(fpin);
for (i=0; i < subspectra; i++)
	{
	fclose(fpout[i]);
	}

/* Adjust td1 and td1s to reflect the new datasize */
td1  /= subspectra;
td1s /= subspectra;

for (i=0; i < subspectra; i++)
{	
	DATASET (name,startexpno+i,procno,disk,user)
	STOREPAR1("TD",td1)
	STOREPAR1S("TD",td1s)
}

	
	
	
	
(void)sprintf(	text,
		"Finished; Subsp = %d, Expno = %d-%d, Rows/inc = %d",
		subspectra, startexpno, startexpno+subspectra-1, rowsperinc);

Show_status(text);

QUIT