/*** ^^A -*-C++-*- ******************************************************/
/*  pshift_no_convdta      06 August 2014                               */
/************************************************************************/
/*  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. Conversion of the data from digital to analog format */
/*  is circumvented in this version by carrying on the group delay at */
/*  the beginning of the FID from the first time increment. */
/*  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)*/
/* 08AUG2014 Circumvent the need to convert to "analog" data by reading */
/*           the data left shift required from GRPDLY. This only works  */
/*           properly, if this parameter is well defined, either in the */
/*           acqus file, or in the approximate lists available for      */
/*           systems using DSPFVS 10, 11 or 12. (LK)                    */
/* 23SEP2014 Replaced use of malloc with calloc in version detection(RA)*/
/************************************************************************/





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

int row[MAXSIZE], temp[MAXSIZE], temp2[256];
int pmod, si1, si2, si3, si12, si13, si12c, drppts, npoints, nzero=0, i4, debug, download, infomsg, twoFiveSixCorr,chunk1;
int aexpno=9999, pexpno=9998, nexpno=expno+1000;
double sw1, sw2, sw3;
double sfo1, sfo2, sfo3;
double grpdly, decim; /*, dz;     Added additional double variables needed.
														grpdly is the Group delay status parameter read out from acqus. It contains
															1) The width of the group delay is given as floor(grpdly) as multiples of the complex dwell time point length
															2) the remaining 1st order phase correction that has to be applied to the raw data after circular shifting as the broken number after the point as multiples of 360deg
															This parameter may be well defined in the acquisition status parameter, for recent hardware. However, for older analog-to-digital converters, it is necessary to read them from a list, which is defined in the 
														decim is the Oversampling dwell time status parameter read out from acqus.
														lk20140806  */
int dspfvs, pnts_in_grpdly, intital_expno;							/*	Added additional integer variables needed.
														dspfvs is the DSP firmware version status parameter read out from acqus.
														pnts_in_grpdly is the number of data points of the group delay (=2*(floor(grpdly)+1)).
														initial_expno is the experiment number of the source data set (used in 3D conversion).
														lk20140806	*/
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=0;
/*download = 1 to display pop up, 0 to skip*/
download=0;
/*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");
	}
} 

/* store current experiment number to intital_expno. lk20140808 */
intital_expno=expno;


