/*** ^^A -*-C++-*- ******************************************************/
/*  pshift2D3D4D      07 August 2013                                    */
/************************************************************************/
/*  Short Description :                                                 */
/*  AU program to reconstruct Zangger-Sterk based                       */
/*  2D, 3D and 4D Pureshift                                             */
/************************************************************************/
/*  Keywords :                                                      	*/
/*  pureshift, Zangger-Sterk, BIRD, DOSY                             	*/
/************************************************************************/
/*  Description/Usage :                                           		*/
/* 	This AU program reconstructs 2D, 3D Pureshift (aqseq 312), not 4D   */
/*  Pureshift (aqseq 4123) data                                         */
/*	data recorded with Zangger-Sterk based pulse sequences.   			*/
/*	Depending on the dimensionality, 2D, 3D or 4D processing is   		*/
/*	employed. First, data is converted to analog form into temporary	*/
/*	expno 9999. Further, the user is promted to enter an unused expno   */
/*  for the output data (default is expno+1). A template for the output */
/*	data is then generated using the macros RSER (for 2D) and 			*/
/*	RSER2D (for 3D, 4D), respectively. Consequently, the data is read   */
/*	in (complete ser file for 2D data, 13 plane ser file for 3D data). 	*/
/*	After concatenation of the data to a single fid (2D) and 			*/
/*	ser file (3D,4D), respectively, the status parameters in the output	*/
/*	dataset are adjusted.			      				                */
/************************************************************************/ 
/*  Author(s) :                                                 		*/
/*  Name          : Ralph Adams            				                */
/*  Organisation  : University of Manchester, UK                        */
/*  Email         : ralph.adams@manchester.ac.uk                        */
/************************************************************************/
/* Version History, original 2D macro written by Aitor Morreno@Bruker   */
/* 01JUN2011 modified to include debugging and correction for non-256   */
/*           block sized data. Also added comments.  	(R. Adams)		*/
/* 01NOV2011 modified to incorporate TS2 and TS3 compatibility using    */
/*           precompiler checking with #ifdef and #ifndef and version   */
/*           checking using getxwinvers() to correct RSER2D  (R. Adams) */
/* 21NOV2011 updates td in acquisition as well as status files (RA)     */
/* 27JUL2012 modified macro to cope with 2D->1D files larger than       */
/*           defined MAXSIZE that resulted in a seqmentation fault (RA) */
/* 14SEP2012 modified macro to correct calculation of sw/sw1 (GAM)      */
/* 07AUG2013 4D->3D loop added but not operational yet(RA)              */
/* 02OCT2013 Added annoying pop up to make sure use is acknowledged (RA)*/
/************************************************************************/


/* For TD = 32 and lower */
#define MAXSIZE 655360


int row[MAXSIZE], temp[MAXSIZE], temp2[256];
int pmod, si1, si2, si3, si12, si13, si12c, drppts, npoints, nzero=0, i4, debug, download, infomsg, twoFiveSixCorr;
int aexpno=9999, pexpno=9998, nexpno=expno+1000;
double sw1, sw2, sw3;
double sfo1, sfo2, sfo3;
float l4;
char infile[PATH_MAX], outfile[PATH_MAX], statustext[256];
FILE *fpin, *fpout;

GETCURDATA

/*Flags - set all of them to 0 for an easy life */
/*debug = 1 for on, 0 for off*/
debug=1;
/*download = 1 to display pop up, 0 to skip*/
download=1;
/*Advice about processing: infomsg = 1 for on, 0 for off*/
infomsg=1;

if(download == 1){
i1 = Proc_err(QUESTION_OPT,"By using this macro you agree to acknowledge \n\"http://nmr.chemistry.manchester.ac.uk\" \nin any resulting publication. \n\nClick OK to agree and continue. \nClick cancel to disagree and abort.\n\n This message can be suppressed by setting the variable \n\"download=0;\" in the macro.");
if ( i1 == ERR_CANCEL )
{
ABORT
}
else{
}
}

