source: trunk/readexp.tcl @ 265

Last change on this file since 265 was 265, checked in by toby, 14 years ago

# on 2000/08/17 23:44:52, toby did:
Fix bug setting gamma
  • Property rcs:author set to toby
  • Property rcs:date set to 2000/08/17 23:44:52
  • Property rcs:lines set to +1 -1
  • Property rcs:rev set to 1.21
  • Property rcs:state set to Exp
  • Property svn:keywords set to Author Date Revision Id
File size: 57.7 KB
Line 
1# $Id: readexp.tcl 265 2009-12-04 23:03:10Z toby $
2# Routines to deal with the .EXP "data structure"
3set expmap(Revision) {$Revision: 265 $ $Date: 2009-12-04 23:03:10 +0000 (Fri, 04 Dec 2009) $}
4
5#  The GSAS data is read from an EXP file.
6#   ... reading an EXP file into an array
7proc expload {expfile} {
8    global exparray
9    # $expfile is the path to the data file.
10    if [catch {set fil [open "$expfile" r]}] {
11        tk_dialog .expFileErrorMsg "File Open Error" \
12                "Unable to open file $expfile" error 0 "Exit" ; return -1
13    }
14    set len [gets $fil line]
15    if {[string length $line] != $len} {
16        tk_dialog .expConvErrorMsg "old tcl" \
17                "You are using an old version of Tcl/Tk and your .EXP file has binary characters; run convstod or upgrade" \
18                error 0 "Exit"
19        return -1
20    }
21    catch {
22        unset exparray
23    }
24    if {$len > 160} {
25        set fmt 0
26        # a UNIX-type file
27        set i1 0
28        set i2 79
29        while {$i2 < $len} {
30            set nline [string range $line $i1 $i2]
31            incr i1 80
32            incr i2 80
33            set key [string range $nline 0 11]
34            set exparray($key) [string range $nline 12 end]
35        }
36    } else {
37        set fmt 1
38        while {$len > 0} {
39            set key [string range $line 0 11]
40            set exparray($key) [string range $line 12 end]
41            set len [gets $fil line]
42        }
43    }
44    close $fil
45    return $fmt
46}
47
48proc createexp {expfile title} {
49    global exparray expmap
50    catch {unset exparray}
51    foreach key   {"     VERSION" "      DESCR" "ZZZZZZZZZZZZ" " EXPR NPHAS"} \
52            value {"   6"         ""            "  Last EXP file record" ""} {
53        # truncate long keys & pad short ones
54        set key [string range "$key        " 0 11]
55        set exparray($key) $value
56    }
57    expinfo title set $title
58    exphistory add " created readexp.tcl [lindex $expmap(Revision) 1] [clock format [clock seconds]]"
59    expwrite $expfile
60}
61
62# get information out from an EXP file
63#   creates the following entries in global array expmap
64#     expmap(phaselist)     gives a list of defined phases
65#     expmap(phasetype)     gives the phase type for each defined phase
66#                           =1 nuclear; 2 mag+nuc; 3 mag; 4 macro
67#     expmap(atomlist_$p)   gives a list of defined atoms in phase $p
68#     expmap(htype_$n)      gives the GSAS histogram type for histogram
69#     expmap(powderlist)    gives a list of powder histograms
70#     expmap(phaselist_$n)  gives a list of phases used in histogram $n
71#
72proc mapexp {} {
73    global expmap exparray
74    # clear out the old array
75    set expmap_Revision $expmap(Revision)
76    unset expmap
77    set expmap(Revision) $expmap_Revision
78    # get the defined phases
79    set line [readexp " EXPR NPHAS"]
80#    if {$line == ""} {
81#       set msg "No EXPR NPHAS entry. This is an invalid .EXP file"
82#       tk_dialog .badexp "Error in EXP" $msg error 0 Exit
83#       destroy .
84#    }
85    set expmap(phaselist) {}
86    set expmap(phasetype) {}
87    # loop over phases
88    foreach iph {1 2 3 4 5 6 7 8 9} {
89        set i5s [expr ($iph - 1)*5]
90        set i5e [expr $i5s + 4]
91        set flag [string trim [string range $line $i5s $i5e]]
92        if {$flag == ""} {set flag 0}
93        if $flag {
94            lappend expmap(phaselist) $iph
95            lappend expmap(phasetype) $flag
96        }
97    }
98    # get the list of defined atoms for each phase
99    foreach iph $expmap(phaselist) {
100        set expmap(atomlist_$iph) {}
101        foreach key [array names exparray "CRS$iph  AT*A"] {
102            regexp { AT *([0-9]+)A} $key a num
103            lappend expmap(atomlist_$iph) $num
104        }
105        # note that sometimes an .EXP file contains more atoms than are actually defined
106        # drop the extra ones
107        set expmap(atomlist_$iph) [lsort -integer $expmap(atomlist_$iph)]
108        set natom [phaseinfo $iph natoms]
109        if {$natom != [llength $expmap(atomlist_$iph)]} {
110            set expmap(atomlist_$iph) [lrange $expmap(atomlist_$iph) 0 [expr $natom-1]]
111        }
112    }
113    # now get the histogram types
114    set nhist [string trim [readexp { EXPR  NHST }]]
115    set n 0
116    set expmap(powderlist) {}
117    for {set i 0} {$i < $nhist} {incr i} {
118        set ihist [expr $i + 1]
119        if {[expr $i % 12] == 0} {
120            incr n
121            set line [readexp " EXPR  HTYP$n"]
122            if {$line == ""} {
123                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
124                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
125            }
126            set j 0
127        } else {
128            incr j
129        }
130        set expmap(htype_$ihist) [lindex $line $j]
131        # at least for now, ignore non-powder histograms
132        if {[string range $expmap(htype_$ihist) 0 0] == "P" && \
133                [string range $expmap(htype_$ihist) 3 3] != "*"} {
134            lappend expmap(powderlist) $ihist
135        }
136    }
137
138    # now process powder histograms
139    foreach ihist $expmap(powderlist) {
140        # make a 2 digit key -- hh
141        if {$ihist < 10} {
142            set hh " $ihist"
143        } else {
144            set hh $ihist
145        }
146        set line [readexp "HST $hh NPHAS"]
147        if {$line == ""} {
148            set msg "No NPHAS entry for Histogram $ihist. This is an invalid .EXP file"
149            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
150        }
151        set expmap(phaselist_$ihist) {}
152        # loop over phases
153        foreach iph {1 2 3 4 5 6 7 8 9} {
154            set i5s [expr ($iph - 1)*5]
155            set i5e [expr $i5s + 4]
156            set flag [string trim [string range $line $i5s $i5e]]
157            if {$flag == ""} {set flag 0}
158            if $flag {lappend expmap(phaselist_$ihist) $iph}
159        }
160    }
161}
162
163# return the value for a ISAM key
164proc readexp {key} {
165    global exparray
166    # truncate long keys & pad short ones
167    set key [string range "$key        " 0 11]
168    if [catch {set val $exparray($key)}] {
169        global expgui
170        if $expgui(debug) {puts "Error accessing record $key"}
171        return ""
172    }
173    return $val
174}
175
176# return the number of records matching ISAM key (may contain wildcards)
177proc existsexp {key} {
178    global exparray
179    # key can contain wild cards so don't pad
180    return [llength [array names exparray  $key]]
181}
182
183
184# replace a section of the exparray with $value
185#   replace $char characters starting at character $start (numbered from 1)
186proc setexp {key value start chars} {
187    global exparray
188    # truncate long keys & pad short ones
189    set key [string range "$key        " 0 11]
190    if [catch {set exparray($key)}] {
191        global expgui
192        if $expgui(debug) {puts "Error accessing record $key"}
193        return ""
194    }
195
196    # pad value to $chars
197    set l0 [expr $chars - 1]
198    set value [string range "$value                                           " 0 $l0]
199
200    if {$start == 1} {
201        set ret {}
202        set l1 $chars
203    } else {
204        set l0 [expr $start - 2]
205        set l1 [expr $start + $chars - 1]
206        set ret [string range $exparray($key) 0 $l0]
207    }
208    append ret $value [string range $exparray($key) $l1 end]
209    set exparray($key) $ret
210}
211
212proc makeexprec {key} {
213    global exparray
214    # truncate long keys & pad short ones
215    set key [string range "$key        " 0 11]
216    if [catch {set exparray($key)}] {
217        # set to 68 blanks
218        set exparray($key) [format %68s " "]
219    }
220}
221
222# delete an exp record
223# returns 1 if OK; 0 if not found
224proc delexp {key} {
225    global exparray
226    # truncate long keys & pad short ones
227    set key [string range "$key        " 0 11]
228    if [catch {unset exparray($key)}] {
229        return 0
230    }
231    return 1
232}
233# test an argument if it is a valid number; reform the number to fit
234proc validreal {val length decimal} {
235    upvar $val value
236    if [catch {expr $value}] {return 0}
237    if [catch {
238        set tmp [format "%${length}.${decimal}f" $value]
239        while {[string length $tmp] > $length} {
240            set tmp [format "%${length}.${decimal}E" $value]
241            incr decimal -1
242        }
243        set value $tmp
244    }] {return 0}
245    return 1
246}
247
248# test an argument if it is a valid integer; reform the number into
249# an integer, if appropriate -- be sure to pass the name of the variable not the value
250proc validint {val length} {
251    upvar $val value
252    # FORTRAN type assumption: blank is 0
253    if {$value == ""} {set value 0}
254    set tmp [expr round($value)]
255    if {$tmp != $value} {return 0}
256    if [catch {
257        set value [format "%${length}d" $tmp]
258    }] {return 0}
259    return 1
260}
261
262# process history information
263#    action == last
264#       returns number and value of last record
265#    action == add
266#
267proc exphistory {action "value 0"} {
268    global exparray
269    if {$action == "last"} {
270        set key [lindex [lsort -decreasing [array names exparray *HSTRY*]] 0]
271        if {$key == ""} {return ""}
272        return [list [string trim [string range $key 9 end]] $exparray($key)]
273    } elseif {$action == "add"} {
274        set key [lindex [lsort -decreasing [array names exparray *HSTRY*]] 0]
275        if {$key == ""} {
276            set index 1
277        } else {
278            set index [string trim [string range $key 9 end]]
279            if {$index != "***"} {
280                if {$index < 999} {incr index}
281                set key [format "    HSTRY%3d" $index]
282                set exparray($key) $value
283            }
284        }
285        set key [format "    HSTRY%3d" $index]
286        set exparray($key) $value
287    }
288}
289# get overall info
290#   parm:
291#     print     -- GENLES print option (*)
292#     cycles    -- number of GENLES cycles (*)
293#     title     -- the overall title (*)
294proc expinfo {parm "action get" "value {}"} {
295    switch ${parm}-$action {
296        title-get {
297            return [string trim [readexp "      DESCR"]]
298        }
299        title-set {
300            setexp "      DESCR" "  $value" 2 68
301        }
302
303        cycles-get {
304            return [string trim [cdatget MXCY]]
305        }
306        cycles-set {
307            if ![validint value 1] {return 0}
308            cdatset MXCY [format %4d $value]
309        }
310        print-get {
311            set print [string trim [cdatget PRNT]]
312            if {$print != ""} {return $print}
313            return 0
314        }
315        print-set {
316            if ![validint value 1] {return 0}
317            cdatset PRNT [format %3d $value]
318        }
319        default {
320            set msg "Unsupported expinfo access: parm=$parm action=$action"
321            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
322        }
323    }
324    return 1
325}
326
327proc cdatget {key} {
328    foreach i {1 2 3 4 5 6 7 8 9} {
329        if {[existsexp "  GNLS CDAT$i"] == 0} break
330        set line [readexp "  GNLS CDAT$i"]
331        if {$line == {}} break
332        foreach i1 {2 10 18 26 34 42 50 58 66} \
333                i2 {9 17 25 33 41 49 57 65 73} {
334            set item [string range $line $i1 $i2]
335            if {[string trim $item] == {}} continue
336            if [regexp "${key}(.*)" $item a b] {return $b}
337        }
338    }
339    return {}
340}
341
342proc cdatset {key value} {
343    # round 1 see if we can find the string
344    foreach i {1 2 3 4 5 6 7 8 9} {
345        set line [readexp "  GNLS CDAT$i"]
346        if {$line == {}} break
347        foreach i1 {2 10 18 26 34 42 50 58 66} \
348                i2 {9 17 25 33 41 49 57 65 73} {
349            set item [string range $line $i1 $i2]
350            if {[string trim $item] == {}} continue
351            if [regexp "${key}(.*)" $item a b] {
352                # found it now replace it
353                incr i1
354                setexp "  GNLS CDAT$i" "${key}${value}" $i1 8
355                return
356            }
357        }
358    }
359    # not found, take the 1st blank space, creating a card if needed
360    foreach i {1 2 3 4 5 6 7 8 9} {
361        set line [readexp "  GNLS CDAT$i"]
362        if {$line == {}} {makeexprec "  GNLS CDAT$i"}
363        foreach i1 {2 10 18 26 34 42 50 58 66} \
364                i2 {9 17 25 33 41 49 57 65 73} {
365            set item [string range $line $i1 $i2]
366            if {[string trim $item] == {}} {
367                # found a blank space: now replace it
368                incr i1
369                setexp "  GNLS CDAT$i" "${key}${value}" $i1 8
370                return
371            }
372        }
373    }
374    return {}
375}
376
377# get phase information: phaseinfo phase parm action value
378#   phase: 1 to 9 (as defined)
379#   parm:
380#     name -- phase name
381#     natoms -- number of atoms (*)
382#     a b c alpha beta gamma -- cell parameters (*)
383#     cellref -- refinement flag for the unit cell(*)
384#     celldamp  -- damping for the unit cell refinement (*)
385#     spacegroup -- space group symbol (*)
386#     ODForder -- spherical harmonic order (*)
387#     ODFsym   -- sample symmetry (0-3) (*)
388#     ODFdampA -- damping for angles (*)
389#     ODFdampC -- damping for coefficients (*)
390#     ODFomega -- omega oriention angle (*)
391#     ODFchi -- chi oriention angle (*)
392#     ODFphi -- phi oriention angle (*)
393#     ODFomegaRef -- refinement flag for omega (*)
394#     ODFchiRef -- refinement flag for chi (*)
395#     ODFphiRef -- refinement flag for phi (*)
396#     ODFterms -- a list of the {l m n} values for each ODF term (*)
397#     ODFcoefXXX -- the ODF coefficient for for ODF term XXX (*)
398#     ODFRefcoef -- refinement flag for ODF terms (*)
399#  action: get (default) or set
400#  value: used only with set
401#  * =>  read+write supported
402proc phaseinfo {phase parm "action get" "value {}"} {
403    switch -glob ${parm}-$action {
404
405        name-get {
406            return [string trim [readexp "CRS$phase    PNAM"]]
407        }
408
409        name-set {
410            setexp "CRS$phase    PNAM" " $value" 2 68
411        }
412
413        spacegroup-get {
414            return [string trim [readexp "CRS$phase  SG SYM"]]
415        }
416
417        spacegroup-set {
418            setexp "CRS$phase  SG SYM" " $value" 2 68
419        }
420
421        natoms-get {
422            return [string trim [readexp "CRS$phase   NATOM"]]     
423        }
424
425        natoms-set {
426            if ![validint value 5] {return 0}
427            setexp "CRS$phase   NATOM" $value 1 5
428        }
429
430        a-get {
431           return [string trim [string range [readexp "CRS$phase  ABC"] 0 9]]
432        }
433        b-get {
434           return [string trim [string range [readexp "CRS$phase  ABC"] 10 19]]
435        }
436        c-get {
437           return [string trim [string range [readexp "CRS$phase  ABC"] 20 29]]
438        }
439        alpha-get {
440           return [string trim [string range [readexp "CRS$phase  ANGLES"] 0 9]]
441        }
442        beta-get {
443           return [string trim [string range [readexp "CRS$phase  ANGLES"] 10 19]]
444        }
445        gamma-get {
446           return [string trim [string range [readexp "CRS$phase  ANGLES"] 20 29]]
447        }
448
449        a-set {
450            if ![validreal value 10 6] {return 0}
451            setexp "CRS$phase  ABC" $value 1 10             
452        }
453        b-set {
454            if ![validreal value 10 6] {return 0}
455            setexp "CRS$phase  ABC" $value 11 10           
456        }
457        c-set {
458            if ![validreal value 10 6] {return 0}
459            setexp "CRS$phase  ABC" $value 21 10           
460        }
461        alpha-set {
462            if ![validreal value 10 4] {return 0}
463            setexp "CRS$phase  ANGLES" $value 1 10         
464        }
465        beta-set {
466            if ![validreal value 10 4] {return 0}
467            setexp "CRS$phase  ANGLES" $value 11 10         
468        }
469        gamma-set {
470            if ![validreal value 10 4] {return 0}
471            setexp "CRS$phase  ANGLES" $value 21 10         
472        }
473        cellref-get {
474            if {[string toupper [string range [readexp "CRS$phase  ABC"] 34 34]] == "Y"} {
475                return 1
476            }
477            return 0
478        }
479        cellref-set {
480            if $value {
481                setexp "CRS$phase  ABC" "Y" 35 1
482            } else {
483                setexp "CRS$phase  ABC" "N" 35 1
484            }       
485        }
486        celldamp-get {
487            set val [string range [readexp "CRS$phase  ABC"] 39 39]
488            if {$val == " "} {return 0}
489            return $val
490        }
491        celldamp-set {
492            setexp "CRS$phase  ABC" $value 40 1
493        }
494
495        ODForder-get {
496            set val [string trim [string range [readexp "CRS$phase  ODF"] 0 4]]
497            if {$val == " "} {return 0}
498            return $val
499        }
500        ODForder-set {
501            if ![validint value 5] {return 0}
502            setexp "CRS$phase  ODF" $value 1 5
503        }
504        ODFsym-get {
505            set val [string trim [string range [readexp "CRS$phase  ODF"] 10 14]]
506            if {$val == " "} {return 0}
507            return $val
508        }
509        ODFsym-set {
510            if ![validint value 5] {return 0}
511            setexp "CRS$phase  ODF" $value 11 5
512        }
513        ODFdampA-get {
514            set val [string range [readexp "CRS$phase  ODF"] 24 24]
515            if {$val == " "} {return 0}
516            return $val
517        }
518        ODFdampA-set {
519            setexp "CRS$phase  ODF" $value 25 1
520        }
521        ODFdampC-get {
522            set val [string range [readexp "CRS$phase  ODF"] 29 29]
523            if {$val == " "} {return 0}
524            return $val
525        }
526        ODFdampC-set {
527            setexp "CRS$phase  ODF" $value 30 1
528        }
529        ODFomegaRef-get {
530            if {[string toupper [string range [readexp "CRS$phase  ODF"] 16 16]] == "Y"} {
531                return 1
532            }
533            return 0
534        }
535        ODFomegaRef-set {
536            if $value {
537                setexp "CRS$phase  ODF" "Y" 17 1
538            } else {
539                setexp "CRS$phase  ODF" "N" 17 1
540            }       
541        }
542        ODFchiRef-get {
543            if {[string toupper [string range [readexp "CRS$phase  ODF"] 17 17]] == "Y"} {
544                return 1
545            }
546            return 0
547        }
548        ODFchiRef-set {
549            if $value {
550                setexp "CRS$phase  ODF" "Y" 18 1
551            } else {
552                setexp "CRS$phase  ODF" "N" 18 1
553            }       
554        }
555        ODFphiRef-get {
556            if {[string toupper [string range [readexp "CRS$phase  ODF"] 18 18]] == "Y"} {
557                return 1
558            }
559            return 0
560        }
561        ODFphiRef-set {
562            if $value {
563                setexp "CRS$phase  ODF" "Y" 19 1
564            } else {
565                setexp "CRS$phase  ODF" "N" 19 1
566            }       
567        }
568        ODFcoef*-get {
569            regsub ODFcoef $parm {} term
570            set k [expr ($term+5)/6]
571            if {$k <= 9} {set k " $k"}
572            set j [expr (($term-1) % 6)+1]
573            set lineB [readexp "CRS$phase  ODF${k}B"]
574            set j0 [expr  ($j-1) *10]
575            set j1 [expr $j0 + 9]
576            set val [string trim [string range $lineB $j0 $j1]]
577            if {$val == ""} {return 0.0}
578            return $val
579        }
580        ODFcoef*-set {
581            regsub ODFcoef $parm {} term
582            if ![validreal value 10 3] {return 0}
583            set k [expr ($term+5)/6]
584            if {$k <= 9} {set k " $k"}
585            set j [expr (($term-1) % 6)+1]
586            set col [expr  ($j-1)*10 + 1]
587            setexp "CRS$phase  ODF${k}B" $value $col 10
588        }
589        ODFRefcoef-get {
590            if {[string toupper [string range [readexp "CRS$phase  ODF"] 19 19]] == "Y"} {
591                return 1
592            }
593            return 0
594        }
595        ODFRefcoef-set {
596            if $value {
597                setexp "CRS$phase  ODF" "Y" 20 1
598            } else {
599                setexp "CRS$phase  ODF" "N" 20 1
600            }       
601        }
602        ODFomega-get {
603           return [string trim [string range [readexp "CRS$phase  ODF"] 30 39]]
604        }
605        ODFchi-get {
606           return [string trim [string range [readexp "CRS$phase  ODF"] 40 49]]
607        }
608        ODFphi-get {
609           return [string trim [string range [readexp "CRS$phase  ODF"] 50 59]]
610        }
611        ODFomega-set {
612            if ![validreal value 10 4] {return 0}
613            setexp "CRS$phase  ODF" $value 31 10
614        }
615        ODFchi-set {
616            if ![validreal value 10 4] {return 0}
617            setexp "CRS$phase  ODF" $value 41 10
618        }
619        ODFphi-set {
620            if ![validreal value 10 4] {return 0}
621            setexp "CRS$phase  ODF" $value 51 10
622        }
623
624        ODFterms-get {
625            set vallist {}
626            set val [string trim [string range [readexp "CRS$phase  ODF"] 5 9]]
627            for {set i 1} {$i <= $val} {incr i 6} {
628                set k [expr 1+($i-1)/6]
629                if {$k <= 9} {set k " $k"}
630                set lineA [readexp "CRS$phase  ODF${k}A"]
631                set k 0
632                for {set j $i} {$j <= $val && $j < $i+6} {incr j} {
633                    set j0 [expr ($k)*10]
634                    set j1 [expr $j0 + 9]
635                    lappend vallist [string trim [string range $lineA $j0 $j1]]
636                    incr k
637                }
638            }
639            return $vallist
640        }
641        ODFterms-set {
642            set key "CRS$phase  ODF   "
643            if {![existsexp $key]} {
644                makeexprec $key
645                set oldlen 0
646            } else {
647                set oldlen [string trim [string range [readexp $key] 5 9]]
648            }
649            set len [llength $value]
650            if ![validint len 5] {return 0}
651            setexp $key $len 6 5
652            set j 0
653            set k 0
654            foreach item $value {
655                incr j
656                if {$j % 6 == 1} {
657                    incr k
658                    if {$k <= 9} {set k " $k"}
659                    set col 1
660                    set keyA "CRS$phase  ODF${k}A"
661                    set keyB "CRS$phase  ODF${k}B"
662                    if {![existsexp $keyA]} {
663                        makeexprec $keyA
664                        makeexprec $keyB
665                    }
666                }
667                set col1 [expr $col + 1]
668                foreach n [lrange $item 0 2] {
669                    if ![validint n 3] {return 0}
670                    setexp $keyA $n $col1 3
671                    incr col1 3
672                }
673                incr col 10
674            }
675            for {incr j} {$j <= $oldlen} {incr j} {
676                if {$j % 6 == 1} {
677                    incr k
678                    if {$k <= 9} {set k " $k"}
679                    set col 1
680                    set keyA "CRS$phase  ODF${k}A"
681                    set keyB "CRS$phase  ODF${k}B"
682                    delexp $keyA
683                    delexp $keyB
684                }
685                if {[existsexp $keyA]} {
686                    setexp $keyA "          " $col 10
687                    setexp $keyB "          " $col 10
688                }
689                incr col 10
690            }
691        }
692
693        default {
694            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
695            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
696        }
697    }
698    return 1
699}
700
701
702# get atom information: atominfo phase atom parm action value
703#   phase: 1 to 9 (as defined)
704#   atom: a valid atom number [see expmap(atomlist_$phase)]
705#      Note that atom and phase can be paired lists, but if there are extra
706#      entries in the atoms list, the last phase will be repeated.
707#      so that atominfo 1 {1 2 3} xset 1
708#               will set the xflag for atoms 1-3 in phase 1
709#      but atominfo {1 2 3} {1 1 1} xset 1
710#               will set the xflag for atoms 1 in phase 1-3
711#   parm:
712#     type -- element code
713#     mult -- atom multiplicity
714#     label -- atom label (*)
715#     x y z -- coordinates (*)
716#     frac --  occupancy (*)
717#     temptype -- I or A for Isotropic/Anisotropic (*)
718#     Uiso  -- Isotropic temperature factor (*)
719#     U11  -- Anisotropic temperature factor (*)
720#     U22  -- Anisotropic temperature factor (*)
721#     U33  -- Anisotropic temperature factor (*)
722#     U12  -- Anisotropic temperature factor (*)
723#     U13  -- Anisotropic temperature factor (*)
724#     U23  -- Anisotropic temperature factor (*)
725#     xref/xdamp -- refinement flag/damping value for the coordinates (*)
726#     uref/udamp -- refinement flag/damping value for the temperature factor(s)  (*)
727#     fref/fdamp -- refinement flag/damping value for the occupancy (*)
728#  action: get (default) or set
729#  value: used only with set
730#  * =>  read+write supported
731
732proc atominfo {phaselist atomlist parm "action get" "value {}"} {
733    foreach phase $phaselist atom $atomlist {
734        if {$phase == ""} {set phase [lindex $phaselist end]}
735        if {$atom < 10} {
736            set key "CRS$phase  AT  $atom"
737        } elseif {$atom < 100} {
738            set key "CRS$phase  AT $atom"
739        } else {
740            set key "CRS$phase  AT$atom"
741        }
742        switch -glob ${parm}-$action {
743            type-get {
744                return [string trim [string range [readexp ${key}A] 2 9] ]
745            }
746            mult-get {
747                return [string trim [string range [readexp ${key}A] 58 61] ]
748            }
749            label-get {
750                return [string trim [string range [readexp ${key}A] 50 57] ]
751            }
752            label-set {
753                setexp ${key}A $value 51 8
754            }
755            temptype-get {
756                return [string trim [string range [readexp ${key}B] 62 62] ]
757            }
758            temptype-set {
759                if {$value == "A"} {
760                    setexp ${key}B A 63 1
761                    # copy the Uiso to the diagonal terms
762                    set Uiso [string range [readexp ${key}B] 0 9]
763                    foreach value [CalcAniso $phase $Uiso] \
764                            col {1 11 21 31 41 51} {
765                        validreal value 10 6
766                        setexp ${key}B $value $col 10
767                    }
768                } else {
769                    setexp ${key}B I 63 1
770                    set value 0.0
771                    catch {
772                        # get the trace
773                        set value [expr ( \
774                                [string range [readexp ${key}B] 0 9] + \
775                                [string range [readexp ${key}B] 10 19] + \
776                                [string range [readexp ${key}B] 20 29])/3.]
777                    }
778                    validreal value 10 6
779                    setexp ${key}B $value 1 10
780                    # blank out the remaining terms
781                    set value " "
782                    setexp ${key}B $value 11 10
783                    setexp ${key}B $value 21 10
784                    setexp ${key}B $value 31 10
785                    setexp ${key}B $value 41 10
786                    setexp ${key}B $value 51 10
787                }
788            }
789            x-get {
790                return [string trim [string range [readexp ${key}A] 10 19] ]
791            }
792            x-set {
793                if ![validreal value 10 6] {return 0}
794                setexp ${key}A $value 11 10
795            }
796            y-get {
797                return [string trim [string range [readexp ${key}A] 20 29] ]
798            }
799            y-set {
800                if ![validreal value 10 6] {return 0}
801                setexp ${key}A $value 21 10
802            }
803            z-get {
804                return [string trim [string range [readexp ${key}A] 30 39] ]
805            }
806            z-set {
807                if ![validreal value 10 6] {return 0}
808                setexp ${key}A $value 31 10
809            }
810            frac-get {
811                return [string trim [string range [readexp ${key}A] 40 49] ]
812            }
813            frac-set {
814                if ![validreal value 10 6] {return 0}
815                setexp ${key}A $value 41 10
816            }
817            U*-get {
818                regsub U $parm {} type
819                if {$type == "iso" || $type == "11"} {
820                    return [string trim [string range [readexp ${key}B] 0 9] ]
821                } elseif {$type == "22"} {
822                    return [string trim [string range [readexp ${key}B] 10 19] ]
823                } elseif {$type == "33"} {
824                    return [string trim [string range [readexp ${key}B] 20 29] ]
825                } elseif {$type == "12"} {
826                    return [string trim [string range [readexp ${key}B] 30 39] ]
827                } elseif {$type == "13"} {
828                    return [string trim [string range [readexp ${key}B] 40 49] ]
829                } elseif {$type == "23"} {
830                    return [string trim [string range [readexp ${key}B] 50 59] ]
831                }
832            }
833            U*-set {
834                if ![validreal value 10 6] {return 0}
835                regsub U $parm {} type
836                if {$type == "iso" || $type == "11"} {
837                    setexp ${key}B $value 1 10
838                } elseif {$type == "22"} {
839                    setexp ${key}B $value 11 10
840                } elseif {$type == "33"} {
841                    setexp ${key}B $value 21 10
842                } elseif {$type == "12"} {
843                    setexp ${key}B $value 31 10
844                } elseif {$type == "13"} {
845                    setexp ${key}B $value 41 10
846                } elseif {$type == "23"} {
847                    setexp ${key}B $value 51 10
848                }
849            }
850            xref-get {
851                if {[string toupper [string range [readexp ${key}B] 64 64]] == "X"} {
852                    return 1
853                }
854                return 0
855            }
856            xref-set {
857                if $value {
858                    setexp ${key}B "X" 65 1
859                } else {
860                    setexp ${key}B " " 65 1
861                }           
862            }
863            xdamp-get {
864                set val [string range [readexp ${key}A] 64 64]
865                if {$val == " "} {return 0}
866                return $val
867            }
868            xdamp-set {
869                setexp ${key}A $value 65 1
870            }
871            fref-get {
872                if {[string toupper [string range [readexp ${key}B] 63 63]] == "F"} {
873                    return 1
874                }
875                return 0
876            }
877            fref-set {
878                if $value {
879                    setexp ${key}B "F" 64 1
880                } else {
881                    setexp ${key}B " " 64 1
882                }           
883            }
884            fdamp-get {
885                set val [string range [readexp ${key}A] 63 63]
886                if {$val == " "} {return 0}
887                return $val
888            }
889            fdamp-set {
890                setexp ${key}A $value 64 1
891            }
892
893            uref-get {
894                if {[string toupper [string range [readexp ${key}B] 65 65]] == "U"} {
895                    return 1
896                }
897                return 0
898            }
899            uref-set {
900                if $value {
901                    setexp ${key}B "U" 66 1
902                } else {
903                    setexp ${key}B " " 66 1
904                }           
905            }
906            udamp-get {
907                set val [string range [readexp ${key}A] 65 65]
908                if {$val == " "} {return 0}
909                return $val
910            }
911            udamp-set {
912                setexp ${key}A $value 66 1
913            }
914            default {
915                set msg "Unsupported atominfo access: parm=$parm action=$action"
916                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
917            }
918        }
919    }
920    return 1
921}
922
923# get histogram information: histinfo histlist parm action value
924# histlist is a list of histogram numbers
925# parm:
926#     title
927#     scale (*)
928#     sref/sdamp -- refinement flag/damping value for the scale factor (*)
929#     lam1, lam2 (*)
930#     ttref refinement flag for the 2theta (ED Xray) (*)
931#     wref refinement flag for the wavelength (*)
932#     ratref refinement flag for the wavelength ratio (*)
933#     difc, difa -- TOF calibration constants (*)
934#     dcref,daref -- refinement flag for difc, difa (*)
935#     zero (*)
936#     zref refinement flag for the zero correction (*)
937#     ipola (*)
938#     pola (*)
939#     pref refinement flag for the polarization (*)
940#     kratio (*)
941#     ddamp -- damping value for the diffractometer constants (*)
942#     backtype -- background function number *
943#     backterms -- number of background terms *
944#     bref/bdamp -- refinement flag/damping value for the background (*)
945#     bterm$n -- background term #n (*)
946#     bank -- Bank number
947#     tofangle -- detector angle (TOF only)
948#     foextract  -- Fobs extraction flag (*)
949proc histinfo {histlist parm "action get" "value {}"} {
950    global expgui
951    foreach hist $histlist {
952        if {$hist < 10} {
953            set key "HST  $hist"
954        } else {
955            set key "HST $hist"
956        }
957        switch -glob ${parm}-$action {
958            foextract-get {
959                set line [readexp "${key} EPHAS"]
960                # add a EPHAS if not exists
961                if {$line == {}} {
962                    makeexprec "${key} EPHAS"
963                    # expedt defaults this to "F", but I think "T" is better
964                    setexp "${key} EPHAS" "T" 50 1
965                    if $expgui(debug) {puts "Warning: creating a ${key} EPHAS record"}
966                }
967                if {[string toupper [string range $line 49 49]] == "T"} {
968                    return 1
969                }
970                return 0
971            }
972            foextract-set {
973                if $value {
974                    setexp "${key} EPHAS" "T" 50 1
975                } else {
976                    setexp "${key} EPHAS" "F" 50 1
977                }           
978            }
979            title-get {
980                return [string trim [readexp "${key}  HNAM"] ]
981            }
982            scale-get {
983                return [string trim [string range [readexp ${key}HSCALE] 0 14]]
984            }
985            scale-set {
986                if ![validreal value 15 6] {return 0}
987                setexp ${key}HSCALE $value 1 15
988            }
989            sref-get {
990                if {[string toupper [string range [readexp ${key}HSCALE] 19 19]] == "Y"} {
991                    return 1
992                }
993                return 0
994            }
995            sref-set {
996                if $value {
997                    setexp ${key}HSCALE "Y" 20 1
998                } else {
999                    setexp ${key}HSCALE "N" 20 1
1000                }           
1001            }
1002            sdamp-get {
1003                set val [string range [readexp ${key}HSCALE] 24 24]
1004                if {$val == " "} {return 0}
1005                return $val
1006            }
1007            sdamp-set {
1008                setexp ${key}HSCALE $value 25 1
1009            }
1010
1011            difc-get -
1012            lam1-get {
1013                return [string trim [string range [readexp "${key} ICONS"] 0 9]]
1014            }
1015            difc-set -
1016            lam1-set {
1017                if ![validreal value 10 7] {return 0}
1018                setexp "${key} ICONS" $value 1 10
1019            }
1020            difa-get -
1021            lam2-get {
1022                return [string trim [string range [readexp "${key} ICONS"] 10 19]]
1023            }
1024            difa-set -
1025            lam2-set {
1026                if ![validreal value 10 7] {return 0}
1027                setexp "${key} ICONS" $value 11 10
1028            }
1029            zero-get {
1030                return [string trim [string range [readexp "${key} ICONS"] 20 29]]
1031            }
1032            zero-set {
1033                if ![validreal value 10 5] {return 0}
1034                setexp "${key} ICONS" $value 21 10
1035            }
1036            ipola-get {
1037                return [string trim [string range [readexp "${key} ICONS"] 54 54]]
1038            }
1039            ipola-set {
1040                if ![validint value 1] {return 0}
1041                setexp "${key} ICONS" $value 55 1
1042            }
1043            pola-get {
1044                return [string trim [string range [readexp "${key} ICONS"] 40 49]]
1045            }
1046            pola-set {
1047                if ![validreal value 10 5] {return 0}
1048                setexp "${key} ICONS" $value 41 10
1049            }
1050            kratio-get {
1051                return [string trim [string range [readexp "${key} ICONS"] 55 64]]
1052            }
1053            kratio-set {
1054                if ![validreal value 10 5] {return 0}
1055                setexp "${key} ICONS" $value 56 10
1056            }
1057
1058            wref-get {
1059            #------------------------------------------------------
1060            # col 33: refine flag for lambda, difc, ratio and theta
1061            #------------------------------------------------------
1062                if {[string toupper [string range \
1063                        [readexp "${key} ICONS"] 32 32]] == "L"} {
1064                    return 1
1065                }
1066                return 0
1067            }
1068            wref-set {
1069                if $value {
1070                    setexp "${key} ICONS" "L" 33 1
1071                } else {
1072                    setexp "${key} ICONS" " " 33 1
1073                }           
1074            }
1075            ratref-get {
1076                if {[string toupper [string range \
1077                        [readexp "${key} ICONS"] 32 32]] == "R"} {
1078                    return 1
1079                }
1080                return 0
1081            }
1082            ratref-set {
1083                if $value {
1084                    setexp "${key} ICONS" "R" 33 1
1085                } else {
1086                    setexp "${key} ICONS" " " 33 1
1087                }           
1088            }
1089            dcref-get {
1090                if {[string toupper [string range \
1091                        [readexp "${key} ICONS"] 32 32]] == "C"} {
1092                    return 1
1093                }
1094                return 0
1095            }
1096            dcref-set {
1097                if $value {
1098                    setexp "${key} ICONS" "C" 33 1
1099                } else {
1100                    setexp "${key} ICONS" " " 33 1
1101                }           
1102            }
1103            ttref-get {
1104                if {[string toupper [string range \
1105                        [readexp "${key} ICONS"] 32 32]] == "T"} {
1106                    return 1
1107                }
1108                return 0
1109            }
1110            ttref-set {
1111                if $value {
1112                    setexp "${key} ICONS" "T" 33 1
1113                } else {
1114                    setexp "${key} ICONS" " " 33 1
1115                }           
1116            }
1117
1118
1119            pref-get {
1120            #------------------------------------------------------
1121            # col 34: refine flag for POLA & DIFA
1122            #------------------------------------------------------
1123                if {[string toupper [string range \
1124                        [readexp "${key} ICONS"] 33 33]] == "P"} {
1125                    return 1
1126                }
1127                return 0
1128            }
1129            pref-set {
1130                if $value {
1131                    setexp "${key} ICONS" "P" 34 1
1132                } else {
1133                    setexp "${key} ICONS" " " 34 1
1134                }           
1135            }
1136            daref-get {
1137                if {[string toupper [string range \
1138                        [readexp "${key} ICONS"] 33 33]] == "A"} {
1139                    return 1
1140                }
1141                return 0
1142            }
1143            daref-set {
1144                if $value {
1145                    setexp "${key} ICONS" "A" 34 1
1146                } else {
1147                    setexp "${key} ICONS" " " 34 1
1148                }           
1149            }
1150
1151            zref-get {
1152            #------------------------------------------------------
1153            # col 34: refine flag for zero correction
1154            #------------------------------------------------------
1155                if {[string toupper [string range [readexp "${key} ICONS"] 34 34]] == "Z"} {
1156                    return 1
1157                }
1158                return 0
1159            }
1160            zref-set {
1161                if $value {
1162                    setexp "${key} ICONS" "Z" 35 1
1163                } else {
1164                    setexp "${key} ICONS" " " 35 1
1165                }           
1166            }
1167
1168            ddamp-get {
1169                set val [string range [readexp "${key} ICONS"] 39 39]
1170                if {$val == " "} {return 0}
1171                return $val
1172            }
1173            ddamp-set {
1174                setexp "${key} ICONS" $value 40 1
1175            }
1176
1177            backtype-get {
1178                set val [string trim [string range [readexp "${key}BAKGD "] 0 4]]
1179                if {$val == " "} {return 0}
1180                return $val
1181            }
1182            backtype-set {
1183                if ![validint value 5] {return 0}
1184                setexp "${key}BAKGD " $value 1 5
1185            }
1186            backterms-get {
1187                set val [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1188                if {$val == " "} {return 0}
1189                return $val
1190            }
1191            backterms-set {
1192                # this takes a bit of work -- if terms are added, add lines as needed to the .EXP
1193                set oldval [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1194                if ![validint value 5] {return 0}
1195                if {$oldval < $value} {
1196                    set line1  [expr 2 + ($oldval - 1) / 4]
1197                    set line2  [expr 1 + ($value - 1) / 4]
1198                    for {set i $line1} {$i <= $line2} {incr i} {
1199                        # create a blank entry if needed
1200                        makeexprec ${key}BAKGD$i
1201                    }
1202                    incr oldval
1203                    for {set num $oldval} {$num <= $value} {incr num} {
1204                        set f1 [expr 15*(($num - 1) % 4)]
1205                        set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1206                        set line  [expr 1 + ($num - 1) / 4]
1207                        if {[string trim [string range [readexp ${key}BAKGD$line] $f1 $f2]] == ""} {
1208                            set f1 [expr 15*(($num - 1) % 4)+1]
1209                            setexp ${key}BAKGD$line 0.0 $f1 15                 
1210                        }
1211                    }
1212                }
1213                setexp "${key}BAKGD " $value 6 5
1214
1215            }
1216            bref-get {
1217                if {[string toupper [string range [readexp "${key}BAKGD"] 14 14]] == "Y"} {
1218                    return 1
1219                }
1220                return 0
1221            }
1222            bref-set {
1223                if $value {
1224                    setexp "${key}BAKGD "  "Y" 15 1
1225                } else {
1226                    setexp "${key}BAKGD "  "N" 15 1
1227                }           
1228            }
1229            bdamp-get {
1230                set val [string range [readexp "${key}BAKGD "] 19 19]
1231                if {$val == " "} {return 0}
1232                return $val
1233            }
1234            bdamp-set {
1235                setexp "${key}BAKGD " $value 20 1
1236            }
1237            bterm*-get {
1238                regsub bterm $parm {} num
1239                set f1 [expr 15*(($num - 1) % 4)]
1240                set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1241                set line  [expr 1 + ($num - 1) / 4]
1242                return [string trim [string range [readexp ${key}BAKGD$line] $f1 $f2] ]
1243            }
1244            bterm*-set {
1245                regsub bterm $parm {} num
1246                if ![validreal value 15 6] {return 0}
1247                set f1 [expr 15*(($num - 1) % 4)+1]
1248                set line  [expr 1 + ($num - 1) / 4]
1249                setexp ${key}BAKGD$line $value $f1 15
1250            }
1251            bank-get {
1252                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1253            }
1254            tofangle-get {
1255                return [string trim [string range [readexp "${key}BNKPAR"] 10 19]]
1256            }
1257            default {
1258                set msg "Unsupported histinfo access: parm=$parm action=$action"
1259                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1260            }
1261        }
1262    }
1263    return 1
1264}
1265
1266# read the information that differs by both histogram and phase (profile & phase fraction)
1267# use: hapinfo hist phase parm action value
1268# 
1269#     frac -- phase fraction (*)
1270#     frref/frdamp -- refinement flag/damping value for the phase fraction (*)
1271#     proftype -- profile function number (*)
1272#     profterms -- number of profile terms (*)
1273#     pdamp -- damping value for the profile (*)
1274#     pcut -- cutoff value for the profile (*)
1275#     pterm$n -- profile term #n (*)
1276#     pref$n -- refinement flag value for profile term #n (*)
1277#     extmeth -- Fobs extraction method (*)
1278#     POnaxis -- number of defined M-D preferred axes
1279proc hapinfo {histlist phaselist parm "action get" "value {}"} {
1280    foreach phase $phaselist hist $histlist {
1281        if {$phase == ""} {set phase [lindex $phaselist end]}
1282        if {$hist == ""} {set hist [lindex $histlist end]}
1283        if {$hist < 10} {
1284            set hist " $hist"
1285        }
1286        set key "HAP${phase}${hist}"
1287        switch -glob ${parm}-$action {
1288            extmeth-get {
1289                set i1 [expr ($phase - 1)*5]
1290                set i2 [expr $i1 + 4]
1291                return [string trim [string range [readexp "HST $hist EPHAS"] $i1 $i2]]
1292            }
1293            extmeth-set {
1294                set i1 [expr ($phase - 1)*5 + 1]
1295                if ![validint value 5] {return 0}
1296                setexp "HST $hist EPHAS" $value $i1 5
1297            }
1298            frac-get {
1299                return [string trim [string range [readexp ${key}PHSFR] 0 14]]
1300            }
1301            frac-set {
1302                if ![validreal value 15 6] {return 0}
1303                setexp ${key}PHSFR $value 1 15
1304            }
1305            frref-get {
1306                if {[string toupper [string range [readexp ${key}PHSFR] 19 19]] == "Y"} {
1307                    return 1
1308                }
1309                return 0
1310            }
1311            frref-set {
1312                if $value {
1313                    setexp ${key}PHSFR "Y" 20 1
1314                } else {
1315                    setexp ${key}PHSFR "N" 20 1
1316                }           
1317            }
1318            frdamp-get {
1319                set val [string range [readexp ${key}PHSFR] 24 24]
1320                if {$val == " "} {return 0}
1321                return $val
1322            }
1323            frdamp-set {
1324                setexp ${key}PHSFR $value 25 1
1325            }
1326            proftype-get {
1327                set val [string range [readexp "${key}PRCF "] 0 4]
1328                if {$val == " "} {return 0}
1329                return $val
1330            }
1331            proftype-set {
1332                if ![validint value 5] {return 0}
1333                setexp "${key}PRCF " $value 1 5
1334            }
1335            profterms-get {
1336                set val [string range [readexp "${key}PRCF "] 5 9]
1337                if {$val == " "} {return 0}
1338                return $val
1339            }
1340            profterms-set {
1341                if ![validint value 5] {return 0}
1342                setexp "${key}PRCF " $value 6 5
1343                # now check that all needed entries exist
1344                set lines [expr 1 + ($value - 1) / 4]
1345                for {set i 1} {$i <= $lines} {incr i} {
1346                    makeexprec "${key}PRCF $i"
1347                }
1348            }
1349            pcut-get {
1350                return [string trim [string range [readexp "${key}PRCF "] 10 19]]
1351            }
1352            pcut-set {
1353                if ![validreal value 10 5] {return 0}
1354                setexp "${key}PRCF " $value 11 10
1355            }
1356            pdamp-get {
1357                set val [string range [readexp "${key}PRCF "] 24 24]
1358                if {$val == " "} {return 0}
1359                return $val
1360            }
1361            pdamp-set {
1362                setexp "${key}PRCF   " $value 25 1
1363            }
1364            pterm*-get {
1365                regsub pterm $parm {} num
1366                set f1 [expr 15*(($num - 1) % 4)]
1367                set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1368                set line  [expr 1 + ($num - 1) / 4]
1369                return [string trim [string range [readexp "${key}PRCF $line"] $f1 $f2] ]
1370            }
1371            pterm*-set {
1372                if ![validreal value 15 6] {return 0}
1373                regsub pterm $parm {} num
1374                set f1 [expr 1+ 15*(($num - 1) % 4)]
1375                set line  [expr 1 + ($num - 1) / 4]
1376                setexp "${key}PRCF $line" $value $f1 15
1377            }
1378            pref*-get {
1379                regsub pref $parm {} num
1380                set f [expr 24+$num]
1381                if {[string toupper [string range [readexp "${key}PRCF  "] $f $f]] == "Y"} {
1382                    return 1
1383                }
1384                return 0
1385            }
1386            pref*-set {
1387                regsub pref $parm {} num
1388                set f [expr 25+$num]
1389                if $value {
1390                    setexp ${key}PRCF "Y" $f 1
1391                } else {
1392                    setexp ${key}PRCF "N" $f 1
1393                }           
1394            }
1395            POnaxis-get {
1396                set val [string trim \
1397                        [string range [readexp "${key}NAXIS"] 0 4]]
1398                if {$val == ""} {return 0}
1399                return $val
1400            }
1401            POnaxis-set {
1402                if ![validint value 5] {return 0}
1403                # there should be a NAXIS record, but if not make one
1404                if {![existsexp "${key}NAXIS"]} {
1405                    makeexprec "${key}NAXIS"
1406                }
1407                setexp "${key}NAXIS  " $value 1 5
1408            }
1409            default {
1410                set msg "Unsupported hapinfo access: parm=$parm action=$action"
1411                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1412            }
1413        }
1414    }
1415    return 1
1416}
1417
1418#  get a logical constraint
1419#
1420#  type action
1421#  -----------
1422#  atom get  number        returns a list of constraints.
1423#   "   set  number value  replaces a list of constraints
1424#                          (value is a list of constraints)
1425#   "   add  number value  inserts a new list of constraints
1426#                          (number is ignored)
1427#   "   delete number      deletes a set of constraint entries
1428# Each item in the list of constraints is composed of 4 items:
1429#              phase, atom, variable, multiplier
1430# If variable=UISO atom can be ALL, otherwise atom is a number
1431# legal variable names: FRAC, X, Y, Z, UISO, U11, U22, U33, U12, U23, U13,
1432#                       MX, MY, MZ
1433#
1434#  type action
1435#  -----------
1436#  profileXX get number         returns a list of constraints for term XX=1-36
1437#                               use number=0 to get # of defined
1438#                                  constraints for term XX
1439#   "        set number value   replaces a list of constraints
1440#                               (value is a list of constraints)
1441#   "        add number value   inserts a new list of constraints
1442#                               (number is ignored)
1443#   "        delete number      deletes a set of constraint entries
1444# Each item in the list of constraints is composed of 3 items:
1445#              phase-list, histogram-list, multiplier
1446# Note that phase-list and/or histogram-list can be ALL
1447
1448proc constrinfo {type action number "value {}"} {
1449    switch -glob ${type}-$action {
1450        atom-get {
1451            # does this constraint exist?
1452            set key [format "LNCN%4d%4d" $number 1]
1453            if {![existsexp $key]} {return -1}
1454            set clist {}
1455            for {set i 1} {$i < 999} {incr i} {
1456                set key [format "LNCN%4d%4d" $number $i]
1457                if {![existsexp $key]} break
1458                set line [readexp $key]
1459                set j1 2
1460                set j2 17
1461                set seg [string range $line $j1 $j2]
1462                while {[string trim $seg] != ""} {
1463                    lappend clist [list \
1464                            [string range $seg 0 0] \
1465                            [string trim [string range $seg 1 3]] \
1466                            [string trim [string range $seg 4 7]] \
1467                            [string trim [string range $seg 8 end]]]
1468                    incr j1 16
1469                    incr j2 16
1470                    set seg [string range $line $j1 $j2]
1471                }
1472            }
1473            return $clist
1474        }
1475        atom-set {
1476            # delete records for current constraint
1477            for {set i 1} {$i < 999} {incr i} {
1478                set key [format "LNCN%4d%4d" $number $i]
1479                if {![existsexp $key]} break
1480                delexp $key
1481            }
1482            set line {}
1483            set i 1
1484            foreach tuple $value {
1485                if {[string toupper [lindex $tuple 1]] == "ALL"} {
1486                    set seg [format %1dALL%-4s%8.4f \
1487                            [lindex $tuple 0] \
1488                            [lindex $tuple 2] \
1489                            [lindex $tuple 3]]
1490                } else {
1491                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
1492                }
1493                append line $seg
1494                if {[string length $line] > 50} {
1495                    set key  [format "LNCN%4d%4d" $number $i]
1496                    makeexprec $key
1497                    setexp $key $line 3 68
1498                    set line {}
1499                    incr i
1500                }
1501            }
1502            if {$line != ""} {
1503                set key  [format "LNCN%4d%4d" $number $i]
1504                makeexprec $key
1505                setexp $key $line 3 68
1506            }
1507            return
1508        }
1509        atom-add {
1510            # loop over defined constraints
1511            for {set j 1} {$j < 9999} {incr j} {
1512                set key [format "LNCN%4d%4d" $j 1]
1513                if {![existsexp $key]} break
1514            }
1515            set number $j
1516            # save the constraint
1517            set line {}
1518            set i 1
1519            foreach tuple $value {
1520                if {[string toupper [lindex $tuple 1]] == "ALL"} {
1521                    set seg [format %1dALL%-4s%8.4f \
1522                            [lindex $tuple 0] \
1523                            [lindex $tuple 2] \
1524                            [lindex $tuple 3]]
1525                } else {
1526                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
1527                }
1528                append line $seg
1529                if {[string length $line] > 50} {
1530                    set key  [format "LNCN%4d%4d" $number $i]
1531                    makeexprec $key
1532                    setexp $key $line 3 68
1533                    set line {}
1534                    incr i
1535                }
1536            }
1537            if {$line != ""} {
1538                set key  [format "LNCN%4d%4d" $number $i]
1539                makeexprec $key
1540                setexp $key $line 3 68
1541            }
1542            return
1543        }
1544        atom-delete {
1545            for {set j $number} {$j < 9999} {incr j} {
1546                # delete records for current constraint
1547                for {set i 1} {$i < 999} {incr i} {
1548                    set key [format "LNCN%4d%4d" $j $i]
1549                    if {![existsexp $key]} break
1550                    delexp $key
1551                }
1552                # now copy records, from the next entry, if any
1553                set j1 $j
1554                incr j1
1555                set key1 [format "LNCN%4d%4d" $j1 1]
1556                # if there is no record, there is nothing to copy -- done
1557                if {![existsexp $key1]} return
1558                for {set i 1} {$i < 999} {incr i} {
1559                    set key1 [format "LNCN%4d%4d" $j1 $i]
1560                    if {![existsexp $key1]} break
1561                    set key  [format "LNCN%4d%4d" $j  $i]
1562                    makeexprec $key
1563                    setexp $key [readexp $key1] 1 68
1564                }
1565            }
1566        }
1567        profile*-delete {
1568            regsub profile $type {} term
1569            if {$term < 10} {
1570                set term " $term"
1571            }
1572            set key "LEQV PF$term   "
1573            # return nothing if no term exists
1574            if {![existsexp $key]} {return 0}
1575
1576            # number of constraint terms
1577            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1578            # don't delete a non-existing entry
1579            if {$number > $nterms} {return 0}
1580            set val [expr $nterms - 1]
1581            validint val 5
1582            setexp $key $val 1 5
1583            for {set i1 $number} {$i1 < $nterms} {incr i1} {
1584                set i2 [expr 1 + $i1]
1585                # move the contents of constraint #i2 -> i1
1586                if {$i1 > 9} {
1587                    set k1 [expr ($i1+1)/10]
1588                    set l1 $i1
1589                } else {
1590                    set k1 " "
1591                    set l1 " $i1"
1592                }
1593                set key1 "LEQV PF$term  $k1"
1594                # number of constraint lines for #i1
1595                set n1 [string trim [string range [readexp ${key1}] \
1596                        [expr ($i1%10)*5] [expr 4+(($i1%10)*5)]] ]
1597                if {$i2 > 9} {
1598                    set k2 [expr ($i2+1)/10]
1599                    set l2 $i2
1600                } else {
1601                    set k2 " "
1602                    set l2 " $i2"
1603                }
1604                set key2 "LEQV PF$term  $k2"
1605                # number of constraint lines for #i2
1606                set n2 [string trim [string range [readexp ${key2}] \
1607                        [expr ($i2%10)*5] [expr 4+(($i2%10)*5)]] ]
1608                set val $n2
1609                validint val 5
1610                # move the # of terms
1611                setexp $key1 $val [expr 1+(($i1%10)*5)] 5
1612                # move the terms
1613                for {set j 1} {$j <= $n2} {incr j 1} {
1614                    set key "LEQV PF${term}${l1}$j"
1615                    makeexprec $key
1616                    setexp $key [readexp "LEQV PF${term}${l2}$j"] 1 68
1617                }
1618                # delete any remaining lines
1619                for {set j [expr $n2+1]} {$j <= $n1} {incr j 1} {
1620                    delexp "LEQV PF${term}${l1}$j"
1621                }
1622            }
1623
1624            # clear the last term
1625            if {$nterms > 9} {
1626                set i [expr ($nterms+1)/10]
1627            } else {
1628                set i " "
1629            }
1630            set key "LEQV PF$term  $i"
1631            set cb [expr ($nterms%10)*5]
1632            set ce [expr 4+(($nterms%10)*5)]
1633            set n2 [string trim [string range [readexp ${key}] $cb $ce] ]
1634            incr cb
1635            setexp $key "     " $cb 5
1636            # delete any remaining lines
1637            for {set j 1} {$j <= $n2} {incr j 1} {
1638                delexp "LEQV PF${term}${nterms}$j"
1639            }
1640        }
1641        profile*-set {
1642            regsub profile $type {} term
1643            if {$term < 10} {
1644                set term " $term"
1645            }
1646            set key "LEQV PF$term   "
1647            # get number of constraint terms
1648            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1649            # don't change a non-existing entry
1650            if {$number > $nterms} {return 0}
1651            if {$number > 9} {
1652                set k1 [expr ($number+1)/10]
1653                set l1 $number
1654            } else {
1655                set k1 " "
1656                set l1 " $number"
1657            }
1658            set key1 "LEQV PF$term  $k1"
1659            # old number of constraint lines
1660            set n1 [string trim [string range [readexp ${key1}] \
1661                    [expr ($number%10)*5] [expr 4+(($number%10)*5)]] ]
1662            # number of new constraints
1663            set j2 [llength $value]
1664            # number of new constraint lines
1665            set val [set n2 [expr ($j2 + 2)/3]]
1666            # store the new # of lines
1667            validint val 5
1668            setexp $key1 $val [expr 1+(($number%10)*5)] 5
1669
1670            # loop over the # of lines in the old or new, whichever is greater
1671            set v0 0
1672            for {set j 1} {$j <= [expr ($n1 > $n2) ? $n1 : $n2]} {incr j 1} {
1673                set key "LEQV PF${term}${l1}$j"
1674                # were there more lines in the old?
1675                if {$j > $n2} {
1676                    # this line is not needed
1677                    if {$j % 3 == 1} {
1678                        delexp %key
1679                    }
1680                    continue
1681                }
1682                # are we adding new lines?
1683                if {$j > $n1} {
1684                    makeexprec $key
1685                }
1686                # add the three constraints to the line
1687                foreach s {3 23 43} \
1688                        item [lrange $value $v0 [expr 2+$v0]] {
1689                    if {$item != ""} {
1690                        set val [format %-10s%9.3f \
1691                                [lindex $item 0],[lindex $item 1] \
1692                                [lindex $item 2]]
1693                        setexp $key $val $s 19
1694                    } else {
1695                        setexp $key " " $s 19
1696                    }
1697                }
1698                incr v0 3
1699            }
1700        }
1701        profile*-add {
1702            regsub profile $type {} term
1703            if {$term < 10} {
1704                set term " $term"
1705            }
1706            set key "LEQV PF$term   "
1707            if {![existsexp $key]} {makeexprec $key}
1708            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1709            if {$nterms == ""} {
1710                set nterms 1
1711            } elseif {$nterms >= 99} {
1712                return 0
1713            } else {
1714                incr nterms
1715            }
1716            # store the new # of constraints
1717            set val $nterms
1718            validint val 5
1719            setexp $key $val 1 5
1720
1721            if {$nterms > 9} {
1722                set k1 [expr ($nterms+1)/10]
1723                set l1 $nterms
1724            } else {
1725                set k1 " "
1726                set l1 " $nterms"
1727            }
1728            set key1 "LEQV PF$term  $k1"
1729
1730            # number of new constraints
1731            set j2 [llength $value]
1732            # number of new constraint lines
1733            set val [set n2 [expr ($j2 + 2)/3]]
1734            # store the new # of lines
1735            validint val 5
1736            setexp $key1 $val [expr 1+(($nterms%10)*5)] 5
1737
1738            # loop over the # of lines to be added
1739            set v0 0
1740            for {set j 1} {$j <= $n2} {incr j 1} {
1741                set key "LEQV PF${term}${l1}$j"
1742                makeexprec $key
1743                # add the three constraints to the line
1744                foreach s {3 23 43} \
1745                        item [lrange $value $v0 [expr 2+$v0]] {
1746                    if {$item != ""} {
1747                        set val [format %-10s%9.3f \
1748                                [lindex $item 0],[lindex $item 1] \
1749                                [lindex $item 2]]
1750                        setexp $key $val $s 19
1751                    } else {
1752                        setexp $key " " $s 19
1753                    }
1754                }
1755                incr v0 3
1756            }
1757        }
1758        profile*-get {
1759            regsub profile $type {} term
1760            if {$term < 10} {
1761                set term " $term"
1762            }
1763            if {$number > 9} {
1764                set i [expr ($number+1)/10]
1765            } else {
1766                set i " "
1767            }
1768            set key "LEQV PF$term  $i"
1769            # return nothing if no term exists
1770            if {![existsexp $key]} {return 0}
1771            # number of constraint lines
1772           
1773            set numline [string trim [string range [readexp ${key}] \
1774                    [expr ($number%10)*5] [expr 4+(($number%10)*5)]] ]
1775            if {$number == 0} {return $numline}
1776            set clist {}
1777            if {$number < 10} {
1778                set number " $number"
1779            }
1780            for {set i 1} {$i <= $numline} {incr i} {
1781                set key "LEQV PF${term}${number}$i"
1782                set line [readexp ${key}]
1783                foreach s {1 21 41} e {20 40 60} {
1784                    set seg [string range $line $s $e]
1785                    if {[string trim $seg] == ""} continue
1786                    # parse the string segment
1787                    set parse [regexp { *([0-9AL]+),([0-9AL]+) +([0-9.]+)} \
1788                            $seg junk phase hist mult]
1789                    # was parse successful
1790                    if {!$parse} {continue}
1791                    lappend clist [list $phase $hist $mult]
1792                }
1793            }
1794            return $clist
1795        }
1796        default {
1797            set msg "Unsupported constrinfo access: type=$type action=$action"
1798            tk_dialog .badexp "Error in readexp access" $msg error 0 OK
1799        }
1800
1801    }
1802}
1803
1804# read the default profile information for a histogram
1805# use: profdefinfo hist set# parm action
1806# 
1807#     proftype -- profile function number
1808#     profterms -- number of profile terms
1809#     pdamp -- damping value for the profile (*)
1810#     pcut -- cutoff value for the profile (*)
1811#     pterm$n -- profile term #n
1812#     pref$n -- refinement flag value for profile term #n (*)
1813
1814proc profdefinfo {hist set parm "action get"} {
1815    global expgui
1816    if {$hist < 10} {
1817        set key "HST  $hist"
1818    } else {
1819        set key "HST $hist"
1820    }
1821    switch -glob ${parm}-$action {
1822        proftype-get {
1823            set val [string range [readexp "${key}PRCF$set"] 0 4]
1824            if {$val == " "} {return 0}
1825            return $val
1826        }
1827        profterms-get {
1828            set val [string range [readexp "${key}PRCF$set"] 5 9]
1829            if {$val == " "} {return 0}
1830            return $val
1831        }
1832        pcut-get {
1833            return [string trim [string range [readexp "${key}PRCF$set"] 10 19]]
1834        }
1835        pdamp-get {
1836                set val [string range [readexp "${key}PRCF$set"] 24 24]
1837            if {$val == " "} {return 0}
1838            return $val
1839        }
1840        pterm*-get {
1841            regsub pterm $parm {} num
1842            set f1 [expr 15*(($num - 1) % 4)]
1843            set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1844            set line  [expr 1 + ($num - 1) / 4]
1845            return [string trim [string range [\
1846                        readexp "${key}PRCF${set}$line"] $f1 $f2] ]
1847        }
1848        pref*-get {
1849            regsub pref $parm {} num
1850            set f [expr 24+$num]
1851            if {[string toupper [string range [readexp "${key}PRCF$set"] $f $f]] == "Y"} {
1852                return 1
1853            }
1854            return 0
1855        }
1856        default {
1857            set msg "Unsupported profdefinfo access: parm=$parm action=$action"
1858            tk_dialog .badexp "Code Error" $msg error 0 Exit
1859        }
1860    }
1861}
1862
1863# get March-Dollase preferred orientation information
1864# use MDprefinfo hist phase axis-number parm action value
1865#    ratio    -- ratio of xtallites in PO direction vs random (>1 for more)
1866#    fraction -- fraction in this direction, when more than one axis is used
1867#    h k & l  -- indices of P.O. axis
1868#    ratioref -- flag to vary ratio
1869#    fracref  -- flag to vary fraction
1870#    damp     -- damping value
1871#    type     -- model type (0 = P.O. _|_ to beam, 1 = || to beam)
1872#    new      -- creates a new record with default values (set only)
1873proc MDprefinfo {histlist phaselist axislist parm "action get" "value {}"} {
1874    foreach phase $phaselist hist $histlist axis $axislist {
1875        if {$phase == ""} {set phase [lindex $phaselist end]}
1876        if {$hist == ""} {set hist [lindex $histlist end]}
1877        if {$axis == ""} {set axis [lindex $axislist end]}
1878        if {$hist < 10} {
1879            set hist " $hist"
1880        }
1881        if {$axis > 9} {
1882            set axis "0"
1883        }
1884        set key "HAP${phase}${hist}PREFO${axis}"
1885        switch -glob ${parm}-$action {
1886            ratio-get {
1887                return [string trim [string range [readexp $key] 0 9]]
1888            }
1889            ratio-set {
1890                if ![validreal value 10 6] {return 0}
1891                setexp $key $value 1 10
1892            }
1893            fraction-get {
1894                return [string trim [string range [readexp $key] 10 19]]
1895            }
1896            fraction-set {
1897                if ![validreal value 10 6] {return 0}
1898                setexp $key $value 11 10
1899            }
1900            h-get {
1901                set h [string trim [string range [readexp $key] 20 29]]
1902                # why not allow negative h values?
1903                #               if {$h < 1} {return 0}
1904                return $h
1905            }
1906            h-set {
1907                if ![validreal value 10 2] {return 0}
1908                setexp $key $value 21 10
1909            }
1910            k-get {
1911                set k [string trim [string range [readexp $key] 30 39]]
1912                #               if {$k < 1} {return 0}
1913                return $k
1914            }
1915            k-set {
1916                if ![validreal value 10 2] {return 0}
1917                setexp $key $value 31 10
1918            }
1919            l-get {
1920                set l [string trim [string range [readexp $key] 40 49]]
1921                #if {$l < 1} {return 0}
1922                return $l
1923            }
1924            l-set {
1925                if ![validreal value 10 2] {return 0}
1926                setexp $key $value 41 10
1927            }
1928            ratioref-get {
1929                if {[string toupper \
1930                        [string range [readexp $key] 53 53]] == "Y"} {
1931                    return 1
1932                }
1933                return 0
1934            }
1935            ratioref-set {
1936                if $value {
1937                    setexp $key "Y" 54 1
1938                } else {
1939                    setexp $key "N" 54 1
1940                }
1941            }
1942            fracref-get {
1943                if {[string toupper \
1944                        [string range [readexp $key] 54 54]] == "Y"} {
1945                    return 1
1946                }
1947                return 0
1948            }
1949            fracref-set {
1950                if $value {
1951                    setexp $key "Y" 55 1
1952                } else {
1953                    setexp $key "N" 55 1
1954              }
1955            }
1956            damp-get {
1957                set val [string trim [string range [readexp $key] 59 59]]
1958                if {$val == " "} {return 0}
1959                return $val
1960            }
1961            damp-set {
1962                setexp $key $value 60 1
1963            }
1964            type-get {
1965                set val [string trim [string range [readexp $key] 64 64]]
1966                if {$val == " "} {return 0}
1967                return $val
1968            }
1969            type-set {
1970                # only valid settings are 0 & 1
1971                if {$value != "0" && $value != "1"} {set value "0"}
1972                setexp $key $value 65 1
1973            }
1974            new-set {
1975                makeexprec $key
1976                setexp $key \
1977                        {  1.000000  1.000000  0.000000  0.000000  1.000000   NN    0    0} \
1978                        1 68
1979            }
1980            default {
1981                set msg "Unsupported MDprefinfo access: parm=$parm action=$action"
1982                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1983            }
1984
1985        }
1986
1987    }
1988}
1989
1990# write the .EXP file
1991proc expwrite {expfile} {
1992    global tcl_platform exparray
1993    set blankline \
1994     "                                                                        "
1995    set fp [open ${expfile} w]
1996    set keylist [lsort [array names exparray]]
1997    # reorder the keys so that VERSION comes 1st
1998    set pos [lsearch -exact $keylist {     VERSION}]
1999    set keylist "{     VERSION} [lreplace $keylist $pos $pos]"
2000    if {$tcl_platform(platform) == "windows"} { 
2001        foreach key $keylist {
2002            puts $fp [string range \
2003                    "$key$exparray($key)$blankline" 0 79]
2004        }
2005    } else {
2006        foreach key $keylist {
2007            puts -nonewline $fp [string range \
2008                    "$key$exparray($key)$blankline" 0 79]
2009        }
2010    }
2011    close $fp
2012}
2013
2014# history commands -- delete all but last $keep history records,
2015# renumber if $renumber is true
2016proc DeleteHistory {keep renumber} {
2017    global exparray
2018    foreach y [lrange [lsort -decreasing \
2019            [array names exparray {    HSTRY*}]] $keep end] {
2020        unset exparray($y)
2021    }
2022    if !$renumber return
2023    # renumber
2024    set i 0
2025    foreach y [lsort -increasing \
2026            [array names exparray {    HSTRY*}]] {
2027        set key [format "    HSTRY%3d" [incr i]]
2028        set exparray($key) $exparray($y)
2029        unset exparray($y)
2030    }
2031    # list all history
2032    #    foreach y [lsort -decreasing [array names exparray {    HSTRY*}]] {puts "$y $exparray($y)"}
2033}
2034
2035proc CountHistory {} {
2036    global exparray
2037    return [llength [array names exparray {    HSTRY*}]]
2038}
2039
2040# set the phase flags for histogram $hist to $plist
2041proc SetPhaseFlag {hist plist} {
2042    # make a 2 digit key -- hh
2043    if {$hist < 10} {
2044        set hh " $hist"
2045    } else {
2046        set hh $hist
2047    }
2048    set key "HST $hh NPHAS"
2049    set str {}
2050    foreach iph {1 2 3 4 5 6 7 8 9} {
2051        if {[lsearch $plist $iph] != -1} {
2052            append str {    1}
2053        } else {
2054            append str {    0}     
2055        }
2056    }
2057    setexp $key $str 1 68
2058}
2059
2060# erase atom $atom from phase $phase
2061# update the list of atom types, erasing the record if not needed.
2062proc EraseAtom {atom phase} {
2063    set type [atominfo $phase $atom type]
2064    if {$type == ""} return
2065    if {$atom < 10} {
2066        set key "CRS$phase  AT  $atom"
2067    } elseif {$atom < 100} {
2068        set key "CRS$phase  AT $atom"
2069    } else {
2070        set key "CRS$phase  AT$atom"
2071    }
2072    # delete the records for the atom
2073    global exparray
2074    foreach k [array names exparray ${key}*] {
2075        delexp $k
2076    }
2077    # change the number of atoms in the phase
2078    phaseinfo $phase natoms set [expr [phaseinfo $phase natoms] -1]
2079
2080    # now adjust numbers in "EXPR ATYP" records and delete, if needed.
2081    set natypes [readexp " EXPR  NATYP"]
2082    if {$natypes == ""} return
2083    set j 0
2084    for {set i 1} {$i <= $natypes} {incr i} {
2085        incr j
2086        if {$j <10} {
2087            set key " EXPR ATYP $j"
2088        } else {
2089            set key " EXPR ATYP$j"
2090        }
2091        while {![existsexp $key]} {
2092            incr j
2093            if {$j > 99} {
2094                return
2095            } elseif {$j <10} {
2096                set key " EXPR ATYP $j"
2097            } else {
2098                set key " EXPR ATYP$j"
2099            }
2100        }
2101        set keytype [string trim [string range $exparray($key) 2 9]]
2102        if {$type == $keytype} {
2103            # found the type record
2104            set val [string trim [string range $exparray($key) 10 14]]
2105            incr val -1
2106            # if this is the last reference, remove the record,
2107            # otherwise, decrement the counter
2108            if {$val <= 0} {
2109                incr natypes -1 
2110                validint natypes 5
2111                setexp " EXPR  NATYP" $natypes 1 5
2112                delexp $key
2113            } else {
2114                validint val 5
2115                setexp $key $val 11 5
2116            }
2117            return
2118        }
2119    }
2120}
2121
2122# compute equivalent anisotropic temperature factor for Uequiv
2123proc CalcAniso {phase Uequiv} {
2124    foreach var {a b c alpha beta gamma} {
2125        set $var [phaseinfo $phase $var]
2126    }
2127
2128    set G(1,1) [expr $a * $a]
2129    set G(2,2) [expr $b * $b]
2130    set G(3,3) [expr $c * $c]
2131    set G(1,2) [expr $a * $b * cos($gamma*0.017453292519943)]
2132    set G(2,1) $G(1,2)
2133    set G(1,3) [expr $a * $c * cos($beta *0.017453292519943)]
2134    set G(3,1) $G(1,3)
2135    set G(2,3) [expr $b * $c * cos($alpha*0.017453292519943)]
2136    set G(3,2) $G(2,3)
2137
2138    # Calculate the volume**2
2139    set v2 0.0
2140    foreach i {1 2 3} {
2141        set J [expr ($i%3) + 1]
2142        set K [expr (($i+1)%3) + 1]
2143        set v2 [expr $v2+ $G(1,$i)*($G(2,$J)*$G(3,$K)-$G(3,$J)*$G(2,$K))]
2144    }
2145    if {$v2 > 0} {
2146        set v [expr sqrt($v2)]
2147        foreach i {1 2 3} {
2148            set i1 [expr ($i%3) + 1]
2149            set i2 [expr (($i+1)%3) + 1]
2150            foreach j {1 2 3} {
2151                set j1 [expr ($j%3) + 1]
2152                set j2 [expr (($j+1)%3) + 1]
2153                set C($j,$i) [expr (\
2154                        $G($i1,$j1) * $G($i2,$j2) - \
2155                        $G($i1,$j2)  * $G($i2,$j1)\
2156                        )/ $v]
2157            }
2158        }
2159        set A(1,2) [expr 0.5 * ($C(1,2) + $C(2,1)) / sqrt( $C(1,1)* $C(2,2) )]
2160        set A(1,3) [expr 0.5 * ($C(1,3) + $C(3,1)) / sqrt( $C(1,1)* $C(3,3) )]
2161        set A(2,3) [expr 0.5 * ($C(2,3) + $C(3,2)) / sqrt( $C(2,2)* $C(3,3) )]
2162        foreach i {1 1 2} j {2 3 3} {
2163            set A($i,$j) [expr 0.5 * ($C($i,$j) + $C($j,$i)) / \
2164                    sqrt( $C($i,$i)* $C($j,$j) )]
2165            # clean up roundoff
2166            if {abs($A($i,$j)) < 1e-5} {set A($i,$j) 0.0}
2167        }
2168    } else {
2169        set A(1,2) 0.0
2170        set A(1,3) 0.0
2171        set A(2,3) 0.0
2172    }
2173    return "$Uequiv $Uequiv $Uequiv \
2174            [expr $Uequiv * $A(1,2)] \
2175            [expr $Uequiv * $A(1,3)] \
2176            [expr $Uequiv * $A(2,3)]"
2177}
Note: See TracBrowser for help on using the repository browser.