/*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*) calloc(6,1);
  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==4){(void)printf("Topspin Version: %i Detected Successfully.\n",versn);}
  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);}   
  }
  
  //if(versn==4){(void)printf("Topspin Version: %i Detected Successfully.\n\nIf this script fails use au program sertoint.ptg to convert dataset to integers from floats then run pshift on the created dataset",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;
	}

/* Determine the number of points within group delay of the digital filter */
FETCHPARS("GRPDLY", &grpdly)
if(debug == 1){(void)printf("Read out: grpdly = %f\n", grpdly);} 
if(grpdly<0)
{
	FETCHPARS("DECIM", &decim)
	FETCHPARS("DSPFVS", &dspfvs)
	if(debug == 1){(void)printf("Group delay not well defined in acquisition status parameters. Determine length of group delay according to values listed in grpdly_from_list\nRead out: decim = %f, dspfvs = %d\n", decim, dspfvs);} 
	
	double grpdly_from_list(double *decim, int *dspfvs)
	/* Reads the group delay resulting from the digital filter settings from a list. */
	/* This should be used on hardware, where the variable grpdly is not properly set by the hardware (grpdly=-1) */
	/* For further information see the comment at the end of the macro.  */
	/* If the grpldy corresponding to the input values of decim and dspfvs is not defined in this function, -1 is returned. */
	{
		int vers=*dspfvs, dcim=*decim;
		
		switch(vers)
		{
			case 10:
				switch(dcim)
				{
					case 2: return 44.7500;
					case 3: return 33.5000;
					case 4: return 66.6250;
					case 6: return 59.0833;
					case 8: return 68.5625;
					case 12: return 60.3750;
					case 16: return 69.5313;
					case 24: return 61.0208;
					case 32: return 70.0156;
					case 48: return 61.3438;
					case 64: return 70.2578;
					case 96: return 61.5052;
					case 128: return 70.3789;
					case 192: return 61.5859;
					case 256: return 70.4395;
					case 384: return 61.6263;
					case 512: return 70.4697;
					case 768: return 61.6465;
					case 1024: return 70.4849;
					case 1536: return 61.6566;
					case 2048: return 70.4924;
				}
				break;
			case 11:
				switch(dcim)
				{
					case 2: return 46.0000;
					case 3: return 36.5000;
					case 4: return 48.0000;
					case 6: return 50.1667;
					case 8: return 53.2500;
					case 12: return 69.5000;
					case 16: return 72.2500;
					case 24: return 70.1667;
					case 32: return 72.7500;
					case 48: return 70.5000;
					case 64: return 73.0000;
					case 96: return 70.6667;
					case 128: return 72.5000;
					case 192: return 71.3333;
					case 256: return 72.2500;
					case 384: return 71.6667;
					case 512: return 72.1250;
					case 768: return 71.8333;
					case 1024: return 72.0625;
					case 1536: return 71.9167;
					case 2048: return 72.0313;
				}
				break;
			case 12:
				switch(dcim)
				{
					case 2: return 46.311;
					case 3: return 36.530;
					case 4: return 47.870;
					case 6: return 50.229;
					case 8: return 53.289;
					case 12: return 69.551;
					case 16: return 71.600;
					case 24: return 70.184;
					case 32: return 72.138;
					case 48: return 70.528;
					case 64: return 72.348;
					case 96: return 70.700;
					case 128: return 72.524;
				}
				break;
		}
		return -1;
	}		
	
	grpdly=grpdly_from_list(&decim, &dspfvs);
	if(grpdly<0){STOPMSG("Lacking value for the group delay length valid on this hardware setup. Data rearrangement aborted.")}
}
if(debug == 1){(void)printf("value of grpdly: %f\n", grpdly);}
pnts_in_grpdly=2*(floor(grpdly)+1);
if(debug == 1){(void)printf("length of group delay (in dwell time points): %d\n", pnts_in_grpdly);}

/* 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 */
	
/* Skipping CONVDTA statement and subsequent experiment folder change. No need for aexpno after full macro rearrangement?, lk20140806		
		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);
		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);	 

/* Deleted CONVDTA statement writing to aexpno. Deleting this data later on therefore is not necessary any more. lk20140807  */
	/* 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")}
		
	/* Write data points in Group delay to output FID. Writing the number of complex points during the group delay (=floor(grpdly)) plus the first complex point lk20140807 */	
		i2=0;
		while(i2<pnts_in_grpdly)
		{
			temp[i2]=row[i2];	/* Not convinced at this point, what the influence of the droppoints would be. This would write the group delay, then skipping the droppoints and then appending the points after CNST4*2 real data points. Thus a part of the FID is cut out. lk20140807 */
			i2++;
		}
		if(fwrite(temp, sizeof(int), pnts_in_grpdly, fpout) != pnts_in_grpdly) {STOPMSG("writing to outfile failed")}
		if(debug == 1){(void)printf("Wrote %d points for the group delay./n", i2+1);}
		
	/* Write data blocks after Group delay to output FID. lk20140807 */	
	
		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);
				
			/* shifted the source data points by floor(grpdly)+1 complex data points, as data has not been converted to analog-style data at this stage. lk20140807 */	
				if(debug == 1){(void)printf("%d.%d(%d)  ",i1,i2,i1*(si2+twoFiveSixCorr)+i2+drppts +pnts_in_grpdly);}		/* Added the +2*(floor(grpdly)+1) shift of the source data. lk20140807 */
				temp[i2]=row[i1*(si2+twoFiveSixCorr)+i2+drppts +pnts_in_grpdly]; /*added correction for the case that si2 is not a multiple of 256*/	/* Added the +2*(floor(grpdly)+1) shift of the source data. lk20140807 */
			
				}
			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 +pnts_in_grpdly)         /*status*/		/* Increased the number of resulting data points by the length of the group delay plus one complex data point. lk20140807 */
		STOREPAR("TD",2*npoints*si1 +pnts_in_grpdly)          /*acq*/			/* Increased the number of resulting data points by the length of the group delay plus one complex data point. lk20140807 */
		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 */
