! ------------------------------------------
! TIME INCREMENTS.
!
! This Set of subroutines keeps track of evolution times.
!
! CMR Lab...UIUC...Trent Franks...4-22-04
! 
! Updated to be consistent with ants.csh 
! 
! WTF updated 12-30-04
! 1. AL start parameters added
! 2. Initial start times are calculated.
!  f
! WTF updated 01-05-06
! 1. Constant time evolution added
! 
! WTF updated 051407
! 1. Constant time with soft pulse, Improved J decoupling.
! ------------------------------------------

include "../includes/t1.inc";

!---- t2
    .time   extern   t2_evol     = 100u;
    .time            t2_evol_out = 100u;
    .time   extern   t2_evol_max = 100u;
    .time   extern   t2_evol_max2= 100u;
    .time            t2_evol_x   = 100u;
    .time            cyclerm2 = 1.0u;
    .time            t2by2  = 1.0u;
    .time            t2by2x = 1.0u; 
    .time            t2by2y = 1.0u; 

!---- T2 TPPM Parameters
    .float t2cycles "no. of TPPM cycles during t2" = 8; 
    .long  cyclest2 ;

!---- t2 constant time (echo recoupling)
    .time   extern   t2_ct   = 100u;
    .time            t2_a    = 100u;
    .time            t2aby2  = 1.0u;
    .time            t2aby2x = 1.0u; 
    .time            t2aby2y = 1.0u; 

    .time            t2_b    = 100u;
    .time            t2bby2  = 1.0u;
    .time            t2bby2x = 1.0u; 
    .time            t2bby2y = 1.0u; 

    .long            cyclest2a ;
    .long            cyclest2b ;
    .time            cyclerm2a = 1.0u;
    .time            cyclerm2b = 1.0u;

!---- Add to update section (uncomment/adjust variable names)
!	.update "t2_evol= dwell - (2.0/3.1415)*offset";
!	.update "t2_evol_max = t2_evol + dwell * (AL - 1)";  
!	.update "t2_evol_max2 = dwell * (AL_max)"; 

DEFINE t2_prog(dwell,AL,ALmax,offset,PW)

	t2_evol= dwell - (2.0/3.14159) * offset;

	t2_evol_max = t2_evol + dwell*(AL-1);
     	if(AL<ALmax){t2_evol_max2 = dwell*(ALmax);}
	else{t2_evol_max2 = dwell*(AL);}
	
	t2cycles = (t2_evol)/(PW);
	cyclest2 = floor(t2cycles);
	
	cyclerm2 = t2_evol - cyclest2*PW;
	if(cyclerm2<1.0u){
	cyclest2 = cyclest2 - 1;
	cyclerm2 = t2_evol - cyclest2*PW;
        }
ENDDEF

DEFINE t2_EVOLVE(dwell,ALstart,INDEX,PW_Evol)       

    	t2_evol_out = t2_evol + (dwell*(INDEX+ALstart-1));
	cyclest2 = ceil(t2_evol_out/PW_Evol)-1;
	if (cyclest2 < 0){cyclest2 = 0;}
	cyclerm2 = t2_evol_out - cyclest2*PW_Evol;
	
        if(cyclerm2<1.0u){
	cyclest2 = cyclest2 - 1;
	if (cyclest2 < 0){cyclest2 = 0;}
	cyclerm2 = t2_evol_out - cyclest2*PW_Evol;
        }         !making sure cyclerm2 isn't too short
    
	t2by2 = t2_evol_out / 2.0;
	if(t2by2>(2.0*PW180X+1u)){ 
	    t2by2x = t2by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2x=t2by2;}
	if(t2by2>(2.0*PW180Y+1u)){ 
	    t2by2y = t2by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2y=t2by2;}

        t2_evol_x = t2_evol_max2 + 10u - t2_evol_out;

	if(t2_evol_x<3.0u){t2_evol_x=3u;}
ENDDEF

