1 puts "======================="
2 puts "Test for Circle/Sphere extrema algorithm"
3 puts "Intersection case (two intersection points should be found"
4 puts "======================="
5 puts ""
7 # Make sphere
8 set x0 0.
9 set y0 0.
10 set z0 0.
11 set sph_radius 10.
12 sphere s \$x0 \$y0 \$z0 \$sph_radius
14 # All circles will be to connect the two points - one inside the sphere, another - outside.
15 # Such circle will definitely intersect the initial sphere in two points.
16 # The points to make the circle will taken two spheres - with smaller and bigger radius.
18 # Set the number of iterations (number of pairs of spheres)
19 set nbpairs 2
21 # Set ratio to increase/decrease the radius if additional spheres
22 set ratio_radius 2.
24 # Number of sampling points on the spheres
25 set nbsamples 5
27 # Number of circles rotations
28 set nbsteps 5
29 set angle [expr 180. / \$nbsteps]
31 # Iteration step
32 set iStep 1
34 for {set i 1} {\$i <= \$nbpairs} {incr i} {
35   # Define the radius for spheres
36   set s_in_radius [expr \$sph_radius / (\$i * \$ratio_radius)]
37   set s_out_radius [expr \$sph_radius * (\$i * \$ratio_radius)]
39   # Make the spheres
40   sphere s_in \$x0 \$y0 \$z0 \$s_in_radius
41   sphere s_out \$x0 \$y0 \$z0 \$s_out_radius
43   # Make the sampling of the spheres
45   # spheres are the same, so there is no difference from which one to take the parameters
46   bounds s_in umin umax vmin vmax
47   set du [dval (umax-umin)/\$nbsteps]
48   set dv [dval (vmax-vmin)/\$nbsteps]
50   for {set u1 1} {\$u1 <= \$nbsamples} {incr u1} {
51     for {set v1 1} {\$v1 <= \$nbsamples} {incr v1} {
53       # point on inner sphere
54       svalue s_in [dval umin+\$u1*\$du] [dval vmin+\$v1*\$dv] x1 y1 z1
56       for {set u2 1} {\$u2 <= \$nbsamples} {incr u2} {
57         for {set v2 1} {\$v2 <= \$nbsamples} {incr v2} {
59           # point on outer sphere
60           svalue s_out [dval umin+\$u2*\$du] [dval vmin+\$v2*\$dv] x2 y2 z2
62           # rotation direction
63           set dx [dval x2-x1]
64           set dy [dval y2-y1]
65           set dz [dval z2-z1]
67           # center of the circle
68           set xc [dval (x1+x2)/2.]
69           set yc [dval (y1+y2)/2.]
70           set zc [dval (z1+z2)/2.]
72           # compute radius for circle
73           set cir_radius [expr sqrt(\$dx*\$dx + \$dy*\$dy + \$dz*\$dz) / 2.]
75           # make plane to get its XAxis
76           plane p \$xc \$yc \$zc \$dx \$dy \$dz
78           regexp {XAxis  :([-0-9.+eE]*), ([-0-9.+eE]*), ([-0-9.+eE]*)} [dump p] full dxx dxy dxz
80           circle c \$xc \$yc \$zc \$dxx \$dxy \$dxz \$cir_radius
82           # make rotation
83           copy c c_rotated
84           for {set j 1} {\$j <= \$nbsteps} {incr j} {
85             rotate c_rotated \$xc \$yc \$zc \$dx \$dy \$dz \$angle
87             set log [extrema c_rotated s]
89             # save each circle if necessary
90             # copy c_rotated c_\$iStep
92             set isfound1 0
93             set isfound2 0
95             if {[regexp "ext_1" \$log]} {
96               set ext_dist [lindex [length ext_1] end]
97               checkreal "Step \$iStep, min distance 1" \$ext_dist 0 1.e-7 1.e-7
98               set isfound1 1
99             }
101             if {[regexp "ext_2" \$log]} {
102               set ext_dist [lindex [length ext_2] end]
103               checkreal "Step \$iStep, min distance 2" \$ext_dist 0 1.e-7 1.e-7
104               set isfound2 1
105             }
107             if {[regexp "Extrema 1 is point" \$log]} {
108               puts "Check of Step \$iStep, min distance 1  OK"
109               set isfound1 1
110             }
112             if {[regexp "Extrema 2 is point" \$log]} {
113               puts "Check of Step \$iStep, min distance 2  OK"
114               set isfound2 1
115             }
117             if {\$isfound1 == 0 || \$isfound2 == 0} {
118               puts "Error: Extrema has not detected the intersection case on step \$iStep"
119             }
121             incr iStep
122           }
123         }
124       }
125     }
126   }
127 }