/* Skipping CONVDTA statement and subsequent experiment folder change. No need for aexpno after full macro rearrangement?, lk20140807		
		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")}
			
			if(i1==0)
			{
			/* For first increment write data points in Group delay to output FID. Writing the number of complex points during the group delay (=floor(grpdly)) plus the first complex point lk20140807 */	
				i2=0;
				while(i2<pnts_in_grpdly)
				{
					temp[i2]=row[i2];	/* Not convinced at this point, what the influence of the droppoints would be. This would write the group delay, then skipping the droppoints and then appending the points after CNST4*2 real data points. Thus a part of the FID is cut out. lk20140807 */
					i2++;
				}
				if(fwrite(temp, sizeof(int), pnts_in_grpdly, fpout) != pnts_in_grpdly) {STOPMSG("writing to outfile failed")}
				if(debug == 1){(void)printf("Wrote %d points for the group delay./n", i2+1);}
			}
			
			for (i2=0; i2<2*npoints; i2++){
				sprintf(statustext,"%d of %d (%d.%d)  ",i1,si1,i1,i2);
				Show_status(statustext);
				/* shifted the source data points by floor(grpdly)+1 complex data points, as data has not been converted to analog-style data at this stage. lk20140807 */	
				if(debug == 1){(void)printf("%d.%d(%d)  ",i1,i2,i1*(si2+twoFiveSixCorr)+i2+drppts +pnts_in_grpdly);}
				temp[i2]=row[i2+drppts +pnts_in_grpdly]; /*added correction for the case that si2 is not a multiple of 256*/ 	/* Added the +2*(floor(grpdly)+1) shift of the source data. lk20140807 */
				}
			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);	

/* Deleted CONVDTA statement writing to aexpno. Deleting this data later on therefore is not necessary any more. lk20140807  */		
	/* 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+pnts_in_grpdly)                                           /*status*/					/* Increased the number of resulting data points by the length of the group delay plus one complex data point. lk20140807 */
		STOREPAR("TD",2*npoints*si1+pnts_in_grpdly)                                               /*acq*/					/* Increased the number of resulting data points by the length of the group delay plus one complex data point. lk20140807 */
		/*
		dz=0;
		STOREPARS("GRPDLY", dz)
		*/
		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/sapphire */
	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: %f\nsw2: %f\nsw3: %f\n",si1,si2,si3,sw1,sw2,sw3);}		/* Changed output format for sw1, sw2 and sw3 from %d to %f. lk20140808 */
	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 */
/* Skipping digital to analog conversion, lk20140808 
		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"));
	/* Jump back to intital_expno instead of the obsolete aexpno. lk20140808 */
	REXPNO(intital_expno)
	/* REXPNO(aexpno)			lk20140808*/
#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(versn==4){
	 	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"));
	/* Jump back to intital_expno instead of the obsolete aexpno. lk20140808 */
	REXPNO(intital_expno)
	/* REXPNO(aexpno)			lk20140808*/
