! ------------------------------------------
! 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/acq.inc";

!---- t1
    .time   extern   t1_evol     = 100u;
    .time            t1_evol_out = 100u;
    .time   extern   t1_evol_max = 100u;
    .time   extern   t1_evol_max2= 100u;
    .time            t1_evol_x   = 100u;
    .time            cyclerm1 = 1.0u;
    .time            t1by2  = 1.0u;
    .time            t1by2x = 1.0u; 
    .time            t1by2y = 1.0u; 
    .time   extern   autofix pwXdump= 10u;
    .time   extern   autofix pwYdump= 10u;

!---- t1 TPPM Parameters
    .float t1cycles "no. of TPPM cycles during t1" = 8; 
    .long  cyclest1 ;

!---- t1 constant time (echo recoupling)
    .time   extern   t1_ct       = 100u;
    .time            t1_a        = 100u;
    .time            t1aby2  = 1.0u;
    .time            t1aby2x = 1.0u; 
    .time            t1aby2y = 1.0u; 

    .time            t1_b        = 100u;
    .time            t1bby2  = 1.0u;
    .time            t1bby2x = 1.0u; 
    .time            t1bby2y = 1.0u; 

    .time            t1_c        = 100u;
    .time            t1cby2  = 1.0u;
    .time            t1cby2x = 1.0u; 
    .time            t1cby2y = 1.0u; 

    .long            cyclest1a ;
    .long            cyclest1b ;
    .long            cyclest1c ;
    .time            cyclerm1a = 1.0u;
    .time            cyclerm1b = 1.0u;
    .time            cyclerm1c = 1.0u;
    
    .time            t1LO        = 100u;
    .time            t1LOby2  = 1.0u;
    .time            t1LOby2x = 1.0u; 
    .time            t1LOby2y = 1.0u; 

    .long            cyclest1LO ;
    .time            cyclerm1LO = 1.0u;

    .long            cycles1tau ;
    .time            cyclerm1tau = 1.0u;

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

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

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

	t1_evol_max = t1_evol + dwell*(AL-1);
     	if(AL<ALmax){t1_evol_max2 = dwell*(ALmax);}
	else{t1_evol_max2 = dwell*(AL);}

	t1cycles = (t1_evol)/(PW);
	cyclest1 = floor(t1cycles);

	cyclerm1 = t1_evol - cyclest1*PW;
	if(cyclerm1<1.0u){
	cyclest1 = cyclest1 - 1;
	cyclerm1 = t1_evol - cyclest1*PW;
        }
ENDDEF

DEFINE t1_EVOLVE(dwell,ALstart,INDEX,PW_Evol)       

    	t1_evol_out = t1_evol + (dwell*(INDEX+ALstart-1));
	cyclest1 = floor((t1_evol_out)/(PW_Evol));
	if (cyclest1 < 0){cyclest1 = 0;}
	cyclerm1 = t1_evol_out - cyclest1*PW_Evol;

        if(cyclerm1<1.0u){
	cyclest1 = cyclest1 - 1;
	if (cyclest1 < 0){cyclest1 = 0;}
	cyclerm1 = t1_evol_out - cyclest1*PW_Evol;
        }         !making sure cyclerm1 isn't too short

	t1by2 = t1_evol_out / 2.0;
	if(t1by2>(2.0*PW180X+1u)){ 
	    t1by2x = t1by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2x=t1by2;}
	if(t1by2>(2.0*PW180Y+1u)){ 
	    t1by2y = t1by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2y=t1by2;}

        t1_evol_x = t1_evol_max2 + 10u - t1_evol_out;

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

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

    	t1_evol_out = (t1_evol + (dwell*(INDEX+ALstart-1)))/2.0;
	cyclest1 = ceil(t1_evol_out/PW_Evol)-1;
	if (cyclest1 < 0){cyclest1 = 0;}
	cyclerm1 = t1_evol_out - cyclest1*PW_Evol;
	
        if(cyclerm1<1.0u){
	cyclest1 = cyclest1 - 1;
	if (cyclest1 < 0){cyclest1 = 0;}
	cyclerm1 = t1_evol_out - cyclest1*PW_Evol;
        }         !making sure cyclerm2 isn't too short
    
	t1by2 = t1_evol_out;

	if((t1by2+(PW_homodec/2.0))>(2.0*PW180X+1u )){ 
	    t1by2x = t1by2 - PW90X; ! 4 pi/2's for composite,1 per half
        } 
        else{t1by2x=t1by2;}
	
	
	if((t1by2+(PW_homodec/2.0))>(2.0*PW180Y+1u )){ 
	t1by2y = t1by2 - 2*PW90Y; ! 4 pi/2's for composite,2 per half
        } 
        else{t1by2y=t1by2;}

        t1_evol_x = t1_evol_max2 + 10u - 2.0*t1_evol_out;

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

