/*  cnoesy_3c_sed_pfg_purge+_500.c

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

    3D C13 edited noesy with separation via the carbon of the origination site 
    with purging so that NOEs are from carbon label to non-carbon label.
                       F1      1H
                       F2      13C
                       F3(acq) 1H 

    Uses three channels:
         1)  1H             - carrier @ 3.0 ppm 
         2) 13C             - carrier @ 43 ppm [CA/CB]
       { 3) 15N             - carrier @ 119ppm }


    Set dm = 'nnnn', dmm = 'cccc' [NO 13C decoupling during acquisition].
    Set dm2 = 'nnnn', dmm2 = 'cccc' 

    Must set phase = 1,2 and phase2 = 1,2 for States-TPPI
    acquisition in t1 [H]  and t2 [C].

    Flags
          fsat       'y' for presaturation of H2O
          fscuba     'y' to apply scuba pulse after presaturation of H2O
          f1180      'y' for 180deg linear phase correction in F1
                         otherwise 0deg linear phase correction in F1.
          f2180      'y' for 180deg linear phase correction in F2
          c180_flg   'y' for C180 at t2/2 when recording F1,F3 2D.
          shared     'y' for extended Bax-Logan trick
 
    Standard Settings:
          fsat='n', fscuba='n', f1180='n', f2180='y', c180_flg='n', shared='y'

   Written by L.E.K on Sept 22/93

   If shared[A] == 'y' then uses the extended Bax/Logan trick to save 6 ms
   (See Muhandiram, Xu and Kay  J. Biomol. NMR 3, 463-470 (1993)).

   NOTE: SET pwn=0 dhpwr2=0 or else the first point in the C dimension
         will not be correct since the delays are greater than half of 
         the dwell. Could use pwn and then use backwards lp to replace
         the first real and imag. plane. In either case f2180='y'.

   REF: Lee et. al.  FEBS Letters 350, 87-90 (1994)

BMRB Pulse Sequence Accession Number: 49
*/

#include <standard.h>
#define PI 3.1416

static int  phi1[2]  = {0,2},
            phi2[16] =  {0,0,0,0,0,0,0,0,2,2,2,2,2,2,2,2},
            phi3[4]  = {0,0,2,2},
            phi4[8]  = {0,0,0,0,1,1,1,1},
            rec[16]  = {0,2,0,2,3,1,3,1,
                        2,0,2,0,1,3,1,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             */
            codecseq[MAXSTR],
            shared[MAXSTR],
            c180_flg[MAXSTR];


 int         phase, phase2, ni, 
             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 */
             jhc1,         /* smallest coupling that you wish to purge */
             jhc2,         /* largest coupling that you wish to purge */
             taud,         /* 1/(2jhc1)   */
             taue,         /* 1/(2jhc2)   */
             zeta,         /* Bax/Logan trick       */
             zeta1,
             tauch,        /* taua/2.0              */ 

             gt1,
             gt2,
             gt3,
             gt4,
             gt5,
             gzlvl1, 
             gzlvl2, 
             gzlvl3, 
             gzlvl4, 
             gzlvl5; 

/*  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("shared",shared);
  getstr("c180_flg",c180_flg);

  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");
  jhc1 = getval("jhc1");
  jhc2 = getval("jhc2");
  ni = getval("ni");

  tauch = taua/2.0;

  gt1 = getval("gt1");
  gt2 = getval("gt2");
  gt3 = getval("gt3");
  gt4 = getval("gt4");
  gt5 = getval("gt5");

  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl4 = getval("gzlvl4");
  gzlvl5 = getval("gzlvl5");

/* LOAD PHASE TABLE */

  settable(t1,2,phi1);
  settable(t2,16,phi2);
  settable(t3,4,phi3);
  settable(t4,8,phi4);
  settable(t5,16,rec);

/* 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' || dm2[D] == 'y'))
    {
        printf("incorrect dec2 decoupler flags! Should be 'nnnn' ");
        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 > 45 )
    {
        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( dhpwr > 62 )
    {
        printf("don't fry the probe, DHPWR 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( gt1 > 15e-3 || gt2 > 15e-3 || gt3 > 15e-3 || gt4 > 15e-3 || gt5 > 15e-3 ) 
    {
        printf("gti values < 15e-3\n");
        abort(1);
    } 

    if(jhc1 > jhc2)
    {
        printf("set jhc1 < jhc2\n");
        abort(1);
    }

    if(shared[A] != 'y' && shared[A] != 'n')
    {
        printf("shared must either by y or n\n");
        abort(1);
    }

    if(shared[A] == 'y') {
       if(ni != 1)
           zeta = (tauch + 4.0*pwc + 4.0e-6 - gt1 - 102.0e-6 - 2.0e-6 - PRG_START_DELAY + PRG_STOP_DELAY - 2.0e-6) / ( (double) (ni-1) );
        else zeta = 0.0;
     }

/*  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 = 8.0e-7;
    }

    if(shared[A] == 'n')
        tau1 = tau1/2.0;

/*  Set up f2180  tau2 = t2               */

    tau2 = d3;

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

        if(tau2 < 0.2e-6)  tau2 = 4.0e-7;
                          
        tau2 = tau2/2.0;

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

/* 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(t5,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(t5,2,4);    
    }

   zeta1 = zeta*( (double) t1_counter );
   if(shared[A] == 'y') {
      tau1 = tau1/4.0 - zeta1;
   }
/* BEGIN ACTUAL PULSE SEQUENCE */

status(A);

   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 pulses       */

/* Presaturation Period */

status(B);
   if (fsat[0] == 'y')
   {
        offset(tofps,TODEV);
	delay(2.0e-5);
        rgpulse(d1,zero,2.0e-6,2.0e-6); /* presaturation */
   	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);
   }
   rlpower(tpwr,TODEV);           /* Set transmitter power for hard 1H pulses */
   offset(tof,TODEV);
   txphase(t1);
   decphase(zero);
   dec2phase(zero);
   delay(1.0e-5);