DEFINE t2_HOMO_DEC_EVOLVE(dwell,ALstart,INDEX,PW_Evol,PW_homodec)       

    	t2_evol_out = (t2_evol + (dwell*(INDEX+ALstart-1)))/2.0;
	cyclest2 = ceil(t2_evol_out/PW_Evol)-1;
	if (cyclest2 < 0){cyclest2 = 0;}
	cyclerm2 = t2_evol_out - cyclest2*PW_Evol;
	
        if(cyclerm2<1.0u){
	cyclest2 = cyclest2 - 1;
	if (cyclest2 < 0){cyclest2 = 0;}
	cyclerm2 = t2_evol_out - cyclest2*PW_Evol;
        }         !making sure cyclerm2 isn't too short
    
	t2by2 = t2_evol_out;

	if((t2by2+(PW_homodec/2.0))>(2.0*PW180X+1u )){ 
	    t2by2x = t2by2 - PW180X; ! 2 pi's for composite,1 per half
        } 
        else{t2by2x=t2by2;}
	
	
	if((t2by2+(PW_homodec/2.0))>(2.0*PW180Y+1u )){ 
	t2by2y = t2by2 - PW180Y; ! 2 pi's for composite,1 per half
        } 
        else{t2by2y=t2by2;}

        t2_evol_x = t2_evol_max2 + 10u - 2.0*t2_evol_out;

	if(t2_evol_x<3.0u){t2_evol_x=3u;}
ENDDEF


DEFINE t2_PROG(dwell,AL,ALmax,PW)
        t2_prog(dwell,AL,ALmax,PW);
ENDDEF

DEFINE t2_Y_PI_DEC

    if(t2by2y<t2by2){
    out     time(t2by2y);
    out     pw90Y		ch3:TG|AP(aY90, 0.0);
    out     time(PW180Y)	ch3:TG|AP(aY90,90.0);
    out     pw90Y   		ch3:TG|AP(aY90, 0.0);
    out     time(t2by2y);
    }

    else{
     out time(t2by2);
     out time(t2by2);
    }
    
ENDDEF

DEFINE t2_Y_PI_HOMODEC(PW_homodec)

    if(t2by2y<t2by2){
    out     time(t2by2y);
    if(PW_homodec>1.0u)
      {out     time(PW_homodec);}
    out     pw90Y		ch3:TG|AP(aY90, 0.0);
    out     time(PW180Y)	ch3:TG|AP(aY90,90.0);
    out     pw90Y   		ch3:TG|AP(aY90, 0.0);
    out     time(t2by2y);
    }

    else{
     out time(t2by2);
     out time(t2by2);
    }
    
ENDDEF

DEFINE t2_N_PI_DEC(pw90Y)

    if(t2by2y<t2by2){
    out     time(t2by2y);
    out     pw90Y		ch3:TG|AP(aY90, 0.0);
    out     time(PW180Y)	ch3:TG|AP(aY90,90.0);
    out     pw90Y   		ch3:TG|AP(aY90, 0.0);
    out     time(t2by2y);
    }

    else{
     out time(t2by2);
     out time(t2by2);
    }
    
ENDDEF

DEFINE t2_C_PI_DEC
t2_X_PI_DEC;
ENDDEF

DEFINE t2_X_PI_DEC

    if(t2by2x<t2by2){
    out     time(t2by2x);
    out     pw90X		ch1:TG|AP(aX90, 0.0);
    out     time(PW180X)	ch1:TG|AP(aX90,90.0);
    out     pw90X   		ch1:TG|AP(aX90, 0.0);
    out     time(t2by2x);
    }

    else{
     out time(t2by2);
     out time(t2by2);
    }
ENDDEF

