%% Simulation and Fitting

%NOTE - FOR MORE HELP ABOUT THE SYNTAX FOR ANY COMMAND, TYPE IN MATLAB:
%>> help routine_name
%
%e.g.
%
%>> help gen_sqw
%
% prints in the Matlab command window the help for the gen_sqw routine
%

%=====================

%In this script we reveal the power of Horace to perform quantitative
%analysis of the data. Note that we have provided some fitting routines that can be used as templates

%=== Simulating a pre-prepared S(Q,w) function ========

%2d slice
my_slice=cut_sqw(sqw_file,proj,[-3,0.05,3],[-1.1,-0.9],[-0.1,0.1],[0,4,280]);

%array of 1d cuts
energy_range=[80:20:160];
for i=1:numel(energy_range)
    my_cut(i)=cut_sqw(sqw_file,proj,[-3,0.05,3],[-1.1,-0.9],[-0.1,0.1],[energy_range(i)-10,energy_range(i)+10]);
end

%simulate on sqw objects
sim_slice=sqw_eval(my_slice,@sr122_xsec,[1,0,0,35,-5,15,10,0.1]);
sim_cut=sqw_eval(my_cut,@sr122_xsec,[1,0,0,35,-5,15,10,0.1]);

%repeat on dnd objects
sim_slice_dnd=sqw_eval(d2d(my_slice),@sr122_xsec,[1,0,0,35,-5,15,10,0.1]);
sim_cut_dnd=sqw_eval(d1d(my_cut),@sr122_xsec,[1,0,0,35,-5,15,10,0.1]);

plot(sim_slice); keep_figure;
plot(sim_slice_dnd); keep_figure;

acolor blue
dl(sim_cut(1));
acolor red
pl(sim_cut_dnd(1));
keep_figure;

%note the differences between simulations of notionally the same data. This
%is because dnd just takes the centre point of the integration range,
%whereas sqw takes all of the contributing detector pixels. This is
%imperative if the dispersion varies significantly in a direction
%perpendicular to your cut/slice, as it introduces broadening that the dnd
%simulation fails to capture.

%=== Simulate a peak function with a cut =====

peak_cut=func_eval(my_cut(1),@multigauss,[1,-1,0.1,2,0,0.1,1,1,0.1]);
acolor black
dl(peak_cut);


%======== Fit a single cut with a peak function ==================

%Allow all parameters to be free
pars_in=[0.4,-0.7,0.1,0.5,-0.2,0.1,0.5,0.2,0.1,0.4,0.6,0.1,0.4,1.3,0.1];
[wfit,fitdata]=fit_func(my_cut(1)-0.3,@multigauss,pars_in);

acolor black
plot(my_cut(1)-0.3);
acolor red
pl(wfit);

%Let's keep the widths of all the peaks fixed, but allow heights and
%centres to vary;
pars_in=[0.4,-0.8,0.1,0.5,-0.22,0.1,0.5,0.22,0.1,0.4,0.8,0.1,0.4,1.2,0.1];
pars_free=[1,1,0,1,1,0,1,1,0,1,1,0,1,1,0];
[wfit,fitdata]=fit_func(my_cut(1)-0.35,@multigauss,pars_in,...
    pars_free);

%That wasn't so good, fix with more realistic values now. Instead bind some
%of the positions to follow symmetry (i.e. position of peaks for Q<0 are 
%reflection of those at Q>0)
pars_in=[0.4,-0.8,0.07,0.5,-0.22,0.07,0.5,0.22,0.07,0.4,0.8,0.07,0.4,1.2,0.07];
pars_free=[1,1,0,1,1,0,1,1,0,1,1,0,1,1,0];
pars_bind={{2,11,0,-1},{5,8,0,-1}};%ensures symmetry about x=0
[wfit,fitdata]=fit_func(my_cut(1)-0.35,@multigauss,pars_in,...
    pars_free,pars_bind);

%Repeat the above, but using some of the options
pars_in=[0.4,-0.8,0.07,0.5,-0.22,0.07,0.5,0.22,0.07,0.4,0.8,0.07,0.4,1.2,0.07];
pars_free=[1,1,0,1,1,0,1,1,0,1,1,0,1,1,0];
pars_bind={{2,11,0,-1},{5,8,0,-1}};%ensures symmetry about x=0
[wfit,fitdata]=fit_func(my_cut(1)-0.35,@multigauss,pars_in,...
    pars_free,pars_bind,...
    'list',2);
