/*Pulse sequence code for the acquisition of:
	1) 1H decoupled 19F spectra
	2) 1H decoupled, 13C filtered and decoupled, 19F spectra using selective 180 degree pulses on 19F
	3) 1H decoupled, 13C edited, 19F spectra using selective 180 degree pulses on 19F */


/*-------------------------------------------------------------------
	Developed by Methodology NMR Group
	School of Chemistry
	University of Manchester
	United Kingdom
	October 2016 
---------------------------------------------------------------------*/


/*User's Guide for experimental setup

	1.	bilevel = 'y' selects the application of the 1H adiabatic bi-level decoupling 
		Note that setting dm=’nnn’ is required to avoid conflicts with default decoupling control. The decoupling pulse shapes are automatically created through the pulse sequence. Execute go(‘check’) before running the experiment, when running the experiment for the first time. 
		The parameters that may need optimisation are bandwidth, pp, pplvl, and npoints. The parameter npoints together with the spectral width (sw) determines the duration of the WURST pulses. The length of the high level decoupling pulse is npoints/sw, and the low level pulse will be 2*npoints/sw. The reference 90 pulse length and power levels used by Pbox are set by the parameters pp and pplvl. 
		The macro kp_readbilevelshapes is essential to run the pulse sequence. 
	2.	satsupp = 'n' selects the acquisition of a 1D 19F 
	3.	satsupp = 'n' and bilevel='y' selects the acquisition of a 1D 1H decoupled 19F spectrum 
	4.	satsupp = 'y' selects the filtration of 13C satellites
					Two options within satsupp = 'y' are:
					jfilter = 'n' the short version of the pulse sequence is applied. In this case a 13C decoupling, aiming to suppress only the long range 13C satellites, should be applied using dm2=’nny’ and setting the relevant parameters: dres2, dpwr2, dmf2, dseq2, dmmp=’ccp’. The j1xh parameter should be set to a value close to 1/2(1JFC).
					jfilter = 'y' the extended version of the pulse sequence is applied. In this case a 2D experiment is acquired. No 13C decoupling is applied. The j1xh parameter should be set to a value close to 1/2(1JFC)
	
The parameters preAmpConfig and probeConnect must be set properly to share a single high-band RF amplifier for proton and fluorine pulses. 
In our spectrometer the standard values are probeConnect=’H1 C13 N15’ preAmpConfig=’hln’, which needed to be changed to probeConnect=’F19 H1 C13’ and preAmpConfig=’hhn’.
Other spectrometer configurations may need different settings. Only channels 1 and 2 can be shared for the highband operation or alternatively channel 3 with 4 (but this would be rather unusual). 

*/



/* 
General Parameters

	pw		:	19F 90 degree pulse width
	tpwr    :   19F pulse power
	pwxlvl  :	13C pulse power
	pwx		:	13C 90 degree pulse width
	dpwr	:	decoupler power on the first decoupling channel when bi-level is not used
	dpwr2	:	decoupler power on the second decoupling channel
	j1xh	:	One-bond fluorine-carbon coupling constant
	d1		:	relaxation delay
	d2		:	t1 evolution delay (normally zero)
	gt1		:	duration of the pulsed field gradient
	gzlvl1	:	amplitude of the pulsed field gradient
	gstab	:	gradient recovery delay
	dutyc	: 	sets the time used for decoupling in the time sharing (usually between 0.1 and 0.6) 
	dbcomp  : 	it is only needed when a conventional decoupling is used, dm='nyy', and dutyc is applied during acquisition. Then compensation should be done for the power of the decoupling during evolution where the dutyc is not applied.
	selshapeA:	shape of the selective pulse	
	selpwA	:	duration of the selective pulse 
	selpwrA	:	power of the selective pulse
	
Bi-level parameters
	
	homorof1, homorof2, homorof3: delays for the receiver's switch between on and off (typically 1-2 μs)
	npoints						: npoints together with the spectral width determines the duration of the WURST pulses (see above). 
	bandwidth					: sets the bandwidth of the 1H decoupling
	pp							: sets the 1H reference 90 degree pulse width for the creation of the WURST pulses
	pplvl						: sets the reference power of the 90 degree 1H pulse for the creation of the WURST pulses 
	hidseq, hidpwr, hidmf		: are automatically set and read from the parameters by the pulse sequence using the macro kp_readbilevelshapes. They refer to the shape, the power and the modulation frequency of the high power WURST-2 period. 
	lodseq,hidpwr,lodmf			: are automatically set and read from the parameters by the pulse sequence using the macro kp_readbilevelshapes. They refer to the shape, the power and the modulation frequency of the low power WURST-2 period.
	t1dseq						: are automatically set and read from the parameters by the pulse sequence using the macro kp_readbilevelshapes. They refer to the shape, the power and the modulation frequency of the decoupling applied during the evolution period.

*/

