/* ghmbced.c "H1-C13 Multiplicity edited HMBC"
Version 2
This version uses two adiabatic pi pulses to perform heteronuclear decoupling
during the inital tau period of the standard HMBC sequence

15th Dec. 2005
Written by Dr. A.J. Benie

Echo/Antiecho gradient selection
uses a low-pass J-filter for suppression of one-bond correlations
(default is a third order filter)

Pulse Program Library
Carlsberg Research Center

Copyright (c), 2005, Carlsberg Research Center, All rights reserved
*/


/*
SETUP INFORMATION
Using decoupler channel 1 for X-pulses
Flags:
lpjf   		Order of low pass J-filter (LPJF)
	1: 	1st order LPJF
	2:	2nd order LPJF
	3:	3rd order LPJF (default)
presat	y:	presaturation during relaxation delay
	n:	no presaturation during relaxation delay
edit    y:	up/down selection
	n:	standard HMBC
altgrad y:	alternative gradient scheme
	n:	standard gradient scheme
MRC	y:	experiment schemes as in Nyberg, N.T. & Sørensen O.W.,
		 Magn. Reson. Chem.,
	n:	experiment schemes as in 

Pulses, Powers & Frequencies:
pw	90 degrees transmitter
pwx	90 degrees decoupler 1
pwlvl	90 degrees transmitter hard pulse power level
pwxlvl	power level for decoupler 1 high power
tpwr    transmitter high power level
satpwr	power level for saturation of solvent signal using transmitter 
	during relaxation delay 
satdly	presaturation delay duration (use d1 for total relaxation delay )
satfrq	presaturation offset
hpdseq  composite pulse decoupling program for heteronuclear decoupling (GARP, CHIRP, etc)
hpdpwr	Power level for composite pulse decoupling
hpdmf	dmf for cpd.
hpdres	angle resolution for cpd.

Coupling constants
Jmin	smallest 1J(XH) for LPJF setting
Jmax	largest 1J(XH) for LPJF setting
Jxh	Optimum 1J(XH) for editing of spectra

Delays
DELTA 	Delay for the evolution of long range coupling constants

Gradient settings
gt1 	length of the gradients used ca. 1msec
gzlvl1  max size of the low pass filter gradients.
gzlvl2  gradient for Echo/AntiEcho selection


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)

Other:
This pulseprogram will NOT work on a GEMINI 2000
*/

#include <standard.h>

