function [t1data]=t1_mn(pfgnmrdata,thresh,Specrange,T1range,FitType, Opts)
%line added to test giy
%   [T1]=t1_mn(pfgnmrdata,thresh,Specrange)
%   DOSY (diffusion-ordered spectroscopy) fitting of
%   PFG-NMR diffusion data (aka DOSY data)
% 
%   -------------------------------INPUT--------------------------------------
%   pfgnmrdata      Data structure containing the PFGNMR experiment.
%                   containing the following members:
%                       filename: the name and path of the original file
%                       np: number of (real) data points in each spectrum
%                       wp: width of spectrum (ppm)
%                       sp: start of spectrum (ppm)
%                       dosyconstant: gamma.^2*delts^2*DELTAprime
%                       Gzlvl: vector of gradient amplitudes (T/m)
%                       ngrad: number of gradient levels
%                       Ppmscale: scale for plotting the spectrum
%                       SPECTRA: spectra (processed)
%   ---------------------------INPUT - OPTIONAL-------------------------------
%   thresh =        Noise threshold. Any data point below won't be used
%                   the fitting. Expressed as % of highest peak. Default 5
%                   (%).
%   Specrange       The spectral range (in ppm) in which the score fitting
%                   will be performed. if set to [0] defaults will be used.
%                   DEFAULT is [sp wp+sp];
%   FitType         The type of Fit e.g. T1 or T2
%
%   --------------------------------------OUTPUT------------------------------
%   t1data         Structure containg the data obtained after t1/t2fit with
%                  the follwing elements:
%
%
%                  fit_time:   The time it took to run the fitting%
%                  wp: width of spectrum (ppm)
%                  sp: start of spectrum (ppm)
%                  Ppmscale: scale for plotting the spectrum
%                  filename: the name and path of the original file
%
%
%   Example:
%
%   See also: dosy_mn, score_mn, decra_mn, mcr_mn, varianimport,
%             brukerimport, jeolimport, peakpick_mn, dosyplot_mn,
%             dosyresidual, dosyplot_gui, scoreplot_mn, decraplot_mn,
%             mcrplot_mn
%
%   References:
%
%
%
%   This is a part of the GNAT
%   Copyright  2017  <Mathias Nilsson>%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%   mathias.nilsson@manchester.ac.uk

% MN      12 July 2017          Fixed correct integrals.
% MN      5  February 2020      Made a better attempt at estimating
% MJS     20 April 2020
            %332 - fitting procedure chnaged to use fitnlm
            %334-362 Old calculations for estimating errors deleted
            %345-347    Read coefficents directly from output
            %fitd.Coefficients
            %366 t1data.FITTED_CALC read directly from output of fitd
            %370 t1data.FITTED changed to read from fitd..Coefficients
            %376,379,382 table read from fitd.Coefficients
% MJS     Jan 2022
            % Biexponential fitting added for relaxation data
            % T1 matrix at end of file changed to allow for ROSY plot of
            % biexponential data


if nargin==0
    disp(' ')
    disp(' T1/T2 fitting')
    disp(' ')
    disp(' Type <<help t1_mn>> for more info')
    disp('  ')
    return
elseif nargin<1
    error(' t1_mn needs a pfgnmrdata stucture as input')
elseif nargin >6
    error(' Too many inputs')
end


if FitType==1
    disp('T1 (Inversion recovery 3 parameter fit)')
    
elseif FitType==2
    disp('T1 (Saturation recovery parameter fit) ')
    
    elseif FitType==3
    disp('T2 (2 parameter fit) ')
    
else
    error('Unknown FitType')
end

t_start=cputime;

%Defaults

t1data.filename=pfgnmrdata.filename;
% User options
if (nargin >1)
    th=thresh;
else
    th=0.1;
