interfaces/xmaxima/Tkmaxima/Matrix.tcl

   1 # -*-mode: tcl; fill-column: 75; tab-width: 8; coding: iso-latin-1-unix -*-
   2 #
   3 #       $Id: Matrix.tcl,v 1.4 2006-10-01 23:58:29 villate Exp $
   4 #
   5 ###### Matrix.tcl ######
   6 ############################################################
   7 # Netmath       Copyright (C) 1998 William F. Schelter     #
   8 # For distribution under GNU public License.  See COPYING. #
   9 ############################################################
  10
  11 # In this file a matrix is represented by a list of M*N entries together
  12 # with an integer N giving the number of columns: {1 0 0 1} 2  would give
  13 # the two by two identity
  14
  15 proc comment {args } {
  16 }
  17 global mee
  18 set mee " } \] \[ expr { "
  19
  20 proc mkMultLeftExpr { mat n prefix { constant "" } } {
  21     #create a function body that does MAT (prefix1,prefix2,..) + constant
  22     global mee
  23     set all ""
  24
  25     set vars ""
  26     for { set i 0} { $i < $n} {incr i} { append vars " $prefix$i" }
  27     set j 0
  28     set k 0
  29
  30     foreach v $mat {
  31         if { $j == 0 } {
  32             set ro ""
  33             # append ans ""
  34             set op ""
  35         }
  36         append ro " $op $v*\$$prefix$j"
  37         set op "+"
  38         if { $j == [expr {$n -1}] } {
  39             append ans " "
  40             if { "[lindex $constant $k]" != "" } {
  41                 append ro " + [lindex $constant $k] "
  42             }
  43             incr k
  44             append ans [concat \[ expr [list $ro] \]]
  45             set j -1
  46         }
  47         incr j
  48     }
  49     # puts [list $vars $ans]
  50     return [list $vars $ans]
  51 }
  52
  53 proc mkMultLeftFun { mat n name { constant ""} } {
  54     set expr [mkMultLeftExpr $mat $n _a $constant]
  55     set bod1 [string trim [lindex $expr 1] " "]
  56     #    set bod "return \"$bod1\""
  57     set bod [concat list [lindex $expr 1]]
  58     proc $name [lindex $expr 0] $bod
  59 }
  60
  61 proc rotationMatrix { th ph } {
  62     return [list \
  63                 [expr {cos($ph)*cos($th)}] [expr {- cos($ph)*sin($th)}] [expr {sin($ph)}] \
  64                 [expr {sin($th)}] [expr {cos($th)}] 0.0 \
  65                 [expr {- sin($ph)*cos($th)}] [expr {sin($ph)*sin($th)}] [expr {cos($ph)}]]
  66 }
  67
  68 proc rotationMatrix { thx thy thz } {
  69     return \
  70         [list  \
  71              [expr { cos($thy)*cos($thz) } ] \
  72              [expr { cos($thy)*sin($thz) } ] \
  73              [expr { sin($thy) } ] \
  74              [expr { sin($thx)*sin($thy)*cos($thz)-cos($thx)*sin($thz) } ] \
  75              [expr { sin($thx)*sin($thy)*sin($thz)+cos($thx)*cos($thz) } ] \
  76              [expr { -sin($thx)*cos($thy) } ] \
  77              [expr { -sin($thx)*sin($thz)-cos($thx)*sin($thy)*cos($thz) } ] \
  78              [expr { sin($thx)*cos($thz)-cos($thx)*sin($thy)*sin($thz) } ] \
  79              [expr { cos($thx)*cos($thy) } ] ]
  80 }
  81
  82 # cross [a,b,c] [d,e,f] == [B*F-C*E,C*D-A*F,A*E-B*D]
  83 # cross_product([a,b,c],[d,e,f]):=[B*F-C*E,C*D-A*F,A*E-B*D]
  84 # cross_product(u,v):=sublis([a=u[1],b=u[2],c=u[3],d=v[1],e=v[2],f=v[3]],[B*F-C*E,C*D-A*F,A*E-B*D]);
  85 # the rotation by azimuth th, and elevation ph
  86 # MATRIX([COS(TH),SIN(TH),0],[-COS(PH)*SIN(TH),COS(PH)*COS(TH),SIN(PH)],
  87 #           [SIN(PH)*SIN(TH),-SIN(PH)*COS(TH),COS(PH)]);
  88
  89 proc rotationMatrix { th ph {ignore {} } } {
  90     return \
  91         [list \
  92              [      expr {cos($th)   } ]\
  93              [expr {sin($th)   } ]\
  94              0 \
  95              [expr {-cos($ph)*sin($th)   } ]\
  96              [expr {cos($ph)*cos($th)   } ]\
  97              [expr {sin($ph)   } ]\
  98              [expr {sin($ph)*sin($th)   } ]\
  99              [expr {-sin($ph)*cos($th)   } ]\
 100              [expr {cos($ph)   } ]]
 101 }
 102
 103 proc setMatFromList {name lis n} {
 104     set i 1
 105     set j 1
 106     foreach v $lis {
 107         uplevel 1 set [set name]($i,$j) $v
 108         if { $j == $n } {set j 1; incr i} else { incr j}
 109     }
 110 }
 111
 112 proc matRef { mat cols i j } {
 113     [lindex $mat [expr {$i*$cols + $j}]]
 114 }
 115
 116 proc matTranspose { mat cols } {
 117     set j 0
 118     set m [expr {[llength $mat ] / $cols}]
 119     while { $j < $cols} {
 120         set i 0
 121         while { $i < $m } {
 122             append ans " [lindex $mat [expr {$i*$cols + $j}]]"
 123             incr i
 124         }
 125         incr j
 126     }
 127     return $ans
 128 }
 129
 130
 131 proc matMul { mat1 cols1 mat2 cols2 } {
 132     mkMultLeftFun $mat1 $cols1 __tem
 133     set tr [matTranspose $mat2 $cols2]
 134     set rows1 [expr {[llength $mat1] / $cols1}]
 135     #puts "tr=$tr"
 136     set upto [llength $tr]
 137     set j 0
 138     set ans ""
 139     set i 0
 140     while { $j < $cols2  } {
 141         append ans " [eval __tem [lrange $tr $i [expr {$i+$cols1 -1}]]]"
 142         incr i $cols1
 143         incr j
 144     }
 145     #   return $ans
 146     # puts "matTranspose $ans $rows1"
 147     return [matTranspose $ans $rows1]
 148 }
 149
 150
 151
 152 proc invMat3 { mat } {
 153     setMatFromList xx $mat 3
 154     set det [expr { double($xx(1,1))*($xx(2,2)*$xx(3,3)-$xx(2,3)*$xx(3,2))-$xx(1,2)* \
 155                         ($xx(2,1)*$xx(3,3)-$xx(2,3)*$xx(3,1))+$xx(1,3)*($xx(2,1)*$xx(3,2)\
 156                                                                             -$xx(2,2)*$xx(3,1)) }]
 157
 158     return [list   [expr { ($xx(2,2)*$xx(3,3)-$xx(2,3)*$xx(3,2))/$det}] \
 159                 [expr { ($xx(1,3)*$xx(3,2)-$xx(1,2)*$xx(3,3))/$det}] \
 160                 [expr { ($xx(1,2)*$xx(2,3)-$xx(1,3)*$xx(2,2))/$det}] \
 161                 \
 162                 [expr { ($xx(2,3)*$xx(3,1)-$xx(2,1)*$xx(3,3))/$det}] \
 163                 [expr { ($xx(1,1)*$xx(3,3)-$xx(1,3)*$xx(3,1))/$det}] \
 164                 [expr { ($xx(1,3)*$xx(2,1)-$xx(1,1)*$xx(2,3))/$det}] \
 165                 \
 166                 [expr { ($xx(2,1)*$xx(3,2)-$xx(2,2)*$xx(3,1))/$det}] \
 167                 [expr { ($xx(1,2)*$xx(3,1)-$xx(1,1)*$xx(3,2))/$det}] \
 168                 [expr { ($xx(1,1)*$xx(2,2)-$xx(1,2)*$xx(2,1))/$det}]]
 169 }
 170
 171
 172 proc vectorOp { a op b} {
 173     set i [llength $a]
 174     set k 0
 175     set ans [expr [list [lindex $a 0]  $op [lindex $b 0]]]
 176     while { [incr k] < $i } {
 177         lappend ans  [expr  [list [lindex $a $k] $op [lindex $b $k]]]
 178     }
 179     return $ans
 180 }
 181 ## endsource matrix.tcl