/*

Real time variables

	v1-v6		: reserved for the the bi-level
	v21,v22		: reserved for the bilevel loop counters
	v31-v37		: reserved for phase cycling

Note: using il=’y’ is not compatible with this pulse sequence code, but interleaved 2D acquisition can be done by acquiring a set of 2D experiments and adding them together after acquisition, or alternatively all initval functions must be replaced by the function F_initval. 

*/



#include <standard.h>


static int	ph1[16] = {0,0,0,0, 1,1,1,1, 2,2,2,2, 3,3,3,3},		//v31; 1st 90 obs
		ph2[16] = {0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0},		//v32; zero; mqfilter 180 obs
		ph3[16] = {0,0,0,0, 2,2,2,2, 0,0,0,0, 2,2,2,2},		//v33; jfilter 180 obs
		ph4[16] = {0,2,0,2, 0,2,0,2, 0,2,0,2, 0,2,0,2},		//v34; 1st 90 on carbon
		ph5[16] = {0,0,2,2, 0,0,2,2, 0,0,2,2, 0,0,2,2},		//v35; 2nd 90 on carbon
		ph6[16] = {0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0},		//v36; zero 180 on carbon; in jfilter
		ph7[16] = {0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0},		//v37; zero last 90 on carbon; in jfilter

	      ph8mq[16] = {0,0,0,0, 3,3,3,3, 2,2,2,2, 1,1,1,1},		//oph if jfilter='n'; 
		ph8[16] = {0,0,0,0, 1,1,1,1, 2,2,2,2, 3,3,3,3},		//oph if jfilter='y';
	      ph8sp[16] = {0,1,2,3, 0,1,2,3, 0,1,2,3, 0,1,2,3};		//oph if jfilter='y';  



		
