/*
SPENms @ Zhiyong @ 20141010
contact: zhiyongxmu@gmail.com or zhiyong.zhang@weizmann.ac.il;
*/      
        
#include <standard.h>
#include <string.h>
#include "group.h"
#include <math.h>
#include <Pbox_psg.h>
#include "sgl.c"
        
static int ph90[8] ={0, 2, 1, 3, 0, 2, 1, 3},
           ph180[8]={0, 0, 0, 0, 2, 2, 2, 2},
	   phoph[8]={0, 2, 3, 1, 0, 2, 3, 1};


pulsesequence()
{

        //Delay variables
        double tref,temin,
               te_delay1, te_delay2, tr_delay, te_delay,
               del1=0, del2=0, del3=0, seqtime, 		//before and after diffusion gradients  
               busy1, busy2;                                    //time spent on rf pulses etc. in TE periods      
	double roff1,roff2;

	double Ta,Tp,BW,Npe,Rvalue=getval("Rvalue"),cpcor=getval("cpcor"),piflag=getval("piflag");
	FILE *f1, *f2;
	char buf[256], *p;                                     //for update the Rvalue into pwbw (pulsewidth to bandwidth product) in WURST180;
	char cmd[MAXSTR];

	shape SPENchirp180;

	SLICE_SELECT_GRADIENT_T chirp180_grad;
	double gamma,Gchirp,tchirp,chirptof;
	double gcrush3=getval("gcrush3"),tcrush3=getval("tcrush3");

	double  freq90[MAXNSLICE];
	int     shapelist90;

        // Diffusion variables
        #define MAXDIR 1024           			// Will anybody do more than 1024 directions or b-values? 
        int    	diff_in_one = 0;
        double 	tmp, tmp_ss2;
        double 	roarr[MAXDIR], pearr[MAXDIR], slarr[MAXDIR];
        int    	nbval,               			// Total number of bvalues*directions 
               	nbro, nbpe, nbsl;    			// bvalues*directions along RO, PE, and SL 
        double 	bro[MAXDIR], bpe[MAXDIR], bsl[MAXDIR], 	// b-values along RO, PE, SL 
         	brs[MAXDIR], brp[MAXDIR], bsp[MAXDIR], 	// the cross-terms 
	 	btrace[MAXDIR],                       	// and the trace 
	 	max_bval=0,
         	dcrush, dgss2,       			// "delta" for crusher and gss2 gradients 
         	Dro, Dcrush, Dgss2;  			// "DELTA" for readout, crusher and gss2 gradients 

	int i;						//loop variable

	int vpe          = v21;   
	int vpe_ctr      = v22;
	int vgro         = v23;   
  	int vms_slices   = v25;   			// Number of slices
  	int vms_ctr      = v26;   			// Slice loop counter
  	int vnseg        = v27;   			// Number of segments
  	int vnseg_ctr    = v28;   			// Segment loop counter


  	get_parameters();
  	euler_test();
	
	if (tep < 0) { 					// adjust by reducing gpropdelay by that amount
          gpropdelay += tep;
          tep = 0;
        }
	setacqmode(WACQ|NZ);


	if (!strcmp(orient,"oblique")) {
	  if ((phi != 90) || (psi != 90) || (theta != 90)) {
	    //oblique slice - this should take care of most cases
	    epiro_grad.slewRate /= 3; // = gmax/trise 
	    epipe_grad.slewRate /= 3; 
	  }
	}

        //preproduce readout gradient and blip gradients to get Ta

	init_readout(&ro_grad,"SPEN180_ro",lro,np,sw);
	init_readout_refocus(&ror_grad,"SPEN180_ror");
	init_phase(&pe_grad,"SPEN180_pe",lpe,nv);
	init_phase(&per_grad,"SPEN180_per",lpe,nv);
	init_readout(&nav_grad,"nav",lro,np,sw);
	init_epi(&epi_grad);
	fract_ky=nv/2;
	etl=nv;
	strcpy(rampsamp,"n");
	strcpy(ky_order,"l");
	calc_epi(&epi_grad,&ro_grad,&pe_grad,&ror_grad,&per_grad,&nav_grad,NOWRITE); 

        //According to Ta, calculate the chirp 180 pulse related parameter (BW, Tp, Gchirp)
	Ta=ro_grad.duration;
	Tp=Ta/2;
	BW=Rvalue/Tp;
	gamma=nuc_gamma();
	Gchirp=BW/gamma/lpe;
	chirptof=poffset(ppe,Gchirp);
        //initialize a temp RF struture for calculating the chirp pulse


	
	sprintf(cmd,"%s/wavelib/decoupling/%s", userdir,"WURST180");
	f1=fopen(cmd,"r");
	sprintf(cmd,"%s/wavelib/decoupling/%s", userdir,"WURST180tmp");
	f2=fopen(cmd,"w");

	while (1)
	{	
            if (NULL == fgets(buf,255,f1)) break;
            if ((p = strstr(buf,"pwbw 	= 92"))) sprintf(buf,"pwbw 	= %f	(pulsewidth to bandwidth product)\n",Rvalue*0.8574);
             fprintf(f2,"%s",buf);
        }
        fclose(f1);
	fclose(f2);
             
	init_rf(&p3_rf,"hard",p3/2,90,rof1,rof2); //temp RF struture for calculating the chirp pulse
  	calc_rf(&p3_rf,"",""); 
	sprintf(cmd, "Pbox SPENchirp180.RF -w \"WURST180tmp %f  %f\" -s 4.0 -p %f -l %f -header i -0", Tp,chirptof,p3_rf.powerCoarse+20*log(p3_rf.powerFine/4095)/log(10),p3*cpcor*1e6/2);
	system(cmd);
	SPENchirp180=getRsh("SPENchirp180");
     	if (SPENchirp180.pwr >50)  abort_message("Chirp power: %f maybe too high\n, please check", SPENchirp180.pwr);


  	init_rf(&p2_rf,SPENchirp180.name,Tp,180,rof1,rof2 );    
	init_slice_butterfly(&chirp180_grad,"SPEN180_grad",lpe*10,gcrush,tcrush);                  
	p2_rf.bandwidth=BW;
  	calc_slice(&chirp180_grad,&p2_rf,WRITE,"gss");

	init_rf(&p3_rf,p3pat,p3,180,rof1,rof2);
  	calc_rf(&p3_rf,"tpwr3","tpwr3f"); 

        //According to the Gchirp, update the readout and blip gradients
	init_readout(&ro_grad,"SPEN180_ro",lro,np,sw);
	init_readout_refocus(&ror_grad,"SPEN180_ror");
	init_phase(&pe_grad,"SPEN180_pe",nv/gamma/Gchirp/Tp/2,nv);
	init_phase(&per_grad,"SPEN180_per",nv/gamma/Gchirp/Tp/2,nv);
	init_readout(&nav_grad,"nav",lro,np,sw);
	init_epi(&epi_grad);
	calc_epi(&epi_grad,&ro_grad,&pe_grad,&ror_grad,&per_grad,&nav_grad,NOWRITE); 

  	roff1 = -poffset(pro,ro_grad.roamp);
  	roff2 = -poffset(pro,-ro_grad.roamp);	

	//Slice selective pulse and gradient
  	init_rf(&p1_rf,p2pat,p1,flip1,rof1,rof2);
  	init_slice(&ss_grad,"ss",thk);
  	init_slice_refocus(&ssr_grad,"ssr");  
	calc_rf(&p1_rf,"tpwr1","tpwr1f"); 
  	calc_slice(&ss_grad,&p1_rf,WRITE,"");
	calc_slice_refocus(&ssr_grad,&ss_grad,NOWRITE,"gssr");
  	offsetlist(pss,ss_grad.ssamp,0,freq90,ns,seqcon[1]);
  	shapelist90 = shapelist(p1_rf.pulseName,ss_grad.rfDuration,freq90,ns,ss_grad.rfFraction,seqcon[1]);
      
	tref=calc_sim_gradient(&ssr_grad, &per_grad, &ror_grad, 0, WRITE);
	
	//Some spoil gradient and crush gradient
	trampfixed=0;
	if (gspoil*tspoil!=0){
	   init_generic(&spoil_grad,"spoil",gspoil,tspoil);
	   calc_generic(&spoil_grad,WRITE,"","");
	}
	
	if (gcrush3*tcrush3!=0)
	   trapezoid(&crush_grad,"gcrush",gcrush3,tcrush3,0,WRITE);
	else 
	   trapezoid(&crush_grad,"gcrush",5,0.002,0,WRITE);

        //diffusing gradient
  	if (diff[0] == 'y') {
    	   init_generic(&diff_grad,"diff",gdiff,tdelta);
    	   diff_grad.maxGrad = gmax;
           calc_generic(&diff_grad,NOWRITE,"","");
          // adjust duration, so tdelta is from start ramp up to start ramp down */
            diff_grad.duration += diff_grad.tramp; 
            calc_generic(&diff_grad,WRITE,"","");
         }


	  if (diff[0] == 'y') {  /* CALCULATE B VALUES */
	    /* Get multiplication factors and make sure they have same # elements */
	    /* All this is only necessary because putCmd only work for ix==1      */
	    nbro = (int) getarray("dro",roarr);  nbval = nbro; 
	    nbpe = (int) getarray("dpe",pearr);  if (nbpe > nbval) nbval = nbpe;
	    nbsl = (int) getarray("dsl",slarr);  if (nbsl > nbval) nbval = nbsl;
	    if ((nbro != nbval) && (nbro != 1))
	      abort_message("%s: Number of directions/b-values must be the same for all axes (readout)",seqfil);
	    if ((nbpe != nbval) && (nbpe != 1))
	      abort_message("%s: Number of directions/b-values must be the same for all axes (phase)",seqfil);
	    if ((nbsl != nbval) && (nbsl != 1))
	      abort_message("%s: Number of directions/b-values must be the same for all axes (slice)",seqfil);

	    if (nbro == 1) for (i = 1; i < nbval; i++) roarr[i] = roarr[0];
	    if (nbpe == 1) for (i = 1; i < nbval; i++) pearr[i] = pearr[0];
	    if (nbsl == 1) for (i = 1; i < nbval; i++) slarr[i] = slarr[0];
	  }
	  else {
	    nbval = 1;
	    roarr[0] = 0;
	    pearr[0] = 0;
	    slarr[0] = 0;
	  }

        if (piflag==0 || piflag==-1) {
	   busy1=ss_grad.rfCenterBack+ssr_grad.duration+chirp180_grad.duration/2;
	   busy2=chirp180_grad.duration/2+tep+Ta/2;

           if (diff[0]=='y') {
              if (tDELTA < (diff_grad.duration + GDELAY))        	//Minimum tDELTA 
                 abort_message("ERROR %s: tDELTA is too short, minimum is %.2fms\n",
	             seqfil,(diff_grad.duration + GDELAY)*1000+0.005);
              busy1=busy1+tDELTA+diff_grad.duration+GDELAY;
	      del1=tDELTA-diff_grad.duration;
           }

	   temin=MAX(busy1,busy2)*2;
	   if (minte[0]=='y')  te=temin;
	   if(te<temin)  abort_message("ERROR %s: te is too short, minimum is %fms\n",seqfil, temin*1000); 


	   del2=te/2-busy1;
	   del3=te/2-busy2;
    	   
	// Diffusing b value caculation
 	   double tarr[20],ty[1024];

	   tarr[0]=0;
           tarr[1]=tarr[0]+ss_grad.duration/2;
           tarr[2]=tarr[1]+ssr_grad.duration;
           tarr[3]=tarr[2]+diff_grad.duration;
           tarr[4]=tarr[3]+tDELTA-diff_grad.duration;
	   tarr[5]=tarr[4]+diff_grad.duration;
           tarr[6]=tarr[5]+del2;
           tarr[7]=tarr[6]+chirp180_grad.rfDelayFront;
	   tarr[8]=1000;
	   tarr[9]=tarr[7]+SPENchirp180.pw;
	   tarr[10]=tarr[9]+chirp180_grad.rfDelayBack;
	   tarr[11]=tarr[10]+del3;
	   tarr[12]=tarr[11]+at/2+epi_grad.skip;
	   tarr[13]=tarr[12]+nv*epi_grad.skip*2;
	
	   double Gxarr[13]={0, ror_grad.amp, gdiff, 0 ,-gdiff, 0 ,0, 0 ,0 ,0 ,ro_grad.amp,0};
	   double Gyarr[13]={0, per_grad.amp, gdiff, 0 ,-gdiff, 0 ,gcrush, Gchirp, Gchirp,gcrush ,0 ,0, pe_grad.amp};
	   double Gzarr[13]={ss_grad.amp,-ss_grad.amp*tarr[1]/(tarr[2]-tarr[1]), gdiff, 0 ,-gdiff, 0 ,0 ,0 ,0, 0, 0 ,0,0};

 	   for (i=0;i<nv;i++)
             ty[i]=(gamma*Gchirp*(i-nv/2)*lpe/nv+BW/2)/(BW/Tp);

	   putarray("tarr",tarr,14);
	   putarray("ty",ty,nv);
	   putarray("Gxarr",Gxarr,13);
	   putarray("Gyarr",Gyarr,13);
	   putarray("Gzarr",Gzarr,13);

	}

	else if (piflag==1) {
           busy1=ss_grad.rfCenterBack+ssr_grad.duration+chirp180_grad.duration/2+crush_grad.duration+p3/2;
           busy2=Ta/2+tep+crush_grad.duration+p3/2;

           if (diff[0]=='y') {
              if (tDELTA < (diff_grad.duration + GDELAY))        	//Minimum tDELTA 
                 abort_message("ERROR %s: tDELTA is too short, minimum is %.2fms\n",
	             seqfil,(diff_grad.duration + GDELAY)*1000+0.005);
              if (tDELTA + diff_grad.duration > Ta/2+GDELAY) 
                 busy2=busy2+tDELTA+diff_grad.duration+GDELAY-Ta/2;
	   }
	   temin=MAX(busy1,busy2)*2;
	   if (minte[0]=='y')  te=temin;
	   if(te<temin)  abort_message("ERROR %s: te is too short, minimum is %fms\n",seqfil, temin*1000); 

	   del1=tDELTA-diff_grad.duration;
	   del2=te/2-busy1+Ta/2;
	   del3=te/2-busy2;
    	   
	   te_delay1=te/2-busy1;
	   te_delay2=te/2-busy2-Ta/2;
	}

       seqtime=ss_grad.rfCenterFront+te+Ta/2;
       if (piflag==-1)  seqtime=seqtime+tspoil*2+p3;
       trmin=seqtime*ns; 

       if (mintr[0]=='y') tr=trmin;

       if (tr <trmin)  abort_message("%s: Requested tr too short.  Min tr = %.2f ms\n", seqfil,ceil(trmin*100000)/100.00);

	tr_delay=tr/ns-tspoil*2-p3-ss_grad.rfCenterFront-te-Ta/2;

	sgl_error_check(sglerror);


	settable(t1,8,ph90);
	settable(t2,8,ph180);
	settable(t3,8,phoph);

	getelem(t1,ct,v1);
	getelem(t2,ct,v2);
	getelem(t3,ct,oph);
	


	Npe=(int) nv/2/nseg;

	F_initval(nseg, vnseg);

 	get_wsparameters();
	if (ws[0] == 'y') create_watersuppress();
  	if (fsat[0] == 'y')  create_fatsat();

	putCmd("tr= %f",tr); 
	putCmd("te= %f",te); 
	putCmd("aqtm= %f",at);
	putCmd("fract_ky= %f",nv/2);
	putCmd("etl= %f",nv);
	putCmd("rampsamp= \'n\'");
	putCmd("ky_order= \'l\'");
	putCmd("BW= %f",BW);
	putCmd("Ta= %f",Ta);
	putCmd("Tp= %f",Tp);
	putCmd("Gchirp= %f",Gchirp);


	/*   =======================================
	     ==    Beginning of Pulse Sequence    ==
	     =======================================   */
        status(A);
  	rotate();
	obsoffset(resto); delay(GDELAY);

	if (ws[0]   == 'y') watersuppress();
        if (fsat[0] == 'y') fatsat();


	F_initval(Npe,vpe);
        setacqmode(WACQ|NZ);
 	ifzero(rtonce); grad_advance(gpropdelay); endif(rtonce);


        loop(vnseg,vnseg_ctr); 
             msloop(seqcon[1],ns,vms_slices,vms_ctr); 


      	if (ticks) {
           xgate(ticks);
           grad_advance(gpropdelay);
           delay(4e-6);
        }

	rcvroff();
	txphase(0);


	/* -----------------------------------------
	      Excitation
	   ----------------------------------------- */


	obsoffset(resto);
	obspower(p1_rf.powerCoarse);
      	obspwrf(p1_rf.powerFine);
      	delay(4e-6);
      	obl_shapedgradient(ss_grad.name,ss_grad.duration,0,0,ss_grad.amp,NOWAIT);
      	delay(ss_grad.rfDelayFront);     
	shapedpulselist(shapelist90,p1_rf.rfDuration,v1,rof1,rof2,seqcon[1],vms_ctr);
      	delay(ss_grad.rfDelayBack);


	pe3_shaped3gradient(ror_grad.name,per_grad.name,ssr_grad.name,ssr_grad.duration,ror_grad.amp,per_grad.amp,-ssr_grad.amp,0,-per_grad.amp*pe_grad.m0/per_grad.m0/nseg,0,zero,vnseg_ctr,zero,WAIT);

        delay(GDELAY);

        if (diff[0] == 'y') {
           obl_shapedgradient(diff_grad.name,diff_grad.duration,diff_grad.amp*dro,diff_grad.amp*dpe,diff_grad.amp*dsl,WAIT);
           delay(del1);
	   obl_shapedgradient(diff_grad.name,diff_grad.duration,-diff_grad.amp*dro,-diff_grad.amp*dpe,-diff_grad.amp*dsl,WAIT);
	}

	delay(del2);
	obsoffset(resto);
	obspower(SPENchirp180.pwr);
	obspwrf(SPENchirp180.pwrf);					
      	delay(4e-6);
	obl_shaped3gradient("",chirp180_grad.name,"",chirp180_grad.duration,0,chirp180_grad.amp,0,NOWAIT); 
      	delay(chirp180_grad.rfDelayFront);
	shaped_pulse(SPENchirp180.name,SPENchirp180.pw,v2,rof1,rof2);
      	delay(chirp180_grad.rfDelayBack);

	delay(del3);
	delay(GDELAY);
	obl_shaped3gradient(ro_grad.name,pe_grad.name,"",ro_grad.duration,ro_grad.roamp,pe_grad.amp,0.0,NOWAIT);
	delay(tep);

	
	nowait_loop(Npe,vpe,vpe_ctr);

    	roff = roff1;  
	delay(epi_grad.skip-alfa); 
	startacq(alfa);
	acquire(np,1/sw);
	delay((epi_grad.dwell[(int) (np/2-1)]-1/sw));
	if (aqtm > at) sample(aqtm-at);
	endacq();
	delay(epi_grad.skip - 1/sw - (aqtm-at));

    	roff = roff2;   
	delay(epi_grad.skip-alfa);
	startacq(alfa);
	acquire(np,1/sw);
	delay((epi_grad.dwell[(int) (np/2-1)]-1/sw));
	if (aqtm > at) sample(aqtm-at);
	endacq();
	delay(epi_grad.skip - 1/sw - (aqtm-at));

	nowait_endloop(vpe_ctr);


	if(piflag==-1) {

	obspower(SPENchirp180.pwr);
	obspwrf(SPENchirp180.pwrf);					
      	delay(4e-6);
	obl_shaped3gradient("",chirp180_grad.name,"",chirp180_grad.duration,0,-chirp180_grad.amp,0,NOWAIT); 
      	delay(chirp180_grad.rfDelayFront);
	shaped_pulse(SPENchirp180.name,SPENchirp180.pw,v2,rof1,rof2);
      	delay(chirp180_grad.rfDelayBack);

/*
	delay(GDELAY);
	obl_shaped3gradient(crush_grad.name,crush_grad.name,crush_grad.name,crush_grad.duration,crush_grad.amp,crush_grad.amp,crush_grad.amp,WAIT);
	obspower(p3_rf.powerCoarse);
      	obspwrf(p3_rf.powerFine);
	shaped_pulse(p3_rf.pulseName,p3,zero,rof1,rof2);
	obl_shaped3gradient(crush_grad.name,crush_grad.name,crush_grad.name,crush_grad.duration,crush_grad.amp,crush_grad.amp,crush_grad.amp,WAIT);
*/
	}

	delay(tr_delay);

    	   	endmsloop(seqcon[1],vms_ctr);   		 /* end multislice loop */
  	endloop(vnseg_ctr); 
	loop_check();               			 	/* end segments loop */
	g_setExpTime(tr*nt*nseg*arraydim+tr*ssc);
}