DEFINE t1_PROG(dwell,AL,ALmax,PW)
        t1_prog(dwell,AL,ALmax,PW);
ENDDEF

DEFINE t1_Y_PI_DEC

    if(t1by2y<t1by2){
    out     time(t1by2y);
    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(t1by2y);
    }

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

DEFINE t1_Y_PI_HOMODEC(PW_homodec)

    if(t1by2y<t1by2){
    out     time(t1by2y);
    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(t1by2y);
    }

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

DEFINE t1_N_PI_DEC(pw90Y)

    if(t1by2y<t1by2){
    out     time(t1by2y);
    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(t1by2y);
    }

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

DEFINE t1_C_PI_DEC
t1_X_PI_DEC;
ENDDEF

DEFINE t1_X_PI_DEC

    if(t1by2x<t1by2){
    out     time(t1by2x);
    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(t1by2x);
    }

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

DEFINE t1_ct_prog(dwell,AL,ALct,pw180CS,PW)

        if(AL< ALct){t1_ct = pw180CS + (dwell*ALct);}
       	if(AL>=ALct){t1_ct = pw180CS + (dwell*AL);  }
        
        t1_evol_out = t1_ct;

        cyclest1 = floor((t1_ct)/(PW));
        cyclest1 = cyclest1 - 1;
        cyclerm1 = t1_ct - cyclest1*PW;

        if(cyclerm1<1.0u){
        cyclest1 = cyclest1 - 1;
        cyclerm1 = t1_evol - cyclest1*PW;
        }
!!!
	t1_evol_max = t1_evol + dwell*(AL-1);
     	if(AL<ALct){t1_evol_max2 = dwell*(ALct);}
	else       {t1_evol_max2 = dwell*(AL);}

