%
%function fitexp1()
%
% Fit diffusion coefficient
%
function fitD1()

    close all

    imsize = 500;

    text_fig = figure('name','notes','position', ...
                      [250,250,imsize,imsize],...
                      'visible','off') ;

    fexp_fig = figure('name','single exponential fit of data','position', ...
                      [250+imsize,250,imsize,imsize],...
                      'visible','off') ;

    flog_fig = figure('name','line fit to log(data)','position', ...
                      [250+2*imsize,250,imsize,imsize],...
                      'visible','off') ;


    % ----- Units conversion --------- %

    sec2ms = 1e3;                       % seconds to milliseconds
    ms2sec = 1./sec2ms;                 % milliseconds to seconds 
    mm2um = 1e3;                        % millimeter to micrometer (micron)
    um2mm = 1./mm2um;                   % micrometer (micron) to millimeter
    cm2mm = 1e1;                        % centimeter to millimeter
    mm2cm = 1./cm2mm;                   % millimeter to centimeter
    cm2um = 1e4;                        % centimeter to micrometer (micron)
    um2cm = 1./cm2um;                   % micrometer (micron) to centimeter 

    % gyromagnetic ratio (Rad/s/G)
    Gam = 4258*2*pi;         

    % ----- User input --------- %

    fprintf('--------- User input ---------\n');

    n = input('number of b-values: [25] ');
    if isempty(n),
        n = 25;
    end

    snr_def = 10;  
    istr = sprintf('snr: (0=inf) [%g] ',snr_def);
    snr = input(istr);
    if isempty(snr),
        snr = snr_def;
    end
    if snr==0,
        nvar = 0;
    else
        nvar = 1./snr;
    end

    % Diffusion coefficient
    D_def = 2;  
    istr = sprintf('Diffusion coefficient (u^2/ms): [%g] ',D_def);
    D = input(istr);
    if isempty(D),
        D = D_def;
    end
    D = D*(um2mm^2/ms2sec);                 % D in (mm^2/sec)
    Dstr = sprintf('\t Diffusion coefficient = %g mm^2/sec',D);

    % Diffusion gradient amplitudes

    Gmin_def = 0.;
    istr = sprintf('Minimum gradient amplitude (G/cm): [%g] ',Gmin_def);
    Gmin = input(istr);
    if isempty(Gmin),
        Gmin = Gmin_def;
    end
    Gmin = Gmin/cm2mm;

    Gmax_def = 4.;
    istr = sprintf('Maximum gradient amplitude (G/cm): [%g] ',Gmax_def);
    Gmax = input(istr);
    if isempty(Gmax),
        Gmax = Gmax_def;
    end
    Gmax = Gmax/cm2mm;
    Gstr = sprintf('\t Minimum/maximum gradient amplitude = %g/%g Gauss/cm',Gmin,Gmax);

    % Diffusion gradient width
    delta_def = 20;
    istr = sprintf('Gradient duration (delta,ms): [%g] ',delta_def);
    delta = input(istr);
    if isempty(delta),
        delta = delta_def;
    end

    % time to play out 180 pulse (in ms)
    rf180_time = 5;

    % Time between diffusion gradients
    Delta_def = delta + rf180_time;
    istr = sprintf('Time between gradients (Delta,ms): [%g] ',Delta_def);
    Delta = input(istr);
    if isempty(Delta),
        Delta = Delta_def;
    end

    delta = delta*ms2sec;
    Delta = Delta*ms2sec;

    Tstr = sprintf(['\t Timing parameters: delta = %g ms, Delta = %g, ' ...
                    'width of rf(180) = %gms'],delta,Delta,rf180_time);

    % ----- Calculate b-values (sec/mm^2) --------- %

    tau = Delta - delta/3.;
    qmin = Gam*Gmin*delta;
    qmax = Gam*Gmax*delta;
    bmin = qmin.^2*tau;
    bmax = qmax.^2*tau;
    bvals = linspace(bmin,bmax,n)';

    Bstr = sprintf('\tb-values (sec/mm^2); (min,max) = (%g,%g)\n',min(bvals),max(bvals));

    Eqn_str = sprintf('\t signal: s(b) = s(0) e^{-b D}');

    if 1
        % ----- text test --------- %

        descr1 = {' ';
                  ' ';
                  ' ';
                  ' ';
                  'Fitting the diffusion coefficient';
                  Eqn_str;
                  ' ';
                  'Physical parameters:';
                  Dstr;
                  ' ';
                  'Experimental parameters:';
                  Gstr;
                  Tstr;
                  Bstr};

        %h = msgbox(descr1,'CreateMode','replace')

        figure(text_fig)
        fontsize = 18;
        fontweight = 'bold';
        text(.1,.9,descr1,'Fontsize',fontsize,'FontWeight',fontweight);
        axis off
    end

    % --------------------------%

    fprintf('\n *** Running analysis with the following parameters:\n');
    fprintf('Physical parameters:\n');
    fprintf('\tDiffusion coefficient = %g (mm^2/ms)\n',D);
    fprintf('\tSNR = %g (mm^2/ms)\n',snr);

    fprintf('Experimental parameters:\n');
    fprintf('\tGmin = %g (G/cm), Gmax = %g (G/cm)\n',Gmin/mm2cm,Gmax/mm2cm);
    fprintf('\tdelta = %g (ms), Delta = %g (ms)\n',delta*sec2ms,Delta*sec2ms);
    fprintf('\tb-values (sec/mm^2); (min,max) = (%g,%g)\n',min(bvals),max(bvals));

    % Create data

    s0 = 2.;
    sb = s0*exp(-bvals*D) + nvar*randn(size(bvals));

    % -------- Exp fit to data ---------%

    fprintf('Fitting data to exponential decay\n');
    figure(fexp_fig)

    fexp = fit(bvals,sb,'exp1')
    %pp = plot(fexp,bvals,sb,'bo','MarkerFaceColor','b')
    plot(fexp,bvals,sb,'bo')
    grid on
    grid minor
    axis tight
    xlabel('b-value (sec/mm^2)')
    ylabel('s(b) = s(0) e^{-b D}')
    title('exponential fit of exp[-s(b)]')

    % -------- Linear fit to -log(data) ---------%

    fprintf('Fitting -log(data) to line\n');
    figure(flog_fig)

    x = [ones(length(bvals),1) bvals];
    flog = lsqlin(x,-log(sb));

    S0_est = exp(-flog(1));
    D_est = flog(2);
    S_est = S0_est*exp(-bvals*D_est);

    fprintf('\tEstimated S0 = %g\n',S0_est);
    fprintf('\tEstimated diffusion coefficient = %g (mm^2/ms)\n',D_est);

    %plot(bvals,sb,'bo','MarkerFaceColor','b')
    plot(bvals,sb,'bo',bvals,S_est,'r')
    hold off
    grid on
    grid minor
    axis tight
    xlabel('b-value (sec/mm^2)')
    ylabel('s(b) = s(0) e^{-b D}')
    title('linear fit of -log[s(b)]')

    %keyboard
    return