 function tbeta = calc_beta(recon_type, kspace, t, fov, target_res, mask, l2_beta, r2s_med, L)
%function tbeta = calc_beta(recon_type, kspace, t, fov, target_res, mask, l2_beta, [r2s_med,] [L])
%
% in:
%       recon_type                      string with type of reconstruction beta is used
%         'dz_recon'                    reconstructing z using linear model
%         'imiter_recon'                reconstructing R2* weighted image
%         'im0iter_recon'               reconstructing m0 (non-R2* weighted)
%       kspace                          kspace trajectory
%       t                               time vector
%       fov                             field of view
%       target_res                      target fwhm [voxels]
%       mask                            image mask
%       l2_beta         [Nb,1]          vector containing log2(beta) to interpolate onto
%       r2s_med                         calc beta using this R2* median (default is 3T avg. from wansapura:99:nrt)
%       L                               number of temp interpolators [must be set if r2s_med is set] (def. L=7)
% out:
%       tbeta                           beta values associated with target resolution

if isvar('r2s_med') & ~isvar('L')
    error 'L has to be set if r2s_med is set'
end

if ~isvar('r2s_med')
    r2s_med = mean([1/39.6e-3 1/55.1e-3]);  % From JMRI 9:531-538 (1999)a
    L = 7;
end

% Size of image
Nd = size(mask);

switch recon_type
case 'dz_recon'
    % Check if target_res is of correct size
    if length(target_res(:)) ~= 2
        error('target_res not of correct size')
    end

    % Build dummy fast_mr objects with the correct k-space trajectory
    z = r2s_med*ones(Nd);
    rho = double(mask);
    if size(t,2) == 1
        A_fast_mr = fast_mr(kspace(:,1),kspace(:,2),fov,Nd(1),2*Nd(1),5,t,-i*z(mask),1,L,1,0,mask);
        A_fast = diff_r2(A_fast_mr,rho(mask));
    else
        str_A_fast = '';
        for ii = 1:size(t,2)
            eval(sprintf('A_fast_mr%i = fast_mr(kspace(:,1),kspace(:,2),fov,Nd(1),2*Nd(1),5,t(:,ii),-i*z(mask),1,L,1,0,mask);',ii));
            eval(sprintf('A_fast%i = diff_r2(A_fast_mr%i,rho(mask));',ii,ii));
            str_A_fast = strcat(str_A_fast,sprintf('A_fast%i,',ii));
        end
        eval(sprintf('A_fast = multi_echo_diffr2(%s);',str_A_fast(1:end-1)));
    end

    % Find the beta value corresponding to the target resolution
    if any(target_res(1:2) == 1)
        warning 'One of the target resolutions equal 1, hence both are set to 1'
        tbeta = zeros(2,1);
    else
        % Build the spatially uniform penalty and calculate the beta value for target_res
        R = Rbuild('tight',mask,8,1); % kappa >= 0
        CC = R.C;

        % Build the range of betas to analyze the data
        [beta1,beta2] = meshgrid(l2_beta);
        Cbeta = 2.^[beta1(:) beta2(:)];
        [psf, fwhm, tbeta] = res2beta(A_fast, CC, @qpwls_psf_st, Cbeta, [0,0], mask, target_res(:));
        if any(isnan(tbeta)), keyboard; end
    end
