/* N15T2_lek_pfg_sel_enh.c
    This pulse sequence will allow one to perform the following experiment:
   
    15N T2 determination with minimal perturbation of water and enhanced
    sensitivity PFG 15N selection.

    2D               F1       15N
                     F2(acq)   1H (NH)

    Experiment setup-
    Uses three channels:
         1)  1H            - carrier (tof) @  4.73ppm [H2O] 
         2) 13C            - carrier (dof) @118.00ppm [middle of CO and Calpha]
         3) 15N            - carrier(dof2) @119.00ppm [centre of amide 15N]

    Set dm = 'nnnn', dmm = 'cccp' 
    Set dm2 = 'nnny', dmm2 = 'cccp'. [15N decoupling during acquisition]

    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' proton purge pulse at the start of sequence
	f1180		'y' for 180 deg linear phase correction (-90,180) in F1
                            (t1 to be started at half dwell time),
			    otherwise 0 deg linear phase correction (0,0).
 	c_flg		'y' to apply C180 at t1/2 if the sample is 15N and 13C 
                            labelled.

	Standard Settings
	fsat='n', fscuba='n', mess_flg='n', f1180='n', c_flg='n'.

    Set ncyc > 0 to vary 15N T2 relaxation (CPMG) period (approximately
    ncyc*14.4ms).

          
       Written by L. E. Kay 4-02-93
       Modified by L. E. Kay, Feb. 1, 1994 from N15T2_lek_pfg.c
       Must sit on water throughout expt.
       Added a c_flg for use with 15N and 13C labelled samples - RM May 6, 1994
   
    REF: Farrow et. al.  Biochemistry   33, 5984-6003 (1994).

BMRB Pulse Sequence Accession Number: 54
*/

#include <standard.h>
#define PI 3.1416


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

static double d2_init = 0.0;

pulsesequence()

{
/* DECLARE VARIABLES */

 char        fscuba[MAXSTR],fsat[MAXSTR],f1180[MAXSTR],shape[MAXSTR],                        c_flg[MAXSTR];


 int	     phase,t1_counter,icosel;

 double      tau1,                 /* t1 delay */
             hscuba,               /* length of 1/2 scuba delay */
             taua,                 /* < 1 / 4J(NH) 2.25 ms      */
             taub,                 /*  1 / 4J(NH) 2.75 ms      */
             pwn,                  /* PW90 for N-nuc            */
             pwc,                  /* PW90 for C-nuc            */
             tsatpwr,              /* low power level for presat */
             dhpwr2,               /* power level for N hard pulses */
             dhpwr,               /* power level for C hard pulses */
             ncyc,                 /* number of times to loop    */
             time_T2,              /* total time for T2 measuring     */
             tofps,                /* water freq */
	     sw1,
             BigT1,
             tpwrsl,               /* power level for water selective pulse */
             pw_sl,

             gt1,
             gt2,
             gt3,
             gt4,
             gt6,
             gt7,
             gt8,
             gzlvl1,
             gzlvl2,
             gzlvl3,
             gzlvl4,
             gzlvl6,
             gzlvl7,
             gzlvl8;

/* LOAD VARIABLES */

  taua = getval("taua");
  taub = getval("taub");
  pwn = getval("pwn");
  pwc = getval("pwc");
  tsatpwr = getval("tsatpwr");
  dhpwr2 = getval("dhpwr2"); 
  dhpwr = getval("dhpwr"); 
  hscuba = getval("hscuba");
  ncyc = getval("ncyc");
  phase = (int) (getval("phase") + 0.5);
  sw1 = getval("sw1");
  tofps = getval("tofps");
  BigT1 = getval("BigT1");
  tpwrsl = getval("tpwrsl");
  pw_sl = getval("pw_sl");


  gt1 = getval("gt1");
  gt2 = getval("gt2");
  gt3 = getval("gt3");
  gt4 = getval("gt4");
  gt6 = getval("gt6");
  gt7 = getval("gt7");
  gt8 = getval("gt8");
  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl4 = getval("gzlvl4");
  gzlvl6 = getval("gzlvl6");
  gzlvl7 = getval("gzlvl7");
  gzlvl8 = getval("gzlvl8");
  

  getstr("fscuba",fscuba); 
  getstr("fsat",fsat); 
  getstr("f1180",f1180);
  getstr("shape",shape);
  getstr("c_flg",c_flg);

/* 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( dpwr > 0 )
    {
	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( dhpwr2 > 59 )
    {
        printf("don't fry the probe, dhpwr2 too large for cpmg !");
        abort(1);
    }

    time_T2 = ncyc*(32.0*pwn + 32.0*450.0e-6);

    if(time_T2 > 0.190)
    {
       printf("net T2 recovery time exceeds 0.190 sec. May be too long \n");
       abort(1);
    }  

   if(ncyc > 12)
   {
      printf("value of ncyc must be less than or equal to 12\n");
      abort(1);
   }

   if(gt1 > 15e-3 || gt2 > 15e-3 || gt3 > 15e-3 || gt4 > 15e-3 || gt6 > 15e-3 || gt7 > 15e-3 || gt8 > 15e-3)
   {
      printf("gti must be less than 15e-3\n");
      abort(1);
   }

/* LOAD VARIABLES */

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

/* Phase incrementation for hypercomplex 2D data */

   if (phase == 2) {
     tsadd(t6,2,4);
     icosel = -1;
   }

   else
     icosel = 1;

/* Set up f1180  */

   tau1 = d2;
   if(f1180[A] == 'y') {
     tau1 += ( 1.0 / (2.0*sw1) );
     if(tau1 < 0.2e-6) tau1 = 0.2e-6;
   }
  
   tau1 = tau1/2.0;

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

   if( ix == 1 ) d2_init = d2;
   t1_counter = (int) ( (d2-d2_init)*sw1 + 0.5 );
   if( t1_counter %2 ) {
      tsadd(t2,2,4);
      tsadd(t7,2,4);
   }


/* 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);

 if(fsat[0] == 'y')
{
  offset(tofps,TODEV);            /* move H carrier to the water */
  rgpulse(d1,zero,2.0e-6,2.0e-6);  /* presat. with transmitter */
  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);
}

  offset(tof,TODEV);

