/*  gtrosy_2D.c - auto-calibrated version of the original sequence

    Nietlspach 2D TROSY

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

    TROSY-based 3D hnca with deuterium decoupling
	F1 	N + JNH/2
	F-acq	HN - JNH/2

    This sequence uses the four channel configuration
         1)  1H             - carrier @ 4.7 ppm [H2O]
         2) 13C             - carrier @ 56 ppm (CA)
         3) 15N             - carrier @ 120 ppm  
         4) 2H		    - carrier @ 4.5 ppm 

    Set dm = 'nnn', dmm = 'ccc' 
    Set dm2 = 'nnn', dmm2 = 'ccc' [NO 15N decoupling during acquisition]
    Set dm3 = 'nnn', dmm2 = 'ccn' 

    Must set phase = 1,2 for States-TPPI
*/


#include <standard.h>
#define DELAY_BLANK 0.0e-6
#include "Pbox_bio.h"

static shape H2Osinc;

static int  phi1[2]  = {0,2},
            phi2[4]  = {0,0,2,2}, 
	    phi3[1]  = {0},
            rec[2]   = {0,2};

static double d2_init=0.0, d3_init=0.0;

pulsesequence()
{

/* DECLARE VARIABLES */

 char       C13refoc[MAXSTR],
            f1180[MAXSTR],    /* Flag to start t1 @ halfdwell             */
            autocal[MAXSTR];  /* auto-calibration flag */

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

 double      tau1,         /*  t1 delay */
             taua,         /*  ~ 1/4JNH =  2.25 ms */
             taub,         /*  ~ 1/4JNH =  2.75 ms */

             sw1,          /* sweep width in f1                    */             

	     compH=getval("compH"),
             pwHs=getval("pwHs"),   /* selective pulse on water      */
             tpwrs,                 /* power for pwHs               */
             tpwrsf_t=getval("tpwrsf_t"),
             phincr_t=getval("phincr_t"),
             tpwrsf_a=getval("tpwrsf_a"),
             phincr_a=getval("phincr_a"),

             compC,       /* C-13 RF calibration parameters */
             pwC,
             pwClvl,

             pwN,          /* PW90 for 15N pulse              */
	     pwNlvl,       /* high dec2 pwr for 15N hard pulses    */

             gstab,

             gt1,
             gt2,
             gt3,
             gt4,
             gt5,
             gt6,
             gt7,
             gt8,
             gt9,
             gt10,

             gzlvl1,
             gzlvl2,
             gzlvl3,
             gzlvl4,
             gzlvl5,
             gzlvl6,
             gzlvl7,
             gzlvl8,
             gzlvl9,
             gzlvl10;
            
/* LOAD VARIABLES */

  getstr("C13refoc",C13refoc);
  getstr("f1180",f1180);
  getstr("autocal",autocal);

  taua   = getval("taua"); 
  taub   = getval("taub"); 
  tpwr = getval("tpwr");
  dpwr = getval("dpwr");
  pwN = getval("pwN");
  pwNlvl = getval("pwNlvl");
  pwC = getval("pwC");
  pwClvl = getval("pwClvl");
  phase = (int) ( getval("phase") + 0.5);
  sw1 = getval("sw1");
  ni = getval("ni");

  gstab = getval("gstab");

  gt1 = getval("gt1");
  if (getval("gt2") > 0) gt2=getval("gt2");
    else gt2=gt1*0.1;
  gt3 = getval("gt3");
  gt5 = getval("gt5");
  gt6 = getval("gt6");
  gt7 = getval("gt7");

  gzlvl1 = getval("gzlvl1");
  gzlvl2 = getval("gzlvl2");
  gzlvl3 = getval("gzlvl3");
  gzlvl5 = getval("gzlvl5");
  gzlvl6 = getval("gzlvl6");
  gzlvl7 = getval("gzlvl7");

/* make small angle phase correction always positive */
  if (phincr_t<0.0) phincr_t = phincr_t +360.0;
  if (phincr_a<0.0) phincr_a = phincr_a +360.0;

  if(autocal[0]=='n')
  {
    if (pwHs > 1e-6) tpwrs = tpwr - 20.0*log10(pwHs/(compH*pw*1.69));
      else tpwrs = 0.0;
    tpwrs = (int) (tpwrs);
  }
  else
  {
    if (FIRST_FID)
    {
      H2Osinc = pbox_Rsh("H2Osinc", "sinc90", pwHs, 0.0, compH*pw, tpwr);
    }
    pwHs = H2Osinc.pw; tpwrs = H2Osinc.pwr-1.0;  /* 1dB correction applied */ 
  }   

/* LOAD PHASE TABLE */

  settable(t1,2,phi1);
  settable(t2,4,phi2);
  settable(t3,1,phi3);
  settable(t6,2,rec);

/* CHECK VALIDITY OF PARAMETER RANGES */

    if((dm[A] == 'y' || dm[B] == 'y' || dm[C] == 'y' ))
    {
        text_error("incorrect dec1 decoupler flags! Should be 'nnn' ");
        psg_abort(1);
    }

    if((dm2[A] == 'y' || dm2[B] == 'y' || dm2[C] == 'y'))
    {
        text_error("incorrect dec2 decoupler flags! Should be 'nnn' ");
        psg_abort(1);
    }

    if( dpwr > 46 )
    {
        text_error("don't fry the probe, DPWR too large!  ");
        psg_abort(1);
    }

    if( dpwr2 > 47 )
    {
        text_error("don't fry the probe, DPWR2 too large!  ");
        psg_abort(1);
    }

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

    if( gt1 > 15e-3 || gt2 > 15e-3 || gt3 > 15e-3 
	|| gt5 > 15e-3 || gt6 > 15e-3 || gt7 > 15e-3 )
    {
       text_error("gti values must be < 15e-3\n");
       psg_abort(1);
    } 

/*  Phase incrementation for hypercomplex 2D data */

   if (phase == 2) { tsadd(t3,2,4); icosel = -1.0; }
      else icosel = 1.0;

/*  Set up f1180  tau1 = t1               */
    tau1 = d2;
    if((f1180[A] == 'y') && (ni > 1.0))
        { tau1 += ( 1.0 / (2.0*sw1) ); if(tau1 < 0.2e-6) tau1 = 0.0; }
    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(t1,2,4);  
      tsadd(t6,2,4);    
    }

