function [t,M,ta]=aPulse(a,M,B1,w,ax,B0,T1,M0,g,dt)
% aPulse <alpha-pulse>
% Usage:: [t,M,ta]=aPulse(a,B1,B0,T1,M0,g)
%
% rotates M about the ax-axis by an angle a using a resonant
% field (B1*cos(w*t+p),B1*sin(w*t+p),B0);
% M is the current 3x1 magnetization vector
% B1,w rotating transverse magnetic field strength vector
%      and frequency
% ax Axis of rotation x-axis [1 0] y-axis [0 1] etc
% B0 is static magnetic field vector
% T1 is the relaxation time T_1
% M0 is the equilibrium magnetization M_0
% g is \gamma the gyromagnetic ratio
% dt time step for output

% revision history:
% 05/13/2026 Mark D. Shattuck <mds> aPulse.m

%% Main
p=atan2(ax(2),ax(1)); % phase angle to set axis of rotation
ta=a/g/B1;            % Time for a alpha-pulse 
ts=0:dt:ta;           % Time points for output 

opts=odeset('reltol',1e-4);  % Higher accuracy

% Magnetic field
B=@(t) [B1*cos(w*t+p);-B1*sin(w*t+p);0]+B0;   

% Solve blochEq
[t,M]=ode45(@(t,M) blochEq(M,B(t),T1,M0,g),ts,M,opts);

