tests/lowalgos/extcs/circ_sph_nointer

   1 puts "======================="
   2 puts "Test for Circle/Sphere extrema algorithm"
   3 puts "No intersection cases - one minimum solution should be found"
   4 puts "======================="
   5 puts ""
   6
   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
  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
  18 set nbstep 5
  19 # Rotation angle
  20 set angle [expr 180. / $nbstep]
  21
  22 # Set the number of Inner/Outer spheres in one direction
  23 set nbpairs 1
  24 # Set the delta for the radius of inner circle
  25 set delta_radius [expr $sph_radius * 0.9 / (2 * $nbpairs)]
  26
  27 # Step for sampling of the circle
  28 set dt [expr [dval 2*pi] / $nbstep]
  29
  30 # Iteration step
  31 set iStep 1
  32
  33 for {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 }