/*** ^^A -*-C++-*- ******************************************************/
/*                         pshift4                      16 April 2020   */
/************************************************************************/
/*  Short Description :                                                 */
/*  AU program to reconstruct interferogram pure shift experiment       */
/*  e.g. Zangger-Sterk, band-selective or PSYCHE                        */
/*  apply on a pseudo 2D data                     */
/*  1D pure shift data will be written to 1000+experiment number 			  */
/************************************************************************/
/*  Keywords :                                                      	*/
/*  pureshift, Zangger-Sterk, PSYCHE                                	*/
/************************************************************************/
/*  Description/Usage :       
                                    		*/
/* 	This AU program reconstructs 1D pure shift data by concatenating  */
/*  the chunks acquired in 2D mode							                      */
/*	Only 2D interferogram experiments are supported        	*/
/*  Data processing is done in double format, if input data is int then */
/*   it is automatically converted to double format (ser_int backup file */
/*   is stored in the raw data folder)                                  */
/*	Only tested for TopSpin 4, does not attempt to handle earlier versions*/
/*  The 1D pure shift data is stored in 1000 + current experiment number */
/*   if that is already in use, then the next available experiment number */
/*  will be used. When repeatedly used in the same 2D dataset it may be */ 
/*   called with argument 'y' to overwrite data in 1000+expno instead of */
/*   creating a new experiment number */
/*	a sensible Gaussian weigthing is automatically applied to suppress */
/*  truncation artefacts in the concatenated 1D pure shift data */
/*  data processing parameters set in the raw 2D experiment are ignored */
/*  the pure shift 1D data is processed using gfp (fp can also be used) */
/************************************************************************/ 
/*  Author(s) :                                                 		*/
/*  Name          : Ralph Adams and Peter Kiraly 				                */
/*  Organisation  : University of Manchester, UK                        */
/*  Email         : ralph.adams@manchester.ac.uk                        */
/************************************************************************/

// change MAXSIZE if RAM is limited; min. value is to store TD2 datapoints
#define MAXSIZE 400000

float taq;
int dtypa, kp_datasize, replace=0; 
double row[MAXSIZE], temp[MAXSIZE], temp2[256];
int pmod, si1, si2, si12, si12c, drppts, npoints,  debug, download, twoFiveSixCorr;
int nexpno=expno+999, maxexpno=expno+2000;
double sw1, sw2;
double sfo1, sfo2;
double grpdly; 
int pnts_in_grpdly;
float l4;
char infile[PATH_MAX], outfile[PATH_MAX], statustext[256];
FILE *fpin, *fpout;


GETCURDATA

//debug = 1 for on, 0 for off
debug=0;
/*download = 1 to display pop up, 0 to skip*/
download=0;

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 }
}

if(debug == 1){(void)printf("Debugging on\nsuppress debugging by setting debug=0 in macro\n");}
if(debug == 1)
{
	if(download == 1){
		(void)printf("Download agreement will be shown during processing\n");
	}
} 



// check for data format
FETCHPARS("DTYPA", &dtypa)
if (dtypa==0)
{
 kp_datasize=sizeof(int);
 if (debug == 1) { (void)printf("integer format detected..\n");}
 if (debug == 1) { (void)printf("converting data to double format..\n");}

 if ( (CPR_exec("kp_sertodb i", WAIT_TERM)) != NORM_TERM )
 {
 STOPMSG("Dataset conversion has failed unexpectedly.. ")	
 }
 else
 {
 if (debug == 1) { (void)printf("data has been converted to double format..\n");}
 STOREPARS("DTYPA",2) 
 STOREPAR("DTYPA",2) 
 FETCHPARS("DTYPA", &dtypa)
 kp_datasize=sizeof(double);
 }
}
if (dtypa==1)
{
kp_datasize=sizeof(float);
if(debug == 1){(void)printf("unsupported float format detected..\n");}
ABORT
}
if (dtypa==2)
{
kp_datasize=sizeof(double);
if(debug == 1){(void)printf("double format detected..\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*) 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);}   
  }
  

	#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
  
FETCHPARS("GRPDLY", &grpdly)
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);}

FETCHPAR("PARMODE",&pmod)
if ( (pmod < 1) || (pmod > 1) )
    STOPMSG("Not a 2D dataset.\nThis macro only supports 2D data input..")

if(debug == 1){(void)printf("Starting conversion of 2D raw data\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=128-(si2%128);
	if(twoFiveSixCorr == 128){twoFiveSixCorr = 0;}
		if(debug == 1){(void)printf("Correction Value for si2 not divisible by 128: %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(debug == 1){(void)printf("Slow version of macro is being used - multiple writes to disk rather than using just RAM\n");}
		
		strcpy(infile, ACQUPATH("ser"));
		    		    		    		    
  if (i_argc > 2)
	{
	if (i_argv[2][0] == 'y')	replace=16; 
	}
	if (replace==0)
	{		    		    		    		    		    		    		    		    		    		    		    		    		    		    
		do
     	nexpno++;
    while ( (access(NEWACQUPATH(nexpno, 0), F_OK) == 0) || (nexpno>maxexpno) );
		if(debug == 1) { GETINT("Enter new expno for 1D pureshift data:  ",nexpno) }
	}
	else nexpno=nexpno+1;	
		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 TS4-TS3\n");}  
		strcpy(outfile,NEWACQUPATH(nexpno, "fid"));
		#endif
	
		#ifndef NEWACQUPATH
		if(debug == 1){(void)printf("ifndef TS2\n");}  
		REXPNO(nexpno)
		strcpy(outfile,ACQUPATH("fid"));
		#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=(int)( floor((double)(l4)) );  //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("'Real+Imag' points to drop: %d\n",drppts);}
	
		npoints=(int)(floor(sw2/sw1) );  //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,kp_datasize,si2+twoFiveSixCorr,fpin)!=si2+twoFiveSixCorr){STOPMSG("Read in of ser file failed!\n")}
			
			if(i1==0)
			{
			// For first increment onlt 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 */
	// the droppoints should not be used in the first chunk of data ! this needs changes in the pulse programs
					i2++;
				}
				if(fwrite(temp, kp_datasize, 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);				
				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]; 
			}
			if(fwrite(temp, kp_datasize, 2*npoints, fpout) != 2*npoints) {STOPMSG("writing to outfile failed")}
			i1++;
			if(debug == 1){(void)printf("\n",i1);}		
		}
		fclose(fpout);
		fclose(fpin);	

		
		REXPNO(nexpno)
		SETCURDATA
		STOREPARS("TD",2*npoints*si1+pnts_in_grpdly)  
		STOREPAR("TD",2*npoints*si1+pnts_in_grpdly)   

		if(debug == 1){(void)printf("Finished conversion\n");}

// apply sensible 1D processing parameters
  STOREPAR("SI",2*2*npoints*si1+pnts_in_grpdly)
  CPR_exec("fp", WAIT_TERM);
 	STOREPAR("TDEFF",0);
	STOREPAR("LB",-0.01);	
	FETCHPARS("AQ",&taq);	
	STOREPAR("GB", 3.141593*0.01*(taq*0.5)*(taq*0.5)/(taq*2.0) );
	STOREPAR("WDW",2);	
	
 	if (debug == 1) (void)printf("setgtc_auto was applied..\n");
 
 	CPR_exec("gfp", WAIT_TERM);
 	VIEWDATA_SAMEWIN

QUIT