!!!
        t1_evol_x = t1_evol_max2 + 10u - t1_evol_out;
	t1by2 = t1_evol_out / 2.0;
	if(t1by2>(2.0*PW180X+1u)){ 
	    t1by2x = t1by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2x=t1by2;}
	if(t1by2>(2.0*PW180Y+1u)){ 
	    t1by2y = t1by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2y=t1by2;}

        t1_evol_x = t1_evol_max2 + 10u - t1_evol_out;


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

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

        if(AL< ALct){t1_evol = pw180CS + (dwell*ALct);}
        if(AL>=ALct){t1_evol = pw180CS + (dwell*AL);  }
        
        t1_ct = 2*tau + pw180CS;
        t1_evol_out = 2*tau;

	!if(t1_evol>t1_ct){error(t1outofbounds); ! This ain't workin
	
        cyclest1 = floor((t1_ct)/(PW));
        cyclest1 = cyclest1 - 1;
        cyclerm1 = t1_ct - cyclest1*PW;

        if(cyclerm1<1.0u){
        cyclest1 = cyclest1 - 1;
        cyclerm1 = t1_evol - cyclest1*PW;
        }
!!!
	t1_evol_max = t1_evol + dwell*(AL-1);
     	if(AL<ALct){t1_evol_max2 = dwell*(ALct);}
	else       {t1_evol_max2 = dwell*(AL);}

!!!
        t1_evol_x = t1_evol_max2 + 10u - t1_evol_out;
	t1by2 = t1_evol_out / 2.0;
	if(t1by2>(2.0*PW180X+1u)){ 
	    t1by2x = t1by2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2x=t1by2;}
	if(t1by2>(2.0*PW180Y+1u)){ 
	    t1by2y = t1by2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1by2y=t1by2;}

        t1_evol_x = t1_evol_max2 + 10u - t1_evol_out;


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

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

        t1_a= (t1_ct-PW180_CS) / 2.0 - ((INDEX+ALstart-1) * (dwell / 2.0));
        t1_b= (t1_ct-PW180_CS) / 2.0 + ((INDEX+ALstart-1) * (dwell / 2.0));
        
        if (t1_a<1.0u){t1_a=1u;}
        if (t1_b<1.0u){t1_b=1u;}

! Proton decoupling in each period
	cyclest1a = ceil(t1_a/PW_Evol)-1;
	cyclerm1a = t1_a - cyclest1a*PW_Evol;
 	cyclest1b = ceil(t1_b/PW_Evol)-1;
	cyclerm1b = t1_b - cyclest1b*PW_Evol;
	
        if(cyclerm1a<1.0u){
	cyclest1a = cyclest1a - 1;
	cyclerm1a = t1_a - cyclest1a*PW_Evol;
        }         
        if(cyclerm1b<1.0u){
	cyclest1b = cyclest1b - 1;
	cyclerm1b = t1_b - cyclest1b*PW_Evol;
        }         !making sure cyclerm1(a,b) isn't too short

        !for heteronuclear J-decoupling
!!!
	t1aby2 = t1_a / 2.0;
	if(t1aby2>(2.0*PW180X+1u)){ 
	    t1aby2x = t1aby2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1aby2x=t1aby2;}
	if(t1aby2>(2.0*PW180Y+1u)){ 
	    t1aby2y = t1aby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1aby2y=t1aby2;}

!!!
        t1bby2y = t1bby2x;
        
	t1bby2 = t1_b / 2.0;
	if(t1bby2>(2.0*PW180X+1u)){ 
	    t1bby2x = t1bby2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1bby2x=t1bby2;}
        
	if(t1bby2>(2.0*PW180Y+1u)){ 
	    t1bby2y = t1bby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1bby2y=t1bby2;}

!!!  Leftover for scalar experiments

	t1LO = t1_ct - 2*t1_a + PW180_CS;
	
	t1LOby2 = t1LO / 2.0;
	if(t1LOby2>(2.0*PW180X+1u)){ 
	    t1LOby2x = t1LOby2 - PW180X ; ! 2 pi's for composite,1 per half
        } 
        else{t1LOby2x=t1LOby2;}
	if(t1LOby2>(2.0*PW180Y+1u)){ 
	    t1LOby2y = t1LOby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1LOby2y=t1LOby2;}
        
        cyclest1LO = ceil(t1LO/PW_Evol)-1;
	cyclerm1LO = t1LO - cyclest1LO*PW_Evol;
	
        if(cyclerm1LO<1.0u){
	cyclest1LO = cyclest1LO - 1;
	cyclerm1LO = t1LO - cyclest1LO*PW_Evol;
        }         

        cycles1tau = ceil(tau/PW_Evol)-1;
	cyclerm1tau = tau - cycles1tau*PW_Evol;
	
        if(cyclerm1tau<1.0u){
	cycles1tau = cycles1tau - 1;
	cyclerm1tau = tau - cycles1tau*PW_Evol;
        }         


!!!
ENDDEF

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

        t1_a= tau1 - ((INDEX + ALstart - 1) * (dwell / 2.0));
        t1_b= tau1 - PW180_CS; ! 
        t1_c= ((INDEX + ALstart - 1) * (dwell / 2.0));
                
        if (t1_a<1.0u){t1_a=1u;}
        if (t1_b<1.0u){t1_b=1u;}
        if (t1_c<1.0u){t1_c=1u;}

! Proton decoupling in each period
	cyclest1a = ceil(t1_a/PW_Evol)-1;
	cyclerm1a = t1_a - cyclest1a*PW_Evol;
 	cyclest1b = ceil(t1_b/PW_Evol)-1;
	cyclerm1b = t1_b - cyclest1b*PW_Evol;
 	cyclest1c = ceil(t1_c/PW_Evol)-1;
	cyclerm1c = t1_c - cyclest1c*PW_Evol;
		
        if(cyclerm1a<1.0u){
	cyclest1a = cyclest1a - 1;
	cyclerm1a = t1_a - cyclest1a*PW_Evol;
        }         
        if(cyclerm1b<1.0u){
	cyclest1b = cyclest1b - 1;
	cyclerm1b = t1_b - cyclest1b*PW_Evol;
        } 
        if(cyclerm1c<1.0u){
	cyclest1c = cyclest1c - 1;
	cyclerm1c = t1_c - cyclest1c*PW_Evol;
        }         !making sure cyclerm1(a,b) isn't too short

        !for heteronuclear J-decoupling
!!!
	t1cby2 = t1_c / 2.0;
	if(t1cby2>(2.0*PW180Y+1u)){ 
	    t1cby2y = t1cby2 - PW180Y ; ! 2 pi's for composite,1 per half
        } 
        else{t1cby2y=t1cby2;}
        
!!!
ENDDEF

DEFINE t1_ct_X(pCS)

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

DEFINE t1_ct_Y(pCS)
        if(t1_a > 1.0u){out time(t1_a);}
        out pw180Y ch3:TG|AP(aY90,pCS);
        if(t1_b > 1.0u){out time(t1_b);}
ENDDEF

DEFINE t1_ct_X_J_DEC(PW_CS)

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

        out time(PW_CS);

        if(t1bby2y>2.0*PW180X){
        out        time(t1bby2x)  ;
        out        pw90X           ch1:TG|AP(aX90,0.0);
        out        pw180X 	   ch1:TG|P(90.0);
        out        pw90X           ch1:TG|P(0.0);
        out        time(t1bby2x)  ;
        }
        else {out time(t1_b);}
      }

        else{
        out time(t1_a);
        out pw180Y;
        out time(t1_b);
        }