end
if nargin>2
    if Specrange==0
        %do nothing
    else
        if length(Specrange)~=2
            error('T1/T2: Specrange should have excatly 2 elements')
        end
        if Specrange(1)<pfgnmrdata.sp
            disp('T1/T2: Specrange(1) is too low. The minumum will be used')
            Specrange(1)=pfgnmrdata.sp;
        end
        if Specrange(2)>(pfgnmrdata.wp+pfgnmrdata.sp)
            disp('T1/T2: Specrange(2) is too high. The maximum will be used')
            Specrange(2)=pfgnmrdata.wp+pfgnmrdata.sp;
        end
        for k=1:length(pfgnmrdata.Ppmscale)
            if (pfgnmrdata.Ppmscale(k)>Specrange(1))
                begin=k-1;
                k1=begin;
                break;
            end
        end
        
        for k=begin:length(pfgnmrdata.Ppmscale)
            if (pfgnmrdata.Ppmscale(k)>=Specrange(2))
                endrange=k;
                break;
            end
        end
        %make a new stucture
        pfgnmrdata.sp=pfgnmrdata.Ppmscale(k1);
        pfgnmrdata.wp=pfgnmrdata.Ppmscale(endrange)-pfgnmrdata.Ppmscale(k1);
        pfgnmrdata.Ppmscale=pfgnmrdata.Ppmscale(k1:endrange);
        pfgnmrdata.SPECTRA=pfgnmrdata.SPECTRA(k1:endrange,:);
        pfgnmrdata.np=length(pfgnmrdata.Ppmscale);
    end
end
if nargin>3
    if (length(T1range) ~= 3)
        if T1range==0
            %use defaults
            T1range=[0 10 256];
        else
            error('T1/T2: Diffrange is a vector of size 3')
        end
    end
end
AutoPlotLimits=0;
if T1range(1)==0 && T1range(2)==0
    %automatically find the plot limits
    AutoPlotLimits=1;
end
fn1=T1range(3);

[nspec,~]=size(pfgnmrdata.SPECTRA);

T1=zeros(nspec,fn1);
%Opts(1)=0;
if Opts(1)==0 %Use peak-picking
    [peak]=peakpick_mn(pfgnmrdata.SPECTRA(:,pfgnmrdata.flipnr),th);
    hp = waitbar(0,'T1/T2: Fitting after peak picking');
    [npeaks]=length(peak);
    disp(['T1/T2: Fitting on: ' num2str(npeaks) ' peaks'])
elseif (Opts(1)==1)
    peak(pfgnmrdata.np)=struct('max',0,'start',0,'stop',0);
    %npeaks=nspec;
    pmp=0;
    for k=1:pfgnmrdata.np
        if pfgnmrdata.SPECTRA(k,pfgnmrdata.flipnr)>max(pfgnmrdata.SPECTRA(:,pfgnmrdata.flipnr))*th/100
            pmp=pmp+1;
            peak(pmp).max=k;
            peak(pmp).start=k;
            peak(pmp).stop=k;
        end
    end
    [npeaks]=pmp;
    disp('T1/T2: Fitting each frequency');
    hp = waitbar(0,'T1/T2: Fitting each frequency');
elseif (Opts(1)==2) %use integrals.
    
    peak=pfgnmrdata.integral_peaks;
    [npeaks]=numel(peak);
    for k=1:npeaks        
        [~, Index]=min(abs(pfgnmrdata.Ppmscale - peak(k).max));
        peak(k).max=Index;
        
        [~, Index]=min(abs(pfgnmrdata.Ppmscale - peak(k).start));
        peak(k).start=Index;
        
        [~, Index]=min(abs(pfgnmrdata.Ppmscale - peak(k).stop));
        peak(k).stop=Index;
        
        
    end
    disp(['T1/T2: Fitting on: ' num2str(npeaks) ' integral regions'])
    hp = waitbar(0,'T1/T2: Fitting integral regions');

else
    error('T1/T2: Illegal (Opts(1)')
end


