/*  c4Dnoesy_d2o_pfg_500.c

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

    4D C13 edited noesy  
                       F1      1H
                       F2      13C
                       F3      13C
                       F4      1H

    Uses two channels:
         1)  1H             - carrier (tof) @ 3.0 ppm 
         2) 13C             - carrier (dof) @ 43 ppm (CA/CB)


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

    Must set phase = 1,2 phase2 = 1,2 and phase3 = 1,2 for States-TPPI
    acquisition in t1 [H]  and t2 [C] and t3 [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
         f3180    'y' for 180deg linear phase correction in F3
         c180_flg 'y' for c180 at t2/2 and t3/2 when recording F1,F4 2D.

    Standard Settings
         fsat='n', fscuba='n', f1180='n', f2180='y', f3180='y', c180_flg='n'

    Set f1180 = 'n' f2180 = 'y' and f3180 = 'y' for (0,0) (-90,180) (-90,180)   
    Note: Zero order phase correct may not be exactly -90. ACTUALLY SEDUCE
          APPEARS NOT TO MAKE A DIFFERENCE TO ZERO ORDER PHASE CORRECTION.
          (-90,180) is correct. (LEK, DEC 30, 1992).

   Written by Lewis Kay 11/11/93 
   Modified to eliminate small acq. dim. quad. images (LK, 11/23/93)

   Ref: Vuister et. al.  J. Magn. Reson. B 101, 210 (1993)

*/

#include <standard.h>
#define PI 3.1416

static int  phi1[1]  = {0},
            phi2[2]  = {0,2},
            phi3[4]  = {0,0,2,2},
            phi4[4]  = {0,0,1,1},
            rec[4]   = {0,2,1,3};
                    
static double d2_init=0.0, d3_init=0.0, d4_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             */
            f3180[MAXSTR],    /* Flag to start t3 @ halfdwell             */
            codecseq[MAXSTR],
            c180_flg[MAXSTR];

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

 double      tau1,         /*  t1 delay */
             tau2,         /*  t2 delay */
             tau3,         /*  t3 delay */
             taua,         /*  ~ 1/2JHC =  3.4 ms */
             mix,          /* mixing time in seconds */
             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  */
             sw1,          /* sweep width in f1                    */          
             sw2,          /* sweep width in f2                    */           
             sw3,          /* sweep width in f3                    */           
             tofps,        /* tof for presat                       */ 
             dressed,      /* decoupler resolution for seduce decoupling */
             gt0,
             gt1,
             gt2,
             gt3,
             gt4,
             gzlvl0, 
             gzlvl1, 
             gzlvl2, 
             gzlvl3, 
             gzlvl4; 

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


/* LOAD VARIABLES */

  getstr("fsat",fsat);
  getstr("f1180",f1180);
  getstr("f2180",f2180);
  getstr("f3180",f3180);
  getstr("fscuba",fscuba);
  getstr("codecseq",codecseq);
  getstr("c180_flg",c180_flg);

  tofps  = getval("tofps");
  taua   = getval("taua"); 
  mix = getval("mix");
  pwc = getval("pwc");
  pwcodec = getval("pwcodec");
  tpwr = getval("tpwr");
  tsatpwr = getval("tsatpwr");
  dhpwr = getval("dhpwr");
  dpwr = getval("dpwr");
  dpwrco = getval("dpwrco");
  phase = (int) ( getval("phase") + 0.5);
  phase2 = (int) ( getval("phase2") + 0.5);
  phase3 = (int) ( getval("phase3") + 0.5);
  sw1 = getval("sw1");
  sw2 = getval("sw2");
  sw3 = getval("sw3");
  dressed = getval("dressed");
  gt0 = getval("gt0");
  gt1 = getval("gt1");
  gt2 = getval("gt2");
  gt3 = getval("gt3");
  gt4 = getval("gt4");
  gzlvl0 = getval("gzlvl0");
  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl4 = getval("gzlvl4");

/* LOAD PHASE TABLE */

  settable(t1,1,phi1);
  settable(t2,2,phi2);
  settable(t3,4,phi3);
  settable(t4,4,phi4);
  settable(t5,4,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 > 48 )
    {
        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( pw > 200.0e-6 )
    {
        printf("dont fry the probe, pw too high ! ");
        abort(1);
    } 

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

    if( pwcodec < 275.0e-6 )
    {
        printf("dont fry the probe, pwcodec too high ! ");
        abort(1);
    } 

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

/*  Phase incrementation for hypercomplex 2D data */

    if (phase == 2)
      tsadd(t1,3,4);
    if (phase2 == 2)
      tsadd(t2,1,4);
    if (phase3 == 2)
      tsadd(t3,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 = 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*pwc - 2.0*pw - PRG_START_DELAY - PRG_STOP_DELAY - 2.0*POWER_DELAY - 4.0e-6 ); 
     else
        tau2 = ( tau2 - 4.0/PI*pwc - 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;

/*  Set up f3180  tau3 = t3               */

    tau3 = d4;

    if(f3180[A] == 'y') 
        tau3 += ( 1.0 / (2.0*sw3) - 4.0/PI*pwc - 2.0*pw - PRG_START_DELAY - PRG_STOP_DELAY - 2.0*POWER_DELAY - 4.0e-6 ); 
     else
        tau3 = ( tau3 - 4.0/PI*pwc - 2.0*pw - PRG_START_DELAY - PRG_STOP_DELAY - 2.0*POWER_DELAY - 4.0e-6 ); 

        if(tau3 < 0.2e-6)  tau3 = 4.0e-7;
                          
        tau3 = tau3/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(t5,2,4);      /* H */
    }

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

   if( ix == 1) d4_init = d4 ;
   t3_counter = (int) ( (d4-d4_init)*sw3 + 0.5 );
   if(t3_counter % 2) {
      tsadd(t3,2,4);  
      tsadd(t5,2,4);     /* C */ 
    }

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

/* Presaturation Period */

status(B);
   if (fsat[0] == 'y')
   {
        offset(tofps,TODEV);
	delay(2.0e-5);
        rgpulse(d1,zero,2.0e-6,2.0e-6);
   	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(t2);
   dec2phase(zero);
   delay(1.0e-5);

/* Begin Pulses */

status(C);

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

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

   decrgpulse(pwc,zero,0.0,0.0);

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

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

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

   txphase(zero); decphase(t2);
   delay(taua - PRG_START_DELAY + PRG_STOP_DELAY + 4.0*pwc + 4.0e-6 - gt1 - 4.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,zero,0.0,0.0);

   delay(tau2);

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

   delay(2.0e-6);

  }

  else {
    delay(PRG_START_DELAY - PRG_STOP_DELAY);
    simpulse(2.0*pw,2.0*pwc,zero,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(2.0e-6);

   decphase(one);

   delay(taua + tau1 - gt1 - 4.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);

   rgpulse(pw,zero,0.0,0.0);

   delay(mix/2.0);

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

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

   delay(mix/2.0 - gt2 - gt3 - 4.0e-6 - 500.0e-6);

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

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

   txphase(zero); decphase(t3);

   delay(taua - PRG_START_DELAY + PRG_STOP_DELAY - gt4 - 4.0e-6);

   decrgpulse(pwc,t3,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(tau3);

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

   delay(tau3);

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

   delay(2.0e-6);

  }

  else {
    delay(PRG_START_DELAY - PRG_STOP_DELAY);
    simpulse(2.0*pw,2.0*pwc,zero,zero,2.0e-6,2.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(taua - gt4 - 4.0e-6 - POWER_DELAY);

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

/* BEGIN ACQUISITION */

status(D);
setreceiver(t5);

}