# Library functions for the scripting library (scrib)
#
# File: scrib/lib/math.tcl
#
# $Revision$ on $Date$
# 
# Copyright (c) 1993-2003 ICEM CFD Engineering
#
#--------------------------------------------------------------------
#
# --> if the name of the first proc already exists in interp, return
#
if {"" != [info commands scrib_lib_math_is_integer]} { return }
#
# --> procs in this file depend on the following scrib/lib/ files--source them
#     I guess we could do this by proc, too, not sure waht's best
#
if 0 {
    # none required yet
    set required_files ""
    scrib_lib_setup_source_files $required_files
}
#

# Most of this came from Michigan --jka
# Most it exists in med_batch/geometry_funcs.tcl
# in procs like ic_geo_vec_* --jka
#
#
# Math stuff---------------------------------------

proc scrib_lib_math_is_int {l} {
  if {[string trim $l] == ""} {return 0}
  foreach k $l {
    if ![regexp {^ *[+-]?[0-9]+ *$} $k] {return 0}
  }
  return 1
}

## is the value <em>val</em> an integer?

proc scrib_lib_math_is_integer {val} {
    # why not just use Tcl commands, like this? -jka
    return [string is integer -strict $val]
    #--------------- old-->
    if {$l == ""} { return 0 }
    foreach k $l {
	if ![regexp {^ *[+-]?[0-9]+ *$} $k] {return 0}
    }
    return 1
}

## is the value <em>val</em> a real value?

proc scrib_lib_math_is_real {val} {
    return [string is double -strict $val]
   
   
}
proc scrib_lib_math_is_real {l} {
    ## return value for 1 & 1.0 are 1
    ## return [string is double -strict $val]
    if {[string trim $l] == ""} {return 0}
      foreach k $l {
        if {![expr [regexp {^ *[+-]?[0-9]*(\.[0-9]*)+([eE][+-]?[0-9]+)? *$} $k] || \
                   [regexp {^ *[+-]?[0-9]*(\.[0-9]*)?([eE][+-]?[0-9]+)+ *$} $k]] || \
                   ($k == ".")} {return 0}
      }
      return 1
}
## is the value <em>val</em> an integer or real value?

proc scrib_lib_math_is_int_real {l} {
  if {[string trim $l] == ""} {return 0}
  foreach k $l {
    if {![scrib_lib_math_is_int $k] && ![scrib_lib_math_is_real $k]} {return 0}
  }
  return 1
}

## is the value <em>val</em> an integer or real value? jason's

proc scrib_lib_math_is_int_real_jka {val} {
    return [expr {[scrib_lib_math_is_integer $val] || [scrib_lib_math_is_real $val]}]
}

## returns the given value <em>val</em> as a real value.

proc scrib_lib_math_make_real {l} {
    
    set r ""
    foreach s $l {
	if [scrib_lib_math_is_integer $s] {set s $s.0}
	lappend r $s
    }
    return $r
}

## is the given string <em>val</em> valid coords?

proc scrib_lib_math_is_point {val} {
    
    if {[llength $val] != 3} { return 0 }
    #mess "DEBUG: scrib_lib_math_is_point val: $val\n"
    foreach c $val {
	if ![scrib_lib_math_is_int_real $c] { 
	    #mess "DEBUG: scrib_lib_math_is_point c: $c\n"	    
	    return 0 
	}
    }
    return 1
}

proc scrib_lib_math_std_unit_vector {dir} {
    
    switch -- $dir {
	"+x" - "+X" - "x+" - "X+" { set v "1 0 0" }
	"-x" - "-X" -" x-" - "X-" { set v "-1 0 0" }
	"+y" - "+Y" - "y+" - "Y+" { set v "0 1 0" }
	"-y" - "-Y" - "y-" - "Y-" { set v "0 -1 0" }
	"+z" - "+Z" - "z+" - "Z+" { set v "0 0 1" }
	"-z" - "-Z" - "z-" - "Z-" { set v "0 0 -1" }
	"NO" - "no"               { set v "0 0 0" }
	default {
	    mess "Warning: direction $dir is not supported for standard unit vector\n"  red
	    set v $dir
	}
    }
    return $v
}


## needs doc