if(debug == 1){(void)printf("Debugging on\nsuppress debugging by setting debug=0 in macro\n");}
if(debug == 1){
	if(infomsg == 1){
		(void)printf("Information messages will be shown during processing on\n");
	}
}
if(debug == 1){
	if(download == 1){
		(void)printf("Download agreement will be shown during processing\n");
	}
}


/*Get the version of Topspin that is running so that the appropriate syntax and macros are used*/
char curversion[80];
getxwinvers(curversion);
if(debug == 1){(void)printf("Topspin Version: %s\n",curversion);}
const char* from = curversion;
  char *vers = (char*) malloc(6);
  strncpy(vers, from+7, 1);
  if(debug == 1){(void)printf("Topspin Version: %s\n",vers);}
  int versn=atoi(vers);
  
  if(debug == 1){
  (void)printf("Version detection char to int confirmation:\n");
  if(versn==3){(void)printf("Topspin Version: %i Detected Successfully.\n",versn);}
  if(versn==2){(void)printf("Topspin Version: %i Detected Successfully.\n",versn);}
  if(versn==1){(void)printf("Topspin Version: %i Detected Successfully.\n",versn);}   
  }

 
	#ifdef NEWACQUPATH
	if(debug == 1){(void)printf("NEWACQUPATH is defined and will be used.\n");}
	#endif
	
	#ifndef NEWACQUPATH
	if(debug == 1){(void)printf("NEWACQUPATH is not defined, ACQUPATH will be used.\n");}
	#endif
	
	#ifdef ACQUPATH
	if(debug == 1){(void)printf("ACQUPATH is defined and will be used.\n");}
	#endif
	
	#ifndef ACQUPATH
	if(debug == 1){(void)printf("ACQUPATH is not defined.\n");}
	#endif
  
/* Function to round double to integer */
int round_to_int(double d) {
	return d<0?d-.5:d+.5;
	}

int round_to_int_f(float f) {
	return f<0?f-.5:f+.5;
	}
	
/* Check experiment type, 2D or 3D */
FETCHPAR("PARMODE",&pmod)
if ( (pmod < 1) || (pmod > 2) )
    STOPMSG("Not a 2D or 3D dataset.\nTry changing pparmod in the processing parameters.")