ENDDEF

DEFINE t1_ct_Y_J_DEC(PW_CS)

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

         out time(PW_CS);

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


DEFINE t1_tau_Y_J_DEC(PW_CS)

         out        time(t1_a);
         out time(PW_CS);
         out        time(t1_a);


        if(t1LOby2y>2.0*PW180Y){
        out        time(t1LOby2y)  ;
        out        pw90Y           ch3:TG|AP(aY90,0.0);
        out        pw180Y 	    ch3:TG|P(90.0);
        out        pw90Y           ch3:TG|P(0.0);
        out        time(t1LOby2y)  ;
        }
        if(t1LOby2y<=2.0*PW180Y){
          if (t1LOby2y > 3u){out time(2*t1LOby2y);} 
        }
ENDDEF

DEFINE t1_CTUC_Y_J_DEC(PW_CS)

         out        time(t1_a);
         out time(PW_CS);
         out        time(tau1);

        if(t1cby2y>PW180Y){
        out        pw90Y           ch3:TG|AP(aY90,0.0);
        out        pw180Y 	    ch3:TG|P(90.0);
        out        pw90Y           ch3:TG|P(0.0);
        out        time(t1cby2y)  ;
        out        time(t1cby2y)  ;
        }
        if(t1cby2y<=PW180Y){
          if (t1cby2y > 3u){out time(2*t1cby2y);} 
        }
ENDDEF

DEFINE C_XMIX(PHASE)
  if(Xmix1flag == 0){
     out pwXdump ch1:TG|AP(aX90,PHASE);
    }
ENDDEF

DEFINE t1_C_XMIX(PHASE)

    if(Xmix1flag == 1){
    out  pw90X  	ch1:TG|AP(aX90,p270+PHASE);
    out	 time(ct1)	ch1:TG|AP(aXcp1,PHASE);
    }

    if(t1by2x<t1by2){

	 if(Xmix1flag == 1){
	 out	pw90X	 ch1:TG|AP(aX90,p270+PHASE);
     } 
	 else{out	pw90X	 ;}

     if(Xmix1flag == 1){
	    out	time(t1by2x)  ;
    	out 	pw90X	        ch1:TG|AP(aX90,PHASE);
    	out 	time(PW180X)	ch1:TG|P(p90+PHASE);
    	out 	pw90X	        ch1:TG|P(PHASE);
    	out	time(t1by2x)  ;
       }
      else{
	    out	time(t1by2x)  ;
    	out 	pw90X	        ch1:TG|AP(aX90,p00);
    	out 	time(PW180X)	ch1:TG|P(p90);
    	out 	pw90X	        ch1:TG|P(p00);
    	out	time(t1by2x)  ;
      }

	 if(Xmix1flag == 1){
	 out	pw90X	ch1:TG|AP(aX90,p270+PHASE); ! polarization along -z, flip up to Xmix
     }
	 else{out	pw90X	 ;}
     }
      
    else{
     out time(t1by2);
     out time(t1by2);
    }
ENDDEF

DEFINE H_XMIX(Xmix1flag)
   if(Xmix1flag == 0){
     out pwXdump ;
    }
ENDDEF

DEFINE H_XMIX2(Xmix1flag)
    if(Xmix1flag == 1){
    out	 pw90X	    ch2:TG ;
    out	 time(ct1)	ch2:TG ;
    }
ENDDEF

DEFINE N_XMIX(Xmix1flag)
    if(Xmix1flag == 0){
     out pwXdump ;
    }
ENDDEF

DEFINE t1_N_XMIX(Xmix1flag)
    if(Xmix1flag == 1){
    out	 pw90Y	    ch3:TG ;
    out	 time(ct1)	ch3:TG ;
    }
ENDDEF


! ERROR CODES FOR t1.inc
define t1outofbounds 0x1453
define OOBT1_ERROR_CODE USER_ERROR_BASE + OOBT1_ERROR
comment "ERROR "OOBT1_ERROR_CODE "Max t1 > CT period!";