/*  noesyc_pfg_h2o_NC+_purge_500.c

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

    3D noesy from unlabeled to labeled (15N,13C). Note the corresponding
     expt from labeled to unlabeled is called cnoesy_3c_sed_pfg_purge+_500.c
     (This one is better with suppression levels like the CN noesy).
                       F1      1H
                       F2      13C and 15N
                       F3(acq) 1H

    NOTE: cnoesy_3c_sed_pfg(1).c 3D C13 edited NOESYs with separation via the
          carbon of the origination site. (D2O experiments)

    Uses three channels:
         1)  1H             - carrier @ water  
         2) 13C             - carrier @ 43 ppm [CA/CB]
         3) 15N             - carrier @ 119ppm [centre of amide N]


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

    Must set phase = 1,2 and phase2 = 1,2 for States-TPPI acquisition in
    t1 [H]  and t2 [C and N].
  
    Flags
         fsat       'y' for presaturation of H2O
         fscuba     'y' to apply scuba pulse after presaturation of H2O
         mess_flg   'y' for 1H (H2O) purging pulse before relaxation delay
         shp_flg    'y' to selective excitation & suppression of H2O
         f1180      'y' for 180deg linear phase correction in F1
                        otherwise 0deg linear phase correction
         f2180      'y' for 180deg linear phase correction in F2
         c180_flg   'y' when recording F1,F3 2D.

    Standard Settings:
       fsat='n', fscuba='n', mess_flg='y', shp_flg='n', f1180='y', f2180='y'
       c180_flg='n'
         

    Set f1180 = 'y' and f2180 = 'y' for (-90,180) in F1 and (-90,180) in  F2.   
    Note: Zero order phase correct may not be exactly -90 in F2 due to seduce.

    Modified by L.E. Kay on Aug. 15, 1994 from noesyc_pfg_h2o_NC+.c 
    Will work on the 600

    Modified by L. E. Kay on March 6, 1995 to include 15N purging

BMRB Pulse Sequence Accession Number: 50
*/

#include <standard.h>
#define PI 3.1416

static int  phi1[8]  = {0,0,0,0,2,2,2,2},
            phi2[16] = {0,0,0,0,0,0,0,0,2,2,2,2,2,2,2,2},
            phi3[8]  = {0,0,0,0,2,2,2,2},
            rec[16]  = {0,3,2,1,2,1,0,3,2,1,0,3,0,3,2,1},
            phi5[4]  = {0,1,2,3},
            phi6[2]  = {0,2},
            phi7[4]  = {1,1,3,3};
                    
static double d2_init=0.0, d3_init=0.0;
            