if FitType==1
    disp('T1 (Inversion recovery 3 parameter fit)')
    t1data.FITSTATS=zeros(npeaks,6);
elseif FitType==2
    disp('T1 (Saturation recovery parameter fit) ')
    t1data.FITSTATS=zeros(npeaks,6);
elseif FitType==3
    disp('T2 (2 parameter fit) ')
    t1data.FITSTATS=zeros(npeaks,4);
else
    error('Unknown FitType')
end


if Opts(2)==0 || Opts(2)>1
    disp('ROSY: Using monoexponential fitting')
    Opts(3)=1;
elseif  Opts(2)==1 && Opts(3)>1
    disp(['ROSY: Trying to fit ' num2str(Opts(3)) ' componets per peak'])
else
    disp('ROSY: Using monoexponential fitting only')
    Opts(3)=1;
end
if Opts(4)==0
    Opts(4)=100;
end

%opts=optimset('lsqcurvefit');

opts=statset('Display','off');
opts=statset(opts,'TolFun',1e-6);
opts=statset(opts,'TolX',1e-6);
opts=statset(opts,'Jacobian','on');
opts=statset(opts,'MaxFunEvals',10000);


for i=1:npeaks

    disp(['T1/T2: Peak: ' num2str(i) ])
    waitbar(i/npeaks);
    if Opts(1)==0 %Use peak-picking
        ydata=pfgnmrdata.SPECTRA(peak(i).max,:);
    elseif (Opts(1)==1) %use all data points
        ydata=pfgnmrdata.SPECTRA(peak(i).max,:);
    elseif (Opts(1)==2) %use integrals
        ydata=sum(pfgnmrdata.SPECTRA(peak(i).start:peak(i).stop,:));
    else
        error('T1/T2: Illegal (Opts(1)')
    end
    

   % ydata=pfgnmrdata.SPECTRA(peak(i).max,:);
    %ydata=[-35.56 -34.57 -32.51 -29.52 -23.03 -12.58 3.59 22.64 35.94 39.50 40.22 40.12] %should give T1 1.081
%     y_norm=norm(ydata); % added by GDP 31x2018
%     %ydata
%     ydata=ydata/y_norm;
    y_norm=1;
    
    if FitType==1 %T1 Inversion recovery
        xdata=pfgnmrdata.d2;
        %fitparm(2)=0.1 ;%should do a better guess from logfit
        %estimating T1 from a simple log
        logy=log(ydata(end)-ydata);
        tmp_slope=(logy(2)-logy(end-1))/(xdata(2)-xdata(end-1));
        fitparm(2)=abs(1/tmp_slope);
        
        M0=ydata(end);
        fitparm(1)=-M0;
        fitparm(3)=M0;
        model=@t1;
    elseif FitType==2 %T1 Saturation recovery
        xdata=pfgnmrdata.d2;
        fitparm(2)=1 ;%should do a better guess from logfit
        %estimating T1 from a simple log
        logy=log(ydata);
       % xdata;
        tmp_slope=(logy(2)-logy(end-1))/(xdata(2)-xdata(end-1));
        fitparm(2)=abs(1/tmp_slope);
        
        M0=ydata(end);
        fitparm(1)=0;
        fitparm(3)=M0;
        model=@t1;        
    elseif FitType==3 %T2
        if strcmp(pfgnmrdata.type,'Bruker')
            xdata=pfgnmrdata.d2;
            disp('Bruker data - using d2 as arrayed times')
        elseif strcmp(pfgnmrdata.type,'Varian')
            xdata=pfgnmrdata.bigtau;
            disp('Varian data - using bigtau as arrayed times')
        elseif strcmp(pfgnmrdata.type,'Jeol')
            xdata=pfgnmrdata.d2;
            disp('JEOL data - using array in Y dimension')
        else
            error('Can not determine data type')
        end
        %fitparm(2)=1 ;%should do a better guess from logfit
        logy=log(ydata);
        tmp_slope=(logy(2)-logy(end-1))/(xdata(2)-xdata(end-1));
        fitparm(2)=abs(1/tmp_slope);
        M0=ydata(1);
        fitparm(1)=M0;
        % fitparm(3)=0;
        fitparm=fitparm(1:2);
        model=@t2;

    else
        error('Unknown FitType')
    end

    try
        fitd=fitnlm(xdata,ydata,model,fitparm,'Options',opts);
    catch ME
        disp('Error Message:')
        disp(ME.message)
        disp(['Ignoring Peak: ' num2str(i)])
    end

   residual(i,:)=fitd.Residuals{:,1}; %#ok<AGROW>
   resnorm=fitd.SSE;


    sderrt1(i)=fitd.Coefficients{2,2};
    fittedt1_mono(i)=fitd.Coefficients{2,1};
    fittedamp(i)=fitd.Coefficients{1,1};
    sderramp(i)=fitd.Coefficients{1,2};
    nmult=1;

