/*  gnoesyChsqc
*/

#include <standard.h>

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        f2180[MAXSTR];

 int         
             t1_counter,   /* used for states tppi in t1           */ 
             t2_counter;   /* used for states tppi in t2           */ 

double    

            pwC10 = getval("pwC10"),  /* 180  sinc pulse on CO(174ppm) */
            rf7,	    /* fine power for the pwC10 ("offC10") pulse */
            rf0,            /* full fine power */
            compC = getval("compC"), /*  C13 amplifier compression */
             tau1,tau1a,tau1b,         /*  t1 delay */
             tau2,tau2a,tau2b,         /*  t2 delay */
             jch,          /*  CH coupling constant */
             pwN,          /* PW90 for 15N pulse              */
             pwC,          /* PW90 for c nucleus @ pwClvl         */
             pwC180,       /* PW180 for c nucleus in INEPT transfers */
             pwClvl,        /* power level for 13C pulses on dec1  */
             pwNlvl,       /* high dec2 pwr for 15N hard pulses    */
             mix,          /* noesy mix time                       */
             sw1,          /* spectral width in t1 (H)             */
             sw2,          /* spectral width in t2 (C) (3D only)   */
             gt0,
             gt1,
             gt2,
             gt3,
             gt4,
             gt5,
             gt6,
             gt7,
             gzlvl0, 
             gzlvl1, 
             gzlvl2, 
             gzlvl3, 
             gzlvl4, 
             gzlvl5,
             gzlvl6,
             gzlvl7; 


/* LOAD VARIABLES */

  getstr("f2180",f2180);
  mix  = getval("mix");
  sw1  = getval("sw1");
  sw2  = getval("sw2");
  jch = getval("jch"); 
  pwC = getval("pwC");
  pwC180 = getval("pwC180");  
  pwN = getval("pwN");
  pwClvl = getval("pwClvl");
  pwNlvl = getval("pwNlvl");
  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[C] == 'y' ))
    {
        printf("incorrect 13C decoupler flags! dm='nnnn' or 'nnny' only  ");
        abort(1);
    }

    if((dm2[A] == 'y' || dm2[C] == 'y' ))
    {
        printf("incorrect 15N decoupler flags! No decoupling in relax or mix periods  ");
        abort(1);
    }

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

    if( dpwr2 > 49 )
    {
        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( pwN > 200.0e-6 )
    {
        printf("dont fry the probe, pwN too high ! ");
        abort(1);
    } 

    if( pwC > 200.0e-6 )
    {
        printf("dont fry the probe, pwC too high ! ");
        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);
     }

/*  Phase incrementation for hypercomplex 2D data */

    if (phase1 == 2)
      tsadd(t1,1,4);
    if (phase2 == 2)
      tsadd(t2,1,4);

/*  Set up f1180  tau1 = t1               */
   
    tau1 = d2;
    tau1a = 0.0; tau1b = 0.0;
    tau1a = d2 - (4.0/PI*pw + 4.0*pwC + 6.0e-6);
    tau1b = d2 - (4.0/PI*pw + 10.934*pwC + 6.0e-6);

    if(tau1b > 0.0) {
	tau1 = tau1b;
     }else if(tau1a >0.0) {
	tau1 = tau1a;
     }else{
        tau1=d2;
     }

    tau1 = tau1/2.0;

/*  Set up f2180  tau2 = t2               */

    tau2 = d3; 
    tau2a=0.0;
    tau2b=0.0;
     tau2a = d3 - (4.0/PI*pwC + 2*pw +4.0e-6);
     tau2b = d3 - (4.0/PI*pwC + 2*pwN +WFG2_START_DELAY +4.0e-6);
    if(f2180[A] == 'y')
     {
      tau2 +=  1.0 / (2.0*sw2); 
      tau2a +=  1.0 / (2.0*sw2); 
      tau2b +=  1.0 / (2.0*sw2); 
     }
     if(tau2b >0.0){
        tau2=tau2b;
       }else if(tau2a >0.0){
	tau2 = tau2a;
       }else{
        tau2 = d3; 
       }

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


    /* maximum fine power for pwC pulses */
	rf0 = 4095.0;

    /* "offC10": 180 degree one-lobe sinc pulse on CO, null at Ca 139ppm away */
	rf7 = (compC*4095.0*pwC*2.0*1.65)/pwC10;/* needs 1.65 times more*/
	rf7 = (int) (rf7 + 0.5);	/* power than a square pulse */

/* BEGIN ACTUAL PULSE SEQUENCE */


