%% 1A Design of RF pulse for 1H (3.0T) with 5mm slice thickness at isocenter

for tbw = [4,16]

samples = 512;

rf = (pi/2)*wsinc(tbw,samples);

dz  = 0.3;               % slice thickness [cm]
gam = 4.257;             % gyromagnetic ratio for 1H [kHz/G]

Gz = 8;                  % slice select gradient amplitude [G/cm]

df = Gz*gam*dz;          % BW in kHz

pulseduration = tbw/df;  % RF pulse duration [ms]

rfs = rfscaleg(rf, pulseduration); % Scaled to Gauss

dt = pulseduration/samples;
t = (1:length(rfs))*dt;  % in msec

b1    = [rfs zeros(1,samples/2)];                   % in Gauss
g     = Gz*[ones(1,samples) -ones(1,samples/2)];     % in G/cm
t_all = (1:length(g))*dt;  % in msec
 
figure; 
subplot(311); plot(t_all,b1); grid on;  title(['TBW: ' int2str(tbw)]);
xlabel('Time (msec)'); ylabel('B_1 (Gauss)');
subplot(312); plot(t_all,g); grid on; 
xlabel('Time (msec)'); ylabel('G_z (Gauss/cm)');

x = (-4:.01:4);          % in cm
f = 0;                  % center frequency in Hz
dt = pulseduration/samples/1e3;
t = (1:length(b1))*dt;  % in sec

% Bloch Simulation
[mx,my,mz] = bloch(b1(:),g(:),t(:),1,.2,f(:),x(:),0); 

% Transverse Magnetization
mxy = mx+1i*my;

subplot(313); plot(-x,abs(mxy)); 
grid on; xlabel('z position (cm)');

end

%% 1B Design of RF pulse for 1H (3.0T) with 3mm slice thickness at z=+10mm

tbw = 16;

samples = 512;

rf = (pi/2)*wsinc(tbw,samples);

dz  = 0.3;               % slice thickness [cm]
gam = 4.257;             % gyromagnetic ratio for 1H [kHz/G]

Gz = 8;                  % slice select gradient amplitude [G/cm]

df = Gz*gam*dz;          % BW in kHz

pulseduration = tbw/df;  % RF pulse duration [ms]

rfs = rfscaleg(rf, pulseduration); % Scaled to Gauss

dt = pulseduration/samples;
t = (1:length(rfs))*dt;  % in msec

b1    = [rfs zeros(1,samples/2)];                   % in Gauss
g     = Gz*[ones(1,samples) -ones(1,samples/2)];     % in G/cm
t_all = (1:length(g))*dt;  % in msec
 
figure; 
subplot(311); plot(t_all,b1); grid on;  title(['TBW: ' int2str(tbw)]);
xlabel('Time (msec)'); ylabel('B_1 (Gauss)');

subplot(312); plot(t_all,g); grid on; 
xlabel('Time (msec)'); ylabel('G_z (Gauss/cm)');

x = (-4:.1:4);          % in cm

% since we are off isocenter we are not as f=0 (in rotating frame)
z = 3;                % slice position [cm]
f = gam*Gz*z*1e3;           % in Hz
dt = pulseduration/samples/1e3;
t = (1:length(b1))*dt;  % in sec

% Bloch Simulation
[mx,my,mz] = bloch(b1(:),g(:),t(:),1,.2,f(:),x(:),0); 

% Transverse Magnetization
mxy = mx+1i*my;

subplot(313); plot(-x,abs(mxy)); 
grid on; xlabel('z position (cm)');