/* hsqc_pfg_sel_500.c        - Varian UNITYplus, UNITY

    This pulse sequence will allow one to perform the following experiment:

    2D HSQC  Heteronuclear Single Quantum Coherence with enhanced sensitivity
             PFG and minimal perturbation of water.
           F1      15N
           F2(acq)  1H (NH)
 
    This sequence uses the standard three channel configuration
         1)  1H             - carrier (tof) @ 4.7 ppm [H2O]
         2) 13C             - carrier (dof) @ 117 ppm [centre of CO and Ca]
         3) 15N             - carrier (dof2)@ 119 ppm [centre of amide 15N]
 
    Set dm = 'nnnn', dmm = 'cccc'
    Set dm2 = 'nnny', dmm2 = 'cccp' [15N decoupling during acquisition]
                      dseq2 = 'waltz16'
 
    Must set phase = 1,2 for States-TPPI acquisition in t1 [15N] with
    gradient selection of 15N magnetization. [The fids must be manipulated
    (add/subtract) with 'grad_sort_nd' program (or equivalent) prior to
    regular processing.]
 
    Flags
        fsat            'y' for presaturation of H2O
        fscuba          'y' to apply scuba pulse after presaturation of H2O
        mess_flg        'y' for Messerlie type purging pulse
        f1180           'y' for 180 deg linear phase correction in F1
                            otherwise 0 deg linear phase correction
        c180_flg        'y' for 15N and 13C labelled samples
 
        Standard Settings
        fsat='n',fscuba='n',mess_flg='n',f1180='n',c180_flg='n'
          

        Written by L. E. Kay July 7 1992

    REF: L. E. Kay, P. Keifer, and T. Saarinen
                                J. Am. Chem. Soc. 114, 10663-10665 (1992).

        Modified by L. E. Kay Dec 15, 1993 to encorporate selective pulses
                    on water 

        REF: Zhang et. al.  J. Biomol. NMR 4, 845 (1994) 

BMRB Pulse Sequence Accession Number: 61
*/

#include <standard.h>
#define PI 3.1416

static double d2_init = 0.0;

static int phi1[1] = {0},
           phi2[1] = {1},
           phi3[8] = {0,0,1,1,2,2,3,3}, 
           phi4[2] = {0,2},
           phi5[1] = {0},
           phi7[1] = {0},
           rec[4] = {0,0,2,2};