/* Begin Pulses */

status(C);

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

/* ensure that magnetization originates on H and not 13C */

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

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

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

   txphase(t3); decphase(t2);
   delay(taua - PRG_START_DELAY + PRG_STOP_DELAY + 4.0*pwc + 4.0e-6 - gt1 - 102.0e-6 + 2*pwn);

   decrgpulse(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,t3,0.0e-6,0.0);

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

   delay(tau2);

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

   delay(2.0e-6);

  }

  else
     simpulse(2*pw,2*pwc,t3,zero,2.0e-6,2.0e-6);

   decrgpulse(pwc,zero,0.0,0.0);

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

   txphase(zero); decphase(one);

   delay(taua + tau1 - gt1 - 102.0e-6);

   decrgpulse(pwc,one,0.0,0.0);
   decrgpulse(2*pwc,zero,2.0e-6,0.0);
   decrgpulse(pwc,one,2.0e-6,0.0);
  
   delay(tau1);

  }

  else {

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

   txphase(zero); decphase(one);

   delay(tauch + zeta1 - gt1 - 102.0e-6 + tau1);

   decrgpulse(pwc,one,0.0,0.0);
   decrgpulse(2.0*pwc,zero,2.0e-6,0.0);
   decrgpulse(pwc,one,2.0e-6,0.0);

   decphase(t2);

   delay(tau1);

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

   delay(tauch - zeta1 - gt1 - 102.0e-6 + 4.0*pwc + 4.0e-6 - PRG_START_DELAY +   PRG_STOP_DELAY + 2.0*pwn);

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

   decrgpulse(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,t3,0.0e-6,0.0);

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

   delay(tau2);

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

   delay(2.0e-6);

  }

  else
     simpulse(2*pw,2*pwc,t3,zero,2.0e-6,2.0e-6);

   decrgpulse(pwc,zero,0.0,0.0);

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

   txphase(zero); decphase(one);

   delay(tauch + zeta1 - gt2 - 102.0e-6 + tau1);

   decrgpulse(pwc,one,0.0,0.0);
   decrgpulse(2.0*pwc,zero,2.0e-6,0.0);
   decrgpulse(pwc,one,2.0e-6,0.0);

   decphase(zero);

   delay(tau1);

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

   delay(tauch - zeta1 - gt2 - 102.0e-6 + 4.0*pwc + 4.0e-6);

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

  }

   rgpulse(pw,zero,0.0,0.0);

   delay(mix/2.0);

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

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

   delay(mix/2.0 - gt3 - gt4 - 4.0e-6 - 502.0e-6 - POWER_DELAY);

   rlpower(dpwr2,DO2DEV);
 
   rgpulse(pw,t4,2.0e-6,0.0);

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

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

   simpulse(2.0*pw,pwc,zero,zero,0.0,0.0);

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

   delay(taue - gt5 - 102.0e-6);

   decrgpulse(pwc,zero,0.0,0.0);

   delay(taud - taue - pwc - POWER_DELAY);

   rlpower(dpwr,DODEV);  /* Set power for decoupling */
    
/*   rcvron();  */          /* Turn on receiver to warm up before acq */ 

/* BEGIN ACQUISITION */

status(D);
setreceiver(t5);

}