proc scrib_lib_math_vector_from_points {p0 p1} {
    
    foreach p {p0 p1} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_vector_from_points: Wrong parameter $p = \{$point\} \n" red
	    return 0
	} else {
	    set $p $point
	}
    }
    
    foreach ax {0 1 2} {
	set v_$ax [expr [lindex $p1 $ax]-[lindex $p0 $ax]]
    }
    return "$v_0 $v_1 $v_2"
}

## needs doc

proc scrib_lib_math_center_point {args} {
    
    set list ""
    foreach p $args {
	set point [scrib_lib_math_make_real $p]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_center_point: Wrong parameter \{$point\} \n" red
	    return 0  
	} else {
	    lappend list $point
	}
    }
    
    set sum "0.0 0.0 0.0"
    foreach p $list {
	set sum [scrib_lib_math_vector_sum $sum $p]        
    }  
    return [scrib_lib_math_vector_mult $sum [expr 1.0/[llength $list]]]
}

## needs doc

proc scrib_lib_math_vector_sum {a b} {
    
    foreach p {a b} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_vector_sum: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set $p $point
	}
    }
    
    foreach ax {0 1 2} {
	set a_$ax [lindex $a $ax]
	set b_$ax [lindex $b $ax]
    }
    set v_0 [expr $a_0 + $b_0]
    set v_1 [expr $a_1 + $b_1]
    set v_2 [expr $a_2 + $b_2]
    return "$v_0 $v_1 $v_2"
}

## needs doc

proc scrib_lib_math_vector_vmult {a b} {
    
    foreach p {a b} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_vector_vmult: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set $p $point
	}
    }
    
    foreach ax {0 1 2} {
	set a_$ax [lindex $a $ax]
	set b_$ax [lindex $b $ax]
    }
    set v_0 [expr $a_1*$b_2 - $a_2*$b_1]
    set v_1 [expr $a_2*$b_0 - $a_0*$b_2]
    set v_2 [expr $a_0*$b_1 - $a_1*$b_0]
    return "$v_0 $v_1 $v_2"
}

## needs doc

proc scrib_lib_math_vector_smult {a b} {
    
    foreach p {a b} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_vector_smult: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set $p $point
	}
    }
    
    foreach ax {0 1 2} {
	set a_$ax [lindex $a $ax]
	set b_$ax [lindex $b $ax]
    }
    return [expr $a_0*$b_0 + $a_1*$b_1 + $a_2*$b_2]
}

## needs doc

proc scrib_lib_math_angle_between_vectors {a b} {
    
    #global rad
    
    foreach p {a b} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_angle_between_vectors: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set $p $point
	}
    }
    set ab [expr [scrib_lib_math_vector_length $a]*[scrib_lib_math_vector_length $b]]
    if {$ab == 0} {return 0}

    set arg [expr [scrib_lib_math_vector_smult $a $b]/ $ab ]
    if {$arg > 1} {set arg 1.0}
    if {$arg < -1} {set arg -1.0}
    return [expr acos($arg)]
}

## needs doc

proc scrib_lib_math_vector_mult {v k} {
    
    set v [scrib_lib_math_make_real $v]
    if ![scrib_lib_math_is_point $v] {
	mess "scrib_lib_math_vector_mult: Wrong parameter v = \{$v\} \n" red
	return 0  
    }
    if ![scrib_lib_math_is_int_real $k] {
	mess "scrib_lib_math_mult: Wrong parameter k = $k \n" red
	return 0  
    }
    
    set vx [lindex $v 0]
    set vy [lindex $v 1]
    set vz [lindex $v 2]
    return "[expr $k*$vx] [expr $k*$vy] [expr $k*$vz]"
}

## needs doc

proc scrib_lib_math_vector_length {v} {
    
    set v [scrib_lib_math_make_real $v]
    if ![scrib_lib_math_is_point $v] {
	mess "scrib_lib_math_vector_length: Wrong parameter \{$v\} \n" red
	return 0  
    }
    
    set vx [lindex $v 0]
    set vy [lindex $v 1]
    set vz [lindex $v 2]
    return [expr sqrt($vx*$vx + $vy*$vy + $vz*$vz)]
}

## needs doc