% [wfit,fitdata]=fit_func(my_cut(1)-0.35,@multigauss,pars_in,...
%     pars_free,pars_bind,...
%     'list',2,'keep',[-1,1]);

%============ Make dispersion plots ==============================

lattice=[alatt,angdeg];
rlp=[0,0,0; 0,0,1; 0,0,0; 1,0,0; 0,0,0; 1,1,0; 0,0,0; 1,1,1];%reciprocal lattice points to draw dispersion between 
pars=[1,0.05,0.05,35,-5,15,10,0.1];%input parameters
ecent=[0,0.1,200];%energy grid
fwhh=2;%energy broadening term
disp2sqw_plot(lattice,rlp,@sr122_disp,pars,ecent,fwhh);

%============== Fit a single cut with an S(Q,w) model ================

%Use the sr122 x-sec here, and fit only a small section...
my_new_cut=cut_sqw(sqw_file,proj,[0.5,0.05,1.5],[-1.1,-0.9],[-0.1,0.1],[100,120]);

pars=[1,0,0,35,-5,15,10,0.1];
pfree=[1,0,0,1,1,1,1,1];

[wfit,fitdata]=fit_sqw(my_new_cut,@sr122_xsec,pars,pfree,'list',1);

acolor black
plot(my_new_cut);
acolor red
pl(wfit)

%==== FM spin-waves function ===========
test_slice=sqw_eval(my_slice,@Fe_FM_spinwaves_FF,[250 0 2.4 5 1]);
plot(test_slice);

pars=[350,0,2.4,5,1];
pfree=[1,1,1,1,1];
[wfit,fitdata]=fit_sqw(my_cut(1)-0.2,@Fe_FM_spinwaves,pars,pfree,'list',1);


%Now add a linear background a linear background to your function
bgpars=[0.08,0.001];
bgfree=[1,0];
[wfit,fitdata]=fit_sqw(my_new_cut,@Fe_FM_spinwaves,pars,pfree,...
    @linear_bg,bgpars,bgfree,'list',1);%notice that Horace treats the foreground and background separately


%Add the Fe form factor to your cross-section:
pars=[250,0,2.4,10,5];
pfree=[1,0,1,0,1];
bgpars=[0.08,0.001];
bgfree=[1,0];
[wfit,fitdata]=fit_sqw(my_new_cut,@Fe_FM_spinwaves_FF,pars,pfree,...
    @linear_bg,bgpars,bgfree,'list',1);


%======Fit multiple cuts simultaneously with a single S(Q,w) model=========

%We will use the array of 1d cuts we made earlier

%To begin just use the same input as above, i.e. single parameter set and a
%single set of parameters for the background functions
pars=[250,0,2.4,10,5];
pfree=[1,0,1,0,1];
bgpars=[0.08,0.001];
bgfree=[1,0];
[wfit,fitdata]=multifit_sqw(my_cut,@Fe_FM_spinwaves_FF,pars,pfree,...
    @linear_bg,bgpars,bgfree,'list',1);

for i=1:numel(my_cut)
    acolor black
    plot(my_cut(i));
    acolor red
    pl(wfit(i));
    %lx -2 2
    keep_figure;
end

%These fits are pretty good, but let's do something a bit more
%sophisticated with the backgrounds
bgfunc={@linear_bg,@linear_bg,@linear_bg,@quadratic_bg,@quadratic_bg};%different background func for different cuts
bgpars={[0.37,0],[0.2,0],[0.14,0],[0.08,0,0],[0.03,0,0]};%use different initial guesses and different free/fixed parameters for the background
bgfree={[1,1],[1,1],[1,1],[1,1,1],[1,1,1]};

[wfit,fitdata]=multifit_sqw(my_cut,@Fe_FM_spinwaves_FF,pars,pfree,...
    bgfunc,bgpars,bgfree,'list',1);

for i=1:numel(my_cut)
    acolor black
    plot(my_cut(i));
    acolor red
    pl(wfit(i));
    %lx -2 2
    keep_figure;
end

%Look at a slice with our fit parameters:
test_slice=sqw_eval(my_slice,@Fe_FM_spinwaves_FF,fitdata.p);
plot(test_slice); keep_figure;
plot(my_slice); keep_figure;