DEFINE t2_ct_prog(dwell,AL,ALct,pw180CT,PW)

        if(AL< ALct){t2_ct = pw180CT + (dwell*ALct);}
       	if(AL>=ALct){t2_ct = pw180CT + (dwell*AL);  }
       	
       	t2_evol_out = t2_ct;
       	
     	cyclest2 = floor((t2_ct)/(PW));
        cyclest2 = cyclest2 - 1;
        cyclerm2 = t2_ct - cyclest2*PW2;

        if(cyclerm2<1.0u){
        cyclest2 = cyclest2 - 1;
        cyclerm2 = t2_evol - cyclest2*PW;
        }
!!!
	t2_evol_max = t2_evol + dwell*(AL-1);
     	if(AL<ALct){t2_evol_max2 = dwell*(ALct);}
	else       {t2_evol_max2 = dwell*(AL);}

!!!
        t2_evol_x = t2_evol_max2 + 10u - t2_evol_out;
	t2by2 = t2_evol_out / 2.0;
	if(t2by2>(2.0*PW180X+1u)){ 
	    t2by2x = t2by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2x=t2by2;}
	if(t2by2>(2.0*PW180Y+1u)){ 
	    t2by2y = t2by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2y=t2by2;}

        t2_evol_x = t2_evol_max2 + 10u - t2_evol_out;


	if(t2_evol_x<3.0u){t2_evol_x=3u;}
ENDDEF

