/*** ^^A -*-C++-*- ********************************************************/
/*  pshift4f                      25 March 2022  			  */
/**************************************************************************/
/*  Short Description :                                                   */
/*  AU program to reconstruct interferogram pure shift data               */
/*  e.g. Zangger-Sterk, band-selective, PSYCHE or BIRD                    */
/*  apply on a pseudo 2D data                     			  */
/*  1D pure shift data will be written to 1000 + experiment number 	  */
/**************************************************************************/
/*  Keywords :                                                      	  */
/*  pure shift, band-selective, 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          : Manchester Methodology Group 			  */
/*  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  (RA)         */
/* 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)  */
/* 16APR2020 Modified to be compatible with TopSpin 4 (PK)                */
/* 25MAR2022 Modified to read parameters SW_h instead of SW and increased */ 
/*	     SI in 1D processing parameters by a factor of 4 (CM,GAM)     */
/**************************************************************************/

// 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_h",&sw1)
	FETCHPARS("SW_h",&sw2)
	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*8*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