-
Notifications
You must be signed in to change notification settings - Fork 2
/
ellipse_metric.m
executable file
·118 lines (87 loc) · 4.59 KB
/
ellipse_metric.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
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
111
112
113
114
115
116
function [g_sqrt, g_down11, g_down22, g_down12] = ellipse_metric(ellipse_parameters,n_points,uni)
R1 = ellipse_parameters(1);
R2 = ellipse_parameters(2);
P1 = ellipse_parameters(3);
P2 = ellipse_parameters(4);
Q1 = ellipse_parameters(5);
Q2 = ellipse_parameters(6);
N_intervals_x = 3;
N_intervals_y = 3;
x1_t_plus = P1/2+Q1;
x1_t_minus = P1/2-Q1;
x2_t_plus = P2/2+Q2;
x2_t_minus = P2/2-Q2;
b_x1 = [0 x1_t_minus x1_t_plus P1];
b_x2 = [0 x2_t_minus x2_t_plus P2];
ksi1 = linspace(-uni,uni,n_points);
ksi2 = linspace(-uni,uni,n_points);
x1 = zeros(N_intervals_x,n_points);
x2 = zeros(N_intervals_y,n_points);
for k1 = 1:N_intervals_x
x1(k1,:) = ( (b_x1(k1+1)-b_x1(k1))*ksi1+(b_x1(k1+1)+b_x1(k1)) )/2;
end
for k2 = 1:N_intervals_y
x2(k2,:) = ( (b_x2(k2+1)-b_x2(k2))*ksi2+(b_x2(k2+1)+b_x2(k2)) )/2;
end
g_sqrt = zeros(N_intervals_x,N_intervals_y,n_points,n_points);
g_sqrt(1,1,:,:) = ones(n_points, n_points);
g_sqrt(1,3,:,:) = ones(n_points, n_points);
g_sqrt(3,1,:,:) = ones(n_points, n_points);
g_sqrt(3,3,:,:) = ones(n_points, n_points);
g_down11 = zeros(N_intervals_x,N_intervals_y,n_points,n_points);
g_down11(1,1,:,:) = ones(n_points, n_points);
g_down11(1,3,:,:) = ones(n_points, n_points);
g_down11(3,1,:,:) = ones(n_points, n_points);
g_down11(3,3,:,:) = ones(n_points, n_points);
g_down22 = zeros(N_intervals_x,N_intervals_y,n_points,n_points);
g_down22(1,1,:,:) = ones(n_points, n_points);
g_down22(1,3,:,:) = ones(n_points, n_points);
g_down22(3,1,:,:) = ones(n_points, n_points);
g_down22(3,3,:,:) = ones(n_points, n_points);
g_down12 = zeros(N_intervals_x,N_intervals_y,n_points,n_points);
g_down12(1,1,:,:) = zeros(n_points, n_points);
g_down12(1,3,:,:) = zeros(n_points, n_points);
g_down12(3,1,:,:) = zeros(n_points, n_points);
g_down12(3,3,:,:) = zeros(n_points, n_points);
for nx = 1:n_points
for ny=1:n_points
g_sqrt(1,2,nx,ny) = (P1/2-(R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5)/x1_t_minus;
g_sqrt(2,1,nx,ny) = (P2/2-(R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5)/x2_t_minus;
g_sqrt(2,2,nx,ny) = ( 1/(x1_t_plus-x1_t_minus)*(x2_t_plus-x2_t_minus) )*...
(4*(R1^2-(x1(2,nx)-P1/2)^2)^0.5 * (R2^2-(x2(2,ny)-P2/2)^2)^0.5 - ...
(2*x1(2,nx)-x1_t_plus-x1_t_minus)*(2*x2(2,ny)-x2_t_plus-x2_t_minus)*...
(x1(2,nx)-P1/2)*(x2(2,ny)-P2/2)/( (R1^2-(x1(2,nx)-P1/2)^2)^0.5*(R2^2-(x2(2,ny)-P2/2)^2)^0.5 ));
g_sqrt(2,3,nx,ny) = ( (R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5 - P2/2 )/(x2_t_plus - P2);
g_sqrt(3,2,nx,ny) = ( (R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5 - P1/2 )/(x1_t_plus - P1);
g_down11(1,2,nx,ny) = (P1/2 - (R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5)^2/x1_t_minus^2;
g_down11(2,1,nx,ny) = 1 + (( (R2/R1)*(x1(2,nx)-P1/2))^2/(R1^2-(x1(2,nx)-P1/2)^2))*...
(x2(1,ny)/x2_t_minus)^2;
g_down11(2,2,nx,ny) = (1/(x1_t_plus-x1_t_minus)^2)*4*(R1/R2)^2*(R2^2-(x2(2,ny)-P2/2)^2) + ...
( ((2*x2(2,ny)-x2_t_plus-x2_t_minus)/(x2_t_plus-x2_t_minus))*(R2/R1)*(x1(2,nx)-P1/2) )^2/...
(R1^2-(x1(2,nx)-P1/2)^2);
g_down11(2,3,nx,ny) = 1 + ( (R2/R1)*((x2(3,ny)-P2)/(x2_t_plus-P2))*(x1(2,nx)-P1/2) )^2/...
(R1^2-(x1(2,nx)-P1/2)^2);
g_down11(3,2,nx,ny) = ( ((R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5 - P1/2)/(x1_t_plus-P1) )^2;
g_down22(1,2,nx,ny) = ( R1/R2*(x2(2,ny)-P2/2)*x1(1,nx)/x1_t_minus )^2/(R2^2-(x2(2,ny)-P2/2)^2)+1;
g_down22(2,1,nx,ny) = (P2/2 - (R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5)^2/x2_t_minus^2;
g_down22(2,2,nx,ny) = ((2*x1(2,nx)-x1_t_plus-x1_t_minus)/(x1_t_plus-x1_t_minus))^2*...
((R1/R2)*(x2(2,ny)-P2/2))^2/(R2^2-(x2(2,ny)-P2/2)^2)+...
4*(R2/R1)^2*(R1^2-(x1(2,nx)-P1/2)^2)/(x2_t_plus-x2_t_minus)^2;
g_down22(2,3,nx,ny) = ((R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5 - P2/2)^2/(x2_t_plus-P2)^2;
g_down22(3,2,nx,ny) = ( ((x1(3,nx)-P1)/(x1_t_plus-P1))*((R1/R2)*(x2(2,ny)-P2/2)) )^2/...
(R2^2-(x2(2,ny)-P2/2)^2)+1;
g_down12(1,2,nx,ny) = (P1/2 - (R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5)/x1_t_minus *...
((R1/R2)*(x2(2,ny)-P2/2)/(R2^2-(x2(2,ny)-P2/2)^2)^0.5)*x1(1,nx)/x1_t_minus;
g_down12(2,1,nx,ny) = (x2(1,ny)/x2_t_minus)*(R2/R1)*((x1(2,nx)-P1/2)/(R1^2-(x1(2,nx)-P1/2)^2)^0.5)*...
(P2/2-(R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5)/x2_t_minus;
g_down12(2,2,nx,ny) = (2*(R1/R2)^2/(x1_t_plus-x1_t_minus)^2) * ...
(2*x1(2,nx)-x1_t_plus-x1_t_minus)*(-(x2(2,ny)-P2/2)) + ...
((2*x2(2,ny)-x2_t_plus-x2_t_minus)/(x2_t_plus-x2_t_minus)^2) * ...
2*(R2/R1)^2*(-(x1(2,nx)-P1/2));
g_down12(2,3,nx,ny) = (((x2(3,ny)-P2)/(x2_t_plus-P2)^2)*(R2/R1)*(-(x1(2,nx)-P1/2))/...
(R1^2-(x1(2,nx)-P1/2)^2)^0.5)*((R2/R1)*(R1^2-(x1(2,nx)-P1/2)^2)^0.5-P2/2);
g_down12(3,2,nx,ny) = (((R1/R2)*(R2^2-(x2(2,ny)-P2/2)^2)^0.5 - P1/2)*(x1(3,nx)-P1)/...
(x1_t_plus-P1)^2)*(R1/R2)*(-(x2(2,ny)-P2/2))/(R2^2-(x2(2,ny)-P2/2)^2)^0.5;
end
end
g_sqrt = abs(g_sqrt);