/*  noesyc_pfg_h2o_NC+_500.c

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

    3D 13C and 15N edited noesy with separation via the carbon or nitrogen of
    the destination site recorded on a water sample 
                            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)

    This sequence uses the standard three channel configuration
         1)  1H             - carrier (tof) @ 4.7 ppm [H2O]
         2) 13C             - carrier (dof) @  67 ppm
         3) 15N             - carrier (dof2)@ 119 ppm [centre of amide 15N]
 
    Set dm = 'nnny', dmm = 'cccp'   [13C decoupling during acquisition]
                     dseq = 'garp1'
    Set dm2 = 'nnny', dmm2 = 'cccp' [15N decoupling during acquisition]
                      dseq2 = 'waltz16'
 
    Must set phase = 1,2 and phase2 = 1,2 for States-TPPI acquisition in t1[1H]
    and t2 [C and N] respectively. 
 
    Flags
        mess_flg        'y' for Messerle type purging at the beginning of seq.
        shp_flg         'y' for sel. excitation and dephasing of water
        fsat            'y' for presaturation of H2O
        fscuba          'y' for apply scuba pulse after presaturation of H2O
        f1180           'y' for 180 deg linear phase correction in F1
        f2180           'y' for 180 deg linear phase correction in F2
        c180_flg        'y' for C180 at t2/2 when checking F1/F3 2D
 
    Standard Settings
        mess_flg='y',shp_flg='n',fsat='n',fscuba='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 to allow for simultaneous N, C acquisition
                July 19, 1993.
    Modified by L.E. Kay on Aug. 15, 1994 to make unity+ friendly 
    Will work on the 600
  
    REF: Pascal et. al.   J. Magn. Reson.  B 103, 197-201 (1994). 

*/

#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              */
             pwndec,      /* PW90 for 15N decoupling during t1   */
             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,
             gt0,
             gt1,
             gt2,
             gt3,
             gt4,
             gt5,
             gt6,
             gt7,
             gzlvl0, 
             gzlvl1, 
             gzlvl2, 
             gzlvl3, 
             gzlvl4, 
             gzlvl5,
             gzlvl6,
             gzlvl7; 

/*  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");
  pwndec = getval("pwndec");
  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");
  gt0 = getval("gt0");
  gt1 = getval("gt1");
  gt2 = getval("gt2");
  gt3 = getval("gt3");
  gt4 = getval("gt4");
  gt5 = getval("gt5");
  gt6 = getval("gt6");
  gt7 = getval("gt7");
  gzlvl0 = getval("gzlvl0");
  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl4 = getval("gzlvl4");
  gzlvl5 = getval("gzlvl5");
  gzlvl6 = getval("gzlvl6");
  gzlvl7 = getval("gzlvl7");

/* 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! Should be 'nnn' ");
        abort(1);
    }

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

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

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

    if( dpwr2 > 47 )
    {
        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( pwndec > 600.0e-6 || pwndec < 100.0e-6 )
    {
        printf("dont fry the probe, pwndec is wrong ! ");
        abort(1);
    } 

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

    if( pwcodec < 260.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( 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,1,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) - 4.0/PI*pw - 4.0*pwc - 6.0e-6 - 2*POWER_DELAY                   - PRG_START_DELAY - PRG_STOP_DELAY );
    }

    else tau1 = tau1 - (4.0/PI*pw + 4.0*pwc + 6.0e-6 + 2*POWER_DELAY                   + PRG_START_DELAY + PRG_STOP_DELAY );

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

    tau1 = tau1/2.0;

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

/* BEGIN ACTUAL PULSE SEQUENCE */

/* Receiver off time */

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(dpwr2,DO2DEV);       /* Set Dec2 to low power for decoupling */
   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);


/*   decrgpulse(pwc,zero,0.0,0.0);  

   delay(2.0e-6);
   rgradient('z',gzlvl0);
   delay(1.5*gt0);
   rgradient('z',0.0);
   delay(150.0e-6);
*/
 
   rgpulse(pw,t1,0.0,0.0);                    /* 90 deg 1H pulse */

   rlpower(dpwr2,DO2DEV);
   dec2prgon("waltz16",pwndec,dres2);
   dec2on();

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

   dec2off();
   dec2prgoff();
   rlpower(dhpwr2,DO2DEV);

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

}