proc scrib_lib_math_distance {a b} {
    
    foreach p {a b} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_distance: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set $p $point
	}
    }
    
    foreach ax {0 1 2} {
	set d_$ax [expr [lindex $a $ax] - [lindex $b $ax]]
    }
    return [expr sqrt($d_0*$d_0 + $d_1*$d_1 + $d_2*$d_2)]
}

## needs doc

proc scrib_lib_math_distance_to_line {p1 p2 p3} {
    
    foreach p {p1 p2 p3} {
	set point [scrib_lib_math_make_real [subst $$p]]
	if ![scrib_lib_math_is_point $point] {
	    mess "scrib_lib_math_distance_to_line: Wrong parameter $p = \{$point\} \n" red
	    return 0  
	} else {
	    set ${p}_x [lindex $point 0]
	    set ${p}_y [lindex $point 1]
	    set ${p}_z [lindex $point 2]
	}
    }
    
    set ax [expr $p2_x - $p1_x]
    set ay [expr $p2_y - $p1_y]
    set az [expr $p2_z - $p1_z]
    set bx [expr $p3_x - $p2_x]
    set by [expr $p3_y - $p2_y]
    set bz [expr $p3_z - $p2_z]
    
    set d1 [expr sqrt((($ay*$bz)-($az*$by))*(($ay*$bz)-($az*$by))+(($ax*$bz)-($az*$bx))*(($ax*$bz)-($az*$bx))+(($ax*$by)-($ay*$bx))*(($ax*$by)-($ay*$bx)))]
    set d2 [expr sqrt($ax*$ax+$ay*$ay+$az*$az)]
    if !$d2 {
	mess "scrib_lib_math_distance_to_line: Wrong line definition: \{$p1\} \{$p2\} \n" red
	return 0  
    }
    
    return [expr $d1/$d2]
}

## needs doc

proc scrib_lib_math_is_point_on_segment {p1 p2 p3 tol} {
    
    set h [scrib_lib_math_distance_to_line $p1 $p2 $p3]
    set d1 [scrib_lib_math_distance $p1 $p3]
    set d2 [scrib_lib_math_distance $p2 $p3]
    set l [expr [scrib_lib_math_distance $p1 $p2]+$tol]
    if {($h <= $tol) && ($d1 <= $l) && ($d2 <= $l)} {
	return 1
    } else {
	return 0
    }
}

## needs doc

proc scrib_lib_math_vector_nrm {v} {
    
    set v [scrib_lib_math_make_real $v]
    if ![scrib_lib_math_is_point $v] {
	mess "scrib_lib_math_vector_nrm: Wrong parameter \{$v\} \n" red
	return 0  
    }
    
    set l [scrib_lib_math_vector_length $v]
    if {$l == 0} {
	mess "scrib_lib_math_vector_nrm: zero vector $v cannot be normalized \n" red
	return 0  
    }
    return [scrib_lib_math_vector_mult $v [expr 1/$l]]
}

proc scrib_lib_math_vector_sub {a b} {

  foreach p {a b} {
    set point [scrib_lib_math_make_real [subst $$p]]
    if ![scrib_lib_math_is_point $point] {
      mess "scrib_lib_math_vector_sub: Wrong parameter $p = \{$point\} \n" red
      return 0
    } else {
      set $p $point
    }
  }

  foreach ax {0 1 2} {
    set a_$ax [lindex $a $ax]
    set b_$ax [lindex $b $ax]
  }
  set v_0 [expr $a_0 - $b_0]
  set v_1 [expr $a_1 - $b_1]
  set v_2 [expr $a_2 - $b_2]
  return "$v_0 $v_1 $v_2"
}

proc scrib_lib_math_transform_cs {x y z center p0} {

  foreach p {x y z center p0} {
    set point [scrib_lib_math_make_real [subst $$p]]
    if ![scrib_lib_math_is_point $point] {
      mess "scrib_lib_math_transform_cs: Wrong parameter $p = \{$point\} \n" red
      return 0
    } else {
      set $p $point
    }
  }

  set m [scrib_lib_math_vector_sub $p0 $center]
  foreach ax {x y z} {lappend T [subst $$ax]}
  return [scrib_lib_math_mult_matrixes $T $m]
}