%% MJS adding multiexponential fitting

if Opts(3)>1
    
    nmult=Opts(3);
        %opts_m=optimset(opts,'Display','off'); %#ok<NASGU>
        opts_m=optimset(opts,'Jacobian','on');
        opts_m_fmin=optimset('fminsearch');
        opts_m_fmin=optimset(opts_m_fmin,'Display','Off');
        opts_m_fmin=optimset(opts_m_fmin,'MaxFunEvals',10000);
        
        resnorm_mono=resnorm;
        %fittedt1_mono= fittedt1_mono(i);
        residual_mono=residual(i);
        sderr_mono=sderrt1(i);
        
        Fit=Opts(6);%
        
    breakflag=0;
    
     
   if FitType==3 && Fit==0 % Monte Carlo minimisation
       model=@multi_t2;
          
         
              nMc=2;
            u_limit=100*ones(1,nmult*2); %not allowing a diffusion coefficient over 100
            u_limit=u_limit(:);
            l_limit=zeros(1,nmult*2); %non-negativty
            nMc_real= round(sqrt(nMc).^nmult);
            CM_best=0; %#ok<NASGU>
            fitd_best=0;
            sderr_best=0;
            resnorm_best=0;
            amps=fittedamp(i);
            ex=fittedt1_mono(i);
            
           %Generate random amplitues and exponent start values
           gamp=abs(normrnd(amps,0.05,[nMc,nmult]));
           gexp=abs(normrnd(ex,0.05,[nMc,nmult]));
           
            nn=1;
            for t=1:nmult
                 gmultifit(:,nn)=gamp(:,t);
                 gmultifit(:,nn+1)=gexp(:,t);
                   nn=nn+2;
            end 
            
            for ii=1:nMc_real
               x0=gmultifit(ii,:);
              
               fitd_m=fitnlm(xdata,ydata,model,x0,'Options',opts);
               
                    
                if ii==1
                        resnorm_best=fitd.SSE;
                        residual_best=fitd_m.Residuals;
                        fitd_best=fitd_m;
                elseif fitd_m.SSE<resnorm_best
                    fitd_best=fitd_m;
                    resnorm_best=fitd_m.RMSE;
                    sderr_best=fitd_m.SSE;
                    residual_best=fitd_m.Residuals;
                    assignin('base','fitd',fitd_best)
                end
            end
            
            if resnorm_best<resnorm_mono && isreal(sderr_best)==1
                for mm=1:nmult*2
                    if resnorm_best>fitd_m.Coefficients{mm,2}*0.2
                        %if stdev is more than 30% we reject
                        %the fit and revert to the monoexponetial values
                        %fitd=fitd_mono;
                        fittedt1(i)=fittedt1_mono(i);
                        sderrt1(i)=sderr_mono;
                        residual(i)=residual_mono;
                        breakflag=1;
                        nmult=nmult-1;
                    
                    end
                end
                
            end
            
            if breakflag==0
                % the multifit is accepted and we use this result
                % we dont use either mono or bi fit if tri is successful
                
                bm=1;
                
                fitd_m=fitd_best;
                % MJS adding loop for biexponentials
                for b=1:size(fitd_m.Coefficients)/2
                fittedamp(b,i)=fitd_m.Coefficients{bm,1};
                sderramp(b,i)=fitd_m.Coefficients{bm,2}; 
                fittedt1(b,i)=fitd_m.Coefficients{bm+1,1};
                sderrt1(b,i)=fitd_m.Coefficients{bm+1,2};
                bm=bm+2;
                end
               nmult_fit(i)=length(fittedt1(:,i));
            end  

