% This is the solution to M219 2016's HW-1 #5. Previously #4 
%
% 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');
%