-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGetPerigee.asv
110 lines (99 loc) · 3.02 KB
/
GetPerigee.asv
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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
function Rho = GetPerigee(E, h, r0)
%GetPerigee (P,h) Íàõîæäåíèÿ íà÷àëüíûõ óñëîâèé äëÿ ðåøåíèÿ óðàâíåíèÿ äâèæåíèÿ.
% Âõîäíûå ïàðàìåòðû:
% E - çíà÷åíèå èíòåãðàëà ýíåðãèè,
% h - çíà÷åíèå èíòåãðàëà ïëîùàäåé.
% r0 - çíà÷åíèå ïàðàìåòðà ìîäåëè ïîòåíöèàëà
% Âûõîäíûå ïàðàìåòðû:
% Rho - hfcc,
%'r^6*(4*E^2) + r^4*(4*E^2*r0^2 - 4*r0^2 - 4*E*h^2) + r^2*(h^4 - 4*E*h^2*r0^2) + h^4*r0^2'
p = [4*E^2, 4*E^2*r0^2 - 4*r0^2 - 4*E*h^2, h^4 - 4*E*h^2*r0^2, h^4*r0^2];
%p = [4/r0^2, 4*E - 4*r0^2 - 4*(E^(1/2)/r0^(1/2))*h^2, h^4 - 4*E^(1/2)*h^2, h^4*r0^2];
rp = roots(p);
[i tt] = size(rp);
if (i > 0)
r = [];
for j = 1:i
if (rp(j) > 0)
a = sqrt(rp(j));
r = [r a];
end;
end;
end;
% r(1) = sqrt(r(1));
% r(2) = sqrt(r(2));
% r(3) = sqrt(r(3));
% fprintf('Çíà÷åíèÿ êîðíåé, ïîëó÷åííûå ìåòîäîì roots():\n');
% fprintf('\t%u -> %u\n', r(1), polyval(p, r(1)));
% fprintf('\t%u -> %u\n', r(2), polyval(p, r(2)));
% fprintf('\t%u -> %u\n', r(3), polyval(p, r(3)));
% % [r1,r2,r3] = getEquationZeros(E, h, r0);
% % fprintf('Çíà÷åíèÿ êîðíåé, ïîëó÷åííûå ìåòîäîì Âèåòà-Êàðäàíî:\n');
% % fprintf('\t%u -> %u\n', r1, polyval(p, r1));
% % fprintf('\t%u -> %u\n', r2, polyval(p, r2));
% % fprintf('\t%u -> %u\n', r3, polyval(p, r3));
% Rho = getMinPositiveValue(r(1), r(2), r(3));
Rho = min(r);
% function [z1, z2, z3] = getEquationZeros (E, h, r0)
% a = (4*E^2*r0^2 - 4*r0^2 - 4*E*h^2)/(4*E^2);
% b = (h^4 - 4*E*h^2*r0^2)/(4*E^2);
% c = (h^4*r0^2)/(4*E^2);
%
% [t1, t2, t3] = doCubic(a, b, c);
%
% %z1 = sqrt(t1);
% %z2 = sqrt(t2);
% %z3 = sqrt(t3);
% z1 = t1;
% z2 = t2;
% z3 = t3;
% end
%
% function [x1, x2, x3] = doCubic (a, b, c) %Âû÷èñëåíèå êîðíåé êóáè÷åñêîãî óðàâíåíèÿ ìåòîäîì Âèåòà-Êàðäàíî
% Q = (a^2-3*b)/9;
% R = (2*a^3 - 9*a*b + 27*c)/54;
%
% R2 = R*R;
% Q3 = Q*Q*Q;
% if (R2 >= Q3)
% t = acos(R/sqrt(Q3))/3;
% x1 = -2*sqrt(Q)*cos(t) - a/3;
% x2 = -2*sqrt(Q)*cos(t + (2*pi/3)) - a/3;
% x3 = -2*sqrt(Q)*cos(t - (2*pi/3)) - a/3;
% else
% Aa = -1 * sign(R) * ((abs(R) + sqrt(R2 - Q3)) ^ (1/3));
% if (Aa ~= 0)
% Bb = Q/Aa;
% else
% Bb = 0;
% end
%
% x1 = (Aa + Bb) - a/3;
%
% if (Aa == Bb)
% x2 = -Aa - a/3;
% x3 = x2;
% else
% x2 = -1; %Ò.ê. çíà÷åíèå êîìïëåêñíûõ êîðíåé íàñ íå èíòåðåñóþò, òî ìû âîçâðàùàåì -1.
% x3 = -1;
% end
% end
% end
%
% function result = getMinPositiveValue (v1, v2, v3)
%
% if (v1 <= 0)
% v1 = max ([v1, v2, v3]);
% end
%
% if (v2 <= 0)
% v2 = max ([v1, v2, v3]);
% end
%
% if (v3 <= 0)
% v3 = max ([v1, v2, v3]);
% end
%
% result = min ([v1, v2, v3]);
% end
end