pulsesequence()
{	
double	j1xh   = getval("j1xh"), 
	pwxlvl = getval("pwxlvl"),
	pw     = getval("pw"),
	pwx    = getval("pwx"),
 	selpwrA = getval("selpwrA"),
	selpwA = getval("selpwA"),
	gt1 = getval("gt1"),		
	gzlvl1 = getval("gzlvl1"),
	gstab = getval("gstab"),
	tau    = 1.0/(2.0*(getval("j1xh")));

int 	kpph = getval("kpph");

char	lkgate_flg[MAXSTR], 
	selshapeA[MAXSTR],
	jfilter[MAXSTR];
 
   getstr("jfilter",jfilter); 
   getstr("lkgate_flg",lkgate_flg);

/* variables for bilevel */

double	dutyc = getval("dutyc"),
	homorof1 = getval("homorof1"),
	homorof2 = getval("homorof2"),
	homorof3 = getval("homorof3"),
	dectime = dutyc/sw,
	acqtime = (1.0-dutyc)/sw -homorof1 - homorof2 - homorof3 - POWER_DELAY, 

	npoints = getval("npoints"),	//defines the length of the decoupler pulse; npoints/sw

	                  //jhf = getval("jhf"),
	pw_hidec = 0.01,  //0.2 / jhf / 2.0,	//pulse length in the high level part

        pw_lodec = 0.002,  //0.2 / jhf,		//pulse length in the low level part		

	bandwidth = getval("bandwidth"), 
	pwr_dec = getval("pplvl"), //58
	pw90_dec = getval("pp")*1000000, //51	//19F when tuned to 1H with HFC probe; 


	histepsize = 5.0, lostepsize = 10.0;

        npoints = (double)( (int)(npoints));

char	bilevel[MAXSTR], satsupp[MAXSTR],
	cmd[MAXSTR],
	hidseq[MAXSTR], lodseq[MAXSTR], t1dseq[MAXSTR];

	getstr("bilevel",bilevel);
	getstr("satsupp",satsupp);
	getstr("hidseq",hidseq);
	getstr("lodseq",lodseq);
	getstr("t1dseq",t1dseq);
	getstr("selshapeA",selshapeA);



if ((npoints/sw*1000000 / histepsize) > 1000)
{
 histepsize = histepsize*(npoints/sw*1000000 / histepsize / 1000);
 printf("Default stepsize (5us) was increased to %.2f us\n",histepsize);
}
if ((npoints/sw*1000000 / lostepsize) > 500)
{
 lostepsize = lostepsize*(npoints/sw*1000000 / lostepsize / 500);
 printf("Default stepsize (10us) was increased to %.2f us\n",lostepsize);
}

if (ix==1)
{
printf("Pulse length for bilevel high decoupling was set to %.3f ms\n",npoints/sw*1000.0);
printf("Pulse length for bilevel low decoupling was set to %.3f ms\n",npoints/sw*2.0*1000.0);
}

if (dectime<0.00001) 
{
abort_message("dectime is too short..");
psg_abort(1);
}
if (acqtime<0.000001) 
{
abort_message("dectime is too short..");
psg_abort(1);
}
if ( npoints/sw < 0.001) 
{
printf("Set npoints more than %0.f or change sw..\n",sw*0.001);
abort_message("npoints is too small..");
psg_abort(1);
}
if ( npoints/sw*20 > at) 
{
abort_message("Reduce npoints or increase at..");
psg_abort(1);
}


if (ix==1)
{
//make decoupler shapes
	sprintf(cmd,"Pbox kp_highdec.DEC -u %s -w \"wurst2i %.1f/%.7f\" -sucyc t5 -s %.1f -dcyc %.2f -p %.0f -l %.1f", userdir, bandwidth, npoints/sw , histepsize, dutyc, pwr_dec, pw90_dec);
	system(cmd);  
	sprintf(cmd,"Pbox kp_lowdec.DEC -u %s -w \"WURST2 %.1f/%.7f\" -sucyc m16 -s %.1f -dcyc %.2f -p %.0f -l %.1f", userdir, bandwidth, npoints/sw*2.0, lostepsize, dutyc, pwr_dec, pw90_dec);
	system(cmd);
	sprintf(cmd,"Pbox kp_t1dec.DEC -u %s -w \"WURST2 %.1f/%.7f\" -sucyc m16 -s %.1f -p %.0f -l %.1f", userdir, bandwidth, npoints/sw*2.0, lostepsize, pwr_dec, pw90_dec);
	system(cmd);
   putCmd("kp_readbilevelshapes\n");
}


double	hidpwr = getval("hidpwr"),
	hidmf = getval("hidmf"),
	hidres = getval("hidres"),
	lodpwr = getval("lodpwr"),
	lodmf = getval("lodmf"),
	lodres = getval("lodres"),
	t1dpwr = getval("t1dpwr"),
	t1dmf = getval("t1dmf"),
	t1dres = getval("t1dres");


if ((hidpwr>56) || (lodpwr>50)) 
{
abort_message("Set npoints is too small; too much power on decoupler. Increase npoints or dutyc..");
psg_abort(1);
}



//phase cycle
//v5: low dec cycle
//v1: high dec cycle
	assign(two,v3);
	mult(v3,v3,v3); mult(v3,v3,v3);
	assign(ct,v1);
	modn(v1,v3,v1);
	assign(one,v2);	add(v1,v2,v1);
	initval(npoints,v6); mult(v6,v1,v1);
	sub(v3,v1,v5);
	initval(np/2-16,v4);
	add(v4,v5,v5);
	

  assign(ct,v17);

if (kpph==0)
{
   settable(t1,16,ph1);
   settable(t2,16,ph2);
   settable(t3,16,ph3);
   settable(t4,16,ph4);
   settable(t5,16,ph5);
   settable(t6,16,ph6);
   settable(t7,16,ph7);

	if (satsupp[0]=='n')
	{
	settable(t8,16,ph8sp);
	}
	if (jfilter[0]=='n') 
	{
	settable(t8,16,ph8mq);
	}
	if ((jfilter[0]=='y') && (satsupp[0]=='y'))
	{
	settable(t8,16,ph8);
	}

}
else
{ 
   settable(t1,kpph,ph1);
   settable(t2,kpph,ph2);
   settable(t3,kpph,ph3);
   settable(t4,kpph,ph4);
   settable(t5,kpph,ph5);
   settable(t6,kpph,ph6);
   settable(t7,kpph,ph7);
	if (satsupp[0]=='n')
	{
	settable(t8,kpph,ph8sp);
	}
	if (jfilter[0]=='n') 
	{
	settable(t8,kpph,ph8mq);
	}
	if ((jfilter[0]=='y') && (satsupp[0]=='y'))
	{
	settable(t8,kpph,ph8);
	}

} 
  getelem(t1, v17, v31);
  getelem(t2, v17, v32);
  getelem(t3, v17, v33);
  getelem(t4, v17, v34);
  getelem(t5, v17, v35);
  getelem(t6, v17, v36);
  getelem(t7, v17, v37);
  getelem(t8, v17, oph);


  

/* equilibrium period */

   status(A);

   obspower(tpwr);
   decpower(dpwr - getval("dbcomp") );
   dec2power(pwxlvl);

   delay(0.02);
   if ( (lkgate_flg[0] == 'y') || (lkgate_flg[0] == 'k') )  lk_sample(); // turn on lock sampling
   delay(d1);
   if ( (lkgate_flg[0] == 'y') || (lkgate_flg[0] == 'k') )  lk_hold(); // turn off lock sampling
   delay(0.02);


if (satsupp[0]=='y')
{

  if  (jfilter[0]=='n')
  {   
status(B);
  delay(0.05);
  rgpulse(pw,v31,rof1,rof1);
   
   obspower(selpwrA);
   delay(tau-gt1-gstab);
   zgradpulse(gzlvl1,gt1);
   delay(gstab);
   shaped_pulse(selshapeA,selpwA,v32,rof1,rof1);
   dec2rgpulse(pwx,v34,rof1,rof1);
   obspower(selpwrA);
   zgradpulse(gzlvl1,gt1);
   delay(gstab);
   delay(tau-pwx-gt1-gstab);
   dec2rgpulse(pwx,v35,rof1,rof2);
  } 
  else
  {
	if (bilevel[0]=='y')
	{  
	  decpower(t1dpwr);
	}
status(B);
	if (bilevel[0]=='y')
	{  
	  decprgon(t1dseq, 1.0/t1dmf, t1dres);
	}
   delay(0.05);
   rgpulse(pw,v31,rof1,rof1);
   
   obspower(selpwrA);
   delay(tau-gt1-gstab);
   zgradpulse(gzlvl1,gt1);
   delay(gstab);
   shaped_pulse(selshapeA,selpwA,v32,rof1,rof1);
   dec2rgpulse(pwx,v34,rof1,rof1);
   obspower(selpwrA);
   zgradpulse(gzlvl1,gt1);
   delay(gstab);
   delay(tau-pwx-gt1-gstab);
   dec2rgpulse(pwx,v35,rof1,rof1);
   obspower(selpwrA);
   delay(d2/2.0);
   shaped_pulse(selshapeA,selpwA,v33,rof1,rof1);
   dec2rgpulse(2.0*pwx,v36,rof1,rof1);
   delay(d2/2.0);
   obspower(selpwrA);
   dec2rgpulse(pwx,v37,rof1,rof2);

	if (bilevel[0]=='y')
	{  
	   decoff(); decpower(-16.0); decblank();
	   decprgoff();
	}
  }
//end of satsupp flag
}
else
{

status(B);
delay(d2);
   rgpulse(pw,oph,rof1,rof2);
}




if (bilevel[0]=='y')
{

//extra delay can be essential for Inova; Hoult
//delay(1.0/(getval("fb")*1.3))



   /* --- observe period --- */

dec2power(dpwr2);
status(C);


//high level part
	decpower(hidpwr);
decprgon(hidseq, 1.0/hidmf, hidres);

	startacq(alfa);
  loop(v1,v21);
	decoff(); decpower(-16.0); decblank();
	delay(homorof3);
	rcvron();  

	delay(homorof1);
	acquire(2.0, acqtime);
	rcvroff(); 
	delay(homorof2);	
	decpower(hidpwr);	decon(); decunblank();
	delay(dectime); 
  endloop(v21);
	decoff(); decpower(-16.0); decblank();
decprgoff();

//low level part
	decpower(lodpwr);
decprgon(lodseq, 1.0/lodmf, lodres);


  loop(v5,v22);
	decoff(); decpower(-16.0); decblank();
	delay(homorof3);
	rcvron();  

	delay(homorof1);
	acquire(2.0, acqtime);
	rcvroff(); 
	delay(homorof2);	
	decpower(lodpwr);	decon(); decunblank();
	delay(dectime); 
  endloop(v22);


	decoff(); decpower(-16.0); decblank();
	delay(homorof3);
decprgoff();

 if (lkgate_flg[0] == 'y') lk_sample();

}
else
{
  /* --- observe period --- */

decpower(dpwr);
dec2power(dpwr2);
status(C);

 if (lkgate_flg[0] == 'y') lk_sample();

}


   

}