proc scrib_lib_math_transform_cs_back {x y z center p_n} {

  foreach p {x y z center p_n} {
    set point [scrib_lib_math_make_real [subst $$p]]
    if ![scrib_lib_math_is_point $point] {
      mess "scrib_lib_math_transform_cs_back: Wrong parameter $p = \{$point\} \n" red
      return 0
    } else {
      set $p $point
    }
  }

  foreach ax {x y z} {lappend T [subst $$ax]}
  return [scrib_lib_math_vector_sum [scrib_lib_math_mult_matrixes [scrib_lib_math_transp_matrix $T] $p_n] $center]
}



proc scrib_lib_math_mult_matrixes {M1 M2} {
  set n1 [llength $M1]
  set n2 [llength [lindex $M1 0]]
  for {set i 1} {$i < $n1} {incr i} {
    if {[llength [lindex $M1 0]] != $n2} {
      mess "Matrix is not rectangular \n" red
      scrib_lib_math_pmatrix $M1
      #va_err "Math error"
      return
    }
  }
  set n3 [llength $M2]
  if {$n3 != $n2} {
    mess "Columns number in matrix " red
    scrib_lib_math_pmatrix $M1
    mess "and rows number in matrix " red
    scrib_lib_math_pmatrix $M2
    mess "are not equal \n" red
    #va_err "Math error"
    return
  }
  set n4 [llength [lindex $M2 0]]
  for {set i 1} {$i < $n3} {incr i} {
    if {[llength [lindex $M2 0]] != $n4} {
      mess "Matrix is not rectangular \n" red
      scrib_lib_math_pmatrix $M2
      #va_err "Math error"
      return
    }
  }
  set M3 ""
  for {set i 0} {$i < $n1} {incr i} {
    set row ""
    for {set j 0} {$j < $n4} {incr j} {
      set c 0
      for {set k 0} {$k < $n2} {incr k} {
        set a [lindex [lindex $M1 $i] $k]
        set b [lindex [lindex $M2 $k] $j]
        set c [expr $c+$a*$b]
      }
      lappend row $c
    }
    lappend M3 "$row"
  }
  return $M3
}

proc scrib_lib_math_transp_matrix {M} {
  set n1 [llength $M]
  set n2 [llength [lindex $M 0]]
  set M1 ""
  for {set i 0} {$i < $n2} {incr i} {
    set row ""
    for {set j 0} {$j < $n1} {incr j} {
      lappend row [lindex [lindex $M $j] $i]
    }
    lappend M1 "$row"
  }
  return $M1
}

#================================================= 

proc _old_val_math_init {} {
    global pi rad
    
    set pi 3.1415927
    set rad [expr $pi/180.0]
}
#================================================= 
# Init 
#================================================= 

# i don't like the scope...--jka
#_old_val_math_init
#
# <------------ end stuff from Michigan

## returns the distance from coords <em>source</em> to <em>target</em>.

proc scrib_lib_math_dist_get {source target} {
    set dist2 0
    if {("" == $source) || ("" == $target)} { return [scrib_lib_math_dbl_max] }
    if {![scrib_lib_math_is_point $source] || ![scrib_lib_math_is_point $target]} {
	mess "invalid parameters given to scrib_lib_math_dist_get:\n \
		\t\{$source\} \{$target\}\n"
	return [scrib_lib_math_dbl_max]
    }
    foreach i {0 1 2} {
	set p [expr {[lindex $target $i] - [lindex $source $i]} ]
	set dist2 [expr {$dist2 + [expr {$p * $p}]}]
    }
    set dist [expr {sqrt($dist2)}]
    #mess "DEBUG: scrib_lib_math_get_dist: $dist from $source to $target\n" blue
    return $dist
}

## get a standard value. <em>mode</em> can be one of <i>{pi, rad, dbl_max}</i>

proc scrib_lib_math_get_val {{mode ""}} {
    set scrib_pi 3.14159265358979323846
    set dbl_max  99999999999999999999 ;# or some such val
    
    switch $mode {
	pi  { return $scrib_pi }
	rad { return [expr $scrib_pi/180.0] }
	dbl_max - flt_max - max { return $dbl_max }
	default {
	    mess "don't know about val $mode!\n" red
	}
    }
    return "__undefined__"
}


proc scrib_lib_math_pi {} { return [scrib_lib_math_get_val pi] }
proc scrib_lib_math_dbl_max {} { return [scrib_lib_math_get_val dbl_max] }