/* 2D */ 
if (pmod == 1){ 
/*Get status parameters */
    if(debug == 1){(void)printf("Starting 2D conversion\n");}
	FETCHPAR1S("TD",&si1)
	FETCHPAR("TD",&si2)
	FETCHPAR1S("SW",&sw1)
	FETCHPAR1S("SFO1",&sfo1)
	sw1=sw1*sfo1;
	FETCHPARS("SW",&sw2)
	FETCHPARS("SFO1",&sfo2)
	sw2=sw2*sfo2;
	si12=si1*si2;

	twoFiveSixCorr=256-(si2%256);
	if(twoFiveSixCorr == 256){twoFiveSixCorr = 0;}
		if(debug == 1){(void)printf("Correction Value for si2 not divisible by 256: %d\n",twoFiveSixCorr);}
		if(debug == 1){(void)printf("Full si2: %d\n",si2+twoFiveSixCorr);}
	si12c=si1*(si2+twoFiveSixCorr);
	if(debug == 1){(void)printf("File size (si12c) calculated as: %d\n",4*si12c);}
	
	if(si12c <=MAXSIZE){
	/* Digitial to analog conversion, define Topspin input ("ser") and output ("fid") data files */
		CONVDTA(aexpno)
		REXPNO(aexpno)
		strcpy(infile, ACQUPATH("ser"));
		GETINT("Enter new expno for 1D pureshift data:  ",nexpno) 
		RSER(1,nexpno,1)

		/*use the precompiler to check if the appropriate Topspin 2 or 3 macros exist*/
		#ifdef NEWACQUPATH
		if(debug == 1){(void)printf("ifdef TS3\n");}  
		strcpy(outfile,NEWACQUPATH(nexpno, "fid"));
		#endif
	
		#ifndef NEWACQUPATH
		if(debug == 1){(void)printf("ifndef TS2\n");}  
		REXPNO(nexpno)
		strcpy(outfile,ACQUPATH("fid"));
		REXPNO(aexpno)
		#endif
	


	/* Define the start position and the number of points to extract for concatenation of the FID */
		FETCHPAR("CNST 4",&l4)
		if(debug == 1){(void)printf("Value of cnst4 passed to macro: %f\n",l4);}
			
		drppts=round_to_int_f(l4);
		if(debug == 1){(void)printf("Value of drppts passed on in macro: %d\n",drppts);}
		drppts=(drppts*2);
	
		/*drppts=2;*/
		if(debug == 1){(void)printf("'Complex' points to drop: %d\n",drppts);}
	
		npoints=round_to_int(sw2/sw1);
	
		/*read and store the entire 2D SER file*/
		if((fpin = fopen(infile,"r")) == NULL){STOPMSG("Open of infile failed!\n")}
		if(fread(row,sizeof(int),si12c,fpin)!=si12c){STOPMSG("Read in of ser file failed!\n")}
		if(debug == 1){(void)printf("Full file size should be: %d\n",si1*(si2+twoFiveSixCorr));}
		fclose(fpin);	 
	
	/* Delete the temporary expno 9999 */
		DELETEEXPNO(name, aexpno, disk, user)

	/* extract and concatenate the chunked FID's to single FID */
		if((fpout=fopen(outfile,"w")) == NULL){STOPMSG("Open of outfile failed!\n")}
		i1=0;
		while(i1 < si1){
			if(debug == 1){(void)printf("%d \n",i1);}
			for (i2=0; i2<2*npoints; i2++){
				/*temp[i2]=row[i1*si2+i2+drppts];*/
				sprintf(statustext,"%d of %d (%d.%d)  ",i1,si1,i1,i2);
				Show_status(statustext);
				if(debug == 1){(void)printf("%d.%d(%d)  ",i1,i2,i1*(si2+twoFiveSixCorr)+i2+drppts);}
				temp[i2]=row[i1*(si2+twoFiveSixCorr)+i2+drppts]; /*added correction for the case that si2 is not a multiple of 256*/
				}
			if(fwrite(temp, sizeof(int), 2*npoints, fpout) != 2*npoints) {STOPMSG("writing to outfile failed")}
			i1++;
			if(debug == 1){(void)printf("\n",i1);}		
		}
		fclose(fpout);

	/* Set status parameter TD in the output dataset according to constructed FID length */
		REXPNO(nexpno)
		SETCURDATA
		STOREPARS("TD",2*npoints*si1)                                           /*status*/
		STOREPAR("TD",2*npoints*si1)                                               /*acq*/
		if(debug == 1){(void)printf("Finished 2D conversion\n");}
		VIEWDATA_SAMEWIN
	}
	else{
		if(debug == 1){(void)printf("Slow version of macro is being used - multiple writes to disk rather than using just RAM\n");}
		
	/* Digitial to analog conversion, define Topspin input ("ser") and output ("fid") data files */
		CONVDTA(aexpno)
		REXPNO(aexpno)
		strcpy(infile, ACQUPATH("ser"));
		GETINT("Enter new expno for 1D pureshift data:  ",nexpno) 
		RSER(1,nexpno,1)
		
		if(infomsg == 1){
			Proc_err(INFO_OPT,"The size of the data set is large.\nProcessing may be slow.\nUnless there is a specific reason for acquiring such a big data set reduce the\nnumber of chunks (TD 1) or the number of points acquired per chunk (TD 2).\n\n This message can be suppressed by setting the infomsg flag in the macro.");
		}

		/*use the precompiler to check if the appropriate Topspin 2 or 3 macros exist*/
		#ifdef NEWACQUPATH
		if(debug == 1){(void)printf("ifdef TS3\n");}  
		strcpy(outfile,NEWACQUPATH(nexpno, "fid"));
		#endif
	
		#ifndef NEWACQUPATH
		if(debug == 1){(void)printf("ifndef TS2\n");}  
		REXPNO(nexpno)
		strcpy(outfile,ACQUPATH("fid"));
		REXPNO(aexpno)
		#endif
	
		if((fpin = fopen(infile,"r")) == NULL){STOPMSG("Open of infile failed!\n")}

	/* Define the start position and the number of points to extract for concatenation of the FID */
		FETCHPAR("CNST 4",&l4)
		if(debug == 1){(void)printf("Value of cnst4 passed to macro: %f\n",l4);}
			
		drppts=round_to_int_f(l4);
		if(debug == 1){(void)printf("Value of drppts passed on in macro: %d\n",drppts);}
		drppts=(drppts*2);
	
		/*drppts=2;*/
		if(debug == 1){(void)printf("'Complex' points to drop: %d\n",drppts);}
	
		npoints=round_to_int(sw2/sw1);
	
	/* extract and concatenate the chunked FID's to single FID */
		if((fpout=fopen(outfile,"w")) == NULL){STOPMSG("Open of outfile failed!\n")}
		i1=0;
		while(i1 < si1){
			if(debug == 1){(void)printf("%d \n",i1);}

			/*read and store each FID from the SER file*/
			if(fread(row,sizeof(int),si2+twoFiveSixCorr,fpin)!=si2+twoFiveSixCorr){STOPMSG("Read in of ser file failed!\n")}
			
			for (i2=0; i2<2*npoints; i2++){
				sprintf(statustext,"%d of %d (%d.%d)  ",i1,si1,i1,i2);
				Show_status(statustext);
				if(debug == 1){(void)printf("%d.%d(%d)  ",i1,i2,i1*(si2+twoFiveSixCorr)+i2+drppts);}
				temp[i2]=row[i2+drppts]; /*added correction for the case that si2 is not a multiple of 256*/
				}
			if(fwrite(temp, sizeof(int), 2*npoints, fpout) != 2*npoints) {STOPMSG("writing to outfile failed")}
			i1++;
			if(debug == 1){(void)printf("\n",i1);}		
		}
		fclose(fpout);
		fclose(fpin);	
		
	/* Delete the temporary expno 9999 */
		DELETEEXPNO(name, aexpno, disk, user)
		
	/* Set status parameter TD in the output dataset according to constructed FID length */
		REXPNO(nexpno)
		SETCURDATA
		STOREPARS("TD",2*npoints*si1)                                           /*status*/
		STOREPAR("TD",2*npoints*si1)                                               /*acq*/
		if(debug == 1){(void)printf("Finished 2D conversion\n");}
		VIEWDATA_SAMEWIN
	
	}

}