/* BEGIN ACTUAL PULSE SEQUENCE */

status(A);
   obspower(tpwr);     /* Set transmitter power for 1H presaturation */
   decpower(pwClvl);       /* Set Dec1 power for Ca 90 pulses         */
   dec2power(pwNlvl);     /* Set Dec2 power for 15N hard pulses         */
   txphase(zero);
   dec2phase(zero);

   delay(d1);
   rcvroff();

/* Begin Pulses */
   rgpulse(pw,zero,0.0,0.0);                    /* 90 deg 1H pulse */

   delay(0.2e-6);
   zgradpulse(gzlvl5,gt5);
   delay(2.0e-6);

   delay(taua -gt5 -2.2e-6);   /* taua <= 1/4JNH */ 

   sim3pulse(2.0*pw,0.0,2.0*pwN,zero,zero,zero,0.0,0.0);

   delay(taua -gt5 -200.2e-6); 

   delay(0.2e-6);
   zgradpulse(gzlvl5,gt5);
   delay(200.0e-6);
 
   rgpulse(pw,three,4.0e-6,0.0);
 
   /* shaped pulse */
   if (tpwrsf_t<4095.0) {obspower(tpwrs +6.0); obspwrf(tpwrsf_t);}
     else obspower(tpwrs);
   txphase(zero); 
   initval(1.0,v2); obsstepsize(phincr_t); xmtrphase(v2);
   shaped_pulse("H2Osinc",pwHs,zero,2.0e-6,2.0e-6);
   xmtrphase(zero);
   obspower(tpwr); if (tpwrsf_t<4095.0) obspwrf(4095.0);
   /* shaped pulse */

   txphase(zero);
   decphase(zero);
   dec2phase(t1);

   delay(0.2e-6);
   zgradpulse(gzlvl3,gt3);
   delay(200.0e-6);

   dec2rgpulse(pwN,t1,2.0e-6,0.0);
