source: trunk/readexp.tcl @ 381

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

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