/*
Cross-term spatiotemporal encoding (xSPEN) sequence: Single-scan MRI with exceptional resilience to field heterogeneities
Reference: Z. Zhang, A. Seginer, L. Frydman,Magn. Reson. Med. DOI: 10.1002/mrm.26145 
Version: Zhiyong Zhang @20150924 
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, 1, 3, 0, 2, 1, 3};


pulsesequence()
{

        //Delay variables
        double tref,temin,tr_delay,
               te_delay=0,te_delay1=0, te_delay2=0,
               del=0, seqtime;		     
	double roff1,roff2;

        // Some xSPEN input variables

	double Rvalue=getval("Rvalue"),                        // time-bandwidth of the chirped sweep pulse
	       belta=getval("belta"),                          // variable to control the tilt of the flipped region in y-z plane, this affects the multi-slicing gap
                                                               // when belta = 0.5, the gap = thickness is the safty setting for multi-slicing experiments.
               Gysign=getval("Gysign"),                        // the sign for the bipolar Gy gradients, this controls the direction for reading the y axis 
                                                               // Gysign=1 means -Gy to Gy; Gysign=-1 means Gy to -Gy
	       Tf=getval("Tf"),                                // The delay to break the xSPEN's fully refocusing to create function contrast; Default Tf=0; no function contrast. 
               cpcor=getval("cpcor");			       // a power factor for the chirp pulse, default: cpcor=1
	char   Resmod[MAXSTR];                                 // Resolution module: Resmod='SR'  means super-resolution module, which is not supported yet
	       getstr("Resmod",Resmod);                        //                    Resmod='FT'  means the resolution is same as the regular Fourier imaging, it will automatic set Rvalue
                 					       //                    Resmod=''  means the resolution is controlled by Rvalue
        // xSPEN inter calculated variables
	double gamma,chirptof,Ta,Tp,BW,BWy,BWz,Gzenc,Gyenc,Gzd1=0,Gzd2=0;
	FILE *f1, *f2;
	int Npe=nv;
	char buf[256], *p;                                     //for update the Rvalue into pwbw (pulsewidth to bandwidth product) in WURST180;
	char cmd[MAXSTR];
	shape xSPENchirp180;

	SLICE_SELECT_GRADIENT_T Gz_enc,Gy_enc;
	GENERIC_GRADIENT_T Gz_per, Gz_acq, Gz_diff1, Gz_diff2;

	double  freq90[MAXNSLICE],chirpoffset1[MAXNSLICE],chirpoffset2[MAXNSLICE],ppes[MAXNSLICE];
	int     shapelist90,shapechirp1,shapechirp2;

        // EPI-based 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
	int vsegsign	 = v29;

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

        // Generate readout gradient and blip gradients (we don't use the blips), then we can get total acquistion time (Ta)

	init_readout(&epiro_grad,"xSPENgro",lro,np,sw);
	init_readout_refocus(&ror_grad,"xSPENgror");
	init_phase(&epipe_grad,"xSPENgpe",10*lpe,nv);    // the lpe here we apply 10 times to reduce Ta in case of small FOV along this axis
	init_phase(&per_grad,"xSPENper",10*lpe,nv);      // instead of blips, we use Gz to read the y image
	init_readout(&nav_grad,"xSPENgnav",lro,np,sw);
	init_epi(&epi_grad);
	fract_ky=nv/2;
	etl=nv;
	strcpy(ky_order,"l");
	calc_epi(&epi_grad,&epiro_grad,&epipe_grad,&ror_grad,&per_grad,&nav_grad,WRITE); 

	/* *************************************************** xSPEN main calculations start ******************************************************************* */ 
	if (belta<0 || belta>1) belta=0.5;
	if (!strcmp(Resmod,"FT"))
	{  
           Rvalue=nv/nseg*0.25/belta/(1-belta);
	   putCmd("Rvalue= %f",Rvalue); 
	}
		
	Ta=epiro_grad.duration;
	Tp=Ta/4/belta;
	BW=Rvalue/Tp;
	gamma=nuc_gamma();
	
	BWz=BW*(1-belta);
	BWy=BW*belta;

	//Slice selective pulse and gradient
  	init_rf(&p1_rf,p1pat,p1,flip1,rof1,rof2);
	p1=p1_rf.header.bandwidth/BWz;
	putCmd("p1= %f",p1*1e6); 
	init_rf(&p1_rf,p1pat,p1,flip1,rof1,rof2);

  	init_slice(&ss_grad,"ss",thk);
	calc_rf(&p1_rf,"tpwr1","tpwr1f"); 
  	calc_slice(&ss_grad,&p1_rf,WRITE,"");
  	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]);

	//Spatiotemporal encoding parameters

	Gzenc=BWz/gamma/(thk*0.1);
	Gyenc=BWy/gamma/lpe/nseg;
        
	
	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 xSPENchirp180.RF -w \"WURST180tmp %f \" -s 4.0 -p %f -l %f -header i -0", Tp,p3_rf.powerCoarse+20*log(p3_rf.powerFine/4095)/log(10),p3*cpcor*1e6/2);
	system(cmd);
	xSPENchirp180=getRsh("xSPENchirp180");
     	if (xSPENchirp180.pwr >50)  abort_message("Chirp power: %f maybe too high\n, please check", xSPENchirp180.pwr);

	for (i=0;i<nseg;i++)  ppes[i]=(-lpe/2+lpe/nseg/2+i*lpe/nseg+Gysign*ppe)*(1-(i%2)*2);
	offsetlist(ppes,-Gyenc*nseg*nseg,0,chirpoffset1,nseg,'c');
	offsetlist(ppes,Gyenc*nseg*nseg,0,chirpoffset2,nseg,'c');
	shapechirp1=shapelist(xSPENchirp180.name,xSPENchirp180.pw,chirpoffset1,nseg,1,'c');
	shapechirp2=shapelist(xSPENchirp180.name,xSPENchirp180.pw,chirpoffset2,nseg,1,'c');

	 
 	init_rf(&p2_rf,xSPENchirp180.name,Tp,180,rof1,rof2 );      
	init_slice_butterfly(&Gy_enc,"Gy",lpe*10/nseg,gcrush,tcrush);   
	init_slice_butterfly(&Gz_enc,"Gz",thk,gcrush,tcrush);                
	p2_rf.bandwidth=BWy;
	calc_slice(&Gy_enc,&p2_rf,WRITE,"Gyenc");	             
	p2_rf.bandwidth=BWz;
	calc_slice(&Gz_enc,&p2_rf,WRITE,"Gzenc");

	trampfixed=epi_grad.skip;
	init_generic(&Gz_acq,"Gzacq",Gzenc,epiro_grad.duration);
	calc_generic(&Gz_acq,WRITE,"Gzacq","tacq");

	trampfixed=0;
	init_dephase(&Gz_per,"Gzper"); 
	calc_dephase(&Gz_per,WRITE,Gz_acq.m0/2+ss_grad.m0ref,"",""); 
	/* *************************************************** xSPEN calculations finish ******************************************************************* */ 

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

	// spoil gradient if needed
	
	if (gspoil*tspoil!=0){
	   init_generic(&spoil_grad,"spoil",gspoil,tspoil);
	   calc_generic(&spoil_grad,WRITE,"","");
	}

	temin=2*(ss_grad.rfCenterBack+Gy_enc.duration+Gz_acq.duration/2+ror_grad.duration);

	te_delay=ss_grad.rfCenterBack+Gz_acq.duration/2+ror_grad.duration-Gz_per.duration;
        //diffusing gradient  
  	if (diff[0] == 'y') {
    	   init_generic(&diff_grad,"diff",gdiff,tdelta);
    	   diff_grad.maxGrad = gmax;
	   trampfixed=0;
           calc_generic(&diff_grad,NOWRITE,"","");
           diff_grad.duration += diff_grad.tramp; 
           calc_generic(&diff_grad,NOWRITE,"","");
	   init_dephase(&Gz_diff1,"Gzdiff1"); 
	   calc_dephase(&Gz_diff1,NOWRITE,Gz_per.m0/2+diff_grad.m0*dsl,"",""); 
	   init_dephase(&Gz_diff2,"Gzdiff2"); 
	   calc_dephase(&Gz_diff2,NOWRITE,Gz_per.m0/2-diff_grad.m0*dsl,"",""); 
	   tref=calc_sim_gradient(&diff_grad, &Gz_diff2, &Gz_diff1, 0, WRITE);
	   if (tref!=diff_grad.duration) abort_message("please increate the diffusing gradient time");
	   Gzd1=Gz_diff1.amp;
	   Gzd2=Gz_diff2.amp;

	   if ((Gz_per.m0/2+diff_grad.m0*dsl)<0) Gzd1=-Gz_diff1.amp;
	   if ((Gz_per.m0/2-diff_grad.m0*dsl)<0) Gzd2=-Gz_diff2.amp;

	   del=ss_grad.rfCenterBack+Gz_acq.duration/2+ror_grad.duration;
	   if ((tDELTA+tdelta+2*diff_grad.tramp)<del)     
	   	te_delay2=del-(tDELTA+tdelta+2*diff_grad.tramp)-te_delay;
	   else
	   {    
		te_delay1=((tDELTA+tdelta+2*diff_grad.tramp)-del);
	   	te_delay2=-te_delay;
	   	temin=temin+2*te_delay1;
	   }
	   putCmd("Gzd1= %f",Gzd1);
	   putCmd("Gzd2= %f",Gzd2);
         }


	   

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

	te_delay1 +=(te-temin)/2;
	te_delay2 +=(te-temin)/2+te_delay;

	if ((Tf>0) && (Tf<te_delay2))
           te_delay2=te_delay2-Tf;
	else if (Tf>=te_delay2)
	{   
            te_delay1=te_delay1+(Tf-te_delay2);
	    te_delay2=0;
	}

       	seqtime=ss_grad.rfCenterFront+te+Ta/2;
       	if (gspoil*tspoil!=0)  seqtime=seqtime+tspoil*2;
  	if (fsat[0]      == 'y') seqtime += fsatTime;

       	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-trmin)/ns;
	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);
	

 	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("ky_order= \'l\'");
	putCmd("BW= %f",BW);
	putCmd("BWy= %f",BWy);
	putCmd("BWz= %f",BWz);
	putCmd("Ta= %f",Ta);
	putCmd("Gzenc= %f",Gzenc);
	putCmd("Gyenc= %f",Gyenc);
	putCmd("Gro= %f",epiro_grad.roamp);
	putCmd("tror=%f",ror_grad.duration);
	putCmd("trortramp=%f",ror_grad.tramp);
	putCmd("gror=%f",ror_grad.amp); 
	putCmd("te_delay1= %f",te_delay1); 
	putCmd("te_delay2= %f",te_delay2); 
	putCmd("tpertramp= %f",Gz_per.tramp); 
	putCmd("tper= %f",Gz_per.duration); 
	putCmd("Gzper= %f",Gz_per.amp); 

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



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

        setacqmode(WACQ|NZ);
 	ifzero(rtonce); grad_advance(gpropdelay); endif(rtonce);

 	F_initval(epi_grad.etl/2, vpe); 
	F_initval(nseg, vnseg);
        loop(vnseg,vnseg_ctr); 

		mod2(vnseg_ctr,v11);  //0101...
		sub(one,v11,v11);     //1010...
		dbl(v11,v12);         //2020...
		sub(v12,one,vsegsign);//1,-1,1,-1...
	
		msloop(seqcon[1],ns,vms_slices,vms_ctr); 
			sp1on(); delay(4e-6); sp1off();
		      	if (ticks) {
			   xgate(ticks);
			   grad_advance(gpropdelay);
			   delay(4e-6);
			}

			rcvroff();
			txphase(0);


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

			S_position_offset_list(pss,ss_grad.ssamp,ns,resto,OBSch,shapelist90,vms_ctr);
			delay(GDELAY);
			obspower(p1_rf.powerCoarse);
		      	obspwrf(p1_rf.powerFine);
		      	delay(GDELAY);  
	      		obl_shapedgradient(ss_grad.name,ss_grad.duration,0.0,0.0,Gzenc,NOWAIT);
      			delay(ss_grad.rfDelayFront);
			shapedpulse(p1_rf.pulseName,p1_rf.rfDuration,v1,rof1,rof2);
			delay(ss_grad.rfDelayBack);

			delay(te_delay1/2);
			/* ------------- Spatiotemporal encoding---------------------------*/
			obspower(xSPENchirp180.pwr);
			obspwrf(xSPENchirp180.pwrf);	
		      	delay(GDELAY);

			nv=1;
			pe_shaped3gradient("",Gy_enc.name,Gz_enc.name,Gy_enc.duration,0.0,0.0,Gz_enc.amp,-Gysign*Gy_enc.amp,vsegsign,NOWAIT);
			delay(Gy_enc.rfDelayFront);
			shapedpulselist(shapechirp1,xSPENchirp180.pw,v2,rof1,rof2,'c',vnseg_ctr);
		      	delay(Gy_enc.rfDelayBack);
				
			delay(te_delay2/2);
			if (diff[0] != 'y')  				
				obl_shapedgradient(Gz_per.name,Gz_per.duration,0.0,0.0,Gz_per.amp,WAIT);

			else 
			{
				obl_shaped3gradient(diff_grad.name,diff_grad.name,Gz_diff1.name,diff_grad.duration,
				      diff_grad.amp*dro,diff_grad.amp*dpe,Gzd1,WAIT);
				delay(tDELTA-tdelta);
				obl_shaped3gradient(diff_grad.name,diff_grad.name,Gz_diff2.name,diff_grad.duration,
				      -diff_grad.amp*dro,-diff_grad.amp*dpe,Gzd2,WAIT);
			}	
			delay(te_delay2/2);	

			pe_shaped3gradient("",Gy_enc.name,Gz_enc.name,Gy_enc.duration,0.0,0.0,Gz_enc.amp,Gysign*Gy_enc.amp,vsegsign,NOWAIT);
		      	delay(Gy_enc.rfDelayFront);
			shapedpulselist(shapechirp2,xSPENchirp180.pw,v2,rof1,rof2,'c',vnseg_ctr);
		      	delay(Gy_enc.rfDelayBack);
			nv=Npe;

			delay(te_delay1/2);
			delay(GDELAY);
			/* ------------- xSPEN acquisition---------------------------*/
			obl_shapedgradient(ror_grad.name,ror_grad.duration,-ror_grad.amp,0,0,WAIT); 
			obl_shaped3gradient(epiro_grad.name,"",Gz_acq.name,epiro_grad.duration,epiro_grad.roamp,0.0,Gz_acq.amp,NOWAIT);
			delay(tep);
			nowait_loop(epi_grad.etl/2,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);                 /* end pe loop */
			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);
}
