% This is the solution to M219 2017's HW-1 #3.
%
% 2017.01.30 Patrick Magrath & DBE@UCLA

%% HW01 Problem #2A
% Define some constants
dt=1e-6;                 % Time steps in seconds
dt_2 = 1e-10;
gamma=42.57e6;             % Gyromagnetic ratio, [Hz/T]
B0=5.87e-6;                    % External magnetic field strength [T]
%PPM=3.0e-2;                % Fat-Water shift was stated as 2.4ppm, but it is really 3.4ppm

N_Steps=1000;                % Number of time steps
N_Steps_2 = 10000;           % Number of time steps for B0
M0=[1 0 0]';                 % Initial condition [A.U.] 
t_max=dt*N_Steps;            % Maximum time of simulation [s]
t_max_2 = dt_2*N_Steps_2;      % For the B0 simulation
t=linspace(0,t_max,N_Steps); % Vector of time points [s]
t_2 = linspace(0,t_max_2,N_Steps_2); %For the second simulation

% Since we're just looking at precession with respect to a given
% field, we treat the B1 field as if it was the B0 field over the duration
% of the RF pulse. 

M_H20=zeros(4,N_Steps);                 % Initialize the magnetization vectors
M_H20(:,1)=[M0; 0];                     % Define the initial condition
M_part2(:,1) = [M0;0];
dB0_H20=PAM_B0_op(gamma,B0,dt);         % Define the incremental precession
dB0_part2 = PAM_B0_op(gamma,3,dt_2);

% Simulate precession
for n=2:N_Steps
  M_H20(:,n) = dB0_H20*M_H20(:,n-1);
end

for n=2:N_Steps_2
    M_part2(:,n) = dB0_part2*M_part2(:,n-1);
end

M_part2(:,n) = dB0_part2*M_part2(:,n-1);

% Plot the results - Note F/W differences are hard to see on short time scales...
figure; hold;
  plot(t,M_H20(1,:),'linewidth',2);
  plot(t,M_H20(2,:),'linewidth',2);
  %plot(t,M_H20(3,:),'linewidth',2);
  
  xlabel('Time [s]');
  ylabel('Magnetization [a.u.]');
  legend('Mx B1 Field','My B1 Field');
  title('Precession in the B1 Field')
  
  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/','HW2_P2_B1')
  
  figure; hold;
  plot(t_2(1:1000),M_part2(1,1:1000),'linewidth',2);
  plot(t_2(1:1000),M_part2(2,1:1000),'linewidth',2);
  %plot(t,M_H20(3,:),'linewidth',2);
  
  xlabel('Time [s]');
  ylabel('Magnetization [a.u.]');
  legend('Mx B0 Field','My B0 Field');
  title('Precession in the B0 Field');
  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/','HW2_P2_B0')