case 'imiter_recon'
    % Check if target_res is of correct size
    if length(target_res(:)) ~= 1
        error('target_res not of correct size')
    end

    % Build dummy fast_mr objects with the correct k-space trajectory
    % L will be set to 0 in fast_mr( )
    z = zeros(Nd);
    if size(t,2) > 1
        str_A_fast_mr = '';
        for ii = 1:size(t,2)
            eval(sprintf('A_fast_mr%i = fast_mr(kspace(:,1),kspace(:,2),fov,Nd(1),2*Nd(1),5,t(:,ii),-i*z(mask),1,L,1,0,mask);',ii));
            str_A_fast_mr = strcat(str_A_fast_mr,sprintf('A_fast_mr%i,',ii));
        end
        eval(sprintf('A_fast_mr = multi_echo_fastmr(%s);',str_A_fast_mr(1:end-1)));
    else
        A_fast_mr = fast_mr(kspace(:,1),kspace(:,2),fov,Nd(1),2*Nd(1),5,t,-i*z(mask),1,1,1,0,mask);
    end

    % Find the beta value corresponding to the target resolution
    if target_res == 1
        tbeta = 0;
    else
        % Build the spatially uniform penalty and calculate the beta value for target_res
        Rimiter = Rbuild('leak',mask,8,1);
        Cimiter = Rimiter.C;

        % tbeta2 does not include average R2* effects
        [psf, fwhm, tbeta] = res2beta(A_fast_mr, Cimiter, @qpwls_psf, 2.^[l2_beta].', [0,0], mask, target_res);
        if any(isnan(tbeta)), keyboard; end
    end
case 'im0iter_recon'
    % Build dummy fast_mr objects with the correct k-space trajectory
    L = 7;  % Set number of time segments
    z = r2s_med*ones(Nd);
    str_A_fast_mr = '';
    for ii = 1:size(t,2)
        eval(sprintf('A_fast_mr%i = fast_mr(kspace(:,1),kspace(:,2),fov,Nd(1),2*Nd(1),5,t(:,ii),-i*z(mask),1,L,1,0,mask);',ii));
        str_A_fast_mr = strcat(str_A_fast_mr,sprintf('A_fast_mr%i,',ii));
    end
    eval(sprintf('A_fast_mr = multi_echo_fastmr(%s);',str_A_fast_mr(1:end-1)));

    % Find the beta value corresponding to the target resolution
    if target_res == 1
        tbeta = 0;
    else
        % Build the spatially uniform penalty and calculate the beta value for target_res
        Rimiter = Rbuild('leak',mask,8,1);
        Cimiter = Rimiter.C;
        
        % tbeta2 does not include average R2* effects
        [psf, fwhm, tbeta] = res2beta(A_fast_mr, Cimiter, @qpwls_psf, 2.^[l2_beta].', [0,0], mask, target_res);
        if any(isnan(tbeta)), keyboard; end
    end
otherwise
    error('recon_type is unknown');
end

% % Build dummy fast_mr objects with the correct k-space trajectory
% z1 = zeros(N);
% z2 = r2s_med*ones(N);
% rho = double(mask);
% A_fast_mr1 = fast_mr(kspace(:,1),kspace(:,2),fov,N,2*N,5,t,-i*z1(mask),1,L,1,0,mask);
% A_fast_mr2 = fast_mr(kspace(:,1),kspace(:,2),fov,N,2*N,5,t,-i*z2(mask),1,L,1,0,mask);
% A_fast = diff_r2(A_fast_mr2,rho(mask));
% 
% if any(target_res(1:2) == 1)
%     warning 'One of the target resolutions equal 1, hence both are set to 1'
%     tbeta1 = zeros(2,1);
% else
%     % Build the spatially uniform penalty and calculate the beta value for target_res
%     R = Rbuild('tight',mask,8,1); % kappa >= 0
%     CC = R.C;
%     % Build the range of betas to analyze the data
%     [beta1,beta2] = meshgrid(-6:.5:0);
%     Cbeta = 2.^[beta1(:) beta2(:)];
%     [psf1, fwhm1, tbeta1] = res2beta(A_fast, CC, @qpwls_psf_st, Cbeta, [0,0], mask, target_res(1:2));
%     if any(isnan(tbeta1)), keyboard; end
% end
% 
% if target_res(3) == 1
%     tbeta2 = 0;
%     tbeta3 = 0;
% else
%     % Build the spatially uniform penalty and calculate the beta value for target_res
%     Rimiter = Rbuild('leak',mask,8,1);
%     Cimiter = Rimiter.C;
%     % Build the range of betas to analyze the data
%     % tbeta2 does not include average R2* effects
%     [psf2, fwhm2, tbeta2] = res2beta(A_fast_mr0, Cimiter, @qpwls_psf, 2.^[2:.5:10].', [0,0], mask, target_res(3));
%     if any(isnan(tbeta2)), keyboard; end
%     % tbeta3 includes average R2* effects
%     [psf_nor2s, fwhm_nor2s, tbeta_nor2s] = res2beta(A_fast_mr1, Cimiter, @qpwls_psf, 2.^[2:.5:10].', [0,0], mask, target_res(3));
%     if any(isnan(tbeta_nor2s)), keyboard; end
% end
% 
% 
% % Now to make the beta values used later
% tbeta = [tbeta1;tbeta2];
% 