status(A);
   rcvroff();
   decoffset(dof);
   decpower(pwClvl);  /* Set Dec1 power for hard 13C pulses         */
   decpwrf(rf0);
   dec2power(pwNlvl);       /* Set Dec2 to low power for decoupling */
   obspower(tpwr);           /* Set transmitter power for hard 1H pulses */
   delay(d1);
status(B);
   rgpulse(pw,t1,1.0e-6,0.0);                  /* 90 deg 1H pulse */
    if(tau1b > 0.0) {
      delay(tau1);
      sim3pulse(0.0,pwC*158.0/90.0, 0.0,zero,zero,zero, 0.0, 0.0);/* Shaka6*/
      sim3pulse(0.0,pwC*171.2/90.0, pwN - pwC*171.4/90.0,zero,two,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC*342.8/90.0, pwC*342.8/90.0,zero,zero,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC*145.5/90.0, pwN - pwC*171.4/90.0,zero,two,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC*81.2/90.0, 0.0,zero,zero,zero, 0.0, 0.0);/* Shaka6*/
      sim3pulse(0.0,pwC*85.3/90.0, 0.0,zero,two,zero, 0.0, 0.0);/* Shaka6*/
      delay(tau1);
     }else if(tau1a >0.0) {
      delay(tau1);
      sim3pulse(0.0,pwC, pwN - pwC,zero,zero,zero, 0.0, 0.0);/*composite 180*/
      sim3pulse(0.0,2*pwC, 2*pwC,zero,one,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC, pwN - pwC,zero,zero,zero, 0.0, 0.0);
      delay(tau1);
     }else{
        delay(2*tau1);
     }

status(C);
   rgpulse(pw,zero,2.0e-6,0.0);             /*  2nd 1H 90 pulse   */

     delay(mix - 10.0e-3 ); 

   zgradpulse(gzlvl1, gt1);
   delay(10.0e-3 - gt1 - 8.0e-6);

   rgpulse(pw,zero,0.0,2.0e-6);
   zgradpulse(gzlvl2, gt2);       /* g3 in paper   */
   delay(1/(4.0*jch) - gt2 - 2.0e-6 -2*pwC);
      sim3pulse(0.0,pwC, 0.0,zero,zero,zero, 0.0, 0.0);/*composite 180*/
      sim3pulse(2*pw,2*pwC, 0.0,zero,one,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC, 0.0,zero,zero,zero, 0.0, 0.0);
   zgradpulse(gzlvl2, gt2);       /* g4 in paper  */
   delay(1/(4.0*jch) - gt2 - 2.0e-6 -2*pwC);
   rgpulse(pw,one,1.0e-6,2.0e-6);
   zgradpulse(gzlvl3, gt3);
   decrgpulse(pwC,t2,150.0e-6,2.0e-6);   
status(B);
   if( tau2b > 0.0 ) 
    {
        decphase(zero);
	decpwrf(rf7);
	delay(tau2);
	 sim3shaped_pulse("","offC10","",2.0*pw, pwC10,2*pwN, zero, zero, zero, 0.0, 0.0);
        decpwrf(rf0);
        delay(tau2);
    }else if(tau2a > 0.0 ){
	delay(tau2);
        rgpulse(2*pw,zero,2.0e-6,2.0e-6);  
        delay(tau2);
    }else{
	delay(tau2);
     simpulse(2*pw,2*pwC,zero,one,2.0e-6,2.0e-6);  
        delay(tau2);
    }
status(C);
   decrgpulse(pwC,zero,2.0e-6,2.0e-6);
   zgradpulse(gzlvl4, gt4);
   rgpulse(pw,t5,152.0e-6,2.0e-6);
   zgradpulse(gzlvl5, gt5);
   delay(1/(4.0*jch) - gt5 - 2.0e-6 -2*pwC);
      sim3pulse(0.0,pwC, 0.0,zero,zero,zero, 0.0, 0.0);/*composite 180*/
      sim3pulse(2*pw,2*pwC, 0.0,zero,one,zero, 0.0, 0.0);
      sim3pulse(0.0,pwC, 0.0,zero,zero,zero, 0.0, 0.0);
   zgradpulse(gzlvl5, gt5);
   delay(1/(4.0*jch) -  gt5 - 4.0e-6 -2*pwC);
  /*
   decrgpulse(pwC,zero,0.0,0.0);     
   decrgpulse(pwC,t3,2.0e-6,0.0);     
  */
   rgpulse(pw,t5,1.0e-6,rof2);                        /* flip-back pulse  */
   decpower(dpwr);
   dec2power(dpwr2);
   rcvron();
   setreceiver(t4);
 status(D);
}
/****************************End of Sulse Sequence************************/