source: trunk/readexp.tcl @ 387

Last change on this file since 387 was 387, checked in by toby, 13 years ago

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