/* Subroutines here */
LPJF(pass,jmin,jmax,het90,rofa,rofb,zgradlen,zgradlvl,zgradrec,tauadjust)
int pass;
double jmin, jmax, tauadjust,
 het90, rofa, rofb,
 zgradlen, zgradlvl, zgradrec; 
{
/* subroutine to provide various low pass filters. 
   The default is a third order low pass filter

   Parameters:
   pass: flag for the low pass filter order
   jmin,jmax: the value of the minimum and maximum 1Jxh constant
   het90,rofa, rofb: the length of the X nucleus 90deg pulse.
   rofa/b correspond to rof1/2
   zgradlen/lvl/rec: z gradient length, level and the recovery delay
   tauadjust: the ammount by which to shorten the inital tau delay to
   compensate for timing differences
*/

 double tau, taumin,taumax;

 if (jmin == 0.0 || jmax == 0.0)
  {
  jmin = 115;
  jmax = 175;
  }

 switch (pass)
  {

  case 1:
  taumax = 1.0/(jmin + jmax);
  delay(taumax - zgradlen - 2*GRADIENT_DELAY - zgradrec - tauadjust);
  zgradpulse(zgradlvl,zgradlen);
  delay(zgradrec);
  decrgpulse(het90,zero,zgradlen,rofb);
  zgradpulse(zgradlvl*-1,zgradlen);
  break;

  case 2:
  zgradlvl=zgradlvl/3.0;
  taumin = 1.0/(2.0*(jmin + 0.146*(jmax - jmin)));
  taumax = 1.0/(2.0*(jmax - 0.146*(jmax - jmin)));
  delay(taumin - zgradlen - 2*GRADIENT_DELAY - zgradrec - tauadjust);
  zgradpulse(zgradlvl*3.0,zgradlen);
  decrgpulse(het90,zero,zgradrec,rofb);
  zgradpulse(zgradlvl*-2.0,zgradlen);
  delay(taumax - zgradlen - 2*GRADIENT_DELAY - rofa - rofb - het90);
  decrgpulse(het90,zero,rofa,rofb);
  zgradpulse(zgradlvl*-1.0,zgradlen);
  break;

  default: /* low pass 3rd order filter DEFAULT */ 
  zgradlvl=zgradlvl/7;
  taumin = 1.0/(2.0*(jmin + 0.070*(jmax - jmin)));
  tau = 1.0/(jmin + jmax); 
  taumax = 1.0/(2.0*(jmax - 0.070*(jmax - jmin)));
  delay(taumin - zgradlen - 2*GRADIENT_DELAY - zgradrec - tauadjust);
  zgradpulse(zgradlvl*7.0,zgradlen);
  decrgpulse(het90,zero,zgradrec,rofb);
  zgradpulse(zgradlvl*-4.0,zgradlen);
  delay(tau - zgradlen - 2*GRADIENT_DELAY - rofb - rofa - het90);
  decrgpulse(het90,zero,rofa,rofb);
  zgradpulse(zgradlvl*-2.0,zgradlen);
  delay(taumax - zgradlen - 2*GRADIENT_DELAY - rofb - rofa - het90);
  decrgpulse(het90,zero,rofa,rofb);
  zgradpulse(zgradlvl*-1.0,zgradlen);
  break;

  }
 } /* end of LPJF subroutine */
 
/* End of subroutines */

/* The pulsesequence starts here */