#endif	

	/* Read data directly from initial_expno. lk20140808
	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*/
	
	/* shortened the Zero Padding length by pnts_in_grpdly. lk20140808
	nzero=256-(2*npoints*si1)%256;	*/
	
	
	nzero=256-(pnts_in_grpdly+2*npoints*(si1-1))%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 (sapphire) 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(versn==4){
		    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*/

	if(i4<=((si2/2)+1))
		{	
	chunk1=((npoints/2)+((i4-1)*round_to_int((2*sw3)/sw2)));
		}
  else
  	{	
	chunk1=(npoints-((i4-1)*round_to_int((2*sw3)/sw2)));
		}
	
	if(debug == 1){(void)printf("Chunk1 length is %d.",chunk1);}
				

																		
		/* Add Group Delay of first F1 increment to the beginning of the FID. Again at this position, if droppoints are used, the final FID will miss out these points in between the group delay and the actual data. lk20140808 */
		i2=0;
		while(i2<pnts_in_grpdly)
		{
			temp[i2]=row[i2];	/* Not convinced at this point, what the influence of the droppoints would be. This would write the group delay, then skipping the droppoints and then appending the points after CNST4*2 real data points. Thus a part of the FID is cut out. lk20140807 */
			i2++;
		}
		if(fwrite(temp, sizeof(int), pnts_in_grpdly, fpout) != pnts_in_grpdly) {STOPMSG("writing to outfile failed")}
		if(debug == 1){(void)printf("Wrote %d points for the group delay of F2 increment No %d./n", i2+1, i4+1);}
		
		/* Addition of the data chuncks after group delay, lk20140808 */
		i1=0;																					/*Set up the pointer for the loop over the ZS loop*/
		while(i1 < si1)	{	
		
						
		if(i1==0)	{																		/*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*chunk1; 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);
			  ssleep(0.5);
				if(debug == 1){(void)printf("%d.%d.%d ",i4,i1,i2);}
				/* Shifted source data by the lenth of the groups at each beginning of the individual FIDs. lk20140808
				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*/
				temp[i2]=row[i1*(si3+twoFiveSixCorr)+i2+drppts+pnts_in_grpdly]; 		/*transfer data in defined position to temporary file, fid, added correcdtion if si2 is not a multiple of 256*/  /* Shifted data source by pnts_in_grpdly points. lk20140808 */
			}
			if(fwrite(temp, sizeof(int), 2*chunk1, fpout) != 2*chunk1){STOPMSG("writing to outfile failed")} /*Write the data out and check it was written out ok*/																
					
		}
		
		else if (i1>0 && i1<(si1-1)) 	{																		/*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);}
				/* Shifted source data by the lenth of the groups at each beginning of the individual FIDs. lk20140808
				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*/
				temp[i2]=row[i1*(si3+twoFiveSixCorr)+i2+drppts+pnts_in_grpdly]; 		/*transfer data in defined position to temporary file, fid, added correcdtion if si2 is not a multiple of 256*/  /* Shifted data source by pnts_in_grpdly points. lk20140808 */
			}
			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*/
																					/*increment the pointer*/
		}
		
	  else {																		/*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)-(2*chunk1)); 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);}
			/* Shifted source data by the lenth of the groups at each beginning of the individual FIDs. lk20140808
			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*/
			temp[i2]=row[i1*(si3+twoFiveSixCorr)+i2+drppts+pnts_in_grpdly]; 		/*transfer data in defined position to temporary file, fid, added correcdtion if si2 is not a multiple of 256*/  /* Shifted data source by pnts_in_grpdly points. lk20140808 */
	   	}
  	if(fwrite(temp, sizeof(int), ((2*npoints)-(2*chunk1)), fpout) != ((2*npoints)-(2*chunk1))){STOPMSG("writing to outfile failed")} /*Write the data out and check it was written out ok*/
																							/*increment the pointer*/
			}
		
		i1++;
			}

		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 */
	/* aexpno was not created if CONVDTA was not used. No need to delete it afterwards. lk20140808
	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
	/* Adjusted to also include the group delay. lk20140808
	STOREPARS("TD",2*npoints*si1) */																					/*status: zs, si1*/
	STOREPARS("TD",2*npoints*(si1-1)+pnts_in_grpdly) 																								/*status: zs, si1*/		/* Adjusted to also include the group delay. lk20140808 */
	STOREPAR1S("TD",si2) 																							/*status: sapphire, si2*/
	/* Adjusted to also include the group delay. lk20140808
	STOREPAR("TD",2*npoints*si1) 	*/																				/*acq: zs, si1*/
	STOREPAR("TD",2*npoints*(si1-1)+pnts_in_grpdly) 																					/*acq: zs, si1*/	/* Adjusted to also include the group delay. lk20140808 */
	STOREPAR1("TD",si2) 																							/*acq: sapphire, si2*/
	/*
	dz=0;
	STOREPARS("GRPDLY", dz)
	*/
	if(debug == 1){(void)printf("Finished 3D conversion\n");}
	VIEWDATA_SAMEWIN
}


