/* 

noesy C13 hsqc with watergate to be used for aromatic residues only

ADD DETAILS HERE

*/


#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=0.0, pw_iC180=0.0;


    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 (dof is set to 125ppm - i.e. aromatic region */
    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);   
}