%% Multi T1 fit with Globalsearch         
       elseif FitType==1 
           disp('T1 multiexponential fitting currently unsupported reverting to monoexponential fit')
%                     % Initial values for global search taken from the
%                     % monoexponential fit
%                     breakflag=0;
%                     u_limit=100*ones(1,nmult*3); %not allowing a relaxation coefficient over 100
%                     l_limit=-100*ones(1,nmult*3); %non-negativty
%              
%                     amps=fittedamp(i);
%                     ex=fittedt1_mono(i);
%                     gamp=normrnd(amps,1,[1,nmult])/2;
%                     gexp=abs(normrnd(ex,0.1,[1,nmult]));
%                     ampi=normrnd(fitd.Coefficients{3,1},1,[1,nmult])/2;
%                     
%                     multifit=[ampi; gamp; gexp];
%                     
%                     multifit=multifit(:)';
%                     assignin('base','multifit',multifit)
%                     %multifit=[1 1 1 1 1 1];
%                     
%                     model=@multi_t1;
%                     RNCF=@(b) norm(ydata- model(b,xdata));
%                     gs=GlobalSearch('Display','off');
%                     problem = createOptimProblem('fmincon','x0',multifit,...
%                     'objective',RNCF);%,'lb',l_limit,'ub',u_limit);
%                     x0=run(gs,problem);
%  
%                         fitd_m=fitnlm(xdata,ydata,model,x0,'Options',opts);
%                         resnorm_best=fitd.SSE;
%                         residual_m=fitd_m.Residuals;
%                         assignin('base','fitd_m',fitd_m);
% % 
% %                 else
% %                     error('unknown FitType')
%                 
%      
%             %check if we like the multiexponential fit more than the
%             %monoexponential
%             if resnorm_best<resnorm_mono && isreal(sderr_best)==1
%                 for mm=1:nmult*2
%                     if fitd_m.Coefficients{mm,2}>fitd_m.Coefficients{mm,1}*0.2
%                         %if stdev is more than 30% we reject
%                         %the fit and revert to the monoexponetial values
%                         %fitd=fitd_mono;
%                         fittedt1=fittedt1_mono;
%                         sderr=sderr_mono;
%                         residual(:,i)=residual_mono;
%                         breakflag=1;
%                         nmult=nmult-1;
%                     end
%                 end
%             end
%             
%             if breakflag==0
%                 % the multifit is accepted and we use this result
%                 % we dont use either mono or bi fit if tri is successful
%                 
%                 bm=1;
%                 
%                 
%                 % MJS adding loop for biexponentials
%                 for b=1:size(fitd_m.Coefficients)/2
%                 fittedamp(b,i)=fitd_m.Coefficients{bm,1};
%                 sderramp(b,i)=fitd_m.Coefficients{bm,2}; 
%                 fittedt1(b,i)=fitd_m.Coefficients{bm+1,1};
%                 sderrt1(b,i)=fitd_m.Coefficients{bm+1,2};
%                 bm=bm+2;
%                 end
%                nmult_fit(i)=length(fittedt1(:,i));
%             end
            %%  Multi t2 fit with Global search
            elseif FitType==3 && Fit==1
                    % Initial values for global search taken from the
                    % monoexponential fit
                    breakflag=0;
                    u_limit=100*ones(1,nmult*2); %not allowing a relaxation coefficient over 100
                    l_limit=zeros(1,nmult*2); %non-negativty
             
                    amps=fittedamp(i);
                    ex=fittedt1_mono(i);
                    gamp=abs(normrnd(amps,1,[1,nmult]));
                    gexp=abs(normrnd(ex,1,[1,nmult]));
                    
                    multifit=[gamp; gexp];
                    multifit=multifit(:)';
                    
                    
                    model=@multi_t2;
                    RNCF=@(b) norm(ydata- model(b,xdata));
                    gs=GlobalSearch('Display','off');
                    problem = createOptimProblem('fmincon','x0',multifit,...
                    'objective',RNCF,'lb',l_limit,'ub',u_limit);
                    x0=run(gs,problem);

                        fitd_m=fitnlm(xdata,ydata,model,x0,'Options',opts);
                        resnorm_best=fitd.SSE;
                        residual_m=fitd_m.Residuals;
                        assignin('base','fitd_m',fitd_m);
