Add Double Angle scripts and data

qMRLab · Oct 4, 2024 · 60c89e1 · 60c89e1
1 parent 9f94e16
commit 60c89e1
Show file tree

Hide file tree

Showing 8 changed files with 328 additions and 0 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/04-B1/01-DoubleAngle/binder/postBuild b/04-B1/01-DoubleAngle/binder/postBuild
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+git clone https://www.github.com/qMRLab/qMRLab.git ../qMRLab
+cd ../qMRLab
+git checkout f481d86da7272ee0f74abdbe9020cf93be06000f
+cd binder
+
diff --git a/04-B1/01-DoubleAngle/binder/requirements.txt b/04-B1/01-DoubleAngle/binder/requirements.txt
@@ -0,0 +1 @@
+repo2data
diff --git a/04-B1/01-DoubleAngle/content/da_fig1.m b/04-B1/01-DoubleAngle/content/da_fig1.m
@@ -0,0 +1,43 @@
+close all
+clear all
+clc
+
+nom_first_ang = 60;
+B1 = linspace(0.7, 1.3, 61);
+
+true_FA = B1*nom_first_ang;
+true_2FA = 2*B1*nom_first_ang;
+
+TR_range = [linspace(10,100,10), linspace(200,1000,9), linspace(2000,10000,9)];
+
+B1_varTR = zeros(length(TR_range), length(B1), 3);
+
+linfit = zeros(length(TR_range),2,3);
+
+for ii=1:length(TR_range)
+    params_FA.EXC_FA = true_FA;
+    params_FA.T1 = 900;
+    params_FA.TR = TR_range(ii);
+
+    params_2FA.EXC_FA = true_2FA;
+    params_2FA.T1 = 900;
+    params_2FA.TR = TR_range(ii);
+
+    B1_varTR(ii,:,1) = acosd(vfa_t1.analytical_solution(params_2FA)./(2.*vfa_t1.analytical_solution(params_FA)))./nom_first_ang;    
+    linfit(ii,:,1) = polyfit(B1,B1_varTR(ii,:,1),1);
+
+    params_FA.T1 = 1500;
+    params_2FA.T1 = 1500;
+
+    B1_varTR(ii,:,2) = acosd(vfa_t1.analytical_solution(params_2FA)./(2.*vfa_t1.analytical_solution(params_FA)))./nom_first_ang;    
+    linfit(ii,:,2) = polyfit(B1,B1_varTR(ii,:,2),1);
+
+    params_FA.T1 = 4000;
+    params_2FA.T1 = 4000;
+
+    B1_varTR(ii,:,3) = acosd(vfa_t1.analytical_solution(params_2FA)./(2.*vfa_t1.analytical_solution(params_FA)))./nom_first_ang;    
+    linfit(ii,:,3) = polyfit(B1,B1_varTR(ii,:,3),1);    
+
+end
+
+save("da_fig1.mat", "nom_first_ang", "B1", "TR_range", "B1_varTR", "linfit")
diff --git a/04-B1/01-DoubleAngle/content/da_fig2.m b/04-B1/01-DoubleAngle/content/da_fig2.m
@@ -0,0 +1,69 @@
+% clear all
+% close all
+% clc
+% 
+% alpha = deg2rad(60);
+% B1 = 1;
+% T1 = 1000;
+% T2 = 30;
+% TE = 20;
+% TR = 10000;
+% crushFlag = 1;
+% partialDephasingFlag = 1;
+% partialDephasing = 0.3;
+% df = 0;
+% inc = 0;
+% Nex = 100;
+% 
+% 
+% [Msig1,MLong1]=da_blochsim(alpha, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+% 
+% alpha2 = deg2rad(120);
+% 
+% [Msig2,MLong2]=da_blochsim(alpha2, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+
+%%
+
+close all
+clear all
+clc
+
+nom_first_ang = 60;
+B1 = linspace(0.7, 1.3, 61);
+
+TR_range = [linspace(10,100,10), linspace(200,1000,9), linspace(2000,10000,9)];
+
+B1_hard = zeros(length(TR_range), length(B1));
+B1_ideal = zeros(length(TR_range), length(B1));
+
+
+for ii=1:length(TR_range)
+    alpha = deg2rad(nom_first_ang);
+    T1 = 900;
+    T2 = 30;
+    TE = 20;
+    TR = TR_range(ii);
+    crushFlag = 0;
+    partialDephasingFlag = 1;
+    partialDephasing = 0.3;
+    df = 0;
+    inc = 0;
+    Nex = 100;
+
+
+    B1_shortTR = zeros(1,length(B1));
+    for jj=1:length(B1)
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard');
+        B1_hard(ii,jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;
+
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'ideal');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'ideal');
+        B1_ideal(ii,jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;        
+    end
+
+end
+
+save("da_fig2.mat", "nom_first_ang", "B1", "TR_range", "B1_hard", "B1_ideal")
diff --git a/04-B1/01-DoubleAngle/content/da_fig3.m b/04-B1/01-DoubleAngle/content/da_fig3.m
@@ -0,0 +1,74 @@
+% clear all
+% close all
+% clc
+% 
+% alpha = deg2rad(60);
+% B1 = 1;
+% T1 = 1000;
+% T2 = 30;
+% TE = 20;
+% TR = 10000;
+% crushFlag = 1;
+% partialDephasingFlag = 1;
+% partialDephasing = 0.3;
+% df = 0;
+% inc = 0;
+% Nex = 100;
+% 
+% 
+% [Msig1,MLong1]=da_blochsim(alpha, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+% 
+% alpha2 = deg2rad(120);
+% 
+% [Msig2,MLong2]=da_blochsim(alpha2, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+
+%%
+
+close all
+clear all
+clc
+
+nom_first_ang = 60;
+B1 = linspace(0.7, 1.3, 61);
+
+TR_range = [linspace(10,100,10), linspace(200,1000,9), linspace(2000,10000,9)];
+
+B1_hard = zeros(length(TR_range), length(B1));
+B1_ideal = zeros(length(TR_range), length(B1));
+B1_composite = zeros(length(TR_range), length(B1));
+
+
+for ii=1:length(TR_range)
+    alpha = deg2rad(nom_first_ang);
+    T1 = 900;
+    T2 = 30;
+    TE = 20;
+    TR = TR_range(ii);
+    crushFlag = 0;
+    partialDephasingFlag = 1;
+    partialDephasing = 0.3;
+    df = 0;
+    inc = 0;
+    Nex = 100;
+
+
+    B1_shortTR = zeros(1,length(B1));
+    for jj=1:length(B1)
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard');
+        B1_hard(ii,jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;
+
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'ideal');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'ideal');
+        B1_ideal(ii,jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;        
+
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'composite');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'composite');
+        B1_composite(ii,jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;  
+    end
+
+end
+
+save("da_fig3.mat", "nom_first_ang", "B1", "TR_range", "B1_hard", "B1_ideal", "B1_composite")
diff --git a/04-B1/01-DoubleAngle/content/da_fig4.m b/04-B1/01-DoubleAngle/content/da_fig4.m
@@ -0,0 +1,72 @@
+% clear all
+% close all
+% clc
+% 
+% alpha = deg2rad(60);
+% B1 = 1;
+% T1 = 1000;
+% T2 = 30;
+% TE = 20;
+% TR = 10000;
+% crushFlag = 1;
+% partialDephasingFlag = 1;
+% partialDephasing = 0.3;
+% df = 0;
+% inc = 0;
+% Nex = 100;
+% 
+% 
+% [Msig1,MLong1]=da_blochsim(alpha, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc)
+% 
+% 
+% alpha2 = deg2rad(120);
+% 
+% [Msig2,MLong2]=da_blochsim(alpha2, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc)
+% 
+
+%%
+
+close all
+clear all
+clc
+
+nom_first_ang = 60;
+B1 = linspace(0.7, 1.3, 61);
+
+TR_range = [linspace(10,100,10), linspace(200,1000,9), linspace(2000,10000,9)];
+
+for ii=1:length(TR_range)
+    alpha = deg2rad(nom_first_ang);
+    T1 = 1000;
+    T2 = 30;
+    TE = 20;
+    TR = TR_range(ii);
+    crushFlag = 0;
+    partialDephasingFlag = 1;
+    partialDephasing = 0.3;
+    df = 0;
+    inc = 0;
+    Nex = 100;
+
+
+    B1_shortTR = zeros(1,length(B1));
+    for jj=1:length(B1)
+        [sig1, ~] = da_blochsim_composite(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'composite');
+        [sig2, ~] = da_blochsim_composite(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc,'composite');
+        B1_shortTR(jj) = acosd(abs(sig2)./(2*abs(sig1)))./nom_first_ang;
+    end
+
+    plot(B1, B1_shortTR)
+    title('TR = ' + string(TR_range(ii)) + ' milliseconds', 'T1 = 1 second');
+
+    P = polyfit(B1,B1_shortTR,1);
+    yfit = polyval(P,B1);
+    hold on;
+    plot(B1,B1,'r-.');
+    eqn = string(" Linear: y = " + P(1)) + "x + " + string(P(2));
+    text(min(B1),max(B1_shortTR),eqn,"HorizontalAlignment","left","VerticalAlignment","top")
+    disp(B1_shortTR)
+    hold off
+    pause(0.5)
+end
+
diff --git a/04-B1/01-DoubleAngle/content/da_fig5.m b/04-B1/01-DoubleAngle/content/da_fig5.m
@@ -0,0 +1,62 @@
+% clear all
+% close all
+% clc
+% 
+% alpha = deg2rad(60);
+% B1 = 1;
+% T1 = 1000;
+% T2 = 30;
+% TE = 20;
+% TR = 10000;
+% crushFlag = 1;
+% partialDephasingFlag = 1;
+% partialDephasing = 0.3;
+% df = 0;
+% inc = 0;
+% Nex = 100;
+% 
+% 
+% [Msig1,MLong1]=da_blochsim(alpha, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+% 
+% alpha2 = deg2rad(120);
+% 
+% [Msig2,MLong2]=da_blochsim(alpha2, B1, T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'hard')
+% 
+
+%%
+
+close all
+clear all
+clc
+
+nom_first_ang = 60;
+B1 = 1;
+
+TR_range = 5000;
+
+B1_double = zeros(length(TR_range), length(B1));
+
+
+for ii=1:length(TR_range)
+    alpha = deg2rad(nom_first_ang);
+    T1 = 900;
+    T2 = 9000;
+    TE = 20;
+    TR = TR_range(ii);
+    crushFlag = 1;
+    partialDephasingFlag = 0;
+    partialDephasing = 0.3;
+    df = 0;
+    inc = 0;
+    Nex = 1;
+
+    for jj=1:length(B1)
+        [sig1, ~] = da_blochsim(alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'double');
+        [sig2, ~] = da_blochsim(2*alpha, B1(jj), T1, T2, TE, TR, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc, 'double');
+        B1_double(ii,jj) = (acos(abs(sig2)./(8*abs(sig1))))^(1/3)./deg2rad(nom_first_ang);
+    end
+
+end
+
+%save("da_fig5.mat", "nom_first_ang", "B1", "TR_range", "B1_hard", "B1_ideal")