DEFINE t2_ct_tau_prog(dwell,AL,ALct,pw180CS,PW,tau)

        if(AL< ALct){t2_evol_x = pw180CS + (dwell*ALct);}
        if(AL>=ALct){t2_evol_x = pw180CS + (dwell*AL);  }
        
        t2_ct = 2*tau;
        t2_evol_out = 2*tau;

	!if(t2_evol_x>t2_evol_out){error(t2outofbounds);
       	
     	cyclest2 = floor((t2_ct)/(PW));
        cyclest2 = cyclest2 - 1;
        cyclerm2 = t2_ct - cyclest2*PW;

        if(cyclerm2<1.0u){
        cyclest2 = cyclest2 - 1;
        cyclerm2 = t2_evol - cyclest2*PW;
        }
        
        
        cyclest2a = floor((tau)/(PW));
        cyclest2a = cyclest2a - 1;
        cyclerm2a = tau - cyclest2a*PW;

        if(cyclerm2a<1.0u){
        cyclest2a = cyclest2a - 1;
        cyclerm2a = tau - cyclest2a*PW;
        }
!!!
	t2_evol_max = t2_evol + dwell*(AL-1);
     	if(AL<ALct){t2_evol_max2 = dwell*(ALct);}
	else       {t2_evol_max2 = dwell*(AL);}

!!!
        t2_evol_x = t2_evol_max2 + 10u - t2_evol_out;
	t2by2 = t2_evol_out / 2.0;
	if(t2by2>(2.0*PW180X+1u)){ 
	    t2by2x = t2by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2x=t2by2;}
	if(t2by2>(2.0*PW180Y+1u)){ 
	    t2by2y = t2by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t2by2y=t2by2;}

        t2_evol_x = t2_evol_max2 + 10u - t2_evol_out;


	if(t2_evol_x<3.0u){t2_evol_x=3u;}
ENDDEF


DEFINE t2_ct_EVOLVE(dwell,ALstart,INDEX,PW180_CS,PW180_J,PW_Evol)

        t2_a= (t2_ct-PW180_CS) / 2.0 - ((INDEX+ALstart-1) * (dwell / 2.0));
        t2_b= (t2_ct-PW180_CS) / 2.0 + ((INDEX+ALstart-1) * (dwell / 2.0));
        
        if (t2_a<1.0u){t2_a=1u;}
        if (t2_b<1.0u){t2_b=1u;}

! Proton decoupling in each period
	cyclest2a = ceil(t2_a/PW_Evol)-1;
	cyclerm2a = t2_a - cyclest2a*PW_Evol;
 	cyclest2b = ceil(t2_b/PW_Evol)-1;
	cyclerm2b = t2_b - cyclest2b*PW_Evol;
	
        if(cyclerm2a<1.0u){
	cyclest2a = cyclest2a - 1;
	cyclerm2a = t2_a - cyclest2a*PW_Evol;
        }         
        if(cyclerm2b<1.0u){
	cyclest2b = cyclest2b - 1;
	cyclerm2b = t2_b - cyclest2b*PW_Evol;
        }         !making sure cyclerm2(a,b) isn't too short

        !for heteronuclear J-decoupling
!!!
	t2aby2 = t2_a / 2.0;
	if(t2aby2>(2.0*PW180X+1u)){ 
	    t2aby2x = t2aby2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t2aby2x=t2aby2;}
	if(t2aby2>(2.0*PW180Y+1u)){ 
	    t2aby2y = t2aby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t2aby2y=t2aby2;}

!!!
        t2bby2y = t2bby2x;
        
	t2bby2 = t2_b / 2.0;
	if(t2bby2>(2.0*PW180X+1u)){ 
	    t2bby2x = t2bby2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t2bby2x=t2bby2;}
        
	if(t2bby2>(2.0*PW180Y+1u)){ 
	    t2bby2y = t2bby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t2bby2y=t2bby2;}

ENDDEF

DEFINE t2_ct_X(pCS)

        if(t2_a > 1.0u){out time(t2_a);}
        out pw180X ch1:TG|AP(aX90,pCS);
        if(t2_b > 1.0u){out time(t2_b);}
ENDDEF

DEFINE t2_ct_Y(pCS)
        if(t2_a > 1.0u){out time(t2_a);}
        out pw180Y ch3:TG|AP(aY90,pCS);
        if(t2_b > 1.0u){out time(t2_b);}
ENDDEF

DEFINE t2_ct_X_J_DEC(PW_CS)

     if( t2_a > 2*PW180X){
        if(t2aby2y>2.0*PW180X){
        out        time(t2aby2x);
        out        pw90X           ch1:TG|AP(aX90,0.0);
        out        pw180X 	   ch1:TG|P(90.0);
        out        pw90X           ch1:TG|P(0.0);
        out        time(t2aby2x);
        }
        else {
        if (t2_a > 3u){out time(t2_a);} 
        }
        
        out time(PW_CS);

        if(t2bby2y>2.0*PW180X){
        out        time(t2bby2x)  ;
        out        pw90X           ch1:TG|AP(aX90,0.0);
        out        pw180X 	   ch1:TG|P(90.0);
        out        pw90X           ch1:TG|P(0.0);
        out        time(t2bby2x)  ;
        }
        else {out time(t2_b);}
      }
      else{
        out time(t2_a);
        out pw180Y;
        out time(t2_b);
      }
ENDDEF

DEFINE t2_ct_Y_J_DEC(PW_CS)

        if(t2aby2y>2.0*PW180Y){
        out        time(t2aby2y);
        out        pw90Y           ch3:TG|AP(aX90,0.0);
        out        pw180Y 	    ch3:TG|P(90.0);
        out        pw90Y           ch3:TG|P(0.0);
        out        time(t2aby2y);
        }
        if(t2aby2y<=2.0*PW180Y){
        if (t2_a > 3u){out time(t2_a);} 
        }
        
        out time(PW_CS);

        if(t2bby2y>2.0*PW180Y){
        out        time(t2bby2y)  ;
        out        pw90Y           ch3:TG|AP(aX90,0.0);
        out        pw180Y 	    ch3:TG|P(90.0);
        out        pw90Y           ch3:TG|P(0.0);
        out        time(t2bby2y)  ;
        }
        if(t2bby2y<=2.0*PW180Y){out time(t2_b);}
ENDDEF

! ERROR CODES FOR t2.inc
define t2outofbounds 0x1454
define OOBt2_ERROR_CODE USER_ERROR_BASE + OOBt2_ERROR
comment "ERROR "OOBt2_ERROR_CODE "Maximum t2 greater than constant time period requested!";