/* Edited H2BC program

Nils T. Nyberg, Jens Ø. Duus, Ole W. Sørensen
  J. Am. Chem. Soc. 127, 6154-6155 (2005)
  MRC, 2005 

Echo/Antiecho gradient selection
uses a 3rd order low-pass J-filter for suppression of one-bond correlations

Using decoupler channel 1 for X-pulses
Flags:
 presat	y:   presaturation during relaxation delay
        n:   no presaturation during relaxation delay
 edit   y:   up/down spectrum
        n:   standard H2BC
	
Pulses and Powers:
 pw	90 degrees transmitter
 pwx	90 degrees decoupler 1
 pwlvl	90 degrees transmitter hard pulse power level
 pwxlvl	power level for decoupler high power
 tpwr   transmitter high power level
 satpwr	power level for saturation of solvent signal using transmitter 
        during relaxation delay 
 satdly	duration of presaturation
 satfrq	presaturation offset
 hpdseq composite pulse decoupling (CPD) program for heteronuclear decoupling
        during T (GARP, CHIRP, etc) optimise for 1Jxh
 hpdmf  dmf for decoupling during T
 hpdpwr decoupler power for decoupling during T
 	NB: For optimal results make sure that the sequence used here sets the fine
	attentuation internally! See Pandora's Box manuals for more details. The same
	is true of the sequence used for dseq
 hpdres angle resolution for hpdseq
 dseq etc: optimise for nJxh where n>1. i.e. use a lower power decoupling
           so you can acquire more points

Coupling constants
 Jmin smallest 1J(XH) for LPJF setting
 Jmax largest 1J(XH) for LPJF setting
 Jxh optimal 1J(XH) for editing (set to near the average 1J(XH))

Delays
 BigT Constant time evolution period between 18-25 msec

Gradient settings
 gt1     length of the gradients used ca. 1msec
 gzlvl1  max size of gradients.
 gstab  gradient recovery delay

Other settings
 array='edit,phase' (for processing with VNMR )

Processing:
 A: wft2d(1,0,-1,0,1,0,-1,0,0,-1,0,-1,0,-1,0,-1)
 B: wft2d(1,0,-1,0,-1,0,1,0,0,-1,0,-1,0,1,0,1)

Version History:
 Version 1: Initial release 08/02/06 by Dr. A.J. Benie

Testing Information:
vnmr61C:Inova 500:Solaris acquisition host (Works OK)
	edit line 160 to change psg_abort(1) to abort(1)
vnmrJ1.1D:Inova 800:Solaris acquisition host (Works OK)
vnmrJ1.1D:n/a:Linux data station (compiles and displays OK)

Disclaimer:
 This pulse program has been tested and found to be working on our
 spectrometers, but adjustments in pulse program syntax and/or
 parameter sets may be necessary for other hardware configurations
 or systems. Please pay attention to the requirements of your system,
 when implementing this pulse program. The Carlsberg Laboratory or
 the Instrument Center does not make any warranties with respect
 to this software. In no event will Carlsberg  Laboratory or the
 Instrument Center be liable for any loss, damages of any nature
 arising or relating to the use or performance of this software.

*/

#include <standard.h>

