/* gnoesyNfhsqcA.c    "3D noesy Fast N15 hsqc using 3919 watergate suppression"

   program coded by Marco Tonelli and Klaas Hallenga, NMRFAM, January 2004

   for regular NOESY N15-HSQC minimum step cycling is 4
   for TROSY NOESY N15-HSQC minimum step cycling is 8
   
   Modified by Eriks Kupce :
   the TROSY phase cycle corrected to pick the correct component and to use the 
   sensitivity enhanced version. The minimum phase cycle for TROSY reduced to 4 steps.
   Corrected the d2 timing in the TROSY version. Added soft watergate option via
   wtg3919 flag.   
   Use f1coef = '1 0 -1 0 0 1 0 1' for TROSY

   C13 offset (dof) should be at 100ppm (jumped in sequence to 70ppm during t1)
   N15 offset in amides (120ppm)
   1H  offset on H2O

   tau1 timing corrected for regular experiment (4*pwN/PI correction added)
   tau1 timing corrected for TROSY experiment
   f1180 flag added for starting t1 at half-dwell
   C13shape flag added for chosing betweem adiabatic or composite 13C refocussing pulse in t1 
   (Marco Tonelli and Klaas Hallenga, NMRFAM, Univ. Wisconsin).

   Auto-calibrated version, E.Kupce, 27.08.2002.

   The autocal and checkofs flags are generated automatically in Pbox_bio.h
   If these flags do not exist in the parameter set, they are automatically 
   set to 'y' - yes. In order to change their default values, create the flag(s) 
   in your parameter set and change them as required. 
   The available options for the checkofs flag are: 'y' (yes) and 'n' (no). 
   The offset (tof, dof, dof2 and dof3) checks can be switched off individually 
   by setting the corresponding argument to zero (0.0).
   For the autocal flag the available options are: 'y' (yes - by default), 
   'q' (quiet mode - suppress Pbox output), 'r' (read from file, no new shapes 
   are created), 's' (semi-automatic mode - allows access to user defined 
   parameters) and 'n' (no - use full manual setup, equivalent to the original code).

   Timing fixed in H1 t1 dimension. Marco@NMRAFAM (Nov.2004).


Notes for noesy N15 hsqc:

   -mix is the mixing time
 
   -rp1 (zero order phase correction for H1 indirect dimension) should be 135 degrees

   WATERGATE OPTIONS:
   -tauWG is the delay for the 3919 watergate pulse train (wtg3919='y')
   -for using watergate with shaped pulses set wtg3919='n', set all the parameters 
    (important are pw, tpwr, pwHs) and fine powers and 
    small angle phase corrections for up and down flipback pulses (phincr_d, tpwrsf_d 
    and phincr_u, tpwrsf_u respectively).

   -to control water magnetization during H1 indirect evolution weak +/- gradients can
    be used. gt6 is the fixed length for these gradients, gzlvl6 is the power level.
    gstab is the recovery delay after the second gradient.
    gt6 should be not shorter than 500us and the power level, gzlvl6, should be low (400).
    Note. Using these gradients is not recomended on coldprobes as it can result in 
    distorion of baseline in H1 indirect dimension. To turn off these gradients set 
    gzlvl6 to zero.

*/


#include <standard.h>
#include "Pbox_bio.h"               /* Pbox Bio Pack Pulse Shaping Utilities */

static int phx[1]={0}, phy[1]={1}, ph_y[1]={3},

           phi2[2]  = {0,2},
           phi6[4]  = {0,0,2,2},
	   phi5[8]  = {0,0,0,0, 2,2,2,2},
           phi3[16] = {0,0,0,0, 0,0,0,0, 2,2,2,2, 2,2,2,2},
           phi4[32] = {0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
                       2,2,2,2, 2,2,2,2, 2,2,2,2, 2,2,2,2},

           rec[16]  = {2,0,0,2, 0,2,2,0, 0,2,2,0, 2,0,0,2};

static double   d2_init = 0.0, d3_init = 0.0;
static double   H1ofs=4.7, C13ofs=0.0, N15ofs=120.0, H2ofs=0.0;

static shape H2Osinc, iC180;

