-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathintegral_curvature_sine.m
67 lines (46 loc) · 1.49 KB
/
integral_curvature_sine.m
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
function K_space = integral_curvature_sine(x_max)
xxx = 1:x_max;
disp(['Resolution is: ' num2str(1/x_max)])
%%
%Intializing.
min_amp = 3;
max_amp = 20;
iter_amp = 1;
min_freq = 1;
max_freq = 10;
iter_freq = 1;
%%
% Defining necessary functions.
function curv_k = int_sine_curv(F,A,x_i)
nominator = 4*pi*pi*F*F*A*sin(2*pi*x_i*F);
denominator = (1 + 2*pi*F*A*cos(2*pi*x_i*F))^(3/2);
curv_k = abs(nominator/denominator);
end
%%
% Readying summation of integral.
amps = min_amp:iter_amp:max_amp;
freqs = min_freq:iter_freq:max_freq;
length_x = length(xxx);
length_amps = length(amps);
length_freqs = length(freqs);
K_space = zeros(length_amps,length_freqs);
for i=1:length_amps
this_amp = amps(i);
for j=1:length_freqs
this_freq = freqs(j);
for k=1:length_x
x_k = xxx(i)/x_max;
K_ijk = int_sine_curv(this_freq,this_amp,x_k);
K_space(i,j) = K_space(i,j) + K_ijk;
end
end
end
disp('Row axis is amplitude.')
disp('Coloumn axis is frequency.')
figure(98)
bar3(K_space)
axis([0 Inf 0 Inf 1 max(K_space(:))])
title('Sine [Numeric integral] curvature per amplitude and frequency.')
xlabel('Frequencies')
ylabel('Amplitudes')
end