1puts "======================="
2puts "Test for Circle/Sphere extrema algorithm"
3puts "No intersection cases - one minimum solution should be found"
4puts "======================="
5puts ""
6
7# Make sphere
8set x0 0.
9set y0 0.
10set z0 0.
11set sph_radius 10.
12sphere s $x0 $y0 $z0 $sph_radius
13
14# The circles will be made on the distance from the surface
15# as intersection of pairs of inner and outer spheres with the plane
16
17# Set the number of iterations
18set nbstep 5
19# Rotation angle
20set angle [expr 180. / $nbstep]
21
22# Set the number of Inner/Outer spheres in one direction
23set nbpairs 1
24# Set the delta for the radius of inner circle
25set delta_radius [expr $sph_radius * 0.9 / (2 * $nbpairs)]
26
27# Step for sampling of the circle
28set dt [expr [dval 2*pi] / $nbstep]
29
30# Iteration step
31set iStep 1
32
33for {set i 1} {$i <= $nbpairs} {incr i} {
34  # Define the inner circle
35  set circ_radius [expr $i * $delta_radius]
36  circle c $x0 $y0 $z0 0 0 1 $circ_radius
37
38  set diff [expr $sph_radius - $circ_radius]
39
40  # Distance between inner sphere on circle and initial sphere
41  set real_dist [expr $sph_radius - 2*$circ_radius]
42
43  # Circle will be rotated around the line
44  line rotation_line $x0 $y0 $z0 1 0 0
45
46  # Line rotation
47  for {set j 1} {$j <= $nbstep} {incr j} {
48    rotate rotation_line $x0 $y0 $z0 0 0 1 $angle
49
50    # Get direction for circle's rotation
51    regexp {Axis   :([-0-9.+eE]*), ([-0-9.+eE]*), ([-0-9.+eE]*)} [dump rotation_line] full dx dy dz
52
53    # Circle rotation
54    copy c c_rotated
55    for {set k 1} {$k <= $nbstep} {incr k} {
56      rotate c_rotated 0 0 0 $dx $dy $dz $angle
57
58      # Sampling of the circle
59      for {set n 1} {$n <= $nbstep} {incr n} {
60        cvalue c_rotated $n*$dt x1 y1 z1
61
62        set x1 [dval x1]
63        set y1 [dval y1]
64        set z1 [dval z1]
65
66        # Normalize the vector
67        set dtx [expr ($x1 - $x0) / $circ_radius]
68        set dty [expr ($y1 - $y0) / $circ_radius]
69        set dtz [expr ($z1 - $z0) / $circ_radius]
70
71        # Create inner and outer spheres
72        set iC 1
73
74        repeat 2 {
75          sphere s_to_int $x1 $y1 $z1 $circ_radius
76
77          # Define the point closest to the initial sphere
78          set x_sol [expr $x1 + $iC * $circ_radius * $dtx]
79          set y_sol [expr $y1 + $iC * $circ_radius * $dty]
80          set z_sol [expr $z1 + $iC * $circ_radius * $dtz]
81
82
83          # Intersect the sphere with the plane originated in closes point
84
85          # Make the sampling of the sphere to define section plane's direction
86
87          bounds s_to_int umin umax vmin vmax
88
89          set du [dval (umax-umin)/$nbstep]
90          set dv [dval (vmax-vmin)/$nbstep]
91
92          for {set u 1} {$u <= $nbstep} {incr u} {
93            for {set v 1} {$v <= $nbstep} {incr v} {
94
95              # Get point on surface
96              svalue s_to_int [dval umin+$u*$du] [dval vmin+$v*$dv] xs ys zs
97
98              # Check that it is not the same point
99              set sqdist [dval (xs-$x_sol)*(xs-$x_sol)+(ys-$y_sol)*(ys-$y_sol)+(zs-$z_sol)*(zs-$z_sol)]
100              if {$sqdist < 1.e-16} {
101                # Skip the sampling point
102                continue;
103              }
104
105              # Create the intersection plane
106              plane p_int $x_sol $y_sol $z_sol [dval xs-$x_sol] [dval ys-$y_sol] [dval zs-$z_sol]
107              # Intersect the sphere by plane to obtain the circle
108              foreach c_int [intersect c_inter s_to_int p_int] {
109
110                # Check if the circle contains the point
111                if {![regexp "Point on curve" [proj $c_int $x_sol $y_sol $z_sol]]} {
112                  if {[lindex [length ext_1] end] >= 1.e-7} {
113                    # run extrema - one of the ends of the curve should be the solution
114                    set log [extrema $c_int s 1]
115                    if {[regexp "prm_1_1" $log]} {
116                      # get parameters of the curve
117                      bounds $c_int fp lp
118                      if {[dval prm_1_1-fp] > 1.e-7 && [dval lp-prm_1_1] > 1.e-7} {
119                        puts "Error: Extrema has failed to find the minimal distance on step $iStep"
120                      }
121                    } else {
122                      puts "Error: Extrema has failed to find the minimal distance on step $iStep"
123                    }
124
125                    # save each circle if necessary
126                    # copy $c_int c_$iStep
127
128                    incr iStep
129                    continue
130                  }
131                }
132
133                # Make extrema computation
134                set log [extrema $c_int s]
135
136                # save each circle if necessary
137                # copy $c_int c_$iStep
138
139                if {![regexp "ext_1" $log]} {
140                  puts "Error: Extrema has failed to find the minimal distance on step $iStep"
141                } else {
142                  set ext_dist [lindex [length ext_1] end]
143                  checkreal "Step $iStep, min distance " $ext_dist $real_dist 1.e-7 1.e-7
144                }
145                incr iStep
146              }
147            }
148          }
149
150          # prepare for the outer sphere
151          set x1 [expr $x1 + 2 * $diff * $dtx]
152          set y1 [expr $y1 + 2 * $diff * $dty]
153          set z1 [expr $z1 + 2 * $diff * $dtz]
154
155          set iC -1
156        }
157      }
158    }
159  }
160}
161