bloch.m 1.69 KB
Newer Older
1
2
3
4
5
6
7
8
% Plot the pulse sequence
% create figure
figure(); 

% get pulse sequence B-fields from 0-5 ms
pulseq = zeros(1000,3);
for i = 1:1000
    pulseq(i,:) = B_eff(i*0.005); 
9
10
end

11
12
13
14
15
16
17
18
19
20
% plot RF
subplot(2,1,1);
plot(pulseq(:,1));
ylabel('B1');

% plot gradient
subplot(2,1,2);
plot(pulseq(:,3));
ylabel('Gradient');

21
%% Integrate ODE
22
23
24
25
26
27
28
29
T1 = 1500;
T2 = 50;
t0 = 0;
t1 = 5;
dt = 0.005;
M0 = [0; 0; 1];
[t, M] = ode45(@(t,M)bloch_ode(t, M, T1, T2), linspace(t0, t1, (t1-t0)/dt), M0);

30
%% Plot Results
31
32
33
34
% create figure
figure();hold on;

% plot x, y and z components of Magnetisation
35
36
37
plot(t, M(:,1), 'linewidth', 2);
plot(t, M(:,2), 'linewidth', 2);
plot(t, M(:,3), 'linewidth', 2);
38
39
40
41
42

% add legend and grid
legend({'Mx','My','Mz'});
grid on;

43
%% define the bloch equation
44
45
46
47
48
49
50
51
52
53
function dM = bloch_ode(t, M, T1, T2)
    % get effective B-field at time t
    B   =   B_eff(t);                               

    % cross product of M and B, add T1 and T2 relaxation terms
    dM  =  [M(2)*B(3) - M(3)*B(2) - M(1)/T2;        
            M(3)*B(1) - M(1)*B(3) - M(2)/T2;        
            M(1)*B(2) - M(2)*B(1) - (M(3)-1)/T1];   
end

54
%% define effective B-field
55
function b = B_eff(t)
56
57
    % Do nothing for 0.25 ms
    if t < 0.25             
58
        b = [0, 0, 0];
59
60
61
62
63
    % Sinc RF along x-axis and slice-select gradient on for 1.00 ms
    elseif t < 1.25             
        b = [pi*sinc((t-0.75)*4), 0, pi];
    % Do nothing for 0.25 ms
    elseif t < 1.50             
64
        b = [0, 0, 0];
65
66
67
68
69
70
    % Slice refocusing gradient on for 1.50 ms
    % Half the area of the slice-select gradient lobe
    elseif t < 3.00             
        b = [0, 0, -(1/3)*pi];
    % pulse sequence finished
    else                       
71
72
        b = [0, 0, 0];
    end
73
end