invert_resolution/vel_distribution0.m at master · abuts/invert_resolution · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
function [dv,f_d,v_peaks,dV_scat_max] = vel_distribution0(dvs)
% Calculate sample probabilityh distribution for velocity transfer
% for modelling and recovering
%
%
if ~isnumeric(dvs)
    dv = '5peaks';
    f_d = [];
    v_peaks = [];
    dV_scat_max = [];
    return
end

e_max = 25;

e_transf_const = 5.22725e-6; % sec^2/m^2
dV_scat_max = sqrt(e_max/e_transf_const);

if numel(dvs) ==  1 % step provided
    [dv,dV,NpDV] = adjust_step(-dV_scat_max,dV_scat_max,dvs);
else   % velocity points are provided
    dv = dvs;
    dV = dvs(2) - dvs(1);
end


% energy peaks
e_exc = [0,0.3,2,10,15];
sigma = [10,60,5,80,160];
Ampl  = [5,9,1,5,5];

e_transf_const = 5.22725e-6; % sec^2/m^2
v_peaks  = sqrt(e_exc/e_transf_const);

sig22 = 2*sigma.*sigma;
par = [Ampl',v_peaks',sig22'];

f_d = zeros(size(dv));
for i=1:numel(e_exc)
    f_d = f_d +par(i,1).*(exp(-(dv-par(i,2)).^2./par(i,3))+exp(-(dv+par(i,2)).^2./par(i,3)))./sqrt(pi*par(i,3));
end


%plot(v.*v*(V_char*V_char*e_transf_const),f_d);