%
% 2015.01.14 DBE@UCLA
% Updated 2016.01.08
%% HW01 Problem #5E - Rotating Frame
% Define some constants
clear
clc
close all
alpha=pi/2; % RF flip angle [rad]
theta=pi/4; % RF phase [rad]
B1_max=1e-6; % Maximum B1 amplitude [T]
gamma=42.57e6; % 1H gyromagnetic ratio in [MHz/T]
B0=5e-5; % Main magnetic field [T]
t_RF=alpha/(2*pi*gamma*B1_max); % Shortest RF pulse duration in [s]
N_steps=1e3; % Number of simulation time points [#]
dt=t_RF/N_steps; % Time step duration [s]
t=0:dt:(t_RF-dt); % Define a time vector [s]
% Initialize the bulk magnetization vectors
M0=[0 0 1]'; % Initial condition [A.U.]
M_RotFrm=zeros(4,N_steps); % Initialize the magnetization vectors
M_RotFrm(:,1)=[M0; 0]; % Define the initial condition
% Define the bulk magnetization operators
% These should be inside the for-loop if they are non-constant...
dB1=PAM_B1_op(gamma,B1_max,dt,theta); % Define the incremental RF excitation
% Simulate forced precession etc.
for n=2:N_steps
M_RotFrm(:,n)=dB1*M_RotFrm(:,n-1); % RF pulse acts on spin
end
% Plot the rotating frame results...
figure; hold on;
p(1)=plot(t,M_RotFrm(1,:)');
p(2)=plot(t,M_RotFrm(2,:)');
p(3)=plot(t,M_RotFrm(3,:)');
legend('M_x','M_y','M_z');
xlabel('Time [s]');
ylabel('Components of M [A.U.]');
title('Forced Precession in the Rotating Coordinate Frame');
set(p,'LineWidth',3);
set(gcf,'Color','w');
set(gca,'Color','w','XColor','k','YColor','k');
set(gca,'Color','w','XColor','k','YColor','k','LineWidth',1.25,'Fontsize',11,'Fontweight','bold');
set(get(gca,'Title'),'Color','k','FontSize',16);
set(get(gca,'Xlabel'),'FontSize',14,'fontweight','bold');
set(get(gca,'Ylabel'),'FontSize',14,'fontweight','bold');
grid('on')
print2desktop('/Users/pmagrath/desktop/','HW1_P5_E')
%% HW01 Problem #4E - Laboratory Frame
% The RF pulse is long compared to the period of precession so we only
% simulate a small percentage of the RF pulse with precession.
%t_RF=0.000025*t_RF; % Part of the shortest RF pulse duration in [s]
N_steps=1e3; % Number of simulation time points [#]
dt=t_RF/N_steps; % Time step duration [s]
t=0:dt:(t_RF-dt); % Define a time vector [s]
% By forcing the initial condition to be a later time point we can see the
% components better...
%M0=M_RotFrm(1:3,500); % Initial condition [A.U.]
M0 = [0 0 1]' %M0=M_RotFrm(1:3,);
M_RotFrm=zeros(4,N_steps); % Re-Initialize the magnetization vectors
M_RotFrm(:,1)=[M0; 0]; % Define the initial condition
M_LabFrm=zeros(4,N_steps); % Initialize the magnetization vectors
dB1=PAM_B1_op(gamma,B1_max,dt,theta); % Define the incremental RF excitation
dB0=PAM_B0_op(gamma,B0,dt) % Define the incremental precession
% Simulate forced precession - the dB1 operator is in the rotating frame
for n=2:N_steps
M_RotFrm(:,n)=dB1*M_RotFrm(:,n-1); % RF pulse acts on spin
end
% Simulate free precession - the dB0 operator transforms from the rotating
% the the laboratory frame
M_LabFrm(:,1)=[[1 0 0]'; 0]; % Define the initial condition
for n=2:N_steps
M_LabFrm(:,n)=dB0*M_LabFrm(:,n-1); % Add precession to the rotating frame solution
end
% Add forced precession data back in and scale based on (1-Mz)
M_LabFrm(3,:) = M_RotFrm(3,:);
for n=1:N_steps
M_LabFrm(1,n) = M_LabFrm(1,n).*(1-M_LabFrm(3,n));
M_LabFrm(2,n) = M_LabFrm(2,n).*(1-M_LabFrm(3,n));
end
% Plot the laboratory frame results...
figure; hold on;
q(1)=plot(1e9*t,M_LabFrm(1,:)','--');
q(2)=plot(1e9*t,M_LabFrm(2,:)');
q(3)=plot(1e9*t,M_LabFrm(3,:)');
legend('M_x','M_y','M_z');
xlabel('Time [ns]');
ylabel('Components of M [A.U.]');
title('Forced Precession in the Laboratory Coordinate Frame');
set(q,'LineWidth',3);
set(gcf,'Color','w');
set(gca,'Color','w','XColor','k','YColor','k');
set(gca,'Color','w','XColor','k','YColor','k','LineWidth',1.25,'Fontsize',11,'Fontweight','bold');
set(get(gca,'Title'),'Color','k','FontSize',16);
set(get(gca,'Xlabel'),'FontSize',14,'fontweight','bold');
set(get(gca,'Ylabel'),'FontSize',14,'fontweight','bold');
grid('on')
%% HW01 Problem #4F - Relaxation in the rotating frame...
% Define some constants
alpha=pi/2; % RF flip angle [rad]
theta=0; % RF phase [rad]
B1_max=1e-6; % Maximum B1 amplitude [T]
gamma=42.57e6; % 1H gyromagnetic ratio in [MHz/T]
B0=5e-6;%1.5; % Main magnetic field [T]
%% Simulate for the "long" relaxation times...
T1=1; % T1 time constant [s]
T2=0.1; % T2 time constant [s]
t_RF=alpha/(2*pi*gamma*B1_max); % Shortest RF pulse duration in [s]
N_steps=1e3; % Number of simulation time points [#]
dt=t_RF/N_steps; % Time step duration [s]
t=0:dt:(t_RF-dt); % Define a time vector [s]
% Initialize the bulk magnetization vectors
M0=[0 0 1]'; % Initial condition [A.U.]
M_RotFrm=zeros(4,N_steps); % Initialize the magnetization vectors
M_RotFrm(:,1)=[M0; 0]; % Define the initial condition
% Define the bulk magnetization operators
dB1=PAM_B1_op(gamma,B1_max,dt,theta); % Define the incremental RF excitation
dRlx=PAM_Relax_op(T1,T2,dt);
% Simulate forced precession etc.
for n=2:N_steps
M_RotFrm(:,n)=dRlx*dB1*M_RotFrm(:,n-1); % RF pulse acts on spin
end
%% Simulate for the "short" relaxation times...
T1=0.250; % T1 time constant [s]
T2=0.025; % T2 time constant [s]
% Initialize the bulk magnetization vectors
M0=[0 0 1]'; % Initial condition [A.U.]
M_RotFrm2=zeros(4,N_steps); % Initialize the magnetization vectors
M_RotFrm2(:,1)=[M0; 0]; % Define the initial condition
% Update the relaxation operators
dRlx=PAM_Relax_op(T1,T2,dt);
% Simulate forced precession etc.
for n=2:N_steps
M_RotFrm2(:,n)=dRlx*dB1*M_RotFrm2(:,n-1); % RF pulse acts on spin
end
%% Plot the rotating frame results...
figure; hold on;
p(1)=plot(t,M_RotFrm(1,:)');
p(2)=plot(t,M_RotFrm(2,:)');
p(3)=plot(t,M_RotFrm(3,:)');
p(4)=plot(t,M_RotFrm2(1,:)','--');
p(5)=plot(t,M_RotFrm2(2,:)','--');
p(6)=plot(t,M_RotFrm2(3,:)','--');
legend('Mx (T1 = 1000)','My (T1 = 1000)','Mz (T1 = 1000)','M_x (T1 = 250)','M_y (T1 = 250)','M_z (T1 = 250)');
xlabel('Time [s]');
ylabel('Components of M [A.U.]');
title('Forced Precession + Relaxation in the Rotating Frame');
set(p,'LineWidth',3);
set(gcf,'Color','w');
set(gca,'Color','w','XColor','k','YColor','k');
set(gca,'Color','w','XColor','k','YColor','k','LineWidth',1.25,'Fontsize',11,'Fontweight','bold');
set(get(gca,'Title'),'Color','k','FontSize',16);
set(get(gca,'Xlabel'),'FontSize',14,'fontweight','bold');
set(get(gca,'Ylabel'),'FontSize',14,'fontweight','bold');
grid('on')
print2desktop('/Users/pmagrath/desktop/','HW1_P5_F')
% %% HW01 Problem #4G EXTRA CREDIT
% clear all;
% N=10; % Number of T2 species...
% DUR=3e-3;
% FLIP=90;
% dt=1e-6;
%
% t=0:dt:DUR;
%
% B1_max=rf_hard_pulse('flip',FLIP,'dur',DUR);
%
% b1=B1_max*ones(numel(t),1);
% gr=zeros(size(b1));
% T2=logspace(-4,-1,N);
% % T1=1000e-3;
% T1=10*T2;
% df=0;
% dp=[0 0 0];
% mode=2;
% mx0=0;
% my0=0;
% mz0=1;
%
% for k=1:N
% [mx(:,k),my(:,k),mz(:,k)]=bloch(b1,gr,dt,T1,T2(k),df,dp,mode,mx0,my0,mz0);
%
% if k==1, figure; hold on; end
% C=hsv_colors(N);
% p(N)=plot((mx(:,k).^2+my(:,k).^2).^0.5,mz(:,k));
% set(p(N),'color',C(k,:));
%
% LEG{k}=[num2str(1000*T2(k),'%0.1g'),'ms'];
% end
%
% legend(LEG);
% title('Mz vs. Mxy for various T2/T1 Moieties');
% xlabel('M_x_y');
% ylabel('M_z');
%