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
|
puts "======================="
puts "Test for Circle/Sphere extrema algorithm"
puts "Intersection case (two intersection points 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
# All circles will be to connect the two points - one inside the sphere, another - outside.
# Such circle will definitely intersect the initial sphere in two points.
# The points to make the circle will taken two spheres - with smaller and bigger radius.
# Set the number of iterations (number of pairs of spheres)
set nbpairs 2
# Set ratio to increase/decrease the radius if additional spheres
set ratio_radius 2.
# Number of sampling points on the spheres
set nbsamples 5
# Number of circles rotations
set nbsteps 5
set angle [expr 180. / $nbsteps]
# Iteration step
set iStep 1
for {set i 1} {$i <= $nbpairs} {incr i} {
# Define the radius for spheres
set s_in_radius [expr $sph_radius / ($i * $ratio_radius)]
set s_out_radius [expr $sph_radius * ($i * $ratio_radius)]
# Make the spheres
sphere s_in $x0 $y0 $z0 $s_in_radius
sphere s_out $x0 $y0 $z0 $s_out_radius
# Make the sampling of the spheres
# spheres are the same, so there is no difference from which one to take the parameters
bounds s_in umin umax vmin vmax
set du [dval (umax-umin)/$nbsteps]
set dv [dval (vmax-vmin)/$nbsteps]
for {set u1 1} {$u1 <= $nbsamples} {incr u1} {
for {set v1 1} {$v1 <= $nbsamples} {incr v1} {
# point on inner sphere
svalue s_in [dval umin+$u1*$du] [dval vmin+$v1*$dv] x1 y1 z1
for {set u2 1} {$u2 <= $nbsamples} {incr u2} {
for {set v2 1} {$v2 <= $nbsamples} {incr v2} {
# point on outer sphere
svalue s_out [dval umin+$u2*$du] [dval vmin+$v2*$dv] x2 y2 z2
# rotation direction
set dx [dval x2-x1]
set dy [dval y2-y1]
set dz [dval z2-z1]
# center of the circle
set xc [dval (x1+x2)/2.]
set yc [dval (y1+y2)/2.]
set zc [dval (z1+z2)/2.]
# compute radius for circle
set cir_radius [expr sqrt($dx*$dx + $dy*$dy + $dz*$dz) / 2.]
# make plane to get its XAxis
plane p $xc $yc $zc $dx $dy $dz
regexp {XAxis :([-0-9.+eE]*), ([-0-9.+eE]*), ([-0-9.+eE]*)} [dump p] full dxx dxy dxz
circle c $xc $yc $zc $dxx $dxy $dxz $cir_radius
# make rotation
copy c c_rotated
for {set j 1} {$j <= $nbsteps} {incr j} {
rotate c_rotated $xc $yc $zc $dx $dy $dz $angle
set log [extrema c_rotated s]
# save each circle if necessary
# copy c_rotated c_$iStep
set isfound1 0
set isfound2 0
if {[regexp "ext_1" $log]} {
set ext_dist [lindex [length ext_1] end]
checkreal "Step $iStep, min distance 1" $ext_dist 0 1.e-7 1.e-7
set isfound1 1
}
if {[regexp "ext_2" $log]} {
set ext_dist [lindex [length ext_2] end]
checkreal "Step $iStep, min distance 2" $ext_dist 0 1.e-7 1.e-7
set isfound2 1
}
if {[regexp "Extrema 1 is point" $log]} {
puts "Check of Step $iStep, min distance 1 OK"
set isfound1 1
}
if {[regexp "Extrema 2 is point" $log]} {
puts "Check of Step $iStep, min distance 2 OK"
set isfound2 1
}
if {$isfound1 == 0 || $isfound2 == 0} {
puts "Error: Extrema has not detected the intersection case on step $iStep"
}
incr iStep
}
}
}
}
}
}
|
