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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
puts "======================="
puts "Test for Circle/Sphere extrema algorithm"
puts "No intersection cases - one minimum solution should be found"
puts "======================="
puts ""
# Make sphere
set x0 0.
set y0 0.
set z0 0.
set sph_radius 10.
sphere s $x0 $y0 $z0 $sph_radius
# The circles will be made on the distance from the surface
# as intersection of pairs of inner and outer spheres with the plane
# Set the number of iterations
set nbstep 5
# Rotation angle
set angle [expr 180. / $nbstep]
# Set the number of Inner/Outer spheres in one direction
set nbpairs 1
# Set the delta for the radius of inner circle
set delta_radius [expr $sph_radius * 0.9 / (2 * $nbpairs)]
# Step for sampling of the circle
set dt [expr [dval 2*pi] / $nbstep]
# Iteration step
set iStep 1
for {set i 1} {$i <= $nbpairs} {incr i} {
# Define the inner circle
set circ_radius [expr $i * $delta_radius]
circle c $x0 $y0 $z0 0 0 1 $circ_radius
set diff [expr $sph_radius - $circ_radius]
# Distance between inner sphere on circle and initial sphere
set real_dist [expr $sph_radius - 2*$circ_radius]
# Circle will be rotated around the line
line rotation_line $x0 $y0 $z0 1 0 0
# Line rotation
for {set j 1} {$j <= $nbstep} {incr j} {
rotate rotation_line $x0 $y0 $z0 0 0 1 $angle
# Get direction for circle's rotation
regexp {Axis :([-0-9.+eE]*), ([-0-9.+eE]*), ([-0-9.+eE]*)} [dump rotation_line] full dx dy dz
# Circle rotation
copy c c_rotated
for {set k 1} {$k <= $nbstep} {incr k} {
rotate c_rotated 0 0 0 $dx $dy $dz $angle
# Sampling of the circle
for {set n 1} {$n <= $nbstep} {incr n} {
cvalue c_rotated $n*$dt x1 y1 z1
set x1 [dval x1]
set y1 [dval y1]
set z1 [dval z1]
# Normalize the vector
set dtx [expr ($x1 - $x0) / $circ_radius]
set dty [expr ($y1 - $y0) / $circ_radius]
set dtz [expr ($z1 - $z0) / $circ_radius]
# Create inner and outer spheres
set iC 1
repeat 2 {
sphere s_to_int $x1 $y1 $z1 $circ_radius
# Define the point closest to the initial sphere
set x_sol [expr $x1 + $iC * $circ_radius * $dtx]
set y_sol [expr $y1 + $iC * $circ_radius * $dty]
set z_sol [expr $z1 + $iC * $circ_radius * $dtz]
# Intersect the sphere with the plane originated in closes point
# Make the sampling of the sphere to define section plane's direction
bounds s_to_int umin umax vmin vmax
set du [dval (umax-umin)/$nbstep]
set dv [dval (vmax-vmin)/$nbstep]
for {set u 1} {$u <= $nbstep} {incr u} {
for {set v 1} {$v <= $nbstep} {incr v} {
# Get point on surface
svalue s_to_int [dval umin+$u*$du] [dval vmin+$v*$dv] xs ys zs
# Check that it is not the same point
set sqdist [dval (xs-$x_sol)*(xs-$x_sol)+(ys-$y_sol)*(ys-$y_sol)+(zs-$z_sol)*(zs-$z_sol)]
if {$sqdist < 1.e-16} {
# Skip the sampling point
continue;
}
# Create the intersection plane
plane p_int $x_sol $y_sol $z_sol [dval xs-$x_sol] [dval ys-$y_sol] [dval zs-$z_sol]
# Intersect the sphere by plane to obtain the circle
foreach c_int [intersect c_inter s_to_int p_int] {
# Check if the circle contains the point
if {![regexp "Point on curve" [proj $c_int $x_sol $y_sol $z_sol]]} {
if {[lindex [length ext_1] end] >= 1.e-7} {
# run extrema - one of the ends of the curve should be the solution
set log [extrema $c_int s 1]
if {[regexp "prm_1_1" $log]} {
# get parameters of the curve
bounds $c_int fp lp
if {[dval prm_1_1-fp] > 1.e-7 && [dval lp-prm_1_1] > 1.e-7} {
puts "Error: Extrema has failed to find the minimal distance on step $iStep"
}
} else {
puts "Error: Extrema has failed to find the minimal distance on step $iStep"
}
# save each circle if necessary
# copy $c_int c_$iStep
incr iStep
continue
}
}
# Make extrema computation
set log [extrema $c_int s]
# save each circle if necessary
# copy $c_int c_$iStep
if {![regexp "ext_1" $log]} {
puts "Error: Extrema has failed to find the minimal distance on step $iStep"
} else {
set ext_dist [lindex [length ext_1] end]
checkreal "Step $iStep, min distance " $ext_dist $real_dist 1.e-7 1.e-7
}
incr iStep
}
}
}
# prepare for the outer sphere
set x1 [expr $x1 + 2 * $diff * $dtx]
set y1 [expr $y1 + 2 * $diff * $dty]
set z1 [expr $z1 + 2 * $diff * $dtz]
set iC -1
}
}
}
}
}
|
