% 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); end % plot RF subplot(2,1,1); plot(pulseq(:,1)); ylabel('B1'); % plot gradient subplot(2,1,2); plot(pulseq(:,3)); ylabel('Gradient'); %% Integrate ODE 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); %% Plot Results % create figure figure();hold on; % plot x, y and z components of Magnetisation plot(t, M(:,1), 'linewidth', 2); plot(t, M(:,2), 'linewidth', 2); plot(t, M(:,3), 'linewidth', 2); % add legend and grid legend({'Mx','My','Mz'}); grid on; %% define the bloch equation 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 %% define effective B-field function b = B_eff(t) % Do nothing for 0.25 ms if t < 0.25 b = [0, 0, 0]; % 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 b = [0, 0, 0]; % 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 b = [0, 0, 0]; end end