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