pulsesequence()
{
  int       t1_counter,
            t2_counter,
	    ni2=getval("ni2");

  char	    N15refoc[MAXSTR],	      /* C13 refocussing pulse in middle of t1 */
            wtg3919[MAXSTR],
	    shp_iC180[MAXSTR],	      /* calls sech/tanh pulses from shapelib */
	    f1180[MAXSTR],   		       /* Flag to start t1 @ halfdwell */
	    f2180[MAXSTR];   		       /* Flag to start t2 @ halfdwell */

  double    dofa,
            tauxh=getval("tauxh"), 
            tau1, tau2, corr,
	    tauWG=getval("tauWG"),
	    mix=getval("mix"),
            gzlvl0=getval("gzlvl0"),
            gzlvl3=getval("gzlvl3"),
            gzlvl4=getval("gzlvl4"),
            gzlvl5=getval("gzlvl5"),
            gzlvl6=getval("gzlvl6"),
            gt0=getval("gt0"),
            gt3=getval("gt3"),
            gt4=getval("gt4"),
            gt5=getval("gt5"),
            gt6=getval("gt6"),

            gstab=getval("gstab"),			/* gradient recovery delay */
            JNH = getval("JNH"),
            pwN = getval("pwN"),
            pwNlvl = getval("pwNlvl"),  

            pwHs, tpwrs=0.0, compH=1.0,          /* H1 90 degree pulse length at tpwrs */               
	    tpwrsf_u = getval("tpwrsf_u"),   /* fine power and small angle phase correction */
	    phincr_u = getval("phincr_u"),   /* for up (u) and down (d) H2O flipback pulses */
	    tpwrsf_d = getval("tpwrsf_d"),   /* of watergate pulse */
	    phincr_d = getval("phincr_d"),

  	    pwClw=getval("pwClw"),
  	    pwNlw=getval("pwNlw"),
  	    pwZlw = 0.0,                     /* largest of pwNlw and 2*pwClw */

            sw1 = getval("sw1"),
            sw2 = getval("sw2"),
                               /* temporary Pbox parameters */
            bw, pws, ofs, ppm, nst,  /* bandwidth, pulsewidth, offset, ppm, # steps */
            pwClvl = getval("pwClvl"), 	         /* coarse power for C13 pulse */
            pwC = getval("pwC"),       /* C13 90 degree pulse length at pwClvl */
            compC = getval("compC"),   /* adjustment for C13 amplifier compr-n */

	    pwr_iC180, pw_iC180;


    getstr("N15refoc",N15refoc);
    getstr("wtg3919",wtg3919);
    getstr("f1180",f1180);
    getstr("f2180",f2180);
    
/* check validity of parameter range */

    if((dm[A] == 'y' || dm[B] == 'y' ))
    { text_error("incorrect Dec1 decoupler flags!  "); psg_abort(1); } 

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

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

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

/* INITIALIZE VARIABLES */
    
    if(wtg3919[0] != 'y')      /* selective H20 one-lobe sinc pulse needs 1.69  */
    {                                   /* times more power than a square pulse */
      pwHs = getval("pwHs");            
      compH = getval("compH");
    }
    else 
      pwHs = pw*2.385+14.0e-6+tauWG*2.5; 

    setautocal();                        /* activate auto-calibration flags */ 
        
    if (autocal[0] == 'n') 
    {
    }
    else        /* if autocal = 'y'(yes), 'q'(quiet), r(read), or 's'(semi) */
    {
      strcpy(shp_iC180, "iC180_on");
      if(FIRST_FID)                                            /* call Pbox */
      {
        ppm = getval("dfrq"); ofs = 0.0; bw=0.0;
        if (find("bw_iC180")>0.0) bw=getval("bw_iC180");
        if (bw==0.0) bw=50.0*ppm; 
        iC180 = pbox_Rsh(shp_iC180, "isnob2", bw, ofs, compC*pwC, pwClvl); 

        if(wtg3919[0] != 'y')
          H2Osinc = pbox_Rsh("H2Osinc", "sinc90", pwHs, 0.0, compH*pw, tpwr);
        ofs_check(H1ofs, C13ofs, N15ofs, H2ofs);
      }
      pwr_iC180 = iC180.pwr;  pw_iC180 = iC180.pw;
      if (wtg3919[0] != 'y') 
      { 
        pwHs = H2Osinc.pw; tpwrs = H2Osinc.pwr-1.0;  /* 1dB correction applied */ 
      }
    }


/* calculate 3db lower power hard pulses for simultaneous CN decoupling during
   indirect H1 evoluion */

  if (N15refoc[A] == 'y')
    {
     if (pwNlw==0.0) pwNlw = pwN*exp(3.0*2.303/20.0);
     if (pwClw==0.0) pwClw = pwC*exp(3.0*2.303/20.0);
     if (pwNlw > 2.0*pwClw)
         pwZlw=pwNlw;
      else
         pwZlw=2.0*pwClw;
     if (d2==0.0 && d3==0.0) printf(" pwClw = %.2f ; pwNlw = %.2f\n", pwClw*1e6,pwNlw*1e6);
    }

/* LOAD VARIABLES */

    if(ix == 1) d2_init = d2;
    t1_counter = (int) ( (d2-d2_init)*sw1 + 0.5);

    if(ix == 1) d3_init = d3;
    t2_counter = (int) ( (d3-d3_init)*sw2 + 0.5);

    
/*  Set up f1180  */
   
    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;


/*  Set up f2180  */
   
    tau2 = d3;
    if((f2180[A] == 'y') && (ni2 > 1.0))
	{ tau2 += ( 1.0 / (2.0*sw2) ); if(tau2 < 0.2e-6) tau2 = 0.0; }
    tau2 = tau2/2.0;


/* LOAD PHASE TABLES */


    assign(one,v7); 
    assign(three,v8);

    settable(t1,  1, phy);
    settable(t2,  2, phi2);
    settable(t3, 16, phi3);
    settable(t4, 32, phi4);
    settable(t5,  8, phi5);
    settable(t6,  4, phi6);
    settable(t11, 16, rec); 

    if ( phase1 == 2 ) tsadd(t5, 1, 4);     /* Hypercomplex in t1 */
       
    if ( phase2 == 2 ) tsadd(t2, 1, 4);     /* Hypercomplex in t2 */
                                   
    if(t1_counter %2)          /* calculate modification to phases based on */
    { tsadd(t5,2,4); tsadd(t11,2,4); }   		/* current t1 values */

    if(t2_counter %2)          /* calculate modification to phases based on */
    { tsadd(t2,2,4); tsadd(t11,2,4); }   /* current t2 values */

    if(wtg3919[0] != 'y') 
    { add(one,v7,v7); add(one,v8,v8); }


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

    /* set dofa to 70ppm */
    ppm = getval("dfrq");
    dofa = dof -(55.0*ppm);
         
         
                           /* sequence starts!! */
   status(A);
     
     decoffset(dofa);
     obspower(tpwr);
     dec2power(pwNlvl);
     decpower(pwClvl);

/* small angle phase correction for H2O watergate flipback pulses */
     initval(phincr_d,v2);
     initval(phincr_u,v3);

     initval(135.0,v1);
     obsstepsize(1.0);
     txphase(t5);
     xmtrphase(v1);

     delay(d1);
     
   status(B);

     rgpulse(pw, t5, 0.0, 0.0);
/* H1 EVOLUTION BEGINS */
     txphase(t6); xmtrphase(zero);

     if ((N15refoc[A]=='y') && (tau1 > pwZlw +2.0*pw/PI +3.0*SAPS_DELAY +2.0*POWER_DELAY +2.0e-6))
       {
        decpower(pwClvl-3.0); dec2power(pwNlvl-3.0);

        corr = pwZlw +2.0*pw/PI +SAPS_DELAY +2.0*POWER_DELAY +2.0e-6;
        delay(tau1 -corr -2.0*SAPS_DELAY);

        if (pwNlw > 2.0*pwClw)
          {
           dec2rgpulse(pwNlw -2.0*pwClw,zero,0.0,0.0);
           sim3pulse(0.0,pwClw,pwClw,zero,zero,zero,0.0,0.0);
           decphase(one);
           sim3pulse(0.0,2*pwClw,2*pwClw,zero,one,zero,2.0e-6,0.0);
           decphase(zero);
           sim3pulse(0.0,pwClw,pwClw,zero,zero,zero,2.0e-6,0.0);
           dec2rgpulse(pwNlw -2.0*pwClw,zero,0.0,0.0);
          }
         else
          {
           decrgpulse(2.0*pwClw-pwNlw,zero,0.0,0.0);
           sim3pulse(0.0,pwNlw-pwClw,pwNlw-pwClw,zero,zero,zero,0.0,0.0);
           decphase(one);
           sim3pulse(0.0,2.0*pwClw,2.0*pwClw,zero,one,zero,2.0e-6,0.0);
           decphase(zero);
           sim3pulse(0.0,pwNlw-pwClw,pwNlw-pwClw,zero,zero,zero,2.0e-6,0.0);
           decrgpulse(2.0*pwClw-pwNlw,zero,0.0,0.0);
          }

        decpower(pwClvl); dec2power(pwNlvl);
        delay(tau1 -corr -SAPS_DELAY);
       }    

     else if (tau1 > (2.0*pwC +2.0*pw/PI +3.0*SAPS_DELAY +2.0e-6))
       {
	delay(tau1 -2.0*pwC -2.0*pw/PI -2.0e-6 -3.0*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.0*pw/PI -2.0e-6 -SAPS_DELAY);
       }
     else if (tau1 > 2.0*pw/PI +SAPS_DELAY)
	delay(2.0*tau1 -4.0*pw/PI -2.0*SAPS_DELAY);

/* H1 EVOLUTION ENDS */
     rgpulse(pw, t6, 0.0, 0.0);
/* NOESY MIXING TIME BEGINS */
     txphase(zero);
     decoffset(dof);
     delay(mix -gt0 -gstab);

     decrgpulse(pwC, zero, rof1, rof1);
     zgradpulse(gzlvl0,gt0);
     delay(gstab);

/* NOESY MIXING TIME ENDS */
     rgpulse(pw, zero, rof1, rof1);
/* CHSQC BEGINS */

     zgradpulse(gzlvl5,gt5);
     txphase(zero);
     decpower(pwr_iC180);
     decphase(zero);
     delay(tauxh -gt5 -WFG2_START_DELAY -0.5*pw_iC180);

     /* coupling evol reduced by 140us using stC pulses. Also WFG2_START_DELAY*/
     simshaped_pulse("", shp_iC180, 2.0*pw, pw_iC180, t4, zero, 0.0, 0.0);

     zgradpulse(gzlvl5,gt5);
     txphase(t1);
     decpower(pwClvl);
     decphase(t2);
     delay(tauxh -gt5 -0.5*pw_iC180);
  
     rgpulse(pw, t1, rof1, rof1);

     zgradpulse(gzlvl3,gt3);
     txphase(t4);      
     delay(gstab); 

     decrgpulse(pwC, t2, 2.0e-6, 2.0e-6);
/* C13 EVOLUTION BEGINS */
     if (ni2>1)
       {
        decphase(t3);
        if ((N15refoc[A]=='y') && (tau2 > 2.0*pwC/PI +pwN +SAPS_DELAY +rof1))
          {
          delay(tau2 -2.0*pwC/PI -pwN -SAPS_DELAY -2.0e-6);
          sim3pulse(2.0*pw, 0.0, 2.0*pwN, t4, zero, zero, 0.0, 0.0);
          delay(tau2 -2.0*pwC/PI -pwN -2.0e-6);
          }
        else if (tau2 > pw +2.0*pwC/PI +SAPS_DELAY +2.0e-6)
          {
          delay(tau2 -pw -2.0*pwC/PI -SAPS_DELAY -2.0e-6);
          rgpulse(2.0*pw, t4, 0.0, 0.0);
          delay(tau2 -2.0*pwC/PI -pw -2.0e-6);
          }
        else
          rgpulse(2.0*pw, t4, 2.0e-6, 2.0e-6);
       }
      else
       {
        decphase(zero);
        simpulse(2.0*pw, 2.0*pwC, t4, zero, 2.0e-6, 2.0e-6);
        decphase(t3);
       }
/* C13 EVOLUTION ENDS */
    decrgpulse(pwC, t3, 2.0e-6, 2.0e-6);
       
    zgradpulse(gzlvl4,gt4);
    delay(gstab);
    rgpulse(pw, two, 0.0, 0.0);
     
     zgradpulse(gzlvl6,gt6);
     decpower(pwr_iC180);
     txphase(v7); decphase(zero);
     delay(tauxh -gt6 -pwHs);
     
     if(wtg3919[0] == 'y')
     {     	
       rgpulse(pw*0.231,v7,2.0e-6,2.0e-6);     
       delay(tauWG);
       rgpulse(pw*0.692,v7,2.0e-6,2.0e-6);
       delay(tauWG);
       rgpulse(pw*1.462,v7,2.0e-6,2.0e-6);

       delay(tauWG/2.0 -0.5*pw_iC180 -WFG_START_DELAY);
       decshaped_pulse(shp_iC180, pw_iC180, zero, 0.0, 0.0);
       txphase(v8);
       delay(tauWG/2.0 -0.5*pw_iC180);

       rgpulse(pw*1.462,v8,2.0e-6,2.0e-6);
       delay(tauWG);
       rgpulse(pw*0.692,v8,2.0e-6,2.0e-6);
       delay(tauWG);
       rgpulse(pw*0.231,v8,2.0e-6,2.0e-6); 
     }
     else
     {
       xmtrphase(v2);
       if (tpwrsf_d < 4095.0)
	 {obspower(tpwrs+6.0); obspwrf(tpwrsf_d);}
       else	
          obspower(tpwrs); 
       shaped_pulse("H2Osinc", pwHs, v7, 2.0e-6, 0.0);

       xmtrphase(zero);
       obspower(tpwr); if (tpwrsf_d < 4095.0) obspwrf(4095.0);
       simpulse(2.0*pw, 2.0*pwC, v8, zero, 0.0, 0.0);

       xmtrphase(v3);
       if (tpwrsf_d < 4095.0)
         {obspower(tpwrs+6.0); obspwrf(tpwrsf_u);}
       else
          obspower(tpwrs);
       shaped_pulse("H2Osinc", pwHs, v7, 2.0e-6, 0.0);
       xmtrphase(zero);
       obspower(tpwr); if (tpwrsf_d < 4095.0) obspwrf(4095.0);
     } 
        
     zgradpulse(gzlvl6,gt6);   
     decpower(dpwr);
     delay(tauxh -gt6 -pwHs -POWER_DELAY); 

   status(C);
     setreceiver(t11);   
}