pulsesequence()
{

 /* Variables */
 double	dectest, tau1, tau2, tau3, t1_d,
        t1ini, tau, epsilon,
 	heatadja, heatadjb, heatadjt1,
	heatadjt2, heatadjd1, heatadjd2,
	decstartd, decstopd, gradtime, gradtimerec,
	d2max, maxni, minBigT, BigTtest,
	jgrad1,jgrad2, jgrad3,
        BigT=getval("BigT"),
        tpwr=getval("tpwr"),
	pwxlvl=getval("pwxlvl"),
	gzlvl1=getval("gzlvl1"),
	gt1=getval("gt1"),
        gstab=getval("gstab"),
	Jmin=getval("Jmin"),
	Jmax=getval("Jmax"),
	Jxh=getval("Jxh"),
	sw1=getval("sw1"),
	dpwr=getval("dpwr"),
	hpdmf=getval("hpdmf"),
        hpdres=getval("hpdres"),
	hpdpwr=getval("hpdpwr");

 int	abortflag=0, phase=(int)(getval("phase") + 0.5);
 char 	presat[MAXSTR],edit[MAXSTR],hpdseq[MAXSTR],teststr[MAXSTR];

 /* Initialize Variables */
 getstr("presat",presat);
 getstr("hpdseq",hpdseq);
 getstr("edit",edit);

 if((Jmin==0.0)||(Jmax==0.0))
 {
  Jmin=125;
  Jmax=165;
 }

 tau1=1/(2*(Jmin + 0.070*(Jmax - Jmin)));
 tau2=1/(Jmin + Jmax);
 tau3=1/(2*(Jmax - 0.070*(Jmax - Jmin)));
 tau=1/(2*Jxh);
 decstartd=POWER_DELAY + PRG_START_DELAY;
 decstopd=POWER_DELAY + PRG_STOP_DELAY;
 gradtime=gt1 + 2*GRADIENT_DELAY;
 gradtimerec=gradtime + gstab;
 epsilon=pw;

 /* check that the gradients or the gradient recovery time are not too long */
 if(gradtimerec>=tau2)
 {
  printf("Error the gradient length plus recovery delay is too long at %f \n",
   gradtimerec);
  printf("Gradient duration gt1 is %f s and the recovery delay gstab is %f \n",
   gt1,gstab);
  abortflag=1;
 }

 /* setup t1: set d2 to zero for first increment(except for in setup exps) */
 if((ix==1)&&(ni>1)) d2=0.0;
 t1ini=rof1 + 1e-6;
 t1_d=d2 + 2.0*t1ini;

 /* check that ni and BigT are OK */
 d2max=t1ini+ni/sw1;
 BigTtest=BigT/2 - 2*tau - epsilon - pw - rof1 - t1ini;
 maxni=BigTtest/(1/(2*sw1));
 minBigT=(d2max*0.5 + epsilon +2*tau + pw + rof1 + t1ini)*2;
 if((maxni<ni)||(BigT<minBigT))
 {
  maxni=floor(maxni);
  printf("Error maximum ni is %.0f for current BigT (%f) \n",maxni,BigT);
  printf("minimum BigT for current ni (%.0f) is %f \n",ni,minBigT);
  abortflag=1;
 }

 /* setup decoupling for both BigT (high power) and acqusition (low power) */
 strcpy(teststr,"nny");
 dectest=strcmp(teststr,dm);
 if(dectest!=0.0)
 {
  printf("Incorrect decoupler flags for dm (%s)! Should be 'nny'  \n",dm);
  abortflag=1;
 }

 if(((satdly>d1)||(satpwr>30)||(satdly>2.0))&&(presat[A]=='y'))
 {
  if(satdly>d1)
   printf("Presaturation time, satdly (%f s) must be less than d1 (%f s)\n",
    satdly,d1);
  if(satpwr>30)
   printf("presaturation power (satpwr) is too high with %f dB \n",satpwr);
  if(satdly>2.0)
   printf("presaturation delay (satdly) is too long with %f s \n",satdly);
   /* abort sequence this should be the only abort in the sequence */
  abortflag=1; 
 }

 if(abortflag==1) psg_abort(1);

 /* Sample heating compenstation setup to account for differences in the up/down
 and standard sequences, and also the changes due to t1 incrementation.
 This is particularly important when using very different decoupling powers for
 hpdseq and lpdseq */
 heatadja=2*gradtimerec + 2*pw;
 heatadjb=2*tau2;
 heatadjt1=20e-6;
 heatadjt2=heatadjb - heatadja + 20e-6;
 heatadjd1=heatadjb - heatadja + 20e-6 + d2max;
 heatadjd2=20e-6 + d2max;

 /* Phase Cycle */
 hlv(ct,v3);
 hlv(v3,v3);
 add(ct,two,v2);
 hlv(v2,v2);
 hlv(v2,v2);
 dbl(v2,v2);
 mod4(v2,v2);
 mod4(ct,oph);
 dbl(oph,oph);
 add(ct,one,v1);
 hlv(v1,v1);
 dbl(v1,v1);
 mod4(v1,v1);
 hlv(v3,v4);
 hlv(v4,v4);
 dbl(v4,v4);
 /* Which translates to */
 /* oph = 0 2 0 2 */
 /* v1 = 0 2 2 0 */
 /* v2 = 0 0 2 2 2 2 0 0 */
 /* v3 = [ 0 1 2 3 ]4 */
 /* v4 = [ 0 2 ]16 */
 /* Pulse Sequence */

 /* Calculate modifications to phases for States-TPPI acquisition */

 mod2(id2,v12);
 dbl(v12,v12);
 add(v2,v12,v2);
 add(oph,v12,oph);

 /* determine the J filter gradients */
 if(edit[A]=='n')
 {
  jgrad1=0.07545*gzlvl1;
  jgrad2=0.17605*gzlvl1;
 }
 else
 {
  jgrad1=0.1509*gzlvl1;
  jgrad2=0.3521*gzlvl1;
 }
 jgrad3=0.7545*gzlvl1;

   
 status(A); /* relaxation delay + decoupling heating compenstation */

 if(presat[A]=='y')
 {
  delay(d1-satdly); /* if doing presat shorten d1 by the presat ammount */
  obsoffset(satfrq);
  delay(4e-6);
  obspower(satpwr);
  rgpulse(satdly - 2.0*POWER_DELAY - rof1 - 8e-6,zero,rof1,0.0);
  obspower(tpwr);
  obsoffset(tof);
  delay(4e-6);
 }
 else
 {
  delay(d1);
 }

 status(B); /* main body of the pulseprogram */

 if(edit[A]=='n') delay(heatadjd1 - t1_d);
 else delay(heatadjd2 - t1_d);
 decpower(hpdpwr);
 decprgon(hpdseq,1.0/hpdmf,hpdres);
 decon();
 if(edit[A]=='n') delay(heatadjt1 + t1_d);  
 else delay(heatadjt2 + t1_d);	
 rgpulse(pw,zero,rof1,rof1);
 delay(BigT/2 - t1_d/2.0 - 3*pwx - pw - 2*tau - epsilon - rof1 - t1ini);
 if(edit[A]=='n') delay(tau + 2*pwx + epsilon + t1ini);
 decoff();
 decprgoff();
 decpower(pwxlvl);
 delay(tau - decstopd - rof1);
 decrgpulse(pwx,v2,rof1,0.0);
 if(edit[A]=='y')
 {
  delay(tau + epsilon - rof1 + t1ini);
  decrgpulse(2*pwx,v3,rof1,0.0);
 }	
 delay(t1_d/2.0 - rof1);
 rgpulse(2.0*pw,zero,rof1,0.0);
 delay(t1_d/2.0);
 if(edit[A]=='y') delay(tau - epsilon - gradtimerec - t1ini);
 if(phase==1)
 {
  zgradpulse(gzlvl1,gt1);
 }
 else
 {
  zgradpulse(gzlvl1*-1,gt1);
 }
 if(edit[A]=='n')decrgpulse(2*pwx,v3,gstab,0.0);
 else decrgpulse(pwx,v1,gstab,0.0);
 if(phase==1)
 {
  zgradpulse(gzlvl1*-1,gt1);
 }
 else
 {
  zgradpulse(gzlvl1,gt1);
 }
 if(edit[A]=='n')
 {
  delay(pw*2 + 2*t1ini);
  decrgpulse(pwx,v1,gstab,0.0);
  delay(gradtimerec);
 }
 else delay(gstab);
 delay(tau - gradtimerec - decstartd - decstopd - rof1);
 decpower(hpdpwr);
 decprgon(hpdseq,1.0/hpdmf,hpdres);
 decon();   
 if(edit[A]=='n')
  delay(BigT/2.0 - t1_d/2 - tau - 2*gradtimerec - 3*pw - 3*pwx - 2*t1ini);
 else delay(BigT/2.0 - t1_d/2 - 2*tau - pwx + epsilon - pw + t1ini);
 decoff();
 decprgoff();
 decpower(pwxlvl);

 simpulse(pw,pwx,one,v4,rof1,0.0); /* start of J filter */
 zgradpulse(jgrad1,gt1);
 delay(tau3 - rof1 - pwx/2 - gradtime);
 decrgpulse(pwx,zero,rof1,0.0);
 zgradpulse(jgrad2,gt1);
 delay(tau2 - rof1 - pwx/2 - gradtime);
 simpulse(2.0*pw,pwx,zero,zero,rof1,rof2);
 zgradpulse(jgrad3,gt1); /* end of J filter */
 delay(tau1 - rof2 - gradtime - decstartd);
 /* switch to low-power decoupling, optimise for nJ(C,H)  */
 decpower(dpwr); 
 
 status(C); /* pre-acquistion delay and acquisition */
 
 delay(tau2+tau3-tau1);

}