pulsesequence()
{
/* DECLARE VARIABLES */

 char        fscuba[MAXSTR],f1180[MAXSTR],fsat[MAXSTR],mess_flg[MAXSTR],c180_flg[MAXSTR],shp_sl[MAXSTR];

 int	     phase,icosel,t1_counter;

 double      hscuba,                /* length of 1/2 scuba delay */
             tauxh,                 /* 1 / 4J(NH)                */
             pwx2,                  /* PW90 for N-nuc            */
             pwx1,                  /* PW90 for C-nuc            */
             tsatpwr,               /* low power level for presat*/
             tpwrmess,              /* power level for Messerlie purge */
             dly_pg,                /* Messerlie purge delay         */
             dhpwr2,                /* power level for N hard pulses */
             dhpwr,                 /* power level for C hard pulses */
             jxh,                   /* coupling for NH           */
	     tau1,	      	    /* t1/2 */
	     sw1,
             pw_sl,                 /* pw for selective pulse at twprsl  */
             tpwrsl,

             gzlvl0,               /* level of grad. */
             gt0,                  /* grad time */
             gzlvl1,               
             gt1,                 
             gzlvl2,
             gt2,
             gzlvl3,               /* grad level 30G/cm */
             gt3,                  /* grad time 1.25ms */
             gzlvl4,               /* grad level ~29G/cm UNITYplus */
             gt4,                  /* grad time 0.125ms */
             gzlvl5,
             gt5,
             BigT,               
             BigT1;


/* LOAD VARIABLES */

  jxh = getval("jxh");
  pwx2 = getval("pwx2");
  pwx1 = getval("pwx1");
  tsatpwr = getval("tsatpwr");
  tpwrmess = getval("tpwrmess");
  dly_pg = getval("dly_pg");
  dhpwr2 = getval("dhpwr2"); 
  dhpwr = getval("dhpwr"); 
  hscuba = getval("hscuba");
  phase = (int) (getval("phase") + 0.5);
  sw1 = getval("sw1");
  BigT = getval("BigT");
  BigT1 = getval("BigT1");
  pw_sl = getval("pw_sl");
  tpwrsl = getval("tpwrsl");

  gt0 = getval("gt0");
  gt1 = getval("gt1");
  gt2 = getval("gt2");
  gt3 = getval("gt3");
  gt4 = getval("gt4");
  gt5 = getval("gt5");
  gzlvl0 = getval("gzlvl0");
  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl4 = getval("gzlvl4");
  gzlvl5 = getval("gzlvl5");

  getstr("fscuba",fscuba); 
  getstr("fsat",fsat); 
  getstr("f1180",f1180); 
  getstr("mess_flg",mess_flg);
  getstr("c180_flg",c180_flg);
  getstr("shp_sl",shp_sl);

/* check validity of parameter range */

    if((dm[A] == 'y' || dm[B] == 'y' || dm[C] == 'y' || dm[D] == 'y'))
	{
	printf("incorrect Dec1 decoupler flags!  ");
	abort(1);
    } 

    if((dm2[A] == 'y' || dm2[B] == 'y' || dm2[C] == 'y'))
	{
	printf("incorrect Dec2 decoupler flags!  ");
	abort(1);
    } 

    if( tsatpwr > 8 )
    {
	printf("tsatpwr too large !!!  ");
	abort(1);
    }

    if ( tpwrmess > 57 )
    {
       printf("tpwrmess too high !!!");
       abort(1);
    }

    if (dly_pg > 0.010 )
    {
       printf("dly_pg too long !! ");
       abort(1);
    }

    if( dpwr > 50 )
    {
	printf("don't fry the probe, dpwr too large!  ");
	abort(1);
    }

    if( dpwr2 > 50 )
    {
	printf("don't fry the probe, dpwr2 too large!  ");
	abort(1);
    }

    if(gt0 > 15.0e-3 || gt1 > 15.0e-3 || gt2 > 15.0e-3 || gt3 > 15.0e-3 || gt4 > 15.0e-3 || gt5 > 15.0e-3)
    {
        printf("gti must be less than 15 ms \n");
        abort(1);
    }

    if(mess_flg[A] == 'y') {
        printf("mess_flg should REALLY be set to n for optimum performance\n");
    }


/* LOAD VARIABLES */

  settable(t1, 1, phi1);
  settable(t2, 1, phi2);
  settable(t3, 8, phi3);
  settable(t4, 2, phi4);
  settable(t5, 1, phi5);
  settable(t7, 1, phi7);
  settable(t6, 4, rec);

/* INITIALIZE VARIABLES */

  tauxh = ((jxh != 0.0) ? 1/(4*(jxh)) : 2.25e-3);

/* Phase incrementation for hypercomplex data */

   if ( phase == 1 )     /* Hypercomplex in t1 */
   {
        tsadd(t5, 2, 4); 
        icosel = 1; /*  and reverse the sign of the gradient  */
   } 
   else icosel = -1;


/* calculate modification to phases based on current t1 values
   to achieve States-TPPI acquisition */
 
 
   if(ix == 1)
      d2_init = d2;
      t1_counter = (int) ( (d2-d2_init)*sw1 + 0.5);

      if(t1_counter %2) {
        tsadd(t7,2,4);
        tsadd(t6,2,4);
      }


/* set up so that get (-90,180) phase corrects in F1 if f1180 flag is y */

   tau1 = d2;
   if(f1180[A] == 'y')  tau1 += ( 1.0/(2.0*sw1) );
   tau1 = tau1/2.0;
   

/* BEGIN ACTUAL PULSE SEQUENCE */

status(A);
   rlpower(tsatpwr,TODEV);            /* Set power for presaturation  */
   rlpower(dhpwr,DODEV);              /* Set decoupler1 power to dhpwr */
   rlpower(dhpwr2,DO2DEV);            /* Set decoupler2 power to dhpwr2 */

/* Presaturation Period */

status(B);

/* option for Messerlie purge */

  if(mess_flg[A] == 'y') {
    rlpower(tpwrmess,TODEV);
    rgpulse(dly_pg,zero,20.0e-6,20.0e-6);
    rgpulse(dly_pg/1.62,one,20.0e-6,20.0e-6);
    rlpower(tsatpwr,TODEV);
 }

 if(fsat[0] == 'y')
{
  txphase(zero);
  rgpulse(d1,zero,20.0e-6,20.0e-6);
  rlpower(tpwr,TODEV);                /* Set power for hard pulses  */

    if (fscuba[0] == 'y')            /* Scuba pulse sequence */
    {  
      hsdelay(hscuba);

      rgpulse(pw,zero,1.0e-6,0.0);	/* 90x180y90x */
      rgpulse(2*pw,one,1.0e-6,0.0);
      rgpulse(pw,zero,1.0e-6,0.0);
 
      txphase(zero);
      delay(hscuba);        
    }
}

else  {
 rlpower(tpwr,TODEV);                /* Set power for hard pulses  */
 delay(d1);
}

status(C);

  rcvroff();
  delay(20.0e-6);

  /* eliminate all magnetization originating on 15N */

  dec2rgpulse(pwx2,zero,0.0,0.0);
  delay(2.0e-6);
  rgradient('z',gzlvl0);
  delay(gt0);
  rgradient('z',0.0);
  delay(150.0e-6);

  rgpulse(pw,zero,1.0e-6,0.0);    

  txphase(zero);
  dec2phase(zero);

  delay(2.0e-6);
  rgradient('z',gzlvl1);
  delay(gt1);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(tauxh - gt1 - 4.0e-6);               /* delay=1/4J(XH)   */

  sim3pulse(2*pw,0.0e-6,2*pwx2,zero,zero,zero,0.0,0.0);

  txphase(t2);
  dec2phase(t7);

  delay(tauxh - gt1 - 202.0e-6);               /* delay=1/4J(XH)   */

  delay(2.0e-6);
  rgradient('z',gzlvl1);
  delay(gt1);
  rgradient('z',0.0);
  delay(200.0e-6);

  rgpulse(pw,t2,0.0,0.0);

  /* shaped pulse */
  rlpower(tpwrsl,TODEV);
  shaped_pulse(shp_sl,pw_sl,two,2.0e-6,0.0);
  delay(2.0e-6);
  rlpower(tpwr,TODEV);
  /* shaped pulse */

  delay(2.0e-6);
  rgradient('z',gzlvl2);
  delay(gt2);
  rgradient('z',0.0);
  delay(150.0e-6);

  dec2rgpulse(pwx2,t7,0.0,0.0);

  txphase(t4); dec2phase(t3);

  delay(tau1);                  /* delay=t2/2      */ 

  if(c180_flg[A] == 'y')
     simpulse(2*pw,2*pwx1,t4,zero,0.0,0.0);
  else
     rgpulse(2*pw,t4,0.0,0.0);

  txphase(zero);

  delay(tau1);

  delay(BigT);

  dec2rgpulse(2*pwx2,t3,0.0,0.0);

  dec2phase(t5);

  delay(2.0e-6);
  rgradient('z',gzlvl3);
  delay(gt3);
  rgradient('z',0.0);
  delay(2.0e-6);
 
  if( (c180_flg[A] == 'y') && (pwx1 > pw) )
     delay(BigT + 2*pwx1 - gt3 - 2*GRADIENT_DELAY - 4.0e-6);  
  else 
     delay(BigT + 2*pw - gt3 - 2*GRADIENT_DELAY - 4.0e-6);  

  sim3pulse(pw,0.0e-6,pwx2,zero,zero,t5,0.0,0.0);   /* X read pulse */

  txphase(zero);
  dec2phase(zero);

  delay(2.0e-6);
  rgradient('z',gzlvl5);
  delay(gt5);
  rgradient('z',0.0);
  delay(200.0e-6);

  delay(tauxh - gt5 - 202.0e-6);           /* delay=1/4J (XH)  */

  sim3pulse(2*pw,0.0e-6,2*pwx2,zero,zero,zero,0.0,0.0);  

  txphase(one); dec2phase(one);

  delay(tauxh - gt5 - 202.0e-6);            /* delay=1/4J (XH)  */

  delay(2.0e-6);
  rgradient('z',gzlvl5);
  delay(gt5);
  rgradient('z',0.0);
  delay(200.0e-6);

  sim3pulse(pw,0.0e-6,pwx2,one,zero,one,0.0,0.0);

  dec2phase(zero); txphase(zero);

  delay(2.0e-6);
  rgradient('z',gzlvl5);
  delay(gt5);
  rgradient('z',0.0);
  delay(200.0e-6);

  delay(tauxh - gt5 - 202.0e-6);            /* delay=1/4J (XH)  */

  sim3pulse(2*pw,0.0e-6,2*pwx2,zero,zero,zero,0.0,0.0);

  rlpower(dpwr2,DO2DEV);  /* lower decoupler power for decoupling on decouper channel 2 */
  rlpower(dpwr,DODEV);

  delay(tauxh - 2.0*POWER_DELAY - gt5 - 202.0e-6);  /* delay=1/4J(XH)   */

  delay(2.0e-6);
  rgradient('z',gzlvl5);
  delay(gt5);
  rgradient('z',0.0);
  delay(200.0e-6);

  rgpulse(pw,zero,0.0,0.0);

  delay(BigT1);

  rgpulse(2*pw,zero,0.0,0.0);

  delay(2.0e-6);
  rgradient('z',icosel*gzlvl4);
  delay(gt4);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(BigT1 - gt4 - 4.0e-6 - 2.0*GRADIENT_DELAY);

/* acquire data */

status(D);
     setreceiver(t6);
}