% 
%                 else
%                     error('unknown FitType')
                
     
            %check if we like the multiexponential fit more than the
            %monoexponential
            if resnorm_best<resnorm_mono && isreal(sderr_best)==1
                for mm=1:nmult*2
                    if fitd_m.Coefficients{mm,2}>fitd_m.Coefficients{mm,1}*0.2
                        %if stdev is more than 30% we reject
                        %the fit and revert to the monoexponetial values
                        %fitd=fitd_mono;
                        fittedt1=fittedt1_mono;
                        sderrt1=sderr_mono;
                        residual(i,:)=residual_mono(i,:);
                        breakflag=1;
                        nmult=nmult-1;
                    end
                end
            end
            
            if breakflag==0
                % the multifit is accepted and we use this result
                % we dont use either mono or bi fit if tri is successful
                
                bm=1;
                
                
                % MJS adding loop for biexponentials
                for b=1:size(fitd_m.Coefficients)/2
                fittedamp(b,i)=fitd_m.Coefficients{bm,1};
                sderramp(b,i)=fitd_m.Coefficients{bm,2}; 
                fittedt1(b,i)=fitd_m.Coefficients{bm+1,1};
                sderrt1(b,i)=fitd_m.Coefficients{bm+1,2};
                bm=bm+2;
                end
               nmult_fit(i)=length(fittedt1(:,i));
            end 
            
  end
end 


           
    %SDERR{i}=sderr;
   % nmult_fit(i)=length(sderr)/2;
 
%% fill out table
    t1data.Tau=xdata;
    t1data.X_CALC=linspace(t1data.Tau(1), t1data.Tau(end), 256);
    t1data.FITTED_CALC(i,:)=model(fitd.Coefficients{:,1},t1data.X_CALC);
    t1data.freqs(i)=pfgnmrdata.Ppmscale(peak(i).max);
    t1data.RESIDUAL(i,:)=residual(i,:);
    assignin('base','residuals',t1data.RESIDUAL)
    t1data.ORIGINAL(i,:)=ydata;
    t1data.FITTED(i,:)=model(fitd.Coefficients{:,1},xdata);
    % GDP 31x2018: "y_norm" was multiplied to all M(0) and their error here,
    %              to use for LC & GDP PS quantification
    
    if FitType==1 && nmult==1 %T1 Inversion recovery
        fittedt1=fittedt1_mono;
        nmult_fit=ones(1,npeaks);
        
       t1data.FITSTATS(i,1:6)=[y_norm.*fittedamp(i) y_norm*fitd.Coefficients{1,2} y_norm*fitd.Coefficients{3,1} y_norm*fitd.Coefficients{3,2} fittedt1(i) sderrt1(i)];
    elseif FitType==2 && nmult==1 %T1 Saturation recovery
         fittedt1=fittedt1_mono;
        nmult_fit=ones(1,npeaks);
             t1data.FITSTATS(i,1:6)=[y_norm*fittedamp(i) y_norm*fitd.Coefficients{1,2} y_norm*fitd.Coefficients{3,1} y_norm*fitd.Coefficients{3,2} fittedt1(i) sderrt1(i)];
    elseif FitType==3 && nmult==1 %T2 PROJECT/CPMG
        nmult_fit=ones(1,npeaks);
        fittedt1=fittedt1_mono;
        t1data.FITSTATS(i,1:4)=[y_norm*fittedamp(i) y_norm*sderramp(i) fittedt1(i) sderrt1(i)];
    elseif FitType==3 && nmult>1
       kk=1;
       kkk=1;
        for bm=1:nmult
        t1data.FITSTATS(i,kkk:kkk+3)=[y_norm*fittedamp(kk,i) y_norm*sderramp(kk,i) fittedt1(kk,i) sderrt1(kk,i)];
        kk=kk+1;
        kkk=kkk+4;
        end 
    elseif FitType==1 && nmult>1
        fittedt1=fittedt1_mono;
        nmult_fit=ones(1,npeaks);
        
       t1data.FITSTATS(i,1:6)=[y_norm.*fittedamp(i) y_norm*fitd.Coefficients{1,2} y_norm*fitd.Coefficients{3,1} y_norm*fitd.Coefficients{3,2} fittedt1(i) sderrt1(i)];
    