pulsesequence()
{
/* DECLARE VARIABLES */

 char       fsat[MAXSTR],
	    fscuba[MAXSTR],
            f1180[MAXSTR],    /* Flag to start t1 @ halfdwell             */
            f2180[MAXSTR],    /* Flag to start t2 @ halfdwell             */
            c180_flg[MAXSTR],
            codecseq[MAXSTR],
            mess_flg[MAXSTR],
            shp_flg[MAXSTR],
            shape[MAXSTR];

 int         phase, phase2, 
             t1_counter,   /* used for states tppi in t1           */ 
             t2_counter;   /* used for states tppi in t2           */ 

 double      tau1,         /*  t1 delay */
             tau2,         /*  t2 delay */
             taua,         /*  ~ 1/4JHC =  1.6 ms */
             mix,         /* mixing time in seconds */
             pwn,          /* PW90 for 15N pulse              */
             pwc,          /* PW90 for c nucleus @ dhpwr         */
             pwcodec,      /* PW for C' nucleus @ dpwrco seduce dec  */
             tsatpwr,      /* low level 1H trans.power for presat  */
             dhpwr,        /* power level for 13C pulses on dec1  */
             dpwrco,       /* power level for C' seduce decoupling  */
             dhpwr2,       /* high dec2 pwr for 15N hard pulses    */
             sw1,          /* sweep width in f1                    */             
             sw2,          /* sweep width in f2                    */             
             tofps,        /* tof for presat                       */ 
             dressed,      /* decoupler resolution for seduce decoupling */
             tpwrmess,    /* power level for Messerlie purge */
             dly_pg1,     /* duration of first part of purge */
             dly_wt,
             sh_pwr,
             sh_pwr2,
             pw_sh_ms,
             p2_sh_ms,
             jhc1,       /* smallest coupling that you wish to purge */
             jhc2,       /* largest coupling that you wish to purge */
             taud,       /* 1/(2jhc1)   */    
             taue,       /* 1/(2jhc2)   */    


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

/*  variables commented out are already defined by the system      */


/* LOAD VARIABLES */


  getstr("fsat",fsat);
  getstr("f1180",f1180);
  getstr("f2180",f2180);
  getstr("fscuba",fscuba);
  getstr("codecseq",codecseq);
  getstr("c180_flg",c180_flg);
  getstr("mess_flg",mess_flg);
  getstr("shp_flg",shp_flg);
  getstr("shape",shape);

  tofps  = getval("tofps");
  taua   = getval("taua"); 
  mix = getval("mix");
  pwc = getval("pwc");
  pwcodec = getval("pwcodec");
  pwn = getval("pwn");
  tpwr = getval("tpwr");
  tsatpwr = getval("tsatpwr");
  dhpwr = getval("dhpwr");
  dpwr = getval("dpwr");
  dpwrco = getval("dpwrco");
  dhpwr2 = getval("dhpwr2");
  phase = (int) ( getval("phase") + 0.5);
  phase2 = (int) ( getval("phase2") + 0.5);
  sw1 = getval("sw1");
  sw2 = getval("sw2");
  dressed = getval("dressed");
  tpwrmess = getval("tpwrmess");
  dly_pg1 = getval("dly_pg1");
  dly_wt = getval("dly_wt");
  sh_pwr = getval("sh_pwr");
  sh_pwr2 = getval("sh_pwr2");
  pw_sh_ms = getval("pw_sh_ms");
  p2_sh_ms = getval("p2_sh_ms");
  jhc1 = getval("jhc1");
  jhc2 = getval("jhc2");

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

/* LOAD PHASE TABLE */

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

/* CHECK VALIDITY OF PARAMETER RANGES */


    if((dm[A] == 'y' || dm[B] == 'y' || dm[C] == '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 > 6 )
    {
        printf("TSATPWR too large !!!  ");
        abort(1);
    }

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

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

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

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


    if( pw > 200.0e-6 )
    {
        printf("dont fry the probe, pw too high ! ");
        abort(1);
    } 

    if( pwn > 200.0e-6 )
    {
        printf("dont fry the probe, pwn too high ! ");
        abort(1);
    } 

    if( pwc > 200.0e-6 )
    {
        printf("dont fry the probe, pwc too high ! ");
        abort(1);
    } 

    if( pwcodec < 300.0e-6 )
    {
        printf("dont fry the probe, pwcodec too high ! ");
        abort(1);
    } 
    if ( tpwrmess > 56 )
    {
        printf("dont fry the probe, tpwrmess too high ! ");
        abort(1);
    }
    if ( dly_pg1 > 0.010)
    {
        printf("dont fry the probe, dly_pg1 too long ! ");
        abort(1);
    }

    if(dmf2 > 4000) {
        printf("15N decoupling pulses must be less than 250e-6\n");
        abort(1);
     }

    if( gt0 > 15e-3 || gt1 > 15e-3 || gt2 > 15e-3 || gt3 > 15e-3 || gt4 > 15e-3 || gt5 > 15e-3 || gt6 > 15e-3 || gt7 > 15e-3 ) 
    {
        printf("gti values < 15e-3\n");
        abort(1);
    } 

   if( gzlvl3*gzlvl4 > 0.0 ) 
    {
        printf("gt3 and gt4 must be of opposite sign for optimal water suppression\n");
        abort(1);
     }
     if (pw_sh_ms > 50.0e-6 || p2_sh_ms > 50.0e-6 )
     {
        printf( " pwi_sh_ms too long !! ");
        abort(1);
     }
     if (sh_pwr > 25 || sh_pwr2 > 25 )
     {
        printf( " sh_pwri too high !! " );
        abort(1);
     }


/*  Phase incrementation for hypercomplex 2D data */

    if (phase == 2)
      tsadd(t1,3,4);
    if (phase2 == 2)
      tsadd(t2,1,4);

/*  Set up f1180  tau1 = t1               */
   
    tau1 = d2;
    if(f1180[A] == 'y') {
        tau1 += ( 1.0 / (2.0*sw1) );
    }

    if(tau1 < 0.2e-6) tau1 = 2.0e-7;

/*  Set up f2180  tau2 = t2               */

    tau2 = d3;
    if(f2180[A] == 'y') {
        tau2 += ( 1.0 / (2.0*sw2) - (4.0/PI)*pwn - 2.0*pw - PRG_START_DELAY - PRG_STOP_DELAY - 2.0*POWER_DELAY - 4.0e-6); 
    }

    else tau2 = tau2 - ((4.0/PI)*pwn + 2.0*pw + PRG_START_DELAY + PRG_STOP_DELAY + 2.0*POWER_DELAY + 4.0e-6); 

        if(tau2 < 0.2e-6)  tau2 = 4.0e-7;
        tau2 = tau2/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(t1,2,4);     
      tsadd(t4,2,4);    
    }

   if( ix == 1) d3_init = d3 ;
   t2_counter = (int) ( (d3-d3_init)*sw2 + 0.5 );
   if(t2_counter % 2) {
      tsadd(t2,2,4);  
      tsadd(t4,2,4);    
    }

   if(jhc1 != 0.0) taud = 1.0/(2.0*jhc1);
   if(jhc2 != 0.0) taue = 1.0/(2.0*jhc2);

/* BEGIN ACTUAL PULSE SEQUENCE */

status(A);
   delay(5.0e-6);
   rlpower(tsatpwr,TODEV);      /* Set transmitter power for 1H presaturation */
   rlpower(dhpwr,DODEV);        /* Set Dec1 power for hard 13C pulses         */
   rlpower(dhpwr2,DO2DEV);      /* Set Dec2 to high power for 15N pulses  */
   delay(5.0e-6);

/* Presaturation Period */

status(B);

   if (mess_flg[A] == 'y') {

      offset(tofps,TODEV);
      rlpower(tpwrmess,TODEV);
      txphase(zero);
      rgpulse(dly_pg1,zero,2.0e-6,2.0e-6);
      txphase(one);
      rgpulse(dly_pg1/1.62,one,2.0e-6,2.0e-6);

      rlpower(tsatpwr,TODEV);
   }

   if (fsat[0] == 'y')
   {
        offset(tofps,TODEV);
	delay(2.0e-5);
    	rgpulse(d1,zero,2.0e-6,2.0e-6); /* presat */
   	rlpower(tpwr,TODEV);      /* Set transmitter power for hard 1H pulses */
	delay(2.0e-5);
	if(fscuba[0] == 'y')
	{
		delay(2.2e-2);
		rgpulse(pw,zero,2.0e-6,0.0);
		rgpulse(2*pw,one,2.0e-6,0.0);
		rgpulse(pw,zero,2.0e-6,0.0);
		delay(2.2e-2);
	}
   }

   else
   {
    delay(d1);
   }

  if (shp_flg[A] == 'y') {

     offset(tofps,TODEV);
     txphase(t6);   
     rlpower(sh_pwr,TODEV);   /* set the power for the shaped pulse */
     shaped_pulse(shape,pw_sh_ms*1000.0,t6,2.0e-6,2.0e-6);
 
     delay(2.0e-6);
     rgradient('z',gzlvl6);
     delay(gt6);
     rgradient('z',0.0);
     delay(50.0e-6);

     txphase(t7);
     rlpower(sh_pwr2,TODEV);
     shaped_pulse(shape,p2_sh_ms*1000.0,t7,2.0e-6,2.0e-6);

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

     delay(dly_wt);
  }

   rlpower(tpwr,TODEV);           /* Set transmitter power for hard 1H pulses */
   offset(tof,TODEV);
   txphase(t1);
   decphase(zero);
   dec2phase(zero);

/* Begin Pulses */

status(C);

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

   rgpulse(pw,t1,0.0,0.0);                    /* 90 deg 1H pulse */

   txphase(zero);

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

   delay(taud - gt8 - 102.0e-6);

   sim3pulse(2.0*pw,pwc,2.0*pwn,zero,zero,zero,0.0,0.0);

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

   delay(5.5e-3 - taud - pwn - gt8 - 102.0e-6); 

   dec2rgpulse(pwn,zero,0.0,0.0);

   delay(taud + taue - 5.5e-3);

   decrgpulse(pwc,zero,0.0,0.0);

   delay(taud - taue - pwc - 2.0e-6);

   delay(tau1);

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

   delay(mix - 10.0e-3);

   delay(2.0e-6);
   rgradient('z',gzlvl0);
   delay(gt0);
   rgradient('z',0.0);
   delay(102.0e-6);

   decrgpulse(pwc,zero,0.0,0.0); 

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

   delay(10.0e-3 - gt1 - gt0 - 8.0e-6);
   delay(150.0e-6);

   rgpulse(pw,zero,0.0,0.0);

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

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

   sim3pulse(2*pw,2*pwc,2*pwn,zero,zero,zero,0.0,0.0);

   txphase(one); decphase(t2); dec2phase(t2);

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

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

  rgpulse(pw,one,0.0,0.0);

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

  dec2rgpulse((PI-2.0)/PI*(pwn-pwc),t2,0.0,0.0);
  sim3pulse(0.0e-6,pwc,pwc,zero,t2,t2,0.0,0.0);
  dec2rgpulse((2.0/PI)*(pwn-pwc),t2,0.0,0.0);

   if( c180_flg[A] == 'n' ) {

   delay(2.0e-6);

   /* CO decoupling on */
   rlpower(dpwrco,DODEV);
   decprgon(codecseq,pwcodec,dressed);
   decon();
   /* CO decoupling on */

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

   delay(tau2);

   /* CO decoupling off */
   decoff();
   decprgoff();
   rlpower(dhpwr,DODEV);
   /* CO decoupling off */

   decphase(zero);
   delay(2.0e-6);

  }

  else
   sim3pulse(2*pw,2*pwc,2*pwn,zero,zero,zero,2.0e-6,2.0e-6);

  dec2rgpulse((2.0/PI)*(pwn-pwc),zero,0.0,0.0);
  sim3pulse(0.0e-6,pwc,pwc,zero,zero,zero,0.0,0.0);
  dec2rgpulse((PI-2.0)/PI*(pwn-pwc),zero,0.0,0.0);

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

   rgpulse(pw,t5,2.0e-6,0.0);

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

   txphase(zero);
   delay(taua - gt5 - 4.0e-6);

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

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

   txphase(t5);
   delay(taua - 2*pwn - 2.0e-6 - gt5 - 4.0e-6);

   sim3pulse(0.0e-6,pwc,pwn,t5,zero,zero,0.0,0.0);
   sim3pulse(0.0e-6,pwc,pwn,t5,t3,t3,2.0e-6,0.0);

   rgpulse(pw,t5,0.0,0.0);

   rlpower(dpwr2,DO2DEV);  /* Turn down for decoupling */
   rlpower(dpwr,DODEV);  /* Set power for decoupling */
    

/* BEGIN ACQUISITION */

status(D);
 /*  rcvron();    */
   setreceiver(t4);

}