source: trunk/readexp.tcl @ 373

Last change on this file since 373 was 373, checked in by toby, 11 years ago

# on 2001/02/26 19:37:17, toby did:
Change EPHAS (Fobs extraction) to Y/N except on IRIX

  • Property rcs:author set to toby
  • Property rcs:date set to 2001/02/26 19:37:17
  • Property rcs:lines set to +23 -7
  • Property rcs:rev set to 1.23
  • Property rcs:state set to Exp
  • Property svn:keywords set to Author Date Revision Id
File size: 58.1 KB
Line 
1# $Id: readexp.tcl 373 2009-12-04 23:05:05Z toby $
2# Routines to deal with the .EXP "data structure"
3set expmap(Revision) {$Revision: 373 $ $Date: 2009-12-04 23:05:05 +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    if [catch {
255        set tmp [expr round($value)]
256        if {$tmp != $value} {return 0}
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" "Y" 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                # the flag has changed to "Y/N" in the latest versions
971                # of GSAS
972                if {[string toupper [string range $line 49 49]] == "Y"} {
973                    return 1
974                }
975                return 0
976            }
977            foextract-set {
978                global tcl_platform
979                if {[string match "IRIX*" $tcl_platform(os)]} {
980                    if $value {
981                        setexp "${key} EPHAS" "T" 50 1
982                    } else {
983                        setexp "${key} EPHAS" "F" 50 1
984                    }       
985                } else {
986                # the flag has changed to "Y/N" in the latest versions
987                # of GSAS
988                    if $value {
989                        setexp "${key} EPHAS" "Y" 50 1
990                    } else {
991                        setexp "${key} EPHAS" "N" 50 1
992                    }
993                }
994            }
995            title-get {
996                return [string trim [readexp "${key}  HNAM"] ]
997            }
998            scale-get {
999                return [string trim [string range [readexp ${key}HSCALE] 0 14]]
1000            }
1001            scale-set {
1002                if ![validreal value 15 6] {return 0}
1003                setexp ${key}HSCALE $value 1 15
1004            }
1005            sref-get {
1006                if {[string toupper [string range [readexp ${key}HSCALE] 19 19]] == "Y"} {
1007                    return 1
1008                }
1009                return 0
1010            }
1011            sref-set {
1012                if $value {
1013                    setexp ${key}HSCALE "Y" 20 1
1014                } else {
1015                    setexp ${key}HSCALE "N" 20 1
1016                }           
1017            }
1018            sdamp-get {
1019                set val [string range [readexp ${key}HSCALE] 24 24]
1020                if {$val == " "} {return 0}
1021                return $val
1022            }
1023            sdamp-set {
1024                setexp ${key}HSCALE $value 25 1
1025            }
1026
1027            difc-get -
1028            lam1-get {
1029                return [string trim [string range [readexp "${key} ICONS"] 0 9]]
1030            }
1031            difc-set -
1032            lam1-set {
1033                if ![validreal value 10 7] {return 0}
1034                setexp "${key} ICONS" $value 1 10
1035            }
1036            difa-get -
1037            lam2-get {
1038                return [string trim [string range [readexp "${key} ICONS"] 10 19]]
1039            }
1040            difa-set -
1041            lam2-set {
1042                if ![validreal value 10 7] {return 0}
1043                setexp "${key} ICONS" $value 11 10
1044            }
1045            zero-get {
1046                return [string trim [string range [readexp "${key} ICONS"] 20 29]]
1047            }
1048            zero-set {
1049                if ![validreal value 10 5] {return 0}
1050                setexp "${key} ICONS" $value 21 10
1051            }
1052            ipola-get {
1053                return [string trim [string range [readexp "${key} ICONS"] 54 54]]
1054            }
1055            ipola-set {
1056                if ![validint value 1] {return 0}
1057                setexp "${key} ICONS" $value 55 1
1058            }
1059            pola-get {
1060                return [string trim [string range [readexp "${key} ICONS"] 40 49]]
1061            }
1062            pola-set {
1063                if ![validreal value 10 5] {return 0}
1064                setexp "${key} ICONS" $value 41 10
1065            }
1066            kratio-get {
1067                return [string trim [string range [readexp "${key} ICONS"] 55 64]]
1068            }
1069            kratio-set {
1070                if ![validreal value 10 5] {return 0}
1071                setexp "${key} ICONS" $value 56 10
1072            }
1073
1074            wref-get {
1075            #------------------------------------------------------
1076            # col 33: refine flag for lambda, difc, ratio and theta
1077            #------------------------------------------------------
1078                if {[string toupper [string range \
1079                        [readexp "${key} ICONS"] 32 32]] == "L"} {
1080                    return 1
1081                }
1082                return 0
1083            }
1084            wref-set {
1085                if $value {
1086                    setexp "${key} ICONS" "L" 33 1
1087                } else {
1088                    setexp "${key} ICONS" " " 33 1
1089                }           
1090            }
1091            ratref-get {
1092                if {[string toupper [string range \
1093                        [readexp "${key} ICONS"] 32 32]] == "R"} {
1094                    return 1
1095                }
1096                return 0
1097            }
1098            ratref-set {
1099                if $value {
1100                    setexp "${key} ICONS" "R" 33 1
1101                } else {
1102                    setexp "${key} ICONS" " " 33 1
1103                }           
1104            }
1105            dcref-get {
1106                if {[string toupper [string range \
1107                        [readexp "${key} ICONS"] 32 32]] == "C"} {
1108                    return 1
1109                }
1110                return 0
1111            }
1112            dcref-set {
1113                if $value {
1114                    setexp "${key} ICONS" "C" 33 1
1115                } else {
1116                    setexp "${key} ICONS" " " 33 1
1117                }           
1118            }
1119            ttref-get {
1120                if {[string toupper [string range \
1121                        [readexp "${key} ICONS"] 32 32]] == "T"} {
1122                    return 1
1123                }
1124                return 0
1125            }
1126            ttref-set {
1127                if $value {
1128                    setexp "${key} ICONS" "T" 33 1
1129                } else {
1130                    setexp "${key} ICONS" " " 33 1
1131                }           
1132            }
1133
1134
1135            pref-get {
1136            #------------------------------------------------------
1137            # col 34: refine flag for POLA & DIFA
1138            #------------------------------------------------------
1139                if {[string toupper [string range \
1140                        [readexp "${key} ICONS"] 33 33]] == "P"} {
1141                    return 1
1142                }
1143                return 0
1144            }
1145            pref-set {
1146                if $value {
1147                    setexp "${key} ICONS" "P" 34 1
1148                } else {
1149                    setexp "${key} ICONS" " " 34 1
1150                }           
1151            }
1152            daref-get {
1153                if {[string toupper [string range \
1154                        [readexp "${key} ICONS"] 33 33]] == "A"} {
1155                    return 1
1156                }
1157                return 0
1158            }
1159            daref-set {
1160                if $value {
1161                    setexp "${key} ICONS" "A" 34 1
1162                } else {
1163                    setexp "${key} ICONS" " " 34 1
1164                }           
1165            }
1166
1167            zref-get {
1168            #------------------------------------------------------
1169            # col 34: refine flag for zero correction
1170            #------------------------------------------------------
1171                if {[string toupper [string range [readexp "${key} ICONS"] 34 34]] == "Z"} {
1172                    return 1
1173                }
1174                return 0
1175            }
1176            zref-set {
1177                if $value {
1178                    setexp "${key} ICONS" "Z" 35 1
1179                } else {
1180                    setexp "${key} ICONS" " " 35 1
1181                }           
1182            }
1183
1184            ddamp-get {
1185                set val [string range [readexp "${key} ICONS"] 39 39]
1186                if {$val == " "} {return 0}
1187                return $val
1188            }
1189            ddamp-set {
1190                setexp "${key} ICONS" $value 40 1
1191            }
1192
1193            backtype-get {
1194                set val [string trim [string range [readexp "${key}BAKGD "] 0 4]]
1195                if {$val == " "} {return 0}
1196                return $val
1197            }
1198            backtype-set {
1199                if ![validint value 5] {return 0}
1200                setexp "${key}BAKGD " $value 1 5
1201            }
1202            backterms-get {
1203                set val [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1204                if {$val == " "} {return 0}
1205                return $val
1206            }
1207            backterms-set {
1208                # this takes a bit of work -- if terms are added, add lines as needed to the .EXP
1209                set oldval [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1210                if ![validint value 5] {return 0}
1211                if {$oldval < $value} {
1212                    set line1  [expr 2 + ($oldval - 1) / 4]
1213                    set line2  [expr 1 + ($value - 1) / 4]
1214                    for {set i $line1} {$i <= $line2} {incr i} {
1215                        # create a blank entry if needed
1216                        makeexprec ${key}BAKGD$i
1217                    }
1218                    incr oldval
1219                    for {set num $oldval} {$num <= $value} {incr num} {
1220                        set f1 [expr 15*(($num - 1) % 4)]
1221                        set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1222                        set line  [expr 1 + ($num - 1) / 4]
1223                        if {[string trim [string range [readexp ${key}BAKGD$line] $f1 $f2]] == ""} {
1224                            set f1 [expr 15*(($num - 1) % 4)+1]
1225                            setexp ${key}BAKGD$line 0.0 $f1 15                 
1226                        }
1227                    }
1228                }
1229                setexp "${key}BAKGD " $value 6 5
1230
1231            }
1232            bref-get {
1233                if {[string toupper [string range [readexp "${key}BAKGD"] 14 14]] == "Y"} {
1234                    return 1
1235                }
1236                return 0
1237            }
1238            bref-set {
1239                if $value {
1240                    setexp "${key}BAKGD "  "Y" 15 1
1241                } else {
1242                    setexp "${key}BAKGD "  "N" 15 1
1243                }           
1244            }
1245            bdamp-get {
1246                set val [string range [readexp "${key}BAKGD "] 19 19]
1247                if {$val == " "} {return 0}
1248                return $val
1249            }
1250            bdamp-set {
1251                setexp "${key}BAKGD " $value 20 1
1252            }
1253            bterm*-get {
1254                regsub bterm $parm {} num
1255                set f1 [expr 15*(($num - 1) % 4)]
1256                set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1257                set line  [expr 1 + ($num - 1) / 4]
1258                return [string trim [string range [readexp ${key}BAKGD$line] $f1 $f2] ]
1259            }
1260            bterm*-set {
1261                regsub bterm $parm {} num
1262                if ![validreal value 15 6] {return 0}
1263                set f1 [expr 15*(($num - 1) % 4)+1]
1264                set line  [expr 1 + ($num - 1) / 4]
1265                setexp ${key}BAKGD$line $value $f1 15
1266            }
1267            bank-get {
1268                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1269            }
1270            tofangle-get {
1271                return [string trim [string range [readexp "${key}BNKPAR"] 10 19]]
1272            }
1273            default {
1274                set msg "Unsupported histinfo access: parm=$parm action=$action"
1275                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1276            }
1277        }
1278    }
1279    return 1
1280}
1281
1282# read the information that differs by both histogram and phase (profile & phase fraction)
1283# use: hapinfo hist phase parm action value
1284
1285#     frac -- phase fraction (*)
1286#     frref/frdamp -- refinement flag/damping value for the phase fraction (*)
1287#     proftype -- profile function number (*)
1288#     profterms -- number of profile terms (*)
1289#     pdamp -- damping value for the profile (*)
1290#     pcut -- cutoff value for the profile (*)
1291#     pterm$n -- profile term #n (*)
1292#     pref$n -- refinement flag value for profile term #n (*)
1293#     extmeth -- Fobs extraction method (*)
1294#     POnaxis -- number of defined M-D preferred axes
1295proc hapinfo {histlist phaselist parm "action get" "value {}"} {
1296    foreach phase $phaselist hist $histlist {
1297        if {$phase == ""} {set phase [lindex $phaselist end]}
1298        if {$hist == ""} {set hist [lindex $histlist end]}
1299        if {$hist < 10} {
1300            set hist " $hist"
1301        }
1302        set key "HAP${phase}${hist}"
1303        switch -glob ${parm}-$action {
1304            extmeth-get {
1305                set i1 [expr ($phase - 1)*5]
1306                set i2 [expr $i1 + 4]
1307                return [string trim [string range [readexp "HST $hist EPHAS"] $i1 $i2]]
1308            }
1309            extmeth-set {
1310                set i1 [expr ($phase - 1)*5 + 1]
1311                if ![validint value 5] {return 0}
1312                setexp "HST $hist EPHAS" $value $i1 5
1313            }
1314            frac-get {
1315                return [string trim [string range [readexp ${key}PHSFR] 0 14]]
1316            }
1317            frac-set {
1318                if ![validreal value 15 6] {return 0}
1319                setexp ${key}PHSFR $value 1 15
1320            }
1321            frref-get {
1322                if {[string toupper [string range [readexp ${key}PHSFR] 19 19]] == "Y"} {
1323                    return 1
1324                }
1325                return 0
1326            }
1327            frref-set {
1328                if $value {
1329                    setexp ${key}PHSFR "Y" 20 1
1330                } else {
1331                    setexp ${key}PHSFR "N" 20 1
1332                }           
1333            }
1334            frdamp-get {
1335                set val [string range [readexp ${key}PHSFR] 24 24]
1336                if {$val == " "} {return 0}
1337                return $val
1338            }
1339            frdamp-set {
1340                setexp ${key}PHSFR $value 25 1
1341            }
1342            proftype-get {
1343                set val [string range [readexp "${key}PRCF "] 0 4]
1344                if {$val == " "} {return 0}
1345                return $val
1346            }
1347            proftype-set {
1348                if ![validint value 5] {return 0}
1349                setexp "${key}PRCF " $value 1 5
1350            }
1351            profterms-get {
1352                set val [string range [readexp "${key}PRCF "] 5 9]
1353                if {$val == " "} {return 0}
1354                return $val
1355            }
1356            profterms-set {
1357                if ![validint value 5] {return 0}
1358                setexp "${key}PRCF " $value 6 5
1359                # now check that all needed entries exist
1360                set lines [expr 1 + ($value - 1) / 4]
1361                for {set i 1} {$i <= $lines} {incr i} {
1362                    makeexprec "${key}PRCF $i"
1363                }
1364            }
1365            pcut-get {
1366                return [string trim [string range [readexp "${key}PRCF "] 10 19]]
1367            }
1368            pcut-set {
1369                if ![validreal value 10 5] {return 0}
1370                setexp "${key}PRCF " $value 11 10
1371            }
1372            pdamp-get {
1373                set val [string range [readexp "${key}PRCF "] 24 24]
1374                if {$val == " "} {return 0}
1375                return $val
1376            }
1377            pdamp-set {
1378                setexp "${key}PRCF   " $value 25 1
1379            }
1380            pterm*-get {
1381                regsub pterm $parm {} num
1382                set f1 [expr 15*(($num - 1) % 4)]
1383                set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1384                set line  [expr 1 + ($num - 1) / 4]
1385                return [string trim [string range [readexp "${key}PRCF $line"] $f1 $f2] ]
1386            }
1387            pterm*-set {
1388                if ![validreal value 15 6] {return 0}
1389                regsub pterm $parm {} num
1390                set f1 [expr 1+ 15*(($num - 1) % 4)]
1391                set line  [expr 1 + ($num - 1) / 4]
1392                setexp "${key}PRCF $line" $value $f1 15
1393            }
1394            pref*-get {
1395                regsub pref $parm {} num
1396                set f [expr 24+$num]
1397                if {[string toupper [string range [readexp "${key}PRCF  "] $f $f]] == "Y"} {
1398                    return 1
1399                }
1400                return 0
1401            }
1402            pref*-set {
1403                regsub pref $parm {} num
1404                set f [expr 25+$num]
1405                if $value {
1406                    setexp ${key}PRCF "Y" $f 1
1407                } else {
1408                    setexp ${key}PRCF "N" $f 1
1409                }           
1410            }
1411            POnaxis-get {
1412                set val [string trim \
1413                        [string range [readexp "${key}NAXIS"] 0 4]]
1414                if {$val == ""} {return 0}
1415                return $val
1416            }
1417            POnaxis-set {
1418                if ![validint value 5] {return 0}
1419                # there should be a NAXIS record, but if not make one
1420                if {![existsexp "${key}NAXIS"]} {
1421                    makeexprec "${key}NAXIS"
1422                }
1423                setexp "${key}NAXIS  " $value 1 5
1424            }
1425            default {
1426                set msg "Unsupported hapinfo access: parm=$parm action=$action"
1427                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1428            }
1429        }
1430    }
1431    return 1
1432}
1433
1434#  get a logical constraint
1435#
1436#  type action
1437#  -----------
1438#  atom get  number        returns a list of constraints.
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 4 items:
1445#              phase, atom, variable, multiplier
1446# If variable=UISO atom can be ALL, otherwise atom is a number
1447# legal variable names: FRAC, X, Y, Z, UISO, U11, U22, U33, U12, U23, U13,
1448#                       MX, MY, MZ
1449#
1450#  type action
1451#  -----------
1452#  profileXX get number         returns a list of constraints for term XX=1-36
1453#                               use number=0 to get # of defined
1454#                                  constraints for term XX
1455#   "        set number value   replaces a list of constraints
1456#                               (value is a list of constraints)
1457#   "        add number value   inserts a new list of constraints
1458#                               (number is ignored)
1459#   "        delete number      deletes a set of constraint entries
1460# Each item in the list of constraints is composed of 3 items:
1461#              phase-list, histogram-list, multiplier
1462# Note that phase-list and/or histogram-list can be ALL
1463
1464proc constrinfo {type action number "value {}"} {
1465    switch -glob ${type}-$action {
1466        atom-get {
1467            # does this constraint exist?
1468            set key [format "LNCN%4d%4d" $number 1]
1469            if {![existsexp $key]} {return -1}
1470            set clist {}
1471            for {set i 1} {$i < 999} {incr i} {
1472                set key [format "LNCN%4d%4d" $number $i]
1473                if {![existsexp $key]} break
1474                set line [readexp $key]
1475                set j1 2
1476                set j2 17
1477                set seg [string range $line $j1 $j2]
1478                while {[string trim $seg] != ""} {
1479                    lappend clist [list \
1480                            [string range $seg 0 0] \
1481                            [string trim [string range $seg 1 3]] \
1482                            [string trim [string range $seg 4 7]] \
1483                            [string trim [string range $seg 8 end]]]
1484                    incr j1 16
1485                    incr j2 16
1486                    set seg [string range $line $j1 $j2]
1487                }
1488            }
1489            return $clist
1490        }
1491        atom-set {
1492            # delete records for current constraint
1493            for {set i 1} {$i < 999} {incr i} {
1494                set key [format "LNCN%4d%4d" $number $i]
1495                if {![existsexp $key]} break
1496                delexp $key
1497            }
1498            set line {}
1499            set i 1
1500            foreach tuple $value {
1501                if {[string toupper [lindex $tuple 1]] == "ALL"} {
1502                    set seg [format %1dALL%-4s%8.4f \
1503                            [lindex $tuple 0] \
1504                            [lindex $tuple 2] \
1505                            [lindex $tuple 3]]
1506                } else {
1507                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
1508                }
1509                append line $seg
1510                if {[string length $line] > 50} {
1511                    set key  [format "LNCN%4d%4d" $number $i]
1512                    makeexprec $key
1513                    setexp $key $line 3 68
1514                    set line {}
1515                    incr i
1516                }
1517            }
1518            if {$line != ""} {
1519                set key  [format "LNCN%4d%4d" $number $i]
1520                makeexprec $key
1521                setexp $key $line 3 68
1522            }
1523            return
1524        }
1525        atom-add {
1526            # loop over defined constraints
1527            for {set j 1} {$j < 9999} {incr j} {
1528                set key [format "LNCN%4d%4d" $j 1]
1529                if {![existsexp $key]} break
1530            }
1531            set number $j
1532            # save the constraint
1533            set line {}
1534            set i 1
1535            foreach tuple $value {
1536                if {[string toupper [lindex $tuple 1]] == "ALL"} {
1537                    set seg [format %1dALL%-4s%8.4f \
1538                            [lindex $tuple 0] \
1539                            [lindex $tuple 2] \
1540                            [lindex $tuple 3]]
1541                } else {
1542                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
1543                }
1544                append line $seg
1545                if {[string length $line] > 50} {
1546                    set key  [format "LNCN%4d%4d" $number $i]
1547                    makeexprec $key
1548                    setexp $key $line 3 68
1549                    set line {}
1550                    incr i
1551                }
1552            }
1553            if {$line != ""} {
1554                set key  [format "LNCN%4d%4d" $number $i]
1555                makeexprec $key
1556                setexp $key $line 3 68
1557            }
1558            return
1559        }
1560        atom-delete {
1561            for {set j $number} {$j < 9999} {incr j} {
1562                # delete records for current constraint
1563                for {set i 1} {$i < 999} {incr i} {
1564                    set key [format "LNCN%4d%4d" $j $i]
1565                    if {![existsexp $key]} break
1566                    delexp $key
1567                }
1568                # now copy records, from the next entry, if any
1569                set j1 $j
1570                incr j1
1571                set key1 [format "LNCN%4d%4d" $j1 1]
1572                # if there is no record, there is nothing to copy -- done
1573                if {![existsexp $key1]} return
1574                for {set i 1} {$i < 999} {incr i} {
1575                    set key1 [format "LNCN%4d%4d" $j1 $i]
1576                    if {![existsexp $key1]} break
1577                    set key  [format "LNCN%4d%4d" $j  $i]
1578                    makeexprec $key
1579                    setexp $key [readexp $key1] 1 68
1580                }
1581            }
1582        }
1583        profile*-delete {
1584            regsub profile $type {} term
1585            if {$term < 10} {
1586                set term " $term"
1587            }
1588            set key "LEQV PF$term   "
1589            # return nothing if no term exists
1590            if {![existsexp $key]} {return 0}
1591
1592            # number of constraint terms
1593            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1594            # don't delete a non-existing entry
1595            if {$number > $nterms} {return 0}
1596            set val [expr $nterms - 1]
1597            validint val 5
1598            setexp $key $val 1 5
1599            for {set i1 $number} {$i1 < $nterms} {incr i1} {
1600                set i2 [expr 1 + $i1]
1601                # move the contents of constraint #i2 -> i1
1602                if {$i1 > 9} {
1603                    set k1 [expr ($i1+1)/10]
1604                    set l1 $i1
1605                } else {
1606                    set k1 " "
1607                    set l1 " $i1"
1608                }
1609                set key1 "LEQV PF$term  $k1"
1610                # number of constraint lines for #i1
1611                set n1 [string trim [string range [readexp ${key1}] \
1612                        [expr ($i1%10)*5] [expr 4+(($i1%10)*5)]] ]
1613                if {$i2 > 9} {
1614                    set k2 [expr ($i2+1)/10]
1615                    set l2 $i2
1616                } else {
1617                    set k2 " "
1618                    set l2 " $i2"
1619                }
1620                set key2 "LEQV PF$term  $k2"
1621                # number of constraint lines for #i2
1622                set n2 [string trim [string range [readexp ${key2}] \
1623                        [expr ($i2%10)*5] [expr 4+(($i2%10)*5)]] ]
1624                set val $n2
1625                validint val 5
1626                # move the # of terms
1627                setexp $key1 $val [expr 1+(($i1%10)*5)] 5
1628                # move the terms
1629                for {set j 1} {$j <= $n2} {incr j 1} {
1630                    set key "LEQV PF${term}${l1}$j"
1631                    makeexprec $key
1632                    setexp $key [readexp "LEQV PF${term}${l2}$j"] 1 68
1633                }
1634                # delete any remaining lines
1635                for {set j [expr $n2+1]} {$j <= $n1} {incr j 1} {
1636                    delexp "LEQV PF${term}${l1}$j"
1637                }
1638            }
1639
1640            # clear the last term
1641            if {$nterms > 9} {
1642                set i [expr ($nterms+1)/10]
1643            } else {
1644                set i " "
1645            }
1646            set key "LEQV PF$term  $i"
1647            set cb [expr ($nterms%10)*5]
1648            set ce [expr 4+(($nterms%10)*5)]
1649            set n2 [string trim [string range [readexp ${key}] $cb $ce] ]
1650            incr cb
1651            setexp $key "     " $cb 5
1652            # delete any remaining lines
1653            for {set j 1} {$j <= $n2} {incr j 1} {
1654                delexp "LEQV PF${term}${nterms}$j"
1655            }
1656        }
1657        profile*-set {
1658            regsub profile $type {} term
1659            if {$term < 10} {
1660                set term " $term"
1661            }
1662            set key "LEQV PF$term   "
1663            # get number of constraint terms
1664            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1665            # don't change a non-existing entry
1666            if {$number > $nterms} {return 0}
1667            if {$number > 9} {
1668                set k1 [expr ($number+1)/10]
1669                set l1 $number
1670            } else {
1671                set k1 " "
1672                set l1 " $number"
1673            }
1674            set key1 "LEQV PF$term  $k1"
1675            # old number of constraint lines
1676            set n1 [string trim [string range [readexp ${key1}] \
1677                    [expr ($number%10)*5] [expr 4+(($number%10)*5)]] ]
1678            # number of new constraints
1679            set j2 [llength $value]
1680            # number of new constraint lines
1681            set val [set n2 [expr ($j2 + 2)/3]]
1682            # store the new # of lines
1683            validint val 5
1684            setexp $key1 $val [expr 1+(($number%10)*5)] 5
1685
1686            # loop over the # of lines in the old or new, whichever is greater
1687            set v0 0
1688            for {set j 1} {$j <= [expr ($n1 > $n2) ? $n1 : $n2]} {incr j 1} {
1689                set key "LEQV PF${term}${l1}$j"
1690                # were there more lines in the old?
1691                if {$j > $n2} {
1692                    # this line is not needed
1693                    if {$j % 3 == 1} {
1694                        delexp %key
1695                    }
1696                    continue
1697                }
1698                # are we adding new lines?
1699                if {$j > $n1} {
1700                    makeexprec $key
1701                }
1702                # add the three constraints to the line
1703                foreach s {3 23 43} \
1704                        item [lrange $value $v0 [expr 2+$v0]] {
1705                    if {$item != ""} {
1706                        set val [format %-10s%9.3f \
1707                                [lindex $item 0],[lindex $item 1] \
1708                                [lindex $item 2]]
1709                        setexp $key $val $s 19
1710                    } else {
1711                        setexp $key " " $s 19
1712                    }
1713                }
1714                incr v0 3
1715            }
1716        }
1717        profile*-add {
1718            regsub profile $type {} term
1719            if {$term < 10} {
1720                set term " $term"
1721            }
1722            set key "LEQV PF$term   "
1723            if {![existsexp $key]} {makeexprec $key}
1724            set nterms [string trim [string range [readexp ${key}] 0 4] ]
1725            if {$nterms == ""} {
1726                set nterms 1
1727            } elseif {$nterms >= 99} {
1728                return 0
1729            } else {
1730                incr nterms
1731            }
1732            # store the new # of constraints
1733            set val $nterms
1734            validint val 5
1735            setexp $key $val 1 5
1736
1737            if {$nterms > 9} {
1738                set k1 [expr ($nterms+1)/10]
1739                set l1 $nterms
1740            } else {
1741                set k1 " "
1742                set l1 " $nterms"
1743            }
1744            set key1 "LEQV PF$term  $k1"
1745
1746            # number of new constraints
1747            set j2 [llength $value]
1748            # number of new constraint lines
1749            set val [set n2 [expr ($j2 + 2)/3]]
1750            # store the new # of lines
1751            validint val 5
1752            setexp $key1 $val [expr 1+(($nterms%10)*5)] 5
1753
1754            # loop over the # of lines to be added
1755            set v0 0
1756            for {set j 1} {$j <= $n2} {incr j 1} {
1757                set key "LEQV PF${term}${l1}$j"
1758                makeexprec $key
1759                # add the three constraints to the line
1760                foreach s {3 23 43} \
1761                        item [lrange $value $v0 [expr 2+$v0]] {
1762                    if {$item != ""} {
1763                        set val [format %-10s%9.3f \
1764                                [lindex $item 0],[lindex $item 1] \
1765                                [lindex $item 2]]
1766                        setexp $key $val $s 19
1767                    } else {
1768                        setexp $key " " $s 19
1769                    }
1770                }
1771                incr v0 3
1772            }
1773        }
1774        profile*-get {
1775            regsub profile $type {} term
1776            if {$term < 10} {
1777                set term " $term"
1778            }
1779            if {$number > 9} {
1780                set i [expr ($number+1)/10]
1781            } else {
1782                set i " "
1783            }
1784            set key "LEQV PF$term  $i"
1785            # return nothing if no term exists
1786            if {![existsexp $key]} {return 0}
1787            # number of constraint lines
1788           
1789            set numline [string trim [string range [readexp ${key}] \
1790                    [expr ($number%10)*5] [expr 4+(($number%10)*5)]] ]
1791            if {$number == 0} {return $numline}
1792            set clist {}
1793            if {$number < 10} {
1794                set number " $number"
1795            }
1796            for {set i 1} {$i <= $numline} {incr i} {
1797                set key "LEQV PF${term}${number}$i"
1798                set line [readexp ${key}]
1799                foreach s {1 21 41} e {20 40 60} {
1800                    set seg [string range $line $s $e]
1801                    if {[string trim $seg] == ""} continue
1802                    # parse the string segment
1803                    set parse [regexp { *([0-9AL]+),([0-9AL]+) +([0-9.]+)} \
1804                            $seg junk phase hist mult]
1805                    # was parse successful
1806                    if {!$parse} {continue}
1807                    lappend clist [list $phase $hist $mult]
1808                }
1809            }
1810            return $clist
1811        }
1812        default {
1813            set msg "Unsupported constrinfo access: type=$type action=$action"
1814            tk_dialog .badexp "Error in readexp access" $msg error 0 OK
1815        }
1816
1817    }
1818}
1819
1820# read the default profile information for a histogram
1821# use: profdefinfo hist set# parm action
1822
1823#     proftype -- profile function number
1824#     profterms -- number of profile terms
1825#     pdamp -- damping value for the profile (*)
1826#     pcut -- cutoff value for the profile (*)
1827#     pterm$n -- profile term #n
1828#     pref$n -- refinement flag value for profile term #n (*)
1829
1830proc profdefinfo {hist set parm "action get"} {
1831    global expgui
1832    if {$hist < 10} {
1833        set key "HST  $hist"
1834    } else {
1835        set key "HST $hist"
1836    }
1837    switch -glob ${parm}-$action {
1838        proftype-get {
1839            set val [string range [readexp "${key}PRCF$set"] 0 4]
1840            if {$val == " "} {return 0}
1841            return $val
1842        }
1843        profterms-get {
1844            set val [string range [readexp "${key}PRCF$set"] 5 9]
1845            if {$val == " "} {return 0}
1846            return $val
1847        }
1848        pcut-get {
1849            return [string trim [string range [readexp "${key}PRCF$set"] 10 19]]
1850        }
1851        pdamp-get {
1852                set val [string range [readexp "${key}PRCF$set"] 24 24]
1853            if {$val == " "} {return 0}
1854            return $val
1855        }
1856        pterm*-get {
1857            regsub pterm $parm {} num
1858            set f1 [expr 15*(($num - 1) % 4)]
1859            set f2 [expr 15*(1 + ($num - 1) % 4)-1]
1860            set line  [expr 1 + ($num - 1) / 4]
1861            return [string trim [string range [\
1862                        readexp "${key}PRCF${set}$line"] $f1 $f2] ]
1863        }
1864        pref*-get {
1865            regsub pref $parm {} num
1866            set f [expr 24+$num]
1867            if {[string toupper [string range [readexp "${key}PRCF$set"] $f $f]] == "Y"} {
1868                return 1
1869            }
1870            return 0
1871        }
1872        default {
1873            set msg "Unsupported profdefinfo access: parm=$parm action=$action"
1874            tk_dialog .badexp "Code Error" $msg error 0 Exit
1875        }
1876    }
1877}
1878
1879# get March-Dollase preferred orientation information
1880# use MDprefinfo hist phase axis-number parm action value
1881#    ratio    -- ratio of xtallites in PO direction vs random (>1 for more)
1882#    fraction -- fraction in this direction, when more than one axis is used
1883#    h k & l  -- indices of P.O. axis
1884#    ratioref -- flag to vary ratio
1885#    fracref  -- flag to vary fraction
1886#    damp     -- damping value
1887#    type     -- model type (0 = P.O. _|_ to beam, 1 = || to beam)
1888#    new      -- creates a new record with default values (set only)
1889proc MDprefinfo {histlist phaselist axislist parm "action get" "value {}"} {
1890    foreach phase $phaselist hist $histlist axis $axislist {
1891        if {$phase == ""} {set phase [lindex $phaselist end]}
1892        if {$hist == ""} {set hist [lindex $histlist end]}
1893        if {$axis == ""} {set axis [lindex $axislist end]}
1894        if {$hist < 10} {
1895            set hist " $hist"
1896        }
1897        if {$axis > 9} {
1898            set axis "0"
1899        }
1900        set key "HAP${phase}${hist}PREFO${axis}"
1901        switch -glob ${parm}-$action {
1902            ratio-get {
1903                return [string trim [string range [readexp $key] 0 9]]
1904            }
1905            ratio-set {
1906                if ![validreal value 10 6] {return 0}
1907                setexp $key $value 1 10
1908            }
1909            fraction-get {
1910                return [string trim [string range [readexp $key] 10 19]]
1911            }
1912            fraction-set {
1913                if ![validreal value 10 6] {return 0}
1914                setexp $key $value 11 10
1915            }
1916            h-get {
1917                set h [string trim [string range [readexp $key] 20 29]]
1918                # why not allow negative h values?
1919                #               if {$h < 1} {return 0}
1920                return $h
1921            }
1922            h-set {
1923                if ![validreal value 10 2] {return 0}
1924                setexp $key $value 21 10
1925            }
1926            k-get {
1927                set k [string trim [string range [readexp $key] 30 39]]
1928                #               if {$k < 1} {return 0}
1929                return $k
1930            }
1931            k-set {
1932                if ![validreal value 10 2] {return 0}
1933                setexp $key $value 31 10
1934            }
1935            l-get {
1936                set l [string trim [string range [readexp $key] 40 49]]
1937                #if {$l < 1} {return 0}
1938                return $l
1939            }
1940            l-set {
1941                if ![validreal value 10 2] {return 0}
1942                setexp $key $value 41 10
1943            }
1944            ratioref-get {
1945                if {[string toupper \
1946                        [string range [readexp $key] 53 53]] == "Y"} {
1947                    return 1
1948                }
1949                return 0
1950            }
1951            ratioref-set {
1952                if $value {
1953                    setexp $key "Y" 54 1
1954                } else {
1955                    setexp $key "N" 54 1
1956                }
1957            }
1958            fracref-get {
1959                if {[string toupper \
1960                        [string range [readexp $key] 54 54]] == "Y"} {
1961                    return 1
1962                }
1963                return 0
1964            }
1965            fracref-set {
1966                if $value {
1967                    setexp $key "Y" 55 1
1968                } else {
1969                    setexp $key "N" 55 1
1970              }
1971            }
1972            damp-get {
1973                set val [string trim [string range [readexp $key] 59 59]]
1974                if {$val == " "} {return 0}
1975                return $val
1976            }
1977            damp-set {
1978                setexp $key $value 60 1
1979            }
1980            type-get {
1981                set val [string trim [string range [readexp $key] 64 64]]
1982                if {$val == " "} {return 0}
1983                return $val
1984            }
1985            type-set {
1986                # only valid settings are 0 & 1
1987                if {$value != "0" && $value != "1"} {set value "0"}
1988                setexp $key $value 65 1
1989            }
1990            new-set {
1991                makeexprec $key
1992                setexp $key \
1993                        {  1.000000  1.000000  0.000000  0.000000  1.000000   NN    0    0} \
1994                        1 68
1995            }
1996            default {
1997                set msg "Unsupported MDprefinfo access: parm=$parm action=$action"
1998                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1999            }
2000
2001        }
2002
2003    }
2004}
2005
2006# write the .EXP file
2007proc expwrite {expfile} {
2008    global tcl_platform exparray
2009    set blankline \
2010     "                                                                        "
2011    set fp [open ${expfile} w]
2012    set keylist [lsort [array names exparray]]
2013    # reorder the keys so that VERSION comes 1st
2014    set pos [lsearch -exact $keylist {     VERSION}]
2015    set keylist "{     VERSION} [lreplace $keylist $pos $pos]"
2016    if {$tcl_platform(platform) == "windows"} { 
2017        foreach key $keylist {
2018            puts $fp [string range \
2019                    "$key$exparray($key)$blankline" 0 79]
2020        }
2021    } else {
2022        foreach key $keylist {
2023            puts -nonewline $fp [string range \
2024                    "$key$exparray($key)$blankline" 0 79]
2025        }
2026    }
2027    close $fp
2028}
2029
2030# history commands -- delete all but last $keep history records,
2031# renumber if $renumber is true
2032proc DeleteHistory {keep renumber} {
2033    global exparray
2034    foreach y [lrange [lsort -decreasing \
2035            [array names exparray {    HSTRY*}]] $keep end] {
2036        unset exparray($y)
2037    }
2038    if !$renumber return
2039    # renumber
2040    set i 0
2041    foreach y [lsort -increasing \
2042            [array names exparray {    HSTRY*}]] {
2043        set key [format "    HSTRY%3d" [incr i]]
2044        set exparray($key) $exparray($y)
2045        unset exparray($y)
2046    }
2047    # list all history
2048    #    foreach y [lsort -decreasing [array names exparray {    HSTRY*}]] {puts "$y $exparray($y)"}
2049}
2050
2051proc CountHistory {} {
2052    global exparray
2053    return [llength [array names exparray {    HSTRY*}]]
2054}
2055
2056# set the phase flags for histogram $hist to $plist
2057proc SetPhaseFlag {hist plist} {
2058    # make a 2 digit key -- hh
2059    if {$hist < 10} {
2060        set hh " $hist"
2061    } else {
2062        set hh $hist
2063    }
2064    set key "HST $hh NPHAS"
2065    set str {}
2066    foreach iph {1 2 3 4 5 6 7 8 9} {
2067        if {[lsearch $plist $iph] != -1} {
2068            append str {    1}
2069        } else {
2070            append str {    0}     
2071        }
2072    }
2073    setexp $key $str 1 68
2074}
2075
2076# erase atom $atom from phase $phase
2077# update the list of atom types, erasing the record if not needed.
2078proc EraseAtom {atom phase} {
2079    set type [atominfo $phase $atom type]
2080    if {$type == ""} return
2081    if {$atom < 10} {
2082        set key "CRS$phase  AT  $atom"
2083    } elseif {$atom < 100} {
2084        set key "CRS$phase  AT $atom"
2085    } else {
2086        set key "CRS$phase  AT$atom"
2087    }
2088    # delete the records for the atom
2089    global exparray
2090    foreach k [array names exparray ${key}*] {
2091        delexp $k
2092    }
2093    # change the number of atoms in the phase
2094    phaseinfo $phase natoms set [expr [phaseinfo $phase natoms] -1]
2095
2096    # now adjust numbers in "EXPR ATYP" records and delete, if needed.
2097    set natypes [readexp " EXPR  NATYP"]
2098    if {$natypes == ""} return
2099    set j 0
2100    for {set i 1} {$i <= $natypes} {incr i} {
2101        incr j
2102        if {$j <10} {
2103            set key " EXPR ATYP $j"
2104        } else {
2105            set key " EXPR ATYP$j"
2106        }
2107        while {![existsexp $key]} {
2108            incr j
2109            if {$j > 99} {
2110                return
2111            } elseif {$j <10} {
2112                set key " EXPR ATYP $j"
2113            } else {
2114                set key " EXPR ATYP$j"
2115            }
2116        }
2117        set keytype [string trim [string range $exparray($key) 2 9]]
2118        if {$type == $keytype} {
2119            # found the type record
2120            set val [string trim [string range $exparray($key) 10 14]]
2121            incr val -1
2122            # if this is the last reference, remove the record,
2123            # otherwise, decrement the counter
2124            if {$val <= 0} {
2125                incr natypes -1 
2126                validint natypes 5
2127                setexp " EXPR  NATYP" $natypes 1 5
2128                delexp $key
2129            } else {
2130                validint val 5
2131                setexp $key $val 11 5
2132            }
2133            return
2134        }
2135    }
2136}
2137
2138# compute equivalent anisotropic temperature factor for Uequiv
2139proc CalcAniso {phase Uequiv} {
2140    foreach var {a b c alpha beta gamma} {
2141        set $var [phaseinfo $phase $var]
2142    }
2143
2144    set G(1,1) [expr $a * $a]
2145    set G(2,2) [expr $b * $b]
2146    set G(3,3) [expr $c * $c]
2147    set G(1,2) [expr $a * $b * cos($gamma*0.017453292519943)]
2148    set G(2,1) $G(1,2)
2149    set G(1,3) [expr $a * $c * cos($beta *0.017453292519943)]
2150    set G(3,1) $G(1,3)
2151    set G(2,3) [expr $b * $c * cos($alpha*0.017453292519943)]
2152    set G(3,2) $G(2,3)
2153
2154    # Calculate the volume**2
2155    set v2 0.0
2156    foreach i {1 2 3} {
2157        set J [expr ($i%3) + 1]
2158        set K [expr (($i+1)%3) + 1]
2159        set v2 [expr $v2+ $G(1,$i)*($G(2,$J)*$G(3,$K)-$G(3,$J)*$G(2,$K))]
2160    }
2161    if {$v2 > 0} {
2162        set v [expr sqrt($v2)]
2163        foreach i {1 2 3} {
2164            set i1 [expr ($i%3) + 1]
2165            set i2 [expr (($i+1)%3) + 1]
2166            foreach j {1 2 3} {
2167                set j1 [expr ($j%3) + 1]
2168                set j2 [expr (($j+1)%3) + 1]
2169                set C($j,$i) [expr (\
2170                        $G($i1,$j1) * $G($i2,$j2) - \
2171                        $G($i1,$j2)  * $G($i2,$j1)\
2172                        )/ $v]
2173            }
2174        }
2175        set A(1,2) [expr 0.5 * ($C(1,2) + $C(2,1)) / sqrt( $C(1,1)* $C(2,2) )]
2176        set A(1,3) [expr 0.5 * ($C(1,3) + $C(3,1)) / sqrt( $C(1,1)* $C(3,3) )]
2177        set A(2,3) [expr 0.5 * ($C(2,3) + $C(3,2)) / sqrt( $C(2,2)* $C(3,3) )]
2178        foreach i {1 1 2} j {2 3 3} {
2179            set A($i,$j) [expr 0.5 * ($C($i,$j) + $C($j,$i)) / \
2180                    sqrt( $C($i,$i)* $C($j,$j) )]
2181            # clean up roundoff
2182            if {abs($A($i,$j)) < 1e-5} {set A($i,$j) 0.0}
2183        }
2184    } else {
2185        set A(1,2) 0.0
2186        set A(1,3) 0.0
2187        set A(2,3) 0.0
2188    }
2189    return "$Uequiv $Uequiv $Uequiv \
2190            [expr $Uequiv * $A(1,2)] \
2191            [expr $Uequiv * $A(1,3)] \
2192            [expr $Uequiv * $A(2,3)]"
2193}
Note: See TracBrowser for help on using the repository browser.