%        kk=1;
%        kkk=1;
%         for bm=1:nmult
%         t1data.FITSTATS(i,kkk:kkk+3)=[y_norm*fittedamp(kk,i) y_norm*sderramp(kk,i) fittedt1(kk,i) sderrt1(kk,i)];
%         kk=kk+1;
%         kkk=kkk+4;
%         end  
%     else
     %    error('Unknown FitType')
    end
    
    %d=npeaks;
    %t1data.FITSTATS(i,1:6)=[fittedamp(i) tmp(1) M0 tmp(3) fittedt1(i) sderrt1(i)];
    
end



% figure
% plot(t1data.FITTED(1,:))
% hold on
% plot(t1data.ORIGINAL(1,:))
%tt=max(fittedt1);

if AutoPlotLimits==1
    t1min=0;
    [Val, Idx]=max(fittedt1(:));
    t1max=1.25*(Val+sderrt1(Idx));
else
    t1min=T1range(1);
    t1max=T1range(2);
end
xdiff=linspace(t1min,t1max,fn1);
dd=xdiff(2) - xdiff(1);

for i=1:npeaks
    
    
    fitd=fittedt1(:,i)';
    sderr=sderrt1(:,i)';
    ii=1;
    for jj=1:nmult_fit(i)
        T1(peak(i).start:peak(i).stop,:)=T1(peak(i).start:peak(i).stop,:)+...
            (abs(pfgnmrdata.SPECTRA(peak(i).start:peak(i).stop,1)))...
            *0.5*fitd(ii)*(erf( (xdiff + dd/2 -fitd(jj))/sderr(jj)) -...
            erf( (xdiff - dd/2 -fitd(jj))/sderr(jj)));
    end
end
assignin('base','T1',T1)

endt=cputime-t_start;
h=fix(endt/3600);
m=fix((endt-3600*h)/60);
s=endt-3600*h - 60*m;
fit_time=[h m s];
disp(['T1: Fitting time was: ' num2str(fit_time(1)) ' h  ' num2str(fit_time(2)) ' m and ' num2str(fit_time(3),2) ' s']);

t1data.type='t1/t2';
t1data.fit_time=fit_time;
t1data.T1=T1;
t1data.FitType=FitType;
if FitType==1 %T1
    t1data.Spectrum=pfgnmrdata.SPECTRA(:,end);
elseif FitType==2 %T1
    t1data.Spectrum=pfgnmrdata.SPECTRA(:,end);
elseif FitType==3 %T2
    t1data.Spectrum=pfgnmrdata.SPECTRA(:,1);
