Powder fitting

[mducle]

Adds a method to fit powder data.

To test, run this script with the attached data:
mnf2.mat.zip

cutpow.m:

function [xx, yy, ee] = cutpow(powstruct, qmin, qmax, emin, emax)
    xx = powstruct.x{1}.';
    if nargin < 5, emax = max(xx); end
    if nargin < 4, emin = min(xx); end
    ie0 = min(find(xx > emin));
    ie1 = max(find(xx <= emax));
    qq = powstruct.x{2};
    iq0 = min(find(qq > qmin));
    iq1 = max(find(qq <= qmax));
    yy = sum(powstruct.y(iq0:iq1, ie0:ie1), 1);
    ee = sqrt(sum(powstruct.e(iq0:iq1, ie0:ie1).^2, 1));
    xx = xx(ie0:ie1);
    if nargout == 1
        xx = struct('x', xx, 'y', yy, 'e', ee, 'qmin', qmin, 'qmax', qmax);
    end
end

mnf2dat = load('mnf2.mat');
q0 = [0.7, 1.3, 2.0];
clear cts;
for iq = 1:numel(q0)
    cts(iq) = cutpow(mnf2dat, q0(iq)-0.1, q0(iq)+0.1, 1.0);
end

J1 = -0.05;
J2 = 0.3;
D = 0.05;

mnf2 = spinw;
mnf2.genlattice('lat_const', [4.87 4.87 3.31], 'angle', [90 90 90]*pi/180, 'sym', 'P 42/m n m');
mnf2.addatom('r', [0 0 0], 'S', 2.5, 'label', 'MMn2', 'color', 'b')
mnf2.gencoupling('maxDistance', 5)
mnf2.addmatrix('label', 'J1', 'value', J1, 'color', 'red');
mnf2.addmatrix('label', 'J2', 'value', J2, 'color', 'green');
mnf2.addcoupling('mat', 'J1', 'bond', 1)
mnf2.addcoupling('mat', 'J2', 'bond', 2)
mnf2.addmatrix('label', 'D', 'value', diag([0 0 D]), 'color', 'black');
mnf2.addaniso('D')
mnf2.genmagstr('mode', 'direct', 'S', [0 0; 0 0; 1 -1])

% Resolution (from PyChop)
Ei = 20;
eres = @(en) 512.17*sqrt((Ei-en)^3 * ( 8.26326e-10*(0.169+0.4*(Ei/(Ei-en))^1.5)^2 + 2.81618e-11*(1.169+0.4*(Ei/(Ei-en))^1.5)^2));
% Q-resolution (parameters for MARI)
e2k = @(en) sqrt( en .* (2*1.67492728e-27*1.60217653e-22) )./1.05457168e-34./1e10;
L1 = 11.8;  % Moderator to Sample distance in m
L2 = 2.5;   % Sample to detector distance in m
ws = 0.05;  % Width of sample in m
wm = 0.12;  % Width of moderator in m
wd = 0.025; % Width of detector in m
ki = e2k(25); 
a1 = ws/L1; % Angular width of sample seen from moderator
a2 = wm/L1; % Angular width of moderator seen from sample
a3 = wd/L2; % Angular width of detector seen from sample
a4 = ws/L2; % Angular width of sample seen from detector
dQ = 2.35 * sqrt( (ki*a1)^2/12 + (ki*a2)^2/12 + (ki*a3)^2/12 + (ki*a4)^2/12 );

fit_s = struct();
fit_s.func = @(obj, p) matparser(obj, 'param', p, 'mat', {'J1', 'J2', 'D(3,3)'}, 'init', true);
fit_s.x0 = [J1 J2 D];
fit_s.xmin = [-1 0 -0.2];  % Force J1 to be ferromagnetic to agree with structure
fit_s.xmax = [0 1 0.2];    % Force J2 to be antiferromagnetic to agree with structure
fit_s.hermit = false;
fit_s.formfact = true;
fit_s.nRand = 1e3;
fit_s.fibo = true;
fit_s.optimizer = 'lm3';   % Either: 'lm3' (Levenberg-Marquardt), 'pso' (Particle-Swarm), 'simplex' (Nelder-Mead)
fit_s.maxiter = 100;
fit_s.conv_func = @(spec) sw_instrument(spec, 'dE', eres, 'dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei);
fit_s.imagChk = false;
fit_s.nQ = 10;
fit_s.plot = true;
fit_s.intensity_scale = 15;    % Use a negative number if you don't want this to be fitted, set to zero to have a start value automatically determined
fit_s.background = 'linear';   % Either 'none', 'flat', 'linear' (x*p(1)+p(2)) or 'quadratic' (x.^2*p(1)+x*p(2)+p(3))
fit_s.background_par = [];     % Give either zero, 1, 2, or 3 numbers. If background is not 'none' and this is empty, the start params are auto-guessed 

fit_out = mnf2.fitpow(cts, fit_s)

It should take 1-2min to run the fit.

Note that for some reason, lm and lm2 (there are 3 implementations of the Levenberg-Marquardt algorithm in ndbase) doesn't work but lm3 does work.

[RichardWaiteSTFC]

What would happen is the cut thickness is << d|Q| resolution (I'm aware that the resolution might not be the same for all Q-vectors of equal |Q| - e.g. depends on scattering angle)?

[RichardWaiteSTFC]

Maybe for each Q point in average, you also simulate a number of points inside the resolution ellipsoid? Would that work? or have I mis-understood how the MC works?

[RichardWaiteSTFC]

Or do you do the resolution convolution on the simulated data before the powder averaging?

[RichardWaiteSTFC]

Discussed - @mducle to work on mex then @RichardWaiteSTFC to review and merge.
In a separate PR future plan is to:

1. Support 2D slices (i.e. IX_dataset_2d)
2. Refactor fitpow to be a separate class not method of spinw (and put in swfiles) - input struct will be setting attributes in that class
3. Move energy resolution to do before energy binning of modes (note Q-resolution OK, will ask user to supply dQ)

[mducle]

Closed in favour of #174