function [b] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%function [b] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%Function that computes the complex-valued SPSP RF pulse waveform,
%using parameters specified in structure rfp.
%Inputs:
%kp: k-space trajectory parameter structure, generated by kparameterSPSP.m
%rfp:RF waveform design parameter structure, generated by rfparameterSPSP.m
%gz: z gradient waveform vector
%kz,kf: SPSP k-space (kf-kz space) trajectory vectors
%Output:
%b:complex-value SPSP RF pulse waveform in units of Gauss (g).
%
%
%Chun-yu Yip, 4/1/2009

%Physical parameters
gam = 26751;                       %rad/s/g; gyromagnetic ratio
gambar = gam/2/pi;                 %Hz/g

%Create grid in SPSP space
dresz = rfp.dfovz/rfp.ddimz;       %cm
dresf = rfp.dfovf/rfp.ddimf;       %Hz
z = [-rfp.dfovz/2:dresz:rfp.dfovz/2-dresz];
f = [-rfp.dfovf/2:dresf:rfp.dfovf/2-dresf];
[Z,F] = ndgrid(z,f);

disp('Creating desired pattern for SPSP pulse...');

%phase pattern
TD = rfp.TE-kp.pointtime*kp.npnts/2; %pulse end to acquisition of DC sample
dphs = gam * rfp.alpha * TD * F .* (Z-rfp.zshift); %Eq. 5

%magnitude pattern
if streq(rfp.profileshape,'rect')

    dmag = zeros(size(Z));    
    dmag(find((Z<(rfp.slthickz/2)) & (Z>(-rfp.slthickz/2)))) = 1;

elseif streq(rfp.profileshape,'gaussian') 

    a = -4*log(2)/rfp.slthickz/rfp.slthickz;
    dmag = exp(a*(Z.^2));

else

    disp('Slice profile shape unknown.');

end
    
if rfp.smoothprofile

    kernel = gaussian(z,0,rfp.kernelstdv)';
    dmag=convn_fft(dmag,kernel);     %Smoothing implemented in Fourier domain

end;

dmag = dmag/max(max(dmag));
dmag = sin(rfp.flipangle/180 * pi) * dmag;
d = exp(1i*dphs).*dmag;

%Define matrices and vectors for conjugate gradient
b0 = zeros(length(kz),1);            %Initiation of CG
mask = ones(size(d));                %ROI in f-z space
W = spdiag(mask);                    %weighting based on ROI definition

%total energy penalty
beta = rfp.beta;
C = beta*speye(length(kz));

M = length(Z(:));
precon =speye(length(kz));           %CG Preconditioner
Niter = rfp.Niter;                   %Number of CG iterations

%Define matlab object "A" to represent the gigantic system matrix.
%"A" is of type "Fatrix" from Jeff Fessler's toolbox.
disp('Forming spatial-spectral fast object...');
A=spsp_Af(z,f,kz,kf,kp.pointtime,rfp.dfovf,rfp.dfovz);

%%%%%%%%%%%Conjugate Gradient%%%%%%%%%%%%%%%%%%%
disp('Running preconditioned conjugate gradient...');

%keyboard
[b_all, info] = qpwls_pcg(b0, A, W, d(:), zeros(M,1), sqrt(C), precon, Niter);
b = b_all(:,Niter);

if rfp.hamming
scale = 2;          %This parameter controls how wide the window is.
hamm=(0.54-0.46.*cos(pi*((kz/scale/max(kz))+1)));
b= b.*hamm;         %Intended to smooth the excited slice profile.

end;

m = A*b;            %Small-tip-angle approximation of excitation pattern
mm = reshape(m,size(d));


disp('Displaying SPSP pulse design results...');

t = kp.pointtime*[0:1:length(gz)-1].';

figure;
subplot(2,1,1);
plot(t*1000,gz);
grid;
x1=xlabel('Time (ms)');
y1=ylabel('g/cm');
t1=title('z gradient');
set(t1,'Fontsize',14);
set(x1,'Fontsize',12);
set(y1,'Fontsize',12);

subplot(2,1,2); 
hold on;
plot(t*1000,real(b),'b');
plot(t*1000,imag(b),'r');
grid;
x2=xlabel('Time (ms)');
y2=ylabel('g');
t2=title('SPSP pulse');
l2=legend('Real','Imaginary');
set(t2,'Fontsize',14);
set(x2,'Fontsize',12);
set(y2,'Fontsize',12);
set(l2,'Fontsize',12);

%Display desired pattern of the SPSP recovery pulse.
figure;
subplot(2,1,1);
imagesc(f,z,abs(d));
colormap default;colorbar;
x1 = xlabel('Frequency offset (Hz)');
y1 = ylabel('z (cm)');
t1 = title('Desired SPSP pattern (magnitude)');

set(x1,'Fontsize',12);
set(y1,'Fontsize',12);
set(t1,'Fontsize',12);


subplot(2,1,2);
imagesc(f,z,angle(d).*abs(d),[-pi,pi]);colorbar; 
x2 =xlabel('Frequency offset (Hz)');
y2=ylabel('z (cm)');
t2=title('Desired SPSP pattern (phase)');
set(x2,'Fontsize',12);
set(y2,'Fontsize',12);
set(t2,'Fontsize',12);


figure;
subplot(2,1,1);
imagesc(f,z,abs(mm));
xlabel('frequency (Hz)');
ylabel('z (cm)');
title('Resulting pattern (linear regime prediction)');

subplot(2,1,2);
imagesc(f,z,angle(mm));
xlabel('frequency (Hz)');
ylabel('z (cm)');
title('Resulting phase pattern (linear regime prediction)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