status(C);

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

  dec2rgpulse(pwn,zero,0.0,0.0);

  delay(2.0e-6);
  rgradient('z',gzlvl1);
  delay(gt1);
  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',gzlvl2);
  delay(gt2);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(taua - gt2 - 4.0e-6);                   /* delay < 1/4J(XH)   */

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

  txphase(t1);
  dec2phase(t2);

  delay(taua - gt2 - 202.0e-6);

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

  rgpulse(pw,t1,0.0,0.0);

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

  delay(2.0e-6);
  rgradient('z',gzlvl3);
  delay(gt3);
  rgradient('z',0.0);
  delay(150.0e-6);

  dec2rgpulse(pwn,t2,0.0,0.0);

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

  delay(taub - gt4 - 4.0e-6);

  sim3pulse(2.0*pw,0.0e-6,2.0*pwn,zero,zero,zero,0.0,0.0);

  delay(taub - gt4 - 202.0e-6 + (2.0/PI)*pwn);

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

  initval(ncyc,v4);

if(ncyc > 0) {

 loop(v4,v5);

    initval(3.0,v1);
 
    loop(v1,v2);
      delay(450.0e-6);
      dec2rgpulse(2*pwn,zero,0.0,0.0);
      delay(450.0e-6);
    endloop(v2);

    delay(450.0e-6);
    dec2rgpulse(2*pwn,zero,0.0,0.0);
    delay(450.0e-6 - pw);
    rgpulse(2*pw,zero,0.0,0.0);
    delay(450.0e-6 - pw);
    dec2rgpulse(2*pwn,zero,0.0,0.0);
    delay(450.0e-6);
  
    initval(3.0,v1);
 
    loop(v1,v3);
      delay(450.0e-6);
      dec2rgpulse(2*pwn,zero,0.0,0.0);
      delay(450.0e-6);
    endloop(v3);

    initval(3.0,v1);
 
    loop(v1,v2);
      delay(450.0e-6);
      dec2rgpulse(2*pwn,zero,0.0,0.0);
      delay(450.0e-6);
    endloop(v2);

    delay(450.0e-6);
    dec2rgpulse(2*pwn,zero,0.0,0.0);
    delay(450.0e-6 - pw);
    rgpulse(2*pw,two,0.0,0.0);
    delay(450.0e-6 - pw);
    dec2rgpulse(2*pwn,zero,0.0,0.0);
    delay(450.0e-6);
  
    initval(3.0,v1);
 
    loop(v1,v3);
      delay(450.0e-6);
      dec2rgpulse(2*pwn,zero,0.0,0.0);
      delay(450.0e-6);
    endloop(v3);

  endloop(v5);
 }

  txphase(t5); 
  dec2phase(t4);

  if(c_flg[A] == 'y')
     delay(taub + 2*pwc + (2.0/PI)*pwn);                   
  else
     delay(taub + 2*pw + (2.0/PI)*pwn);                   

  dec2rgpulse(2*pwn,t4,0.0,0.0);

  delay(tau1);

  if(c_flg[A] == 'y')
     simpulse(2*pw,2*pwc,t5,zero,0.0,0.0);
  else
     rgpulse(2*pw,t5,0.0,0.0);

  dec2phase(t6);
  txphase(zero);

  delay(2.0e-6);
  rgradient('z',icosel*gzlvl6);
  delay(gt6);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(taub - gt6 - 4.0e-6 - 2.0*GRADIENT_DELAY);

  delay(tau1);

  sim3pulse(pw,0.0e-6,pwn,zero,zero,t6,0.0,0.0);

  dec2phase(zero);

  delay(2.0e-6);
  rgradient('z',gzlvl7);
  delay(gt7);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(taua - gt7 - 4.0e-6);                /* delay=1/4J (XH)  */

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

  delay(taua - gt7 - 202.0e-6 - 2.0e-6);     /* delay=1/4J (XH)  */

  delay(2.0e-6);
  rgradient('z',gzlvl7);
  delay(gt7);
  rgradient('z',0.0);
  delay(200.0e-6);

  txphase(one); dec2phase(one);

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

  dec2phase(zero); txphase(zero);

  delay(2.0e-6);
  rgradient('z',gzlvl7);
  delay(gt7);
  rgradient('z',0.0);
  delay(2.0e-6);

  delay(taua - gt7 - 4.0e-6);                /* delay=1/4J (XH)  */

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

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

  delay(taua - 2.*POWER_DELAY - gt7 - 202.0e-6 - 2.0e-6);  /* delay=1/4J(XH)   */

  delay(2.0e-6);
  rgradient('z',gzlvl7);
  delay(gt7);
  rgradient('z',0.0);
  delay(200.0e-6);

  rgpulse(pw,two,2.0e-6,0.0);
  txphase(zero);

  delay(BigT1);

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

  delay(2.0e-6);
  rgradient('z',gzlvl8);
  delay(gt8);
  rgradient('z',0.0);
  delay(2.0e-6);

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

/* acquire data */

status(D);
     setreceiver(t7);
}