/* N15 EVOLUTION BEGINS */

   if ( (C13refoc[A]=='y') && ((tau1 -2.0*pwC -2.0e-6 -SAPS_DELAY)>0.0))
     {
      delay(tau1 -2.0*pwC -2.0e-6 -SAPS_DELAY);

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

      delay(tau1 -2.0*pwC -2.0e-6 -SAPS_DELAY);
     }
    else    delay(2.0*tau1);

    dec2phase(t2);
    delay(gt1 +gstab);

    dec2rgpulse(2.0*pwN, t2, 0.0, 0.0);

    zgradpulse(gzlvl1, gt1);          /* 2.0*GRADIENT_DELAY */

    dec2phase(t3);
    delay(gstab -2.0*GRADIENT_DELAY);

/* N15 EVOLUTION ENDS */
   sim3pulse(pw,0.0e-6,pwN,zero,zero,t3,0.0,0.0);

  /* shaped pulse */
  if (tpwrsf_a<4095.0) {obspower(tpwrs +6.0); obspwrf(tpwrsf_a);}
    else obspower(tpwrs);
  initval(1.0,v4); obsstepsize(phincr_a); xmtrphase(v4);
  shaped_pulse("H2Osinc",pwHs,zero,2.0e-6,2.0e-6);
  xmtrphase(zero);
  obspower(tpwr); if (tpwrsf_a<4095.0) obspwrf(4095.0);
  /* shaped pulse */

  delay(0.2e-6);
  zgradpulse(gzlvl6,gt6);
  delay(2.0e-6);
 
  dec2phase(zero);
  delay(taub -POWER_DELAY -WFG_START_DELAY -4.0e-6 -pwHs -WFG_STOP_DELAY
		-2.0*SAPS_DELAY -2.0e-6 -POWER_DELAY -gt6 -2.2e-6);
 
  sim3pulse(2.0*pw,0.0e-6,2.0*pwN,zero,zero,zero,0.0,0.0);
 
  txphase(one);
  dec2phase(zero);
  delay(taub - gt6 - 200.2e-6);
 
  delay(0.2e-6);
  zgradpulse(gzlvl6,gt6);
  delay(200.0e-6);

  sim3pulse(pw,0.0e-6,pwN,one,zero,zero,0.0,0.0);
 
  delay(0.2e-6);
  zgradpulse(gzlvl7,gt7);
  delay(2.0e-6);
 
  txphase(zero);
  dec2phase(zero);
  delay(taub - gt7 - 2.2e-6);
 
  sim3pulse(2.0*pw,0.0e-6,2.0*pwN,zero,zero,zero,0.0,0.0);

  txphase(one); 
  dec2phase(one);
  delay(taub - gt7 - 200.2e-6);
 
  delay(0.2e-6);
  zgradpulse(gzlvl7,gt7);
  delay(200.0e-6);

  sim3pulse(pw,0.0e-6,pwN,one,zero,one,0.0,0.0);
 
  txphase(zero);

  delay(gstab +gt2 -0.5*(pwN -pw) -2.0*pw/PI);
 
  rgpulse(2.0*pw,zero,0.0,0.0);
 
  delay(2.0e-6);
  zgradpulse(icosel*gzlvl2,gt2);
  decpower(dpwr);
  dec2power(dpwr2);
  delay(gstab -2.0e-6 -2.0*GRADIENT_DELAY -2.0*POWER_DELAY);

/* BEGIN ACQUISITION */
status(C);
  setreceiver(t6);
}