QUIT



// INFORMATION ABOUT BRUKER DIGITAL FILTER
// For older versions of the Bruker hardware:
//
// A nice piece was found on the internet on how to calculate the number of points
// semi-automatically. Note that currently matNMR doesn't allow for the necessary
// negative-time apodization.
//
//
//    W. M. Westler and F.  Abildgaard
//    July 16, 1996
//
//    The introduction of digital signal processing by Bruker in their DMX
//    consoles also introduced an unusual feature associated with the data. The
//    stored FID no longer starts at its maximum followed by a decay, but is
//    prepended with an increasing signal that starts from zero at the
//    first data point and rises to a maximum after several tens of data points.
//    On transferring this data to a non-Bruker processing program such as FELIX,
//    which is used at NMRFAM, the Fourier transform leads to an unusable spectrum
//    filled with wiggles. Processing the same data with Bruker's Uxnmr
//    program yields a correct spectrum. Bruker has been rather reluctant
//    to describe what tricks are implemented during their processing protocol.
//
//    They suggest the data should be first transformed in Uxnmr and then inverse
//    transformed, along with a GENFID command, before shipping the data to another
//    program. Bruker now supplies a piece of software to convert the digitally
//    filtered data to the equivalent analog form.
//    We find it unfortunate that the vendor has decided to complicate
//    the simple task of Fourier transformation. We find that the procedure
//    suggested by Bruker is cumbersome, and more so, we are irritated since
//    we are forced to use data that has been treated with an unknown procedure.
//    Since we do not know any details of Bruker's digital filtration procedure
//    or the "magic" conversion routine that is used in Uxnmr, we have been forced
//    into observation and speculation. We have found a very simple, empirical
//    procedure that leads to spectra processed in FELIX that are identical,
//    within the noise limits, to spectra processed with Uxnmr. We deposit
//    this information here in the hope that it can be of some
//    use to the general community.
//    The application of a nonrecursive (or recursive) digital filter to time
//    domain data is accomplished by performing a weighted running average of
//    nearby data points. A problem is encountered at the beginning of
//    the data where, due to causality, there are no prior values. The
//    weighted average of the first few points, therefore, must include data
//    from "negative" time. One naive procedure, probably appropriate to NMR
//    data, is to supply values for negative time points is to pad the data with
//    zeros. Adding zeros (or any other data values) to the beginning of
//    the FID, however, shifts the beginning of the time domain data (FID) to
//    a later positive time. It is well known that a shift in the time
//    domain data is equivalent to the application of a frequency-dependent,
//    linear phase shift in the frequency domain. The 1st order phase shift
//    corresponding to a time shift of a single complex dwell point is 360 degrees
//    across the spectral width. The typical number of prepended points
//    found in DMX digitally filtered data is about 60 data points (see below),
//
//    the corresponding 1st order phase correction is ~21,000 degrees.
//    This large linear phase correction can be applied to the transformed data
//    to obtain a normal spectrum. Another, equivalent approach is to time
//    shift the data back to its original position. This results in the need
//    of only a small linear phase shift on the transformed data.
//    There is a question as what to do with the data preceding the actual
//    FID. The prepended data can be simply eliminated with the addition
//    of an equal number of zeros at the end of the FID (left shift). This
//    procedure, however, introduces "frowns" (some have a preference to refer
//    to these as "smiles") at the edge of the spectrum. If the sweep
//    width is fairly wide this does not generally cause a problem. The
//    (proper) alternative is to retain this data by applying a circular left
//    shift of the data, moving the first 60 or so points (see recommendations
//    below) to the end of the FID. This is identical to a Fourier transformation
//    followed by the large linear phase correction mentioned above. The
//    resulting FID is periodic with the last of the data rising to meet the
//    first data point (in the next period). Fourier transform of this
//    data results in an approximately phased spectrum. Further linear
//    phase corrections of up to 180 degrees are necessary. A zero fill applied
//    after a circular shift of the data will cause a discontinuity and thus
//    introduce sinc wiggles on the peaks. The usual correction for DC
//    offset and apodization of the data, if not done correctly, also results
//    in the frowns at the edges of the spectrum.
//
//    In our previous document on Bruker digital filters, we presented deduced
//    rules for calculating the appropriate number of points to be circular left
//    shifted. However, since then, newer versions of hardware (DQD) and software
//    has introduced a new set of values. Depending on the firmware versions
//    (DSPFVS) and the decimation rate (DECIM), the following lookup table will
//    give the circular shift values needed to correct the DMX data. The values
//    of DECIM and DSPFVS can be found in the acqus file in the directory containing
//    the data.
//
//     DECIM           DSPFVS 10       DSPFVS 11      DSPFVS 12
//
//       2              44.7500         46.0000        46.311
//       3              33.5000         36.5000        36.530
//       4              66.6250         48.0000        47.870
//       6              59.0833         50.1667        50.229
//       8              68.5625         53.2500        53.289
//      12              60.3750         69.5000        69.551
//      16              69.5313         72.2500        71.600
//      24              61.0208         70.1667        70.184
//      32              70.0156         72.7500        72.138
//      48              61.3438         70.5000        70.528
//      64              70.2578         73.0000        72.348
//      96              61.5052         70.6667        70.700
//     128              70.3789         72.5000        72.524
//     192              61.5859         71.3333
//     256              70.4395         72.2500
//     384              61.6263         71.6667
//     512              70.4697         72.1250
//     768              61.6465         71.8333
//    1024              70.4849         72.0625
//    1536              61.6566         71.9167
//    2048              70.4924         72.0313
//
//
//    The number of points obtained from the table are usually not integers.  The appropriate procedure is to circular shift (see protocol for details) by the integer obtained from truncation of the obtained value and then the residual 1st order phase shift that needs to be applied can be obtained by multiplying the decimal portion of the calculated number of points by 360.
//
//    For example,
//
//    If DECIM = 32, and DSPFVS = 10,
//    then #points 70.0156
//
//    The circular shift performed on the data should be 70 complex points and the linear
//    phase correction after Fourier transformation is approximately 0.0156*360= 5.62 degrees.
//
//    Protocol:
//
//       1. Circular shift (rotate) the appropriate number of points in the data indicated by
//       the  DECIM parameter. (see above formulae).
//
//       2. After the circular shift, resize the data to the original size minus
//       the number of shifted points. This will leave only the part of the
//       data that looks like an FID. Baseline correction (BC) and/or apodization
//       (EM etc.) should be applied only on this data, otherwise "In come the frowns."
//
//       Since the first part of the data (the points that are shifted) represents
//       negative time, a correct apodization would also multiply the shifted points
//       by a negative time apodization. The data size is now returned to
//       its original size to reincorporate the shifted points. There may
//       still be a discontinuity between the FID portion and the shifted points
//       if thelast point of the FID portion is not at zero. This will cause
//       sinc wiggles in the peaks.
//
//       3. Applying a zero fill to this data will lead to a discontinuity in the data
//       between the rising portion of the shifted points and the zero padding.
//       To circumvent this problem, the shifted points are returned (by circular
//       shift) to the front of the data, the data is zero filled, and then the
//       first points are shifted again to the end of the zero filled data.
//
//       4) The data can now be Fourier transformed and the residual calculated
//       1st order phase correction can be applied.
//
//
//
// For newer versions of the Bruker hardware:
//
//     For firmware versions of 20 <= DSPFVS <= 23 the GRPDLY parameter directly shows the number of
//     points that need to be shifted.
//
//     Thanks to Bruker for supplying this information!
//
//
//
//26-09-'00