%% mds spin echo simulation for final project
% 5/13/2026

%% Parameters
Ns=100;       % Number of spins
B0=1;         % Static B-field in z-dir
g=50;         % Gyromagnetic ratio
Meq=1;        % Equilbrium magnetization
T1=inf;       % Relaxation time
B1=1/5;       % Excitation strength
M0=[0;0;Meq]; % Initial magnetization in equilibrium 
Te=1;         % Approx time between 90 and 180
noise=1/10;   % Strength of noise in B

%% Calculations
w1=g*B1;       % Excitation frequency
w0=g*B0;       % Larmor frequency 
T0=2*pi/g/B0;  % Larmor period 
dt=T0/100;     % 100 time steps per Larmor period

T90=pi/g/B1/2; % Time for a 90-pulse 

% Te must be a multiple of T0 for phase coherence
Te=fix((T90+Te)/T0)*T0-T90;


xa=[1 0];     % x-axis
ya=[0 1];     % y-axis 
M=0;  % Initialize M

%% Loop over all spins
for n=1:Ns
  disp(n);  % User feedback   
  
  Bv0=[0;0;B0];            % Average B vector
  Bs=Bv0+noise*randn(3,1); % Add noise to each spin
  
  % Set blochEq to solve during free evolution 
  Free=@(t,M) blochEq(M,Bs,T1,Meq,g);
  
  % Apply 90-deg rotation about the x-axis
  % this will rotate (0,0,Meq) into (Meq,0,0)
  % Note: use Bv0 to approximate fast rotation ie T90<<Te
  [t90,M90,ta]=aPulse(pi/2,M0,B1,w0,xa,Bv0,T1,Meq,g,dt);
  
  % Let evolve from last point of M90 for Te time
  ts=0:dt:Te;
  [t1,M1]=ode45(Free,ts,M90(end,:));
  
  
  % Apply 180-deg rotation about the y-axis
  % this will rotate flip spins causing them to reverse
  % Note: use Bv0 to approximate fast rotation ie T90<<Te
  % Start from last point of M1
  [t180,M180,tb]=aPulse(pi,M1(end,:),B1,w0,ya,Bv0,T1,Meq,g,dt);
  
  % Let evolve from last point of M180 for 2*Te time  
  ts=0:dt:2*Te;
  [t2,M2]=ode45(Free,ts,M180(end,:));
  
  % Concatenate all of the times and M's together
  t=[t90;ta+t1;ta+Te+t180;ta+Te+tb+t2];
  M=M+[M90;M1;M180;M2];
  
  tau=ta+Te;       % time from center of 180 to echo
  
  % Nice plot
  MR=rotFr(t,M,g*B0);
  h=plot(t,M/n,...
    t,.2*((t<ta)+(t>ta+Te & t<ta+Te+tb))-1,'k',...
    ta+tb/2+Te+tau*[1 1],[-1 1],'k--');
  axis([0 t(end) -Meq Meq]);
  drawnow;
  
end