Changeset 261 for trunk/readexp.tcl

Timestamp:

Dec 4, 2009 5:03:06 PM (13 years ago)

Author:

toby

Message:

# on 2000/08/08 23:21:55, toby did:
Correct creation of off-diagonal Uij terms for non-orthogonal cells

File:

• trunk/readexp.tcl (modified) (2 diffs, 3 props)

trunk/readexp.tcl

@@ -760 +760 @@
                     setexp ${key}B A 63 1
                     # copy the Uiso to the diagonal terms
-                    set value [string range [readexp ${key}B] 0 9]
-                    setexp ${key}B $value 11 10
-                    setexp ${key}B $value 21 10
-                    set value 0.0
-                    validreal value 10 6
-                    setexp ${key}B $value 31 10
-                    setexp ${key}B $value 41 10
-                    setexp ${key}B $value 51 10
+                    set Uiso [string range [readexp ${key}B] 0 9]
+                    foreach value [CalcAniso $phase $Uiso] \
+                            col {1 11 21 31 41 51} {
+                        validreal value 10 6
+                        setexp ${key}B $value $col 10
+                    }
                 } else {
                     setexp ${key}B I 63 1
@@ -2121 +2119 @@
     }
 }
+
+# compute equivalent anisotropic temperature factor for Uequiv
+proc CalcAniso {phase Uequiv} {
+    foreach var {a b c alpha beta gamma} {
+        set $var [phaseinfo $phase $var]
+    }
+
+    set G(1,1) [expr $a * $a]
+    set G(2,2) [expr $b * $b]
+    set G(3,3) [expr $c * $c]
+    set G(1,2) [expr $a * $b * cos($gamma*0.017453292519943)]
+    set G(2,1) $G(1,2)
+    set G(1,3) [expr $a * $c * cos($beta *0.017453292519943)]
+    set G(3,1) $G(1,3)
+    set G(2,3) [expr $b * $c * cos($alpha*0.017453292519943)]
+    set G(3,2) $G(2,3)
+
+    # Calculate the volume**2
+    set v2 0.0
+    foreach i {1 2 3} {
+        set J [expr ($i%3) + 1]
+        set K [expr (($i+1)%3) + 1]
+        set v2 [expr $v2+ $G(1,$i)*($G(2,$J)*$G(3,$K)-$G(3,$J)*$G(2,$K))]
+    }
+    if {$v2 > 0} {
+        set v [expr sqrt($v2)]
+        foreach i {1 2 3} {
+            set i1 [expr ($i%3) + 1]
+            set i2 [expr (($i+1)%3) + 1]
+            foreach j {1 2 3} {
+                set j1 [expr ($j%3) + 1]
+                set j2 [expr (($j+1)%3) + 1]
+                set C($j,$i) [expr (\
+                        $G($i1,$j1) * $G($i2,$j2) - \
+                        $G($i1,$j2)  * $G($i2,$j1)\
+                        )/ $v]
+            }
+        }
+        set A(1,2) [expr 0.5 * ($C(1,2) + $C(2,1)) / sqrt( $C(1,1)* $C(2,2) )]
+        set A(1,3) [expr 0.5 * ($C(1,3) + $C(3,1)) / sqrt( $C(1,1)* $C(3,3) )]
+        set A(2,3) [expr 0.5 * ($C(2,3) + $C(3,2)) / sqrt( $C(2,2)* $C(3,3) )]
+        foreach i {1 1 2} j {2 3 3} {
+            set A($i,$j) [expr 0.5 * ($C($i,$j) + $C($j,$i)) / \
+                    sqrt( $C($i,$i)* $C($j,$j) )]
+            # clean up roundoff
+            if {abs($A($i,$j)) < 1e-5} {set A($i,$j) 0.0}
+        }
+    } else {
+        set A(1,2) 0.0
+        set A(1,3) 0.0
+        set A(2,3) 0.0
+    }
+    return "$Uequiv $Uequiv $Uequiv \
+            [expr $Uequiv * $A(1,2)] \
+            [expr $Uequiv * $A(1,3)] \
+            [expr $Uequiv * $A(2,3)]"
+}