invert_resolution/I_mk2.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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
function int = I_mk2(k_v,n_t,v_min,v_max,v_index,t_index)
% v_min == v_min[cm/sec]/DV[cm/sec]
% v_max == v_min[cm/sec]/DV[cm/sec]
kv = v_index(k_v);
nt = t_index(n_t);

if kv~=0
    if nt==0
        int = 0;
        %int = (exp(1i*pi*kv)-exp(-1i*pi*kv))/(2i*pi*kv);
        %         fun0 = @(v)(exp(2i*pi*(kv*v)));
        %         int0 = integral(fun0,-0.5,0.5);

    else
        %         fp = 0.5*kv-2*nt*v_min;
        %         fun = @(u)(exp(2i*pi*u).*u./sqrt(u.*u+4*nt*kv*v_min));
        %         int = sin(2*pi*fp)/(2*pi) + integral(fun,-fp,fp)/kv;

        funD = @(v)(cos(2*pi*(kv*v-v_min*nt./v)));
%         funI = @(u)(cos(2*pi*(kv./u-v_min*nt*u))./(u.*u));
%         int1 = 2*integral(funI,2,Inf);
%
%
%         int0 = integral(funD,-0.5,0.5);
        AbsTol = 1.e-10;

        delta = 1.e-3/abs(2*pi*kv);
        u0 = 1/delta;
        if delta<0.5
            fd = @(v)(kv*v-v_min*nt./v);
            np = root_range(fd,delta,0.5);
            if ~isempty(np)
                det = sqrt(0.25*np.^2/(kv*kv)+nt*v_min/kv);
                shi = np*0.25/kv;
                roots = [shi+det,shi-det];
                valid = roots>delta&roots<0.5;
                wp = roots(valid);
                int20 =integral(funD,delta,0.5,'Waypoints',wp,'RelTol',0,'AbsTol',AbsTol);
            else
                int20 = integral(funD,delta,0.5,'RelTol',0,'AbsTol',AbsTol);
            end


            funIS = @(u)(cos(2*pi*v_min*nt*u)./(u).^4);
            max_lim = (1/AbsTol)^(1/3);
            if u0<max_lim
                intC = integral(funIS,u0,max_lim);
            else
                intC  = 0;
            end
            betta = v_min*nt;
            alpha = kv;
            albet =alpha/betta;
            int22 = 2*albet*cos(2*pi*betta/delta)*delta^3+(4*albet/3-(2*pi*alpha)^2)*intC;
            arg = 2*pi*nt*v_min;
            int21 = 2*delta*cos(arg/delta)+1i*arg*(expint(-1i*arg/delta)-expint(1i*arg/delta));
            int   = 2*int20+int22+int21;
        else
            funI = @(u)(cos(2*pi*(kv./u-v_min*nt*u))./(u.*u));
            int = 2*integral(funI,2,Inf);
        end

    end
else % kv==0
    if nt == 0
        int = 1;
    else
        arg = 2*pi*nt*v_min;
        int = cos(2*arg)+1i*arg*(expint(-2i*arg)-expint(2i*arg));

        %
        %         fun0 = @(v)(exp(-2i*pi*v_min*nt./v));
        %         int0 = integral(fun0,-0.5,0.5);

    end
end

function range = root_range(fun,min_lim,max_lim)

min_val = min(fun(min_lim),fun(max_lim));
max_val = max(fun(min_lim),fun(max_lim));

min_val_n = round(min_val);
if min_val_n<min_val
    min_val_n = min_val_n+1;
end
max_val_n = round(max_val);
if max_val_n>max_val
    max_val_n = max_val_n-1;
end
if min_val_n<=max_val_n
    range = min_val_n:1:max_val_n;
else
    range = [];
end