 function [b, d, f, z, mm] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%function [b, d, f, z, mm] = 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);

printm('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

    printm('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
W = diag_sp(mask); % jf version

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

M = length(Z(:));
%precon =speye(length(kz));           %CG Preconditioner
precon = 1; % jf version
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.
printm('Forming spatial-spectral fast object...');
A=spsp_Af(z,f,kz,kf,kp.pointtime,rfp.dfovf,rfp.dfovz);


%%%%%%%%%%%Conjugate Gradient%%%%%%%%%%%%%%%%%%%
printm('Running preconditioned conjugate gradient...');
[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));


if 0 % no, superceded by ir_mri_rf_spsp_plot.m
	printm('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