/* 3D */
if (pmod == 2){
	/* Get required status parameters for all dimensions */
    if(debug == 1){(void)printf("Starting 3D conversion\n");}
	FETCHPAR3S("TD",&si1)   /*  zs  */
	FETCHPARNS(2,"TD",&si2) /*  noe/roe/tocsy */
	FETCHPAR("TD",&si3)     /*  fid */
	FETCHPAR3S("SW",&sw1)   /**/
	FETCHPAR3S("SFO1",&sfo1)   /**/
	FETCHPARNS(2,"SW",&sw2) /**/
	FETCHPARNS(2,"SFO1",&sfo2) /**/
	FETCHPARS("SW",&sw3)    /**/
	FETCHPARS("SFO1",&sfo3)    /**/
	sw1=sw1*sfo1;
	sw2=sw2*sfo2;
	sw3=sw3*sfo3;
	si13=si3*si1;
	
	/*STOPMSG("This version of the macro is not compatible with such a large amount of data, \n reduce the number of chunks (TD 1) or the number of points acquired per chunk (TD 2).")*/
	
	if(debug == 1){(void)printf("si1: %d\nsi2: %d\nsi3: %d\nsw1: %d\nsw2: %d\nsw3: %d\n",si1,si2,si3,sw1,sw2,sw3);}
	if(debug == 1){(void)printf("Data Size for each plane: %d\n",si13);}
	
/*If using Topspin 3, define Topspin input data file */
#ifdef NEWACQUPATH
	if(debug == 1){(void)printf("ifdef TS3\n");}
	strcpy(infile,NEWACQUPATH(pexpno, "ser"));
	if(debug == 1){
		sprintf(statustext,"If using Topspin 3, define Topspin input data file pexpno");
		Show_status(statustext);
	}
#endif
	
/* digital to analog conversion */
	CONVDTA(aexpno)
	REXPNO(aexpno)
	if(debug == 1){(void)printf("Converted data to analog\n");}
	
/*If using Topspin 2, define Topspin input data file */
#ifndef NEWACQUPATH
	if(debug == 1){(void)printf("ifndef TS2\n");}  
 	if(debug == 1){(void)printf("RSER2D \"s13\"\n");}   
	if(debug == 1){
		sprintf(statustext,"If using Topspin 2, define Topspin input data file pexpno (1)");
		Show_status(statustext);
	}
	RSER2D("s13",1,pexpno)	
	if(debug == 1){
		sprintf(statustext,"If using Topspin 2, define Topspin input data file pexpno (2)");
		Show_status(statustext);
	}
	REXPNO(pexpno) 
	strcpy(infile,ACQUPATH("ser"));
	REXPNO(aexpno)
#endif

/* Define Topspin output data file, create 2D output file using RSER2D */ 
	GETINT("Enter new expno for 2D pureshift data:  ",nexpno)
	
	/*RSER2D version compatibility*/
	if(versn==3){
	 	if(debug == 1){(void)printf("RSER2D 23\n");}   
		RSER2D(23,1,nexpno)
	}
	else{
	 	if(debug == 1){(void)printf("RSER2D \"s23\"\n");}
		RSER2D("s23",1,nexpno)
	}

/*Define Topspin input data file, version dependent */
#ifdef NEWACQUPATH
	if(debug == 1){(void)printf("ifdef TS3\n");}  
	strcpy(outfile,NEWACQUPATH(nexpno, "ser"));
#endif	
#ifndef NEWACQUPATH
	if(debug == 1){(void)printf("ifndef TS2\n");}  
	REXPNO(nexpno)
	strcpy(outfile,ACQUPATH("ser"));
	REXPNO(aexpno)
#endif	

	REXPNO(aexpno)

/* Define the start position and the number of points to extract for concatenation of the FID */
	FETCHPAR("CNST 4",&l4)
	if(debug == 1){(void)printf("Value of cnst4 passed to macro: %f\n",l4);}
		
	drppts=round_to_int_f(l4);
	if(debug == 1){(void)printf("Value of drppts passed on in macro: %d\n",drppts);}
	drppts=(drppts*2);

	/*drppts=2;*/
	if(debug == 1){(void)printf("Real + Imaginary points to drop: %d\n",drppts);}

	npoints=round_to_int(sw3/sw1);
	
/* Create array filled with zeros to properly process the direct dimension */
/* The data has to be a multiple of 256 due to the way Bruker saves it in blocks*/
	nzero=256-(2*npoints*si1)%256;
	if (nzero!=256)	{
		for(i3=0; i3<nzero; i3++){
			temp2[i3]=NULL;
		}
	}
	if(debug == 1){(void)printf("Zero Padding: %d\n",nzero);}
	
/* As the data was saved in blocks which are multiples of 256 we have to jump the zeroes at the end of each fid if each one isn't a multiple of 256 */
/*	twoFiveSixCorr=256-(si3%256);*/
/*	if(debug == 1){(void)printf("Correction Value for si3 not divisible by 256: %d\n",twoFiveSixCorr);}*/
/*	if(debug == 1){(void)printf("Remainder for si3/256 : %d\n",si3%256);} */
	
	twoFiveSixCorr=256-(si3%256);
	if(twoFiveSixCorr == 256){twoFiveSixCorr = 0;}
	if(debug == 1){(void)printf("Correction Value for si3 not divisible by 256: %d\n",twoFiveSixCorr);}
	
	
/* extract and concatenate the chunked FID's to single FID's and add them into 2D ser file */
	if ((fpout=fopen(outfile,"w")) == NULL){
		STOPMSG("Open of outfile failed!\n")
	}
	
	/*Put some information into the status bar*/
	sprintf(statustext,"Starting Processing Loop");
	Show_status(statustext);
	ssleep(2); /*Wait for a couple of seconds so that it can be seen*/
	
	for(i4=1; i4<=si2; i4++){ 																/*for pointer i4 go from 1 to si2 (noe) through si2 planes of size si1*si3 */ /*changed "i4=i4++" to "i4++" */
		if(debug == 1){(void)printf("\nLoop position %d.",i4);}
		
		if(versn==3){
		    if(debug == 1){(void)printf("RSER2D 13  (v3)");} 
			if(debug == 1){
				sprintf(statustext,"RSER2D(13,i4,pexpno)");
				Show_status(statustext);
			}	
			RSER2D(13,i4,pexpno)
			}														/*Detect which version of RSER2D to use then... */
		else{
			if(debug == 1){
			sprintf(statustext,"Delete any redundant data");
			Show_status(statustext);
			}
			DELETEEXPNO(name, pexpno, disk, user)													/*Delete any redundant data*/
			if(debug == 1){(void)printf("RSER2D \"s13\" (v1 or 2)");} 
			if(debug == 1){
				sprintf(statustext,"RSER2D(s13,i4,pexpno)");
				Show_status(statustext);
			}
			RSER2D("s13",i4,pexpno)
			}															/*... read 13 plane with increment in si2*/																
		
		if ((fpin = fopen(infile,"r")) == NULL){STOPMSG("Open of ser file failed!\n")}			/*Read the data in then check if the data was read in*/
		if (fread(row,sizeof(int),si13,fpin)!=si13){STOPMSG("Read in of ser file failed!\n")}   /*Read the data into parameter 'row; and check the size of the data*/
		fclose(fpin); 																			/*Close the data file*/
		i1=0;																					/*Set up the pointer for the loop over the ZS loop*/
		while(i1 < si1)	{																		/*Loop over si1 (ZS) with pointer i1 from 0 to si1-1*/
			if(debug == 1){(void)printf("\nLoop position %d.%d\n",i4,i1);}
			for (i2=0; i2<2*npoints; i2++){														/*Loop over 2*npoints for the data transfer to the temporary file*/
				sprintf(statustext,"%d of %d (%d.%d.%d)  ",i4,si2,i4,i1,i2);
				Show_status(statustext);
				if(debug == 1){(void)printf("%d.%d.%d ",i4,i1,i2);}
				temp[i2]=row[i1*(si3+twoFiveSixCorr)+i2+drppts]; 			/*transfer data in defined position to temporary file, fid, added correcdtion if si2 is not a multiple of 256*/
			}
			if(fwrite(temp, sizeof(int), 2*npoints, fpout) != 2*npoints){STOPMSG("writing to outfile failed")} /*Write the data out and check it was written out ok*/
			i1++;																				/*increment the pointer*/
		}
		if (nzero!=256){												
			if(fwrite(temp2, sizeof(int), nzero, fpout) != nzero){STOPMSG("writing of zero's to outfile failed")} /*add some zeros if the file size is not a multiple of 256 then check this worked ok*/
		}
	}
	if(debug == 1){(void)printf("\n");}
	fclose(fpout); 		/*close the output file*/

/* 2D loop end */
	
/* Delete the temporary expno's 9998 and 9999 */
	DELETEEXPNO(name, aexpno, disk, user)
	if(debug == 1){
		sprintf(statustext,"DELETEEXPNO(name, pexpno, disk, user)");
		Show_status(statustext);
	}
	DELETEEXPNO(name, pexpno, disk, user)

/* Set status parameter TD according to constructed ser file length */
	REXPNO(nexpno)
	SETCURDATA
	STOREPARS("TD",2*npoints*si1) 																					/*status: zs, si1*/
	STOREPAR1S("TD",si2) 																							/*status: noe, si2*/
	STOREPAR("TD",2*npoints*si1) 																					/*acq: zs, si1*/
	STOREPAR1("TD",si2) 																							/*acq: noe, si2*/
	if(debug == 1){(void)printf("Finished 3D conversion\n");}
	VIEWDATA_SAMEWIN
}


QUIT