else
    error('Unknown FitType')
end
t1data.Ppmscale=pfgnmrdata.Ppmscale';
t1data.sp=pfgnmrdata.sp;
t1data.wp=pfgnmrdata.wp;
t1data.np=pfgnmrdata.np;
t1data.fn1=fn1;
if (isfield(pfgnmrdata,'NUC1'))
t1data.NUC1=pfgnmrdata.NUC1;
end 
% t1data.dmin=dmin;
% t1data.dmax=dmax;

t1data.T1scale=linspace(t1min,t1max,fn1);
t1data.threshold=th;
t1data.d2=xdata;
%t1data
close(hp);
%t1plot(t1data,5,th);
%dosyplot_gui(t1data,5,th);
% 
% max(max(t1data.T1))
% % min(min(t1data.T1))
% 
% figure
% contour(t1data.T1)
%-----------------------END of t1_mn-------------------------------------

%-----------------Auxillary functions--------------------------------------

    function [y, J] = t1(a,xdata)
        % for fitting a number of pure exponentials
        % this function works for both T1 and T2 but the initial guesses
        % should be different. However, I will use a 2 parameter for T2 as both signal
        %and noise should decay to zero.
        %
        %T1: M(t)=(M(0) - M0)*exp(-t/T1)+M0
        %where M0 is the equilibrium Z magnetization and M(0) is the magnetization
        %at time zero (e.g., immediately after the 180° pulse for an inversion
        %recovery T1 experiment). Notice that this equation will fit inversion
        %recovery data (for which M(0) is approximately equal to -M0) or saturation
        %recovery data (for which M(0) is 0).
        %
        %T2:  M(t) = (M(0) - M(inf))*exp(-t/T2) + M(inf)
        %where M(0) is the magnetization at time zero (i.e., the full
        %magnetization excited by the observe pulse) and M(inf) is the
        %xy-magnetization at infinite time (zero unless the peak is sitting on an
        %offset baseline).        
        
        %xdata=xdata';
        %y=zeros(1,length(xdata));
        y =(a(1)-a(3))*exp(-xdata/a(2)) + a(3);
        if nargin>1
            J=zeros(length(xdata),length(a));
            J(:,1)=exp(-xdata./a(2));
            J(:,2)=((a(1) - a(3))*exp(-xdata/a(2)).*xdata)./a(2).^2;
            J(:,3)=1- exp(-xdata./a(2));
            
        end
    end
    function [y, J] = t2(a,xdata)
        %T2:  M(t) = (M(0) - M(inf))*exp(-t/T2) + M(inf)
        %where M(0) is the magnetization at time zero (i.e., the full
        %magnetization excited by the observe pulse) and M(inf) is the
        %xy-magnetization at infinite time (zero unless the peak is sitting on an
        %offset baseline). So using:
        %T2: M(t) = M(0)*exp(-t/T2)
         %xdata=xdata';
        %y=zeros(1,length(xdata));
        y =a(1)*exp(-xdata/a(2));
        if nargin>1
            J=zeros(length(xdata),length(a));
            J(:,1)=exp(-xdata./a(2));
            J(:,2)=(a(1)*exp(-xdata/a(2)).*xdata)./a(2).^2;
        end
    end
    function [y,J] = multi_t2(a,xdata)
        % for fitting a number of pure exponentials
      
         q=1;
        p=length(a)/2;
       y=zeros(1,length(xdata(1,:)));
        for r=1:p
        y =y+a(q)*exp(-xdata/a(q+1));
        q=q+2;
        end
    end
    function [y,J] = multi_t1(a,xdata)
        % for fitting a number of pure exponentials
      
        q=1;
        p=length(a)/3;
        y=zeros(1,length(xdata(1,:)));
        for r=1:p
        y =y+(((a(q)-a(q+2))*(exp(-xdata/a(q+1))))+a(q+2));
        q=q+3;
        end
end

end