pulsesequence()
{

 /* Variables */
 double
  decstopd, DELTA=getval("DELTA"),
  Epsilon, Epsilonpwx, gradlen,
  gstab=getval("gstab"),
  gt1=getval("gt1"),
  gzlvl1=getval("gzlvl1"),
  gzlvl2=getval("gzlvl2"),
  hpdmf=getval("hpdmf"),
  hpdpwr=getval("hpdpwr"),
  hpdres=getval("hpdres"),
  Jmax=getval("Jmax"),
  Jmin=getval("Jmin"),
  Jxh=getval("Jxh"),
  t1evol, t1min, tau,
  tpwr=getval("tpwr"),
	
  gH=26.752196, /* gyromagnetic ratio of the detected nucleus */
  gX=6.72828; /* gyromagnetic ratio for the indirectly detected nucleus
                 change as needed e.g. for N15 use 2.712621 (ignore the sign) */

 double grad[5]={0.0,0.0,0.0,0.0,0.0};

 char altgrad[MAXSTR],edit[MAXSTR],hpdseq[MAXSTR],MRC[MAXSTR],presat[MAXSTR];

 int lpjf=(int)(getval("lpjf") + 0.5),	
  phase=(int)(getval("phase") + 0.5);

 /* get character strings */
 getstr("altgrad",altgrad);
 getstr("edit",edit);
 getstr("MRC",MRC);
 getstr("presat",presat);
 getstr("hpdseq",hpdseq);

 /* Initialise Variables */
 decstopd=POWER_DELAY + PRG_STOP_DELAY;
 gradlen=gt1 + 2*GRADIENT_DELAY;
 tau=1/(Jxh*2);

 /* make sure that for the first increment d2 is zero
    (except for in setup experiments) */
 if((ix==0)&&(ni>1)) d2=0.0;

 /* check that rof1 is not too long and set up t1 */
/* if (rof1>2.0e-6) rof1=2.0e-6; */
 t1evol=2*rof1 + 1.0e-6;
 t1evol=0.5*(d2 + t1evol);
 t1min=2*rof1 + 1.0e-6;
 Epsilon=t1min+2*pw;
 Epsilonpwx=Epsilon+2*pwx;
 
 /* Start of Real time maths do not mix with 'normal' maths! */
 /* Phase Cycle calculations
    sets up v14 to be used instead of ct so that steady state scans cycle
    sets up the counter v11 for the up/down spectrum selection */
 sub(ct,ssctr,v14);
 mod2(v14,v11);
 /* changes the phase cycle for the up/down experiment(s)
    depending on which scheme is chosen */
 if((edit[A]=='y')&&(MRC[A]!='y')) hlv(v14,v14); 
 /*the rest of the phase cycled phases*/
 dbl(v14,oph);
 hlv(v14,v2);
 hlv(v2,v3);
 dbl(v2,v2);
 dbl(v3,v13);
 add(v2,v13,v2);
 add(v14,one,v1);
 hlv(v1,v1);
 dbl(v1,v1);
 /* Translation of the above
 v1=0 2 2 0
 v2=0 0 2 2 2 2 0 0
 v3={0}4 {1}4 {2}4 {3}4
 oph=0 2
 v13&v14 are only used to calculate phases 1 to 4 and oph
 v11 is used for selecting up/down spectra it is just a series of the form 0 1..
 */
 
 /* Calculate modifications to phases for States-TPPI acquisition
    id2 is a real time variable that contains the increment number */
 mod2(id2,v12);
 dbl(v12,v12);
 add(v2,v12,v2);
 add(oph,v12,oph);

 if(MRC[A]=='y') assign(zero,v11);

 /* Gradient setup
   define a standard basic gradient strength */
 if(altgrad[A]=='y') gzlvl2=gzlvl2/((gH+gX)-((gH*gX)/(gH-gX)));
 else gzlvl2=gzlvl2/(gH+gX);
 /* define the actual strengths and perform the necessary adjustments for
    Echo/AntiEcho the adjustment for the up down spectra is controlled by
    real time if statements */
 if((edit[A]=='n')&&(altgrad[A]=='n'))
  {
  grad[2]=(gH+gX)*gzlvl2;
  grad[3]=(gH-gX)*gzlvl2*-1;
  if(phase==2)
   {
   grad[2]=(gH-gX)*gzlvl2*-1;
   grad[3]=(gH+gX)*gzlvl2;
   }
  }
 else if((edit[A]=='n')&&(altgrad[A]=='y'))
  {
  grad[0]=gX*gzlvl2*-1;
  grad[2]=((gH+gX)-((gH*gX)/(gH-gX)))*gzlvl2;
  grad[3]=(gH-gX)*gzlvl2*-1;
  if(phase==2)
   {
   grad[2]=(gH-gX)*gzlvl2*-1;
   grad[3]=((gH+gX)-((gH*gX)/(gH-gX)))*gzlvl2;
   }
  }
 else if((edit[A]=='y')&&(altgrad[A]=='n'))
  {
  grad[0]=(gH+gX)*gzlvl2*-1;
  grad[4]=(gH-gX)*gzlvl2*-1;
  grad[1]=gH*gzlvl2;
  grad[2]=grad[1]*-1;
  if(phase==2)
   {
    grad[0]=(gH-gX)*gzlvl2*-1;
    grad[4]=(gH+gX)*gzlvl2*-1;
   }
  }
 else if((edit[A]=='y')&&(altgrad[A]=='y'))
  {
  grad[0]=2*((gH+gX)/gH)*gX*gzlvl2;
  grad[4]=2*((gH-gX)/gH)*gX*gzlvl2;
  grad[1]=gX*gzlvl2*-1;
  grad[2]=grad[1]*-1;
  if(phase==2)
   {
   grad[0]=2*((gH-gX)/gH)*gX*gzlvl2;
   grad[4]=2*((gH+gX)/gH)*gX*gzlvl2;
   }
  }
 else
  {
  if((edit[A]!='y')&&(edit[A]!='n'))
   printf("Error edit must be set to 'y' or 'n' or 'y','n' or 'n','y' \n");
  if((altgrad[A]!='y')&&(altgrad[A]!='n'))
   printf("Error altgrad must be set to either 'y' or 'n' \n");
  psg_abort(1);
  }

 /* check that we are not using too much power for the decoupling */
 if(hpdpwr>48)
 {
  printf("Decoupler power too high with %f dB",hpdpwr);
  psg_abort(1);
 }

 /* pulsesequence starts here with pulses */
 
status(A); 

 if(presat[A]=='y')
  {
  if((satdly>d1)||(satpwr>30)||(satdly>2.0))
   {
   if(satdly>d1)
    printf("Presaturation time, satdly (%f) must be shorter than d1 (%f) \n",
     satdly,d1);
   if(satpwr>30)
    printf("Presaturation power, satpwr is too high with %f (max = 30 dB) \n",
     satpwr);
   if(satdly>2.0)
    printf("Presaturation delay, satdly is too long with %f (max = 2 s) \n",
     satdly);
   psg_abort(1);
   }
  else
   {
   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);

 if(edit[A]=='n') delay(tau + Epsilonpwx);
 decpower(hpdpwr);
 decprgon(hpdseq,1.0/hpdmf,hpdres);
 decon();
 if(edit[A]=='y') delay(tau + Epsilonpwx); /* heating compensation */
 rgpulse(pw,zero,rof1,0.0);
 if(edit[A]=='n') delay(tau + Epsilonpwx);
 decoff();
 decprgoff();
 decpower(pwxlvl);
 if(edit[A]=='y')
 {
  ifzero(v11);
  delay(Epsilonpwx);
  endif(v11);
 }
 LPJF(lpjf,Jmin,Jmax,pwx,rof1,0.0,gt1,gzlvl1,gstab,decstopd);
 delay(DELTA - gradlen*2 - gstab);
 if((edit[A]=='n')&&(altgrad[A]=='y')) zgradpulse(grad[0],gt1);
 else if(edit[A]=='y')
  {
  ifzero(v11);
   zgradpulse(grad[0],gt1);
  elsenz(v11);
   zgradpulse(grad[4],gt1);
  endif(v11);
  }
 else delay(gradlen);
 decrgpulse(pwx,v2,gstab,0.0);
 if(edit[A]=='y')
  {
  zgradpulse(grad[1],gt1);
  delay(tau - gradlen - rof1);
  ifzero(v11);
   delay(rof1);
  elsenz(v11);
   delay(Epsilon);
   decrgpulse(pwx*2,v3,rof1,0.0);
  endif(v11);
  }
 delay(t1evol - rof1);
 rgpulse(pw*2,zero,rof1,0.0);
 delay(t1evol - rof1);
 if(edit[A]=='n') delay(rof1);
 if(edit[A]=='y')
  {
  ifzero(v11);
   decrgpulse(pwx*2,v3,rof1,0.0); 
  elsenz(v11);
   delay(rof1);
  endif(v11);
  }
 if((edit[A]=='y')&&(altgrad[A]=='y')) zgradpulse(grad[2],gt1); 
 delay(tau*0.5 - gradlen - gstab);
 if(edit[A]=='y') delay(tau*0.5);
 if(edit[A]=='y')
  {
  ifzero(v11);
   delay(Epsilon);
  endif(v11);
  }
 if((edit[A]=='n')||(altgrad[A]=='n')) zgradpulse(grad[2],gt1); 
 if(edit[A]=='n')
  {
  decrgpulse(pwx*2,v3,gstab,0.0);
  zgradpulse(grad[3],gt1);
  delay(tau*0.5 - gradlen - gstab + Epsilon);
  }
 decrgpulse(pwx,v1,gstab,rof2);
 if(altgrad[A]=='n')
  {
  if(edit[A]=='n') delay(gradlen);
  else
   {
   ifzero(v11);
    zgradpulse(grad[0]*-1,gt1);
   elsenz(v11);
    zgradpulse(grad[4]*-1,gt1);
   endif(v11);
   }
  delay(gstab);
  }
 if(edit[A]=='y')
  {
  ifzero(v11);
  elsenz(v11);
   delay(Epsilonpwx);
  endif(v11);
  }
	
 status(C); /* acquistiion */

} /* end of pulseprogram */