source: trunk/readexp.tcl @ 1175

Last change on this file since 1175 was 1175, checked in by toby, 9 years ago

update stable release to fix bugs & match sandbox: many RB updates; show RB atoms in phase

  • Property svn:keywords set to Author Date Revision Id
File size: 121.0 KB
Line 
1# $Id: readexp.tcl 1175 2011-10-25 04:17:42Z toby $
2# Routines to deal with the .EXP "data structure"
3set expmap(Revision) {$Revision: 1175 $ $Date: 2011-10-25 04:17:42 +0000 (Tue, 25 Oct 2011) $}
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 "ns {}"} {
12    # expfile is the path to the data file.
13    # ns is the namespace to place the output array (default is global)
14    if {$ns != ""} {
15        namespace eval $ns {}
16    }
17    if [catch {set fil [open "$expfile" r]}] {
18        tk_dialog .expFileErrorMsg "File Open Error" \
19                "Unable to open file $expfile" error 0 "Exit" 
20        return -1
21    }
22    fconfigure $fil -translation lf
23    set len [gets $fil line]
24    if {[string length $line] != $len} {
25        tk_dialog .expConvErrorMsg "old tcl" \
26                "You are using an old version of Tcl/Tk and your .EXP file has binary characters; run convstod or upgrade" \
27                error 0 "Exit"
28        return -1
29    }
30    catch {
31        unset ${ns}::exparray
32    }
33    if {$len > 160} {
34        set fmt 0
35        # a UNIX-type file
36        set i1 0
37        set i2 79
38        while {$i2 < $len} {
39            set nline [string range $line $i1 $i2]
40            incr i1 80
41            incr i2 80
42            set key [string range $nline 0 11]
43            set ${ns}::exparray($key) [string range $nline 12 end]
44        }
45    } else {
46        set fmt 1
47        while {$len > 0} {
48            set key [string range $line 0 11]
49            set ${ns}::exparray($key) [string range $line 12 79]
50            if {$len != 81 || [string range $line end end] != "\r"} {set fmt 2}
51            set len [gets $fil line]
52        }
53    }
54    close $fil
55    return $fmt
56}
57
58proc createexp {expfile title} {
59    global exparray expmap
60    catch {unset exparray}
61    foreach key   {"     VERSION" "      DESCR" "ZZZZZZZZZZZZ" " EXPR NPHAS"} \
62            value {"   6"         ""            "  Last EXP file record" ""} {
63        # truncate long keys & pad short ones
64        set key [string range "$key        " 0 11]
65        set exparray($key) $value
66    }
67    expinfo title set $title
68    exphistory add " created readexp.tcl [lindex $expmap(Revision) 1] [clock format [clock seconds] -format %Y-%m-%dT%T]"
69    expwrite $expfile
70}
71
72# get information out from an EXP file
73#   creates the following entries in global array expmap
74#     expmap(phaselist)     gives a list of defined phases
75#     expmap(phasetype)     gives the phase type for each defined phase
76#                           =1 nuclear; 2 mag+nuc; 3 mag; 4 macro
77#     expmap(atomlist_$p)   gives a list of defined atoms in phase $p
78#     expmap(htype_$n)      gives the GSAS histogram type for histogram (all)
79#     expmap(powderlist)    gives a list of powder histograms in use
80#     expmap(phaselist_$n)  gives a list of phases used in histogram $n
81#     expmap(nhst)          the number of GSAS histograms
82#
83proc mapexp {} {
84    global expgui expmap exparray
85    # clear out the old array
86    set expmap_Revision $expmap(Revision)
87    unset expmap
88    set expmap(Revision) $expmap_Revision
89    # apply any updates to the .EXP file
90    updateexpfile
91    # get the defined phases
92    set line [readexp " EXPR NPHAS"]
93#    if {$line == ""} {
94#       set msg "No EXPR NPHAS entry. This is an invalid .EXP file"
95#       tk_dialog .badexp "Error in EXP" $msg error 0 Exit
96#       destroy .
97#    }
98    set expmap(phaselist) {}
99    set expmap(phasetype) {}
100    # loop over phases
101    foreach iph {1 2 3 4 5 6 7 8 9} {
102        set i5s [expr {($iph - 1)*5}]
103        set i5e [expr {$i5s + 4}]
104        set flag [string trim [string range $line $i5s $i5e]]
105        if {$flag == ""} {set flag 0}
106        if $flag {
107            lappend expmap(phaselist) $iph
108            lappend expmap(phasetype) $flag
109        }
110    }
111    # get the list of defined atoms for each phase
112    foreach iph $expmap(phaselist) {
113        set expmap(atomlist_$iph) {}
114        if {[lindex $expmap(phasetype) [expr {$iph - 1}]] != 4} {
115            foreach key [array names exparray "CRS$iph  AT*A"] {
116                regexp { AT *([0-9]+)A} $key a num
117                lappend expmap(atomlist_$iph) $num
118            }
119        } else {
120            foreach key [array names exparray "CRS$iph  AT*"] {
121                scan [string range $key 8 11] %x atm
122                lappend expmap(atomlist_$iph) $atm
123            }
124        }
125        # note that sometimes an .EXP file contains more atoms than are actually defined
126        # drop the extra ones
127        set expmap(atomlist_$iph) [lsort -integer $expmap(atomlist_$iph)]
128        set natom [phaseinfo $iph natoms]
129        if {$natom != [llength $expmap(atomlist_$iph)]} {
130            set expmap(atomlist_$iph) [lrange $expmap(atomlist_$iph) 0 [expr {$natom-1}]]
131        }
132    }
133    # now get the histogram types
134    set expmap(nhst) [string trim [readexp { EXPR  NHST }]]
135    set n 0
136    set expmap(powderlist) {}
137    for {set i 0} {$i < $expmap(nhst)} {incr i} {
138        set ihist [expr {$i + 1}]
139        if {[expr {$i % 12}] == 0} {
140            incr n
141            set line [readexp " EXPR  HTYP$n"]
142            if {$line == ""} {
143                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
144                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
145            }
146            set j 0
147        } else {
148            incr j
149        }
150        set expmap(htype_$ihist) [string range $line [expr 2+5*$j] [expr 5*($j+1)]]
151        # is this a dummy histogram?
152        if {$ihist <=9} {
153            set key "HST  $ihist DUMMY"
154        } else {
155            set key "HST $ihist DUMMY"
156        }
157        # at least for now, ignore non-powder histograms
158        if {[string range $expmap(htype_$ihist) 0 0] == "P" && \
159                [string range $expmap(htype_$ihist) 3 3] != "*"} {
160            if {[existsexp $key]} {
161                set expmap(htype_$ihist) \
162                        [string range $expmap(htype_$ihist) 0 2]D
163            }
164            lappend expmap(powderlist) $ihist
165        }
166    }
167
168    # now process powder histograms
169    foreach ihist $expmap(powderlist) {
170        # make a 2 digit key -- hh
171        if {$ihist < 10} {
172            set hh " $ihist"
173        } else {
174            set hh $ihist
175        }
176        set line [readexp "HST $hh NPHAS"]
177        if {$line == ""} {
178            set msg "No NPHAS entry for Histogram $ihist. This is an invalid .EXP file"
179            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
180        }
181        set expmap(phaselist_$ihist) {}
182        # loop over phases
183        foreach iph {1 2 3 4 5 6 7 8 9} {
184            set i5s [expr {($iph - 1)*5}]
185            set i5e [expr {$i5s + 4}]
186            set flag [string trim [string range $line $i5s $i5e]]
187            if {$flag == ""} {set flag 0}
188            if $flag {lappend expmap(phaselist_$ihist) $iph}
189        }
190    }
191    # load the constrained parameters
192    atom_constraint_load
193    # construct tables of mapped atoms in rigid bodies
194    foreach phase $::expmap(phaselist) {
195        set expmap(rbatoms_$phase) {}
196        foreach bodnum [RigidBodyList] {
197            set natoms [llength [lindex [lindex [lindex [ReadRigidBody $bodnum] 1] 0] 3]]
198            foreach mapnum [RigidBodyMappingList $phase $bodnum] {
199                set atomnum [lindex [ReadRigidBodyMapping $phase $bodnum $mapnum] 0]
200                set st [lsearch $::expmap(atomlist_$phase) $atomnum]
201                set en [expr {$st+$natoms-1}]
202                set atoms [lrange $::expmap(atomlist_$phase) $st $en]
203                set expmap(rbatoms_$phase) [concat $expmap(rbatoms_$phase) $atoms]
204            }
205        }
206    }
207    set expgui(mapstat) 1
208}
209
210# this routine is called to update changes in the .EXP file
211proc updateexpfile {} {
212    global exparray
213    # change "CRSx  ODFxxx" records to "CRSx  OD xxx" records
214    # needed by the June 8, 2005 GSAS release
215    set ODFlist [array names exparray "CRS?  ODF*"]
216    if {[llength $ODFlist] > 0} {
217        catch {incr ::expgui(changed)}
218        foreach key $ODFlist {
219            regsub "ODF" $key "OD " newkey
220            set exparray($newkey) $exparray($key)
221            unset exparray($key)
222        }
223    }
224}
225
226
227# return the value for a ISAM key
228proc readexp {key} {
229    global exparray
230    # truncate long keys & pad short ones
231    set key [string range "$key        " 0 11]
232    if [catch {set val $exparray($key)}] {
233        global expgui
234        if $expgui(debug) {puts "Error accessing record $key"}
235        return ""
236    }
237    return $val
238}
239
240# return the number of records matching ISAM key (may contain wildcards)
241proc existsexp {key} {
242    global exparray
243    # key can contain wild cards so don't pad
244    return [llength [array names exparray  $key]]
245}
246
247
248# replace a section of the exparray with $value
249#   replace $char characters starting at character $start (numbered from 1)
250proc setexp {key value start chars} {
251    global exparray
252    # truncate long keys & pad short ones
253    set key [string range "$key        " 0 11]
254    if [catch {set exparray($key)}] {
255        global expgui
256        if $expgui(debug) {puts "Error accessing record $key"}
257        return ""
258    }
259
260    # pad value to $chars
261    set l0 [expr {$chars - 1}]
262    set value [string range "$value                                           " 0 $l0]
263
264    if {$start == 1} {
265        set ret {}
266        set l1 $chars
267    } else {
268        set l0 [expr {$start - 2}]
269        set l1 [expr {$start + $chars - 1}]
270        set ret [string range $exparray($key) 0 $l0]
271    }
272    append ret $value [string range $exparray($key) $l1 end]
273    set exparray($key) $ret
274}
275
276proc makeexprec {key} {
277    global exparray
278    # truncate long keys & pad short ones
279    set key [string range "$key        " 0 11]
280    if [catch {set exparray($key)}] {
281        # set to 68 blanks
282        set exparray($key) [format %68s " "]
283    }
284}
285
286# delete an exp record
287# returns 1 if OK; 0 if not found
288proc delexp {key} {
289    global exparray
290    # truncate long keys & pad short ones
291    set key [string range "$key        " 0 11]
292    if [catch {unset exparray($key)}] {
293        return 0
294    }
295    return 1
296}
297
298# test an argument if it is a valid number; reform the number to fit
299proc validreal {val length decimal} {
300    upvar $val value
301    # is this a number?
302    if [catch {expr {$value}}] {return 0}
303    if [catch {
304        # how many digits are needed to the left of the decimal?
305        set sign 0
306        if {$value > 0} {
307            set digits [expr {1 + int(log10($value))}]
308        } elseif {$value < 0} {
309            set digits [expr {1 + int(log10(-$value))}]
310            set sign 1
311        } else {
312            set digits 1
313        }
314        if {$digits <= 0} {set digits 1}
315        if {$digits + $sign >= $length} {
316            # the number is much too big -- use exponential notation
317            set decimal [expr {$length - 6 - $sign}]
318            # drop more decimal places, as needed
319            set tmp [format "%${length}.${decimal}E" $value]
320            while {[string length $tmp] > $length && $decimal >= 0} {
321                incr decimal -1
322                set tmp [format "%${length}.${decimal}E" $value]
323            }
324        } elseif {$digits + $sign >= $length - $decimal} {
325            # we will have to trim the number of decimal digits
326            set decimal [expr {$length - $digits - $sign - 1}]
327            set tmp [format "%#.${decimal}f" $value]
328        } elseif {abs($value) < pow(10,2-$decimal) && $length > 6 && $value != 0} {
329            # for small values, switch to exponential notation (2-$decimal -> three sig figs.)
330            set decimal [expr {$length - 6 - $sign}]
331            # drop more decimal places, as needed
332            set tmp [format "%${length}.${decimal}E" $value]
333            while {[string length $tmp] > $length && $decimal >= 0} {
334                incr decimal -1
335                set tmp [format "%${length}.${decimal}E" $value]
336            }
337        } else {
338            # use format as specified
339            set tmp [format "%${length}.${decimal}f" $value]
340        }
341        set value $tmp
342    } errmsg] {return 0}
343    return 1
344}
345
346# test an argument if it is a valid integer; reform the number into
347# an integer, if appropriate -- be sure to pass the name of the variable not the value
348proc validint {val length} {
349    upvar $val value
350    # FORTRAN type assumption: blank is 0
351    if {$value == ""} {set value 0}
352    if [catch {
353        set tmp [expr {round($value)}]
354        if {$tmp != $value} {return 0}
355        set value [format "%${length}d" $tmp]
356    }] {return 0}
357    return 1
358}
359
360# process history information
361#    action == last
362#       returns number and value of last record
363#    action == add
364#
365proc exphistory {action "value 0"} {
366    global exparray
367    if {$action == "last"} {
368        set key [lindex [lsort -decreasing [array names exparray *HSTRY*]] 0]
369        if {$key == ""} {return ""}
370        return [list [string trim [string range $key 9 end]] $exparray($key)]
371    } elseif {$action == "add"} {
372        set key [lindex [lsort -decreasing [array names exparray *HSTRY*]] 0]
373        if {$key == ""} {
374            set index 1
375        } else {
376            set index [string trim [string range $key 9 end]]
377            if {$index != "***"} {
378                if {$index < 999} {incr index}
379                set key [format "    HSTRY%3d" $index]
380                set exparray($key) $value
381            }
382        }
383        set key [format "    HSTRY%3d" $index]
384        set exparray($key) $value
385    }
386}
387# get overall info
388#   parm:
389#     print     -- GENLES print option (*)
390#     cycles    -- number of GENLES cycles (*)
391#     title     -- the overall title (*)
392#     convg     -- convergence criterion: -200 to 200 (*)
393#     marq      -- Marquardt damping factor: 1.0 to 9.99 (*)
394#     mbw       -- LS matrix bandwidth; =0 for full matrix (*)
395proc expinfo {parm "action get" "value {}"} {
396    switch ${parm}-$action {
397        title-get {
398            return [string trim [readexp "      DESCR"]]
399        }
400        title-set {
401            setexp "      DESCR" "  $value" 2 68
402        }
403        cycles-get {
404            return [string trim [cdatget MXCY]]
405        }
406        cycles-set {
407            if ![validint value 1] {return 0}
408            cdatset MXCY [format %4d $value]
409        }
410        cyclesrun-get {
411            set cycle -1
412            regexp {.*cycles run *([0-9]*) } [readexp "  GNLS  RUN"] x cycle
413            return $cycle
414        }
415        print-get {
416            set print [string trim [cdatget PRNT]]
417            if {$print != ""} {return $print}
418            return 0
419        }
420        print-set {
421            if ![validint value 1] {return 0}
422            cdatset PRNT [format %4d $value]
423        }
424        convg-get {
425            set cvg [string trim [cdatget CVRG]]
426            if {$cvg == ""} {return -200}
427            if [catch {expr {$cvg}}] {return -200}
428            return $cvg
429        }
430        convg-set {
431            if ![validint value 1] {return 0}
432            set value [expr {-200>$value?-200:$value}]
433            set value [expr {200<$value?200:$value}]
434            cdatset CVRG [format %4d $value]
435        }
436        marq-get {
437            set mar [string trim [cdatget MARQ]]
438            if {$mar == ""} {return 1.0}
439            if [catch {expr $mar}] {return 1.}
440            return $mar
441        }
442        marq-set {
443            if [catch {
444                set value [expr {1.0>$value?1.0:$value}]
445                set value [expr {9.99<$value?9.99:$value}]
446            }] {return 0}
447            if ![validreal value 4 2] {return 0}
448            cdatset MARQ $value
449        }
450        mbw-get {
451            set mbw [string trim [cdatget MBW]]
452            if {$mbw == ""} {return 0}
453            if [catch {expr $mbw}] {return 0}
454            return $mbw
455        }
456        mbw-set {
457            if ![validint value 1] {return 0}
458            if {$value < 0} {return 0}
459            cdatset MBW [format %5d $value]
460        }
461        default {
462            set msg "Unsupported expinfo access: parm=$parm action=$action"
463            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
464        }
465    }
466    return 1
467}
468
469proc cdatget {key} {
470    foreach i {1 2 3 4 5 6 7 8 9} {
471        if {[existsexp "  GNLS CDAT$i"] == 0} break
472        set line [readexp "  GNLS CDAT$i"]
473        if {$line == {}} break
474        foreach i1 {2 10 18 26 34 42 50 58 66} \
475                i2 {9 17 25 33 41 49 57 65 73} {
476            set item [string range $line $i1 $i2]
477            if {[string trim $item] == {}} continue
478            if [regexp "${key}(.*)" $item a b] {return $b}
479        }
480    }
481    return {}
482}
483
484proc cdatset {key value} {
485    # round 1 see if we can find the string
486    foreach i {1 2 3 4 5 6 7 8 9} {
487        set line [readexp "  GNLS CDAT$i"]
488        if {$line == {}} break
489        foreach i1 {2 10 18 26 34 42 50 58 66} \
490                i2 {9 17 25 33 41 49 57 65 73} {
491            set item [string range $line $i1 $i2]
492            if {[string trim $item] == {}} continue
493            if [regexp "${key}(.*)" $item a b] {
494                # found it now replace it
495                incr i1
496                setexp "  GNLS CDAT$i" "${key}${value}" $i1 8
497                return
498            }
499        }
500    }
501    # not found, take the 1st blank space, creating a card if needed
502    foreach i {1 2 3 4 5 6 7 8 9} {
503        set line [readexp "  GNLS CDAT$i"]
504        if {$line == {}} {makeexprec "  GNLS CDAT$i"}
505        foreach i1 {2 10 18 26 34 42 50 58 66} \
506                i2 {9 17 25 33 41 49 57 65 73} {
507            set item [string range $line $i1 $i2]
508            if {[string trim $item] == {}} {
509                # found a blank space: now replace it
510                incr i1
511                setexp "  GNLS CDAT$i" "${key}${value}" $i1 8
512                return
513            }
514        }
515    }
516    return {}
517}
518
519proc disagldat_get {phase} {
520    set key "  DSGL CDAT$phase"
521    if {[existsexp $key] == 0} {return "{none} {none}"}
522    set line [readexp $key]
523    set i1 2
524    # read atom-atom distance parameter
525    set dist {}
526    set item [string range $line $i1 [expr {$i1+3}]]
527    if {$item == "DMAX"} {
528        set val [string range $line [expr {$i1+4}] [expr {$i1+11}]]
529        set dist [string trim $val]
530        incr i1 13
531    } else {
532        set dist "radii"
533        incr i1 5
534    }
535    # read atom-atom-atom angle parameter
536    set ang {}
537    set item [string range $line $i1 [expr {$i1+3}]]
538    if {$item == "DAGL"} {
539        set val [string range $line [expr {$i1+4}] [expr {$i1+11}]]
540        set ang [string trim $val]
541        incr i1 13
542    } else {
543        set ang "radii"
544        incr i1 5
545    }
546    # note there are two more parameters, NOFO/FORA & ONCR/DFLT, but they are not being processed yet
547    return "$dist $ang"
548}
549
550# get phase information: phaseinfo phase parm action value
551#   phase: 1 to 9 (as defined)
552#   parm:
553#     name -- phase name
554#     natoms -- number of atoms (*)
555#     a b c alpha beta gamma -- cell parameters (*)
556#     cellref -- refinement flag for the unit cell(*)
557#     celldamp  -- damping for the unit cell refinement (*)
558#     spacegroup -- space group symbol (*)
559#     ODForder -- spherical harmonic order (*)
560#     ODFsym   -- sample symmetry (0-3) (*)
561#     ODFdampA -- damping for angles (*)
562#     ODFdampC -- damping for coefficients (*)
563#     ODFomega -- omega oriention angle (*)
564#     ODFchi -- chi oriention angle (*)
565#     ODFphi -- phi oriention angle (*)
566#     ODFomegaRef -- refinement flag for omega (*)
567#     ODFchiRef -- refinement flag for chi (*)
568#     ODFphiRef -- refinement flag for phi (*)
569#     ODFterms -- a list of the {l m n} values for each ODF term (*)
570#     ODFcoefXXX -- the ODF coefficient for for ODF term XXX (*)
571#     ODFRefcoef -- refinement flag for ODF terms (*)
572#     DistCalc   -- returns "radii", "none" or a number (*)
573#                   none: no distance or angle computation for the phase
574#                   radii: computation will be done by sums of radii
575#                          (see AtmTypInfo and DefAtmTypInfo)
576#                   other: a distance specifing the maximum distance
577#     AngCalc    -- returns "radii", "none" or a number (*)
578#                   none: no distance or angle computation for the phase
579#                   radii: computation will be done by sums of radii
580#                          (see AtmTypInfo and DefAtmTypInfo)
581#                   other: a distance specifing the maximum distance
582#  action: get (default) or set
583#  value: used only with set
584#  * =>  read+write supported
585proc phaseinfo {phase parm "action get" "value {}"} {
586    switch -glob ${parm}-$action {
587
588        name-get {
589            return [string trim [readexp "CRS$phase    PNAM"]]
590        }
591
592        name-set {
593            setexp "CRS$phase    PNAM" " $value" 2 68
594        }
595
596        spacegroup-get {
597            return [string trim [readexp "CRS$phase  SG SYM"]]
598        }
599
600        spacegroup-set {
601            setexp "CRS$phase  SG SYM" " $value" 2 68
602        }
603
604        natoms-get {
605            return [string trim [readexp "CRS$phase   NATOM"]]     
606        }
607
608        natoms-set {
609            if ![validint value 5] {return 0}
610            setexp "CRS$phase   NATOM" $value 1 5
611        }
612
613        a-get {
614           return [string trim [string range [readexp "CRS$phase  ABC"] 0 9]]
615        }
616        b-get {
617           return [string trim [string range [readexp "CRS$phase  ABC"] 10 19]]
618        }
619        c-get {
620           return [string trim [string range [readexp "CRS$phase  ABC"] 20 29]]
621        }
622        alpha-get {
623           return [string trim [string range [readexp "CRS$phase  ANGLES"] 0 9]]
624        }
625        beta-get {
626           return [string trim [string range [readexp "CRS$phase  ANGLES"] 10 19]]
627        }
628        gamma-get {
629           return [string trim [string range [readexp "CRS$phase  ANGLES"] 20 29]]
630        }
631
632        a-set {
633            if ![validreal value 10 6] {return 0}
634            setexp "CRS$phase  ABC" $value 1 10             
635        }
636        b-set {
637            if ![validreal value 10 6] {return 0}
638            setexp "CRS$phase  ABC" $value 11 10           
639        }
640        c-set {
641            if ![validreal value 10 6] {return 0}
642            setexp "CRS$phase  ABC" $value 21 10           
643        }
644        alpha-set {
645            if ![validreal value 10 4] {return 0}
646            setexp "CRS$phase  ANGLES" $value 1 10         
647        }
648        beta-set {
649            if ![validreal value 10 4] {return 0}
650            setexp "CRS$phase  ANGLES" $value 11 10         
651        }
652        gamma-set {
653            if ![validreal value 10 4] {return 0}
654            setexp "CRS$phase  ANGLES" $value 21 10         
655        }
656        cellref-get {
657            if {[string toupper [string range [readexp "CRS$phase  ABC"] 34 34]] == "Y"} {
658                return 1
659            }
660            return 0
661        }
662        cellref-set {
663            if $value {
664                setexp "CRS$phase  ABC" "Y" 35 1
665            } else {
666                setexp "CRS$phase  ABC" "N" 35 1
667            }       
668        }
669        celldamp-get {
670            set val [string range [readexp "CRS$phase  ABC"] 39 39]
671            if {$val == " "} {return 0}
672            return $val
673        }
674        celldamp-set {
675            setexp "CRS$phase  ABC" $value 40 1
676        }
677
678        ODForder-get {
679            set val [string trim [string range [readexp "CRS$phase  OD "] 0 4]]
680            if {$val == " "} {return 0}
681            return $val
682        }
683        ODForder-set {
684            if ![validint value 5] {return 0}
685            setexp "CRS$phase  OD " $value 1 5
686        }
687        ODFsym-get {
688            set val [string trim [string range [readexp "CRS$phase  OD "] 10 14]]
689            if {$val == " "} {return 0}
690            return $val
691        }
692        ODFsym-set {
693            if ![validint value 5] {return 0}
694            setexp "CRS$phase  OD " $value 11 5
695        }
696        ODFdampA-get {
697            set val [string range [readexp "CRS$phase  OD "] 24 24]
698            if {$val == " "} {return 0}
699            return $val
700        }
701        ODFdampA-set {
702            setexp "CRS$phase  OD " $value 25 1
703        }
704        ODFdampC-get {
705            set val [string range [readexp "CRS$phase  OD "] 29 29]
706            if {$val == " "} {return 0}
707            return $val
708        }
709        ODFdampC-set {
710            setexp "CRS$phase  OD " $value 30 1
711        }
712        ODFomegaRef-get {
713            if {[string toupper [string range [readexp "CRS$phase  OD "] 16 16]] == "Y"} {
714                return 1
715            }
716            return 0
717        }
718        ODFomegaRef-set {
719            if $value {
720                setexp "CRS$phase  OD " "Y" 17 1
721            } else {
722                setexp "CRS$phase  OD " "N" 17 1
723            }       
724        }
725        ODFchiRef-get {
726            if {[string toupper [string range [readexp "CRS$phase  OD "] 17 17]] == "Y"} {
727                return 1
728            }
729            return 0
730        }
731        ODFchiRef-set {
732            if $value {
733                setexp "CRS$phase  OD " "Y" 18 1
734            } else {
735                setexp "CRS$phase  OD " "N" 18 1
736            }       
737        }
738        ODFphiRef-get {
739            if {[string toupper [string range [readexp "CRS$phase  OD "] 18 18]] == "Y"} {
740                return 1
741            }
742            return 0
743        }
744        ODFphiRef-set {
745            if $value {
746                setexp "CRS$phase  OD " "Y" 19 1
747            } else {
748                setexp "CRS$phase  OD " "N" 19 1
749            }       
750        }
751        ODFcoef*-get {
752            regsub ODFcoef $parm {} term
753            set k [expr {($term+5)/6}]
754            if {$k <= 9} {
755                set k "  $k"
756            } elseif {$k <= 99} {
757                set k " $k"
758            }
759            set j [expr {(($term-1) % 6)+1}]
760            set lineB [readexp "CRS$phase  OD${k}B"]
761            set j0 [expr { ($j-1) *10}]
762            set j1 [expr {$j0 + 9}]
763            set val [string trim [string range $lineB $j0 $j1]]
764            if {$val == ""} {return 0.0}
765            return $val
766        }
767        ODFcoef*-set {
768            regsub ODFcoef $parm {} term
769            if ![validreal value 10 3] {return 0}
770            set k [expr {($term+5)/6}]
771            if {$k <= 9} {
772                set k "  $k"
773            } elseif {$k <= 99} {
774                set k " $k"
775            }
776            set j [expr {(($term-1) % 6)+1}]
777            set col [expr { ($j-1)*10 + 1}]
778            setexp "CRS$phase  OD${k}B" $value $col 10
779        }
780        ODFRefcoef-get {
781            if {[string toupper [string range [readexp "CRS$phase  OD "] 19 19]] == "Y"} {
782                return 1
783            }
784            return 0
785        }
786        ODFRefcoef-set {
787            if $value {
788                setexp "CRS$phase  OD " "Y" 20 1
789            } else {
790                setexp "CRS$phase  OD " "N" 20 1
791            }       
792        }
793        ODFomega-get {
794           return [string trim [string range [readexp "CRS$phase  OD "] 30 39]]
795        }
796        ODFchi-get {
797           return [string trim [string range [readexp "CRS$phase  OD "] 40 49]]
798        }
799        ODFphi-get {
800           return [string trim [string range [readexp "CRS$phase  OD "] 50 59]]
801        }
802        ODFomega-set {
803            if ![validreal value 10 4] {return 0}
804            setexp "CRS$phase  OD " $value 31 10
805        }
806        ODFchi-set {
807            if ![validreal value 10 4] {return 0}
808            setexp "CRS$phase  OD " $value 41 10
809        }
810        ODFphi-set {
811            if ![validreal value 10 4] {return 0}
812            setexp "CRS$phase  OD " $value 51 10
813        }
814
815        ODFterms-get {
816            set vallist {}
817            set val [string trim [string range [readexp "CRS$phase  OD "] 5 9]]
818            for {set i 1} {$i <= $val} {incr i 6} {
819                set k [expr {1+($i-1)/6}]
820                if {$k <= 9} {
821                    set k "  $k"
822                } elseif {$k <= 99} {
823                    set k " $k"
824                }
825                set lineA [readexp "CRS$phase  OD${k}A"]
826                set k 0
827                for {set j $i} {$j <= $val && $j < $i+6} {incr j} {
828                    set j0 [expr {($k)*10}]
829                    set j1 [expr {$j0 + 9}]
830                    lappend vallist [string trim [string range $lineA $j0 $j1]]
831                    incr k
832                }
833            }
834            return $vallist
835        }
836        ODFterms-set {
837            set key "CRS$phase  OD    "
838            if {![existsexp $key]} {
839                makeexprec $key
840                set oldlen 0
841            } else {
842                set oldlen [string trim [string range [readexp $key] 5 9]]
843            }
844            set len [llength $value]
845            if ![validint len 5] {return 0}
846            setexp $key $len 6 5
847            set j 0
848            set k 0
849            foreach item $value {
850                incr j
851                if {$j % 6 == 1} {
852                    incr k
853                    if {$k <= 9} {
854                        set k "  $k"
855                    } elseif {$k <= 99} {
856                        set k " $k"
857                    }
858                    set col 1
859                    set keyA "CRS$phase  OD${k}A"
860                    set keyB "CRS$phase  OD${k}B"
861                    if {![existsexp $keyA]} {
862                        makeexprec $keyA
863                        makeexprec $keyB
864                    }
865                }
866                set col1 [expr {$col + 1}]
867                foreach n [lrange $item 0 2] {
868                    if ![validint n 3] {return 0}
869                    setexp $keyA $n $col1 3
870                    incr col1 3
871                }
872                incr col 10
873            }
874            for {incr j} {$j <= $oldlen} {incr j} {
875                if {$j % 6 == 1} {
876                    incr k
877                    if {$k <= 9} {
878                        set k "  $k"
879                    } elseif {$k <= 99} {
880                        set k " $k"
881                    }
882                    set col 1
883                    set keyA "CRS$phase  OD${k}A"
884                    set keyB "CRS$phase  OD${k}B"
885                    delexp $keyA
886                    delexp $keyB
887                }
888                if {[existsexp $keyA]} {
889                    setexp $keyA "          " $col 10
890                    setexp $keyB "          " $col 10
891                }
892                incr col 10
893            }
894        }
895        DistCalc-get {
896            set val [disagldat_get $phase]
897            return [lindex $val 0]
898        }
899        DistCalc-set {
900            set key "  DSGL CDAT$phase"
901            # for none delete the record & thats all folks
902            if {$value == "none"} {
903                catch {unset ::exparray($key)}
904                return
905            }
906            if {[existsexp $key] == 0} {
907                makeexprec $key
908            }
909            set line [readexp $key]
910            if {[string trim $line] == ""} {
911                # blank set to defaults
912                set line [string replace $line 2 15 "DRAD ARAD NOFO"]
913            }
914            if {$value == "radii"} {
915                if {[string range $line 2 5] == "DMAX"} {
916                    set line [string replace $line 2 13 "DRAD"]
917                } else {
918                    set line [string replace $line 2 5 "DRAD"]
919                }
920            } else {
921                if ![validreal value 8 2] {return 0}
922                if {[string range $line 2 5] == "DMAX"} {
923                    set line [string replace $line 6 13 $value]
924                } else {
925                    set line [string replace $line 2 5 "DMAX"]
926                    set line [string replace $line 6 6 "$value "]
927                }
928            }
929            setexp $key $line 0 68
930        }
931        AngCalc-get {
932            set val [disagldat_get $phase]
933            return [lindex $val 1]
934        }
935        AngCalc-set {
936            set key "  DSGL CDAT$phase"
937            # for none delete the record & thats all folks
938            if {$value == "none"} {
939                catch {unset ::exparray($key)}
940                return
941            }
942            if {[existsexp $key] == 0} {
943                makeexprec $key
944            }
945            set line [readexp $key]
946            if {[string trim $line] == ""} {
947                # blank set to defaults
948                set line [string replace $line 2 15 "DRAD ARAD NOFO"]
949            }
950            if {[string range $line 2 5] == "DMAX"} {
951                set i2 8
952            } else {
953                set i2 0
954            }
955            if {$value == "radii"} {
956                if {[string range $line [expr {$i2+7}] [expr {$i2+10}]] == "DAGL"} {
957                    set line [string replace $line [expr {$i2+7}] [expr {$i2+18}] "ARAD"]
958                } else {
959                    set line [string replace $line [expr {$i2+7}] [expr {$i2+10}] "ARAD"]
960                }
961            } else {
962                if ![validreal value 8 2] {return 0}
963                if {[string range $line [expr {$i2+7}] [expr {$i2+10}]] == "DAGL"} {
964                    set line [string replace $line [expr {$i2+11}] [expr {$i2+18}] $value]
965                } else {
966                    set line [string replace $line [expr {$i2+7}] [expr {$i2+10}] "DAGL"]
967                    set line [string replace $line [expr {$i2+11}] [expr {$i2+11}] "$value "]
968                }
969            }
970            setexp $key $line 0 68
971        }
972        default {
973            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
974            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
975        }
976    }
977    return 1
978}
979
980
981
982# get atom information: atominfo phase atom parm action value
983#   phase: 1 to 9 (as defined)
984#   atom: a valid atom number [see expmap(atomlist_$phase)]
985#      Note that atom and phase can be paired lists, but if there are extra
986#      entries in the atoms list, the last phase will be repeated.
987#      so that atominfo 1 {1 2 3} xset 1
988#               will set the xflag for atoms 1-3 in phase 1
989#      but atominfo {1 2 3} {1 1 1} xset 1
990#               will set the xflag for atoms 1 in phase 1-3
991#   parm:
992#     type -- element code
993#     mult -- atom multiplicity
994#     label -- atom label (*)
995#     x y z -- coordinates (*)
996#     frac --  occupancy (*)
997#     temptype -- I or A for Isotropic/Anisotropic (*)
998#     Uiso  -- Isotropic temperature factor (*)
999#     U11  -- Anisotropic temperature factor (*)
1000#     U22  -- Anisotropic temperature factor (*)
1001#     U33  -- Anisotropic temperature factor (*)
1002#     U12  -- Anisotropic temperature factor (*)
1003#     U13  -- Anisotropic temperature factor (*)
1004#     U23  -- Anisotropic temperature factor (*)
1005#     xref/xdamp -- refinement flag/damping value for the coordinates (*)
1006#     uref/udamp -- refinement flag/damping value for the temperature factor(s)  (*)
1007#     fref/fdamp -- refinement flag/damping value for the occupancy (*)
1008#  action: get (default) or set
1009#  value: used only with set
1010#  * =>  read+write supported
1011proc atominfo {phaselist atomlist parm "action get" "value {}"} {
1012    foreach phase $phaselist atom $atomlist {
1013        if {$phase == ""} {set phase [lindex $phaselist end]}
1014        if {$atom < 10} {
1015            set key "CRS$phase  AT  $atom"
1016        } elseif {$atom < 100} {
1017            set key "CRS$phase  AT $atom"
1018        } else {
1019            set key "CRS$phase  AT$atom"
1020        }
1021        switch -glob ${parm}-$action {
1022            type-get {
1023                return [string trim [string range [readexp ${key}A] 2 9] ]
1024            }
1025            mult-get {
1026                return [string trim [string range [readexp ${key}A] 58 61] ]
1027            }
1028            label-get {
1029                return [string trim [string range [readexp ${key}A] 50 57] ]
1030            }
1031            label-set {
1032                setexp ${key}A $value 51 8
1033            }
1034            temptype-get {
1035                return [string trim [string range [readexp ${key}B] 62 62] ]
1036            }
1037            temptype-set {
1038                if {$value == "A"} {
1039                    setexp ${key}B A 63 1
1040                    # copy the Uiso to the diagonal terms
1041                    set Uiso [string range [readexp ${key}B] 0 9]
1042                    foreach value [CalcAniso $phase $Uiso] \
1043                            col {1 11 21 31 41 51} {
1044                        validreal value 10 6
1045                        setexp ${key}B $value $col 10
1046                    }
1047                } else {
1048                    setexp ${key}B I 63 1
1049                    set value 0.0
1050                    catch {
1051                        # get the trace
1052                        set value [expr {( \
1053                                [string range [readexp ${key}B] 0 9] + \
1054                                [string range [readexp ${key}B] 10 19] + \
1055                                [string range [readexp ${key}B] 20 29])/3.}]
1056                    }
1057                    validreal value 10 6
1058                    setexp ${key}B $value 1 10
1059                    # blank out the remaining terms
1060                    set value " "
1061                    setexp ${key}B $value 11 10
1062                    setexp ${key}B $value 21 10
1063                    setexp ${key}B $value 31 10
1064                    setexp ${key}B $value 41 10
1065                    setexp ${key}B $value 51 10
1066                }
1067            }
1068            x-get {
1069                return [string trim [string range [readexp ${key}A] 10 19] ]
1070            }
1071            x-set {
1072                if ![validreal value 10 6] {return 0}
1073                setexp ${key}A $value 11 10
1074            }
1075            y-get {
1076                return [string trim [string range [readexp ${key}A] 20 29] ]
1077            }
1078            y-set {
1079                if ![validreal value 10 6] {return 0}
1080                setexp ${key}A $value 21 10
1081            }
1082            z-get {
1083                return [string trim [string range [readexp ${key}A] 30 39] ]
1084            }
1085            z-set {
1086                if ![validreal value 10 6] {return 0}
1087                setexp ${key}A $value 31 10
1088            }
1089            frac-get {
1090                return [string trim [string range [readexp ${key}A] 40 49] ]
1091            }
1092            frac-set {
1093                if ![validreal value 10 6] {return 0}
1094                setexp ${key}A $value 41 10
1095            }
1096            U*-get {
1097                regsub U $parm {} type
1098                if {$type == "iso" || $type == "11"} {
1099                    return [string trim [string range [readexp ${key}B] 0 9] ]
1100                } elseif {$type == "22"} {
1101                    return [string trim [string range [readexp ${key}B] 10 19] ]
1102                } elseif {$type == "33"} {
1103                    return [string trim [string range [readexp ${key}B] 20 29] ]
1104                } elseif {$type == "12"} {
1105                    return [string trim [string range [readexp ${key}B] 30 39] ]
1106                } elseif {$type == "13"} {
1107                    return [string trim [string range [readexp ${key}B] 40 49] ]
1108                } elseif {$type == "23"} {
1109                    return [string trim [string range [readexp ${key}B] 50 59] ]
1110                }
1111            }
1112            U*-set {
1113                if ![validreal value 10 6] {return 0}
1114                regsub U $parm {} type
1115                if {$type == "iso" || $type == "11"} {
1116                    setexp ${key}B $value 1 10
1117                } elseif {$type == "22"} {
1118                    setexp ${key}B $value 11 10
1119                } elseif {$type == "33"} {
1120                    setexp ${key}B $value 21 10
1121                } elseif {$type == "12"} {
1122                    setexp ${key}B $value 31 10
1123                } elseif {$type == "13"} {
1124                    setexp ${key}B $value 41 10
1125                } elseif {$type == "23"} {
1126                    setexp ${key}B $value 51 10
1127                }
1128            }
1129            xref-get {
1130                if {[string toupper [string range [readexp ${key}B] 64 64]] == "X"} {
1131                    return 1
1132                }
1133                return 0
1134            }
1135            xref-set {
1136                if $value {
1137                    setexp ${key}B "X" 65 1
1138                } else {
1139                    setexp ${key}B " " 65 1
1140                }           
1141            }
1142            xdamp-get {
1143                set val [string range [readexp ${key}A] 64 64]
1144                if {$val == " "} {return 0}
1145                return $val
1146            }
1147            xdamp-set {
1148                setexp ${key}A $value 65 1
1149            }
1150            fref-get {
1151                if {[string toupper [string range [readexp ${key}B] 63 63]] == "F"} {
1152                    return 1
1153                }
1154                return 0
1155            }
1156            fref-set {
1157                if $value {
1158                    setexp ${key}B "F" 64 1
1159                } else {
1160                    setexp ${key}B " " 64 1
1161                }           
1162            }
1163            fdamp-get {
1164                set val [string range [readexp ${key}A] 63 63]
1165                if {$val == " "} {return 0}
1166                return $val
1167            }
1168            fdamp-set {
1169                setexp ${key}A $value 64 1
1170            }
1171
1172            uref-get {
1173                if {[string toupper [string range [readexp ${key}B] 65 65]] == "U"} {
1174                    return 1
1175                }
1176                return 0
1177            }
1178            uref-set {
1179                if $value {
1180                    setexp ${key}B "U" 66 1
1181                } else {
1182                    setexp ${key}B " " 66 1
1183                }           
1184            }
1185            udamp-get {
1186                set val [string range [readexp ${key}A] 65 65]
1187                if {$val == " "} {return 0}
1188                return $val
1189            }
1190            udamp-set {
1191                setexp ${key}A $value 66 1
1192            }
1193            default {
1194                set msg "Unsupported atominfo access: parm=$parm action=$action"
1195                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1196            }
1197        }
1198    }
1199    return 1
1200}
1201
1202# get macromolecular atom information: mmatominfo phase atom parm action value
1203#   phase: 1 (at present only one mm phase can be defined)
1204#   atom: a valid atom number [see expmap(atomlist_$phase)]
1205#      Note that atoms can be lists
1206#      so that mmatominfo 1 {1 2 3} xset 1
1207#               will set the xflag for atoms 1-3 in phase 1
1208#   parm:
1209#     type -- element code
1210#     frac --  occupancy (*)
1211#     x y z -- coordinates (*)
1212#     Uiso  -- Isotropic temperature factor (*)
1213#     label -- atom label (*)
1214#     residue -- residue label (*)
1215#     group -- group label (*)
1216#     resnum -- residue number (*)
1217#     xref/xdamp -- refinement flag/damping value for the coordinates (*)
1218#     uref/udamp -- refinement flag/damping value for the temperature factor(s)  (*)
1219#     fref/fdamp -- refinement flag/damping value for the occupancy (*)
1220#  action: get (default) or set
1221#  value: used only with set
1222#  * =>  read+write supported
1223proc mmatominfo {phaselist atomlist parm "action get" "value {}"} {
1224    foreach phase $phaselist atom $atomlist {
1225        if {$phase == ""} {set phase [lindex $phaselist end]}
1226        set num [string toupper [format %.4x $atom]]
1227        set key "CRS$phase  AT$num"
1228        switch -glob ${parm}-$action {
1229            type-get {
1230                return [string trim [string range [readexp ${key}] 2 9] ]
1231            }
1232            frac-get {
1233                return [string trim [string range [readexp ${key}] 10 15] ]
1234            }
1235            frac-set {
1236                if ![validreal value 6 4] {return 0}
1237                setexp ${key} $value 11 6
1238            }
1239            x-get {
1240                return [string trim [string range [readexp ${key}] 16 23] ]
1241            }
1242            x-set {
1243                if ![validreal value 8 5] {return 0}
1244                setexp ${key} $value 17 8
1245            }
1246            y-get {
1247                return [string trim [string range [readexp ${key}] 24 31] ]
1248            }
1249            y-set {
1250                if ![validreal value 8 5] {return 0}
1251                setexp ${key} $value 25 8
1252            }
1253            z-get {
1254                return [string trim [string range [readexp ${key}] 32 39] ]
1255            }
1256            z-set {
1257                if ![validreal value 8 5] {return 0}
1258                setexp ${key} $value 33 8
1259            }
1260            Uiso-get {
1261                return [string trim [string range [readexp ${key}] 40 45] ]
1262            }
1263            Uiso-set {
1264                if ![validreal value 6 4] {return 0}
1265                setexp ${key} $value 41 6
1266            }
1267            label-get {
1268                return [string trim [string range [readexp ${key}] 46 50] ]
1269            }
1270            label-set {
1271                setexp ${key} $value 47 5
1272            }
1273            residue-get {
1274                return [string range [readexp ${key}] 51 53]
1275            }
1276            residue-set {
1277                setexp ${key} $value 52 3
1278            }
1279            group-get {
1280                return [string range [readexp ${key}] 54 55]
1281            }
1282            group-set {
1283                setexp ${key} $value 55 2
1284            }
1285            resnum-get {
1286                return [string trim [string range [readexp ${key}] 56 59] ]
1287            }
1288            resnum-set {
1289                if ![validint value 4] {return 0}
1290                setexp "${key} EPHAS" $value 57 4
1291            }
1292            fref-get {
1293                if {[string toupper [string range [readexp $key] 60 60]] == "F"} {
1294                    return 1
1295                }
1296                return 0
1297            }
1298            fref-set {
1299                if $value {
1300                    setexp $key "F" 61 1
1301                } else {
1302                    setexp $key " " 61 1
1303                }           
1304            }
1305            xref-get {
1306                if {[string toupper [string range [readexp $key] 61 61]] == "X"} {
1307                    return 1
1308                }
1309                return 0
1310            }
1311            xref-set {
1312                if $value {
1313                    setexp $key "X" 62 1
1314                } else {
1315                    setexp ${key}B " " 62 1
1316                }           
1317            }
1318            uref-get {
1319                if {[string toupper [string range [readexp $key] 62 62]] == "U"} {
1320                    return 1
1321                }
1322                return 0
1323            }
1324            uref-set {
1325                if $value {
1326                    setexp $key "U" 63 1
1327                } else {
1328                    setexp $key " " 63 1
1329                }           
1330            }
1331
1332            fdamp-get {
1333                set val [string range [readexp ${key}] 63 63]
1334                if {$val == " "} {return 0}
1335                return $val
1336            }
1337            fdamp-set {
1338                setexp ${key} $value 64 1
1339            }
1340            xdamp-get {
1341                set val [string range [readexp ${key}] 64 64]
1342                if {$val == " "} {return 0}
1343                return $val
1344            }
1345            xdamp-set {
1346                setexp ${key} $value 65 1
1347            }
1348
1349            udamp-get {
1350                set val [string range [readexp ${key}] 65 65]
1351                if {$val == " "} {return 0}
1352                return $val
1353            }
1354            udamp-set {
1355                setexp ${key} $value 66 1
1356            }
1357            default {
1358                set msg "Unsupported mmatominfo access: parm=$parm action=$action"
1359                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
1360            }
1361        }
1362    }
1363    return 1
1364}
1365
1366
1367
1368# get histogram information: histinfo histlist parm action value
1369# histlist is a list of histogram numbers
1370# parm:
1371#     title
1372#     file  -- file name of raw data for histogram (*)
1373#     scale (*)
1374#     sref/sdamp -- refinement flag/damping value for the scale factor (*)
1375#     lam1, lam2 (*)
1376#     ttref refinement flag for the 2theta (ED Xray) (*)
1377#     wref refinement flag for the wavelength (*)
1378#     ratref refinement flag for the wavelength ratio (*)
1379#     difc, difa -- TOF calibration constants (*)
1380#     dcref,daref -- refinement flag for difc, difa (*)
1381#     zero (*)
1382#     zref refinement flag for the zero correction (*)
1383#     ipola (*)
1384#     pola (*)
1385#     pref refinement flag for the polarization (*)
1386#     kratio (*)
1387#     ddamp -- damping value for the diffractometer constants (*)
1388#     backtype -- background function number *
1389#     backterms -- number of background terms *
1390#     bref/bdamp -- refinement flag/damping value for the background (*)
1391#     bterm$n -- background term #n (*)
1392#     bank -- Bank number
1393#     tofangle -- detector angle (TOF only)
1394#     foextract  -- Fobs extraction flag (*)
1395#     LBdamp  -- LeBail damping value (*)
1396#     tmin/tmax -- minimum & maximum usable 2theta/TOF/energy
1397#     excl -- excluded regions (*)
1398#     dmin -- minimum d-space for reflection generation (*)
1399#     use  -- use flag; 1 = use; 0 = do not use (*)
1400#     dstart -- dummy histogram starting tmin/emin/2theta (*)
1401#     dstep -- dummy histogram step size tmin/emin/2theta (*)
1402#     dpoints -- dummy histogram number of points (*)
1403#     dtype   -- dummy histogram type (CONST or SLOG)
1404#     abscor1 -- 1st absorption correction (*)
1405#     abscor2 -- 2nd absorption correction (*)
1406#     abstype -- absorption correction type (*)
1407#     absdamp -- damping for absorption refinement (*)
1408#     absref -- refinement damping for absorption refinement (*)
1409#     proftbl -- returns number of profile table terms
1410#     anomff -- returns a list of elements, f' and f"
1411#   parameters transferred from the instrument parameter file:
1412#     ITYP    -- returns the contents of the ITYP record
1413proc histinfo {histlist parm "action get" "value {}"} {
1414    global expgui
1415    foreach hist $histlist {
1416        if {$hist < 10} {
1417            set key "HST  $hist"
1418        } else {
1419            set key "HST $hist"
1420        }
1421        switch -glob ${parm}-$action {
1422            foextract-get {
1423                set line [readexp "${key} EPHAS"]
1424                # add a EPHAS if not exists
1425                if {$line == {}} {
1426                    makeexprec "${key} EPHAS"
1427                    # expedt defaults this to "F", but I think "T" is better
1428                    setexp "${key} EPHAS" "Y" 50 1
1429                    if $expgui(debug) {puts "Warning: creating a ${key} EPHAS record"}
1430                }
1431                if {[string toupper [string range $line 49 49]] == "T"} {
1432                    return 1
1433                }
1434                # the flag has changed to "Y/N" in the latest versions
1435                # of GSAS
1436                if {[string toupper [string range $line 49 49]] == "Y"} {
1437                    return 1
1438                }
1439                return 0
1440            }
1441            foextract-set {
1442                # the flag has changed to "Y/N" in the latest versions
1443                # of GSAS
1444                if $value {
1445                    setexp "${key} EPHAS" "Y" 50 1
1446                } else {
1447                    setexp "${key} EPHAS" "N" 50 1
1448                }
1449            }
1450            LBdamp-get {
1451                set v [string trim [string range [readexp "${key} EPHAS"] 54 54]]
1452                if {$v == ""} {return 0}
1453                return $v
1454            }
1455            LBdamp-set {
1456                if ![validint value 5] {return 0}
1457                setexp "${key} EPHAS" $value 51 5
1458            }
1459            title-get {
1460                return [string trim [readexp "${key}  HNAM"] ]
1461            }
1462            scale-get {
1463                return [string trim [string range [readexp ${key}HSCALE] 0 14]]
1464            }
1465            scale-set {
1466                if ![validreal value 15 6] {return 0}
1467                setexp ${key}HSCALE $value 1 15
1468            }
1469            sref-get {
1470                if {[string toupper [string range [readexp ${key}HSCALE] 19 19]] == "Y"} {
1471                    return 1
1472                }
1473                return 0
1474            }
1475            sref-set {
1476                if $value {
1477                    setexp ${key}HSCALE "Y" 20 1
1478                } else {
1479                    setexp ${key}HSCALE "N" 20 1
1480                }           
1481            }
1482            sdamp-get {
1483                set val [string range [readexp ${key}HSCALE] 24 24]
1484                if {$val == " "} {return 0}
1485                return $val
1486            }
1487            sdamp-set {
1488                setexp ${key}HSCALE $value 25 1
1489            }
1490
1491            difc-get -
1492            lam1-get {
1493                return [string trim [string range [readexp "${key} ICONS"] 0 9]]
1494            }
1495            difc-set -
1496            lam1-set {
1497                if ![validreal value 10 7] {return 0}
1498                setexp "${key} ICONS" $value 1 10
1499                # set the powpref warning (1 = suggested)
1500                catch {
1501                    global expgui
1502                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1503                    set msg "Diffractometer constants" 
1504                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1505                        append expgui(needpowpref_why) "\t$msg were changed\n"
1506                    }
1507                }
1508            }
1509            difa-get -
1510            lam2-get {
1511                return [string trim [string range [readexp "${key} ICONS"] 10 19]]
1512            }
1513            difa-set -
1514            lam2-set {
1515                if ![validreal value 10 7] {return 0}
1516                setexp "${key} ICONS" $value 11 10
1517                # set the powpref warning (1 = suggested)
1518                catch {
1519                    global expgui
1520                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1521                    set msg "Diffractometer constants" 
1522                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1523                        append expgui(needpowpref_why) "\t$msg were changed\n"
1524                    }
1525                }
1526            }
1527            zero-get {
1528                return [string trim [string range [readexp "${key} ICONS"] 20 29]]
1529            }
1530            zero-set {
1531                if ![validreal value 10 5] {return 0}
1532                setexp "${key} ICONS" $value 21 10
1533                # set the powpref warning (1 = suggested)
1534                catch {
1535                    global expgui
1536                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1537                    set msg "Diffractometer constants" 
1538                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1539                        append expgui(needpowpref_why) "\t$msg were changed\n"
1540                    }
1541                }
1542            }
1543            ipola-get {
1544                return [string trim [string range [readexp "${key} ICONS"] 54 54]]
1545            }
1546            ipola-set {
1547                if ![validint value 1] {return 0}
1548                setexp "${key} ICONS" $value 55 1
1549            }
1550            pola-get {
1551                return [string trim [string range [readexp "${key} ICONS"] 40 49]]
1552            }
1553            pola-set {
1554                if ![validreal value 10 5] {return 0}
1555                setexp "${key} ICONS" $value 41 10
1556            }
1557            kratio-get {
1558                set val [string trim [string range [readexp "${key} ICONS"] 55 64]]
1559                if {$val == ""} {set val 0}
1560                # N.B. this code is used w/CW, where Kratio may not be 0.0
1561                set lam2 [string trim [string range [readexp "${key} ICONS"] 10 19]]
1562                if {$lam2 == ""} {set lam2 0}
1563                # Change kratio & flag the change (this is rather kludged)
1564                if {$val == 0 && $lam2 != 0} {
1565                    set val 0.5
1566                    validreal val 10 5
1567                    setexp "${key} ICONS" $val 56 10
1568                    catch {incr ::expgui(changed)}
1569                }
1570                return $val
1571            }
1572            kratio-set {
1573                if ![validreal value 10 5] {return 0}
1574                setexp "${key} ICONS" $value 56 10
1575            }
1576
1577            wref-get {
1578            #------------------------------------------------------
1579            # col 33: refine flag for lambda, difc, ratio and theta
1580            #------------------------------------------------------
1581                if {[string toupper [string range \
1582                        [readexp "${key} ICONS"] 32 32]] == "L"} {
1583                    return 1
1584                }
1585                return 0
1586            }
1587            wref-set {
1588                if $value {
1589                    setexp "${key} ICONS" "L" 33 1
1590                } else {
1591                    setexp "${key} ICONS" " " 33 1
1592                }           
1593            }
1594            ratref-get {
1595                if {[string toupper [string range \
1596                        [readexp "${key} ICONS"] 32 32]] == "R"} {
1597                    return 1
1598                }
1599                return 0
1600            }
1601            ratref-set {
1602                if $value {
1603                    setexp "${key} ICONS" "R" 33 1
1604                } else {
1605                    setexp "${key} ICONS" " " 33 1
1606                }           
1607            }
1608            dcref-get {
1609                if {[string toupper [string range \
1610                        [readexp "${key} ICONS"] 32 32]] == "C"} {
1611                    return 1
1612                }
1613                return 0
1614            }
1615            dcref-set {
1616                if $value {
1617                    setexp "${key} ICONS" "C" 33 1
1618                } else {
1619                    setexp "${key} ICONS" " " 33 1
1620                }           
1621            }
1622            ttref-get {
1623                if {[string toupper [string range \
1624                        [readexp "${key} ICONS"] 32 32]] == "T"} {
1625                    return 1
1626                }
1627                return 0
1628            }
1629            ttref-set {
1630                if $value {
1631                    setexp "${key} ICONS" "T" 33 1
1632                } else {
1633                    setexp "${key} ICONS" " " 33 1
1634                }           
1635            }
1636
1637
1638            pref-get {
1639            #------------------------------------------------------
1640            # col 34: refine flag for POLA & DIFA
1641            #------------------------------------------------------
1642                if {[string toupper [string range \
1643                        [readexp "${key} ICONS"] 33 33]] == "P"} {
1644                    return 1
1645                }
1646                return 0
1647            }
1648            pref-set {
1649                if $value {
1650                    setexp "${key} ICONS" "P" 34 1
1651                } else {
1652                    setexp "${key} ICONS" " " 34 1
1653                }           
1654            }
1655            daref-get {
1656                if {[string toupper [string range \
1657                        [readexp "${key} ICONS"] 33 33]] == "A"} {
1658                    return 1
1659                }
1660                return 0
1661            }
1662            daref-set {
1663                if $value {
1664                    setexp "${key} ICONS" "A" 34 1
1665                } else {
1666                    setexp "${key} ICONS" " " 34 1
1667                }           
1668            }
1669
1670            zref-get {
1671            #------------------------------------------------------
1672            # col 34: refine flag for zero correction
1673            #------------------------------------------------------
1674                if {[string toupper [string range [readexp "${key} ICONS"] 34 34]] == "Z"} {
1675                    return 1
1676                }
1677                return 0
1678            }
1679            zref-set {
1680                if $value {
1681                    setexp "${key} ICONS" "Z" 35 1
1682                } else {
1683                    setexp "${key} ICONS" " " 35 1
1684                }           
1685            }
1686
1687            ddamp-get {
1688                set val [string range [readexp "${key} ICONS"] 39 39]
1689                if {$val == " "} {return 0}
1690                return $val
1691            }
1692            ddamp-set {
1693                setexp "${key} ICONS" $value 40 1
1694            }
1695
1696            backtype-get {
1697                set val [string trim [string range [readexp "${key}BAKGD "] 0 4]]
1698                if {$val == " "} {return 0}
1699                return $val
1700            }
1701            backtype-set {
1702                if ![validint value 5] {return 0}
1703                setexp "${key}BAKGD " $value 1 5
1704            }
1705            backterms-get {
1706                set val [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1707                if {$val == " "} {return 0}
1708                return $val
1709            }
1710            backterms-set {
1711                # this takes a bit of work -- if terms are added, add lines as needed to the .EXP
1712                set oldval [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1713                if ![validint value 5] {return 0}
1714                if {$oldval < $value} {
1715                    set line1  [expr {2 + ($oldval - 1) / 4}]
1716                    set line2  [expr {1 + ($value - 1) / 4}]
1717                    for {set i $line1} {$i <= $line2} {incr i} {
1718                        # create a blank entry if needed
1719                        makeexprec ${key}BAKGD$i
1720                    }
1721                    incr oldval
1722                    for {set num $oldval} {$num <= $value} {incr num} {
1723                        set f1 [expr {15*(($num - 1) % 4)}]
1724                        set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
1725                        set line  [expr {1 + ($num - 1) / 4}]
1726                        if {[string trim [string range [readexp ${key}BAKGD$line] $f1 $f2]] == ""} {
1727                            set f1 [expr {15*(($num - 1) % 4)+1}]
1728                            setexp ${key}BAKGD$line 0.0 $f1 15                 
1729                        }
1730                    }
1731                }
1732                setexp "${key}BAKGD " $value 6 5
1733
1734            }
1735            bref-get {
1736                if {[string toupper [string range [readexp "${key}BAKGD"] 14 14]] == "Y"} {
1737                    return 1
1738                }
1739                return 0
1740            }
1741            bref-set {
1742                if $value {
1743                    setexp "${key}BAKGD "  "Y" 15 1
1744                } else {
1745                    setexp "${key}BAKGD "  "N" 15 1
1746                }
1747            }
1748            bdamp-get {
1749                set val [string range [readexp "${key}BAKGD "] 19 19]
1750                if {$val == " "} {return 0}
1751                return $val
1752            }
1753            bdamp-set {
1754                setexp "${key}BAKGD " $value 20 1
1755            }
1756            bterm*-get {
1757                regsub bterm $parm {} num
1758                set f1 [expr {15*(($num - 1) % 4)}]
1759                set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
1760                set line  [expr {1 + ($num - 1) / 4}]
1761                return [string trim [string range [readexp ${key}BAKGD$line] $f1 $f2] ]
1762            }
1763            bterm*-set {
1764                regsub bterm $parm {} num
1765                if ![validreal value 15 6] {return 0}
1766                set f1 [expr {15*(($num - 1) % 4)+1}]
1767                set line  [expr {1 + ($num - 1) / 4}]
1768                setexp ${key}BAKGD$line $value $f1 15
1769            }
1770            bank-get {
1771                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1772            }
1773            tofangle-get {
1774                return [string trim [string range [readexp "${key}BNKPAR"] 10 19]]
1775            }
1776            tmin-get {
1777                return [string trim [string range [readexp "${key} TRNGE"] 0 9]]
1778            }
1779            tmax-get {
1780                return [string trim [string range [readexp "${key} TRNGE"] 10 19]]
1781            }
1782            excl-get {
1783                set n [string trim [string range [readexp "${key} NEXC"] 0 4]]
1784                set exlist {}
1785                for {set i 1} {$i <= $n} {incr i} {
1786                    set line [readexp [format "${key}EXC%3d" $i]]
1787                    lappend exlist [list \
1788                            [string trim [string range $line  0  9]] \
1789                            [string trim [string range $line 10 19]]]
1790                }
1791                return $exlist
1792            }
1793            excl-set {
1794                set n [llength $value]
1795                if ![validint n 5] {return 0}
1796                setexp "${key} NEXC" $n 1 5
1797                set i 0
1798                foreach p $value {
1799                    incr i
1800                    foreach {r1 r2} $p {}
1801                    validreal r1 10 3
1802                    validreal r2 10 3
1803                    set k [format "${key}EXC%3d" $i]
1804                    if {![existsexp $k]} {
1805                        makeexprec $k
1806                    }
1807                    setexp $k ${r1}${r2} 1 20
1808                }
1809                # set the powpref warning (2 = required)
1810                catch {
1811                    global expgui
1812                    set expgui(needpowpref) 2
1813                    set msg "Excluded regions" 
1814                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1815                        append expgui(needpowpref_why) "\t$msg were changed\n"
1816                    }
1817                }
1818            }
1819            file-get {
1820                return [string trim [readexp "${key}  HFIL"] ]
1821            }
1822            file-set {
1823                setexp "${key}  HFIL" $value 3 65
1824            }
1825            bank-get {
1826                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1827            }
1828            dmin-get {
1829                return [string trim [string range [readexp "${key} NREF"] 5 14]]
1830            }
1831            dmin-set {
1832                if ![validreal value 10 4] {return 0}
1833                setexp "${key} NREF" $value 6 10
1834                # set the powpref warning (2 = required)
1835                catch {
1836                    global expgui
1837                    set expgui(needpowpref) 2
1838                    set msg "Dmin (reflection range)" 
1839                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1840                        append expgui(needpowpref_why) "\t$msg was changed\n"
1841                    }
1842                }
1843            }
1844            use-get {
1845                set k [expr {($hist+11)/12}]
1846                set line [readexp " EXPR  HTYP$k"]
1847                set j [expr {((($hist-1) % 12)+1)*5}]
1848                if {[string range $line $j $j] == "*"} {return 0}
1849                return 1
1850            }
1851            use-set {
1852                set k [expr {($hist+11)/12}]
1853                set line [readexp " EXPR  HTYP$k"]
1854                set j [expr {((($hist-1) % 12)+1)*5+1}]
1855                if {$value} {
1856                    setexp " EXPR  HTYP$k" " " $j 1
1857                } else {
1858                    setexp " EXPR  HTYP$k" "*" $j 1
1859                }
1860                # set the powpref warning (2 = required)
1861                catch {
1862                    global expgui
1863                    set expgui(needpowpref) 2
1864                    set msg "Histogram use flags" 
1865                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1866                        append expgui(needpowpref_why) "\t$msg were changed\n"
1867                    }
1868                }
1869            }
1870            dstart-get {
1871                return [string trim [string range [readexp "${key} DUMMY"] 20 29]]
1872            }
1873            dstart-set {
1874                if ![validreal value 10 3] {return 0}
1875                setexp "${key} DUMMY" $value 21 10
1876                # set the powpref warning (1 = suggested)
1877                catch {
1878                    global expgui
1879                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1880                    set msg "Dummy histogram parameters" 
1881                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1882                        append expgui(needpowpref_why) "\t$msg were changed\n"
1883                    }
1884                }
1885            }
1886            dstep-get {
1887                return [string trim [string range [readexp "${key} DUMMY"] 30 39]]
1888            }
1889            dstep-set {
1890                if ![validreal value 10 3] {return 0}
1891                setexp "${key} DUMMY" $value 31 10
1892                catch {
1893                    global expgui
1894                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1895                    set msg "Dummy histogram parameters" 
1896                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1897                        append expgui(needpowpref_why) "\t$msg were changed\n"
1898                    }
1899                }
1900            }
1901            dpoints-get {
1902                return [string trim [string range [readexp "${key} DUMMY"] 0 9]]
1903            }
1904            dpoints-set {
1905                if ![validint value 10] {return 0}
1906                setexp "${key} DUMMY" $value 1 10
1907                catch {
1908                    global expgui
1909                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1910                    set msg "Dummy histogram parameters" 
1911                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1912                        append expgui(needpowpref_why) "\t$msg were changed\n"
1913                    }
1914                }
1915            }
1916            dtype-get {
1917                return [string trim [string range [readexp "${key} DUMMY"] 10 19]]
1918            }
1919            abscor1-get {
1920                return [string trim [string range [readexp "${key}ABSCOR"] 0 14]]
1921            }
1922            abscor1-set {
1923                if ![validreal value 15 7] {return 0}
1924                setexp "${key}ABSCOR" $value 1 15
1925            }
1926            abscor2-get {
1927                return [string trim [string range [readexp "${key}ABSCOR"] 15 29]]
1928            }
1929            abscor2-set {
1930                # this must have a decimal as the 5th character, so that we end up with a
1931                # decimal point in column 20.
1932                set tmp $value
1933                if ![validreal tmp 12 7] {return 0}
1934                set pos [string first "." $tmp]
1935                while {$pos < 4} {
1936                    set tmp " $tmp"
1937                    set pos [string first "." $tmp]
1938                }
1939                if {$pos == 4} {
1940                    setexp "${key}ABSCOR" $tmp 16 15
1941                    return
1942                }
1943                catch {
1944                    set tmp [format "%12.6E" $value]
1945                    set pos [string first "." $tmp]
1946                    while {$pos < 4} {
1947                        set tmp " $tmp"
1948                        set pos [string first "." $tmp]
1949                    }
1950                    if {$pos == 4} {
1951                        setexp "${key}ABSCOR" $tmp 16 15
1952                        return
1953                    }
1954                }
1955                return 0
1956            }
1957            abstype-get {
1958                set val [string trim [string range [readexp "${key}ABSCOR"] 40 44]]
1959                if {$val == ""} {set val 0}
1960                return $val
1961            }
1962            abstype-set {
1963                if ![validint value 5] {return 0}
1964                setexp "${key}ABSCOR" $value 41 5
1965            }
1966            absdamp-get {
1967                set val [string range [readexp "${key}ABSCOR"] 39 39]
1968                if {$val == " "} {return 0}
1969                return $val
1970            }
1971            absdamp-set {
1972                if ![validint value 5] {return 0}
1973                setexp "${key}ABSCOR" $value 36 5
1974            }
1975            absref-get {
1976                if {[string toupper \
1977                        [string range [readexp "${key}ABSCOR"] 34 34]] == "Y"} {
1978                    return 1
1979                }
1980                return 0
1981            }
1982            absref-set {
1983                if $value {
1984                    setexp "${key}ABSCOR" "    Y" 31 5
1985                } else {
1986                    setexp "${key}ABSCOR" "    N" 31 5
1987                }
1988            }
1989            ITYP-get {
1990                return [string trim [readexp "${key}I ITYP"]]
1991            }
1992            proftbl-get {
1993                set line [readexp "${key}PAB3"]
1994                if {$line == ""} {return 0}
1995                set val [string trim [string range $line 0 4]]
1996                if {$val == ""} {return 0}
1997                return $val
1998            }
1999            anomff-get {
2000                set l {}
2001                foreach i {1 2 3 4 5 6 7 8 9} {
2002                    if {![existsexp "${key}FFANS$i"]} continue
2003                    set line [readexp "${key}FFANS$i"]
2004                    set elem [string trim [string range $line 2 9]]
2005                    set fp [string trim [string range $line 10 19]]
2006                    set fpp [string trim [string range $line 20 29]]
2007                    lappend l [list $elem $fp $fpp]
2008                }
2009                return $l
2010            }
2011            anomff-set {
2012                # match up input against elements in list.
2013                # change elements included, return any elements that are
2014                # not found.
2015                set errorlist {}
2016                foreach triplet $value {
2017                    foreach {e fp fpp} $triplet {}               
2018                    foreach i {1 2 3 4 5 6 7 8 9} {
2019                        if {![existsexp "${key}FFANS$i"]} continue
2020                        # note that the element name is not used or validated
2021                        set elem [string trim [string range \
2022                                                   [readexp "${key}FFANS$i"] 2 9]]
2023                        if {[string match -nocase $e $elem]} { 
2024                            if ![validreal fp 10 3] {return 0}
2025                            setexp "${key}FFANS$i" $fp 11 10
2026                            if ![validreal fpp 10 3] {return 0}
2027                            setexp "${key}FFANS$i" $fpp 21 10
2028                            set e {}
2029                            break
2030                        }
2031                    }
2032                    if {$e != ""} {
2033                        # oops, no match
2034                        lappend errorlist $e
2035                    }
2036                }
2037                if {$errorlist != ""} {return [list 0 $errorlist]}
2038            }
2039            default {
2040                set msg "Unsupported histinfo access: parm=$parm action=$action"
2041                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2042            }
2043        }
2044    }
2045    return 1
2046}
2047
2048# read the information that differs by both histogram and phase (profile & phase fraction)
2049# use: hapinfo hist phase parm action value
2050
2051#     frac -- phase fraction (*)
2052#     frref/frdamp -- refinement flag/damping value for the phase fraction (*)
2053#     proftype -- profile function number (*)
2054#     profterms -- number of profile terms (*)
2055#     pdamp -- damping value for the profile (*)
2056#     pcut -- cutoff value for the profile (*)
2057#     pterm$n -- profile term #n (*)
2058#     pref$n -- refinement flag value for profile term #n (*)
2059#     extmeth -- Fobs extraction method (*)
2060#     POnaxis -- number of defined M-D preferred axes
2061proc hapinfo {histlist phaselist parm "action get" "value {}"} {
2062    foreach phase $phaselist hist $histlist {
2063        if {$phase == ""} {set phase [lindex $phaselist end]}
2064        if {$hist == ""} {set hist [lindex $histlist end]}
2065        if {$hist < 10} {
2066            set hist " $hist"
2067        }
2068        set key "HAP${phase}${hist}"
2069        switch -glob ${parm}-$action {
2070            extmeth-get {
2071                set i1 [expr {($phase - 1)*5}]
2072                set i2 [expr {$i1 + 4}]
2073                return [string trim [string range [readexp "HST $hist EPHAS"] $i1 $i2]]
2074            }
2075            extmeth-set {
2076                set i1 [expr {($phase - 1)*5 + 1}]
2077                if ![validint value 5] {return 0}
2078                setexp "HST $hist EPHAS" $value $i1 5
2079            }
2080            frac-get {
2081                return [string trim [string range [readexp ${key}PHSFR] 0 14]]
2082            }
2083            frac-set {
2084                if ![validreal value 15 6] {return 0}
2085                setexp ${key}PHSFR $value 1 15
2086            }
2087            frref-get {
2088                if {[string toupper [string range [readexp ${key}PHSFR] 19 19]] == "Y"} {
2089                    return 1
2090                }
2091                return 0
2092            }
2093            frref-set {
2094                if $value {
2095                    setexp ${key}PHSFR "Y" 20 1
2096                } else {
2097                    setexp ${key}PHSFR "N" 20 1
2098                }           
2099            }
2100            frdamp-get {
2101                set val [string range [readexp ${key}PHSFR] 24 24]
2102                if {$val == " "} {return 0}
2103                return $val
2104            }
2105            frdamp-set {
2106                setexp ${key}PHSFR $value 25 1
2107            }
2108            proftype-get {
2109                set val [string range [readexp "${key}PRCF "] 0 4]
2110                if {$val == " "} {return 0}
2111                return $val
2112            }
2113            proftype-set {
2114                if ![validint value 5] {return 0}
2115                setexp "${key}PRCF " $value 1 5
2116                # set the powpref warning (1 = suggested)
2117                catch {
2118                    global expgui
2119                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2120                    set msg "Profile parameters" 
2121                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2122                        append expgui(needpowpref_why) "\t$msg were changed\n"
2123                    }
2124                }
2125            }
2126            profterms-get {
2127                set val [string range [readexp "${key}PRCF "] 5 9]
2128                if {$val == " "} {return 0}
2129                return $val
2130            }
2131            profterms-set {
2132                if ![validint value 5] {return 0}
2133                setexp "${key}PRCF " $value 6 5
2134                # now check that all needed entries exist
2135                set lines [expr {1 + ($value - 1) / 4}]
2136                for {set i 1} {$i <= $lines} {incr i} {
2137                    makeexprec "${key}PRCF $i"
2138                }
2139                # set the powpref warning (1 = suggested)
2140                catch {
2141                    global expgui
2142                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2143                    set msg "Profile parameters" 
2144                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2145                        append expgui(needpowpref_why) "\t$msg were changed\n"
2146                    }
2147                }
2148            }
2149            pcut-get {
2150                return [string trim [string range [readexp "${key}PRCF "] 10 19]]
2151            }
2152            pcut-set {
2153                if ![validreal value 10 5] {return 0}
2154                setexp "${key}PRCF " $value 11 10
2155                # set the powpref warning (1 = suggested)
2156                catch {
2157                    global expgui
2158                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2159                    set msg "Profile parameters" 
2160                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2161                        append expgui(needpowpref_why) "\t$msg were changed\n"
2162                    }
2163                }
2164            }
2165            pdamp-get {
2166                set val [string range [readexp "${key}PRCF "] 24 24]
2167                if {$val == " "} {return 0}
2168                return $val
2169            }
2170            pdamp-set {
2171                setexp "${key}PRCF   " $value 25 1
2172            }
2173            pterm*-get {
2174                regsub pterm $parm {} num
2175                set f1 [expr {15*(($num - 1) % 4)}]
2176                set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
2177                set line  [expr {1 + ($num - 1) / 4}]
2178                return [string trim [string range [readexp "${key}PRCF $line"] $f1 $f2] ]
2179            }
2180            pterm*-set {
2181                if ![validreal value 15 6] {return 0}
2182                regsub pterm $parm {} num
2183                set f1 [expr {1+ 15*(($num - 1) % 4)}]
2184                set line  [expr {1 + ($num - 1) / 4}]
2185                setexp "${key}PRCF $line" $value $f1 15
2186                # set the powpref warning (1 = suggested)
2187                catch {
2188                    global expgui
2189                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2190                    set msg "Profile parameters" 
2191                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2192                        append expgui(needpowpref_why) "\t$msg were changed\n"
2193                    }
2194                }
2195            }
2196            pref*-get {
2197                regsub pref $parm {} num
2198                set f [expr {24+$num}]
2199                if {[string toupper [string range [readexp "${key}PRCF  "] $f $f]] == "Y"} {
2200                    return 1
2201                }
2202                return 0
2203            }
2204            pref*-set {
2205                regsub pref $parm {} num
2206                set f [expr {25+$num}]
2207                if $value {
2208                    setexp ${key}PRCF "Y" $f 1
2209                } else {
2210                    setexp ${key}PRCF "N" $f 1
2211                }           
2212            }
2213            POnaxis-get {
2214                set val [string trim \
2215                        [string range [readexp "${key}NAXIS"] 0 4]]
2216                if {$val == ""} {return 0}
2217                return $val
2218            }
2219            POnaxis-set {
2220                if ![validint value 5] {return 0}
2221                # there should be a NAXIS record, but if not make one
2222                if {![existsexp "${key}NAXIS"]} {
2223                    makeexprec "${key}NAXIS"
2224                }
2225                setexp "${key}NAXIS  " $value 1 5
2226            }
2227            default {
2228                set msg "Unsupported hapinfo access: parm=$parm action=$action"
2229                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2230            }
2231        }
2232    }
2233    return 1
2234}
2235
2236#  read fixed constraints
2237
2238proc atom_constraint_read {phase} {
2239    set fix_list ""
2240    foreach k {1 2 3 4 5 6 7 8 9} {
2241        set key [format "LEQV HOLD%1d%2d" $phase $k]
2242        set line [readexp $key]
2243        foreach j {2 10 18 26 34 42 50 58} {
2244            set fix_param [string range $line $j [expr $j+7]]
2245            if {[string trim $fix_param] == ""} {return $fix_list}
2246            lappend fix_list $fix_param
2247        }
2248    }
2249}
2250
2251# load all atom constraints into global array fix_param
2252proc atom_constraint_load { } {
2253    catch {unset ::fix_param}
2254    foreach i $::expmap(phaselist) {
2255        set temp [atom_constraint_read $i]
2256        foreach j $temp {
2257            set atomnum [string trim [string range $j 2 3]]
2258            set param [string trim [string range $j 4 6]]
2259            set ::fix_param($i,$atomnum,$param) 1   
2260        }
2261    }
2262}
2263
2264proc atom_constraint_write {phase fix_list} {
2265    foreach key [array names ::exparray "LEQV HOLD$phase*"] {
2266        delexp $key
2267    }
2268    set k 0
2269    set j 1
2270    set line ""
2271    foreach fix $fix_list {
2272        incr k
2273        append line $fix
2274        if {$k == 8} {
2275            set key [format "LEQV HOLD%1d%2d" $phase $j]
2276            makeexprec $key
2277            setexp $key $line 3 [expr ($k * 8) + 2]
2278            set k 0
2279            incr j
2280            set line ""
2281        }
2282    }
2283    if {$line != ""} {
2284        set key [format "LEQV HOLD%1d%2d" $phase $j]
2285        makeexprec $key
2286        setexp $key $line 3 [expr ($k * 8) + 2]
2287    }   
2288}
2289
2290
2291#  get a logical constraint
2292#
2293#  type action
2294#  -----------
2295#  atom get  number        returns a list of constraints.
2296#   "   set  number value  replaces a list of constraints
2297#                          (value is a list of constraints)
2298#   "   add  number value  inserts a new list of constraints
2299#                          (number is ignored)
2300#   "   delete number      deletes a set of constraint entries
2301# Each item in the list of constraints is composed of 4 items:
2302#              phase, atom, variable, multiplier
2303# If variable=UISO atom can be ALL, otherwise atom is a number
2304# legal variable names: FRAC, X, Y, Z, UISO, U11, U22, U33, U12, U23, U13,
2305#                       MX, MY, MZ
2306#
2307#  type action
2308#  -----------
2309#  profileXX get number         returns a list of constraints for term XX=1-36
2310#                               use number=0 to get # of defined
2311#                                  constraints for term XX
2312#   "        set number value   replaces a list of constraints
2313#                               (value is a list of constraints)
2314#   "        add number value   inserts a new list of constraints
2315#                               (number is ignored)
2316#   "        delete number      deletes a set of constraint entries
2317# Each item in the list of constraints is composed of 3 items:
2318#              phase-list, histogram-list, multiplier
2319# Note that phase-list and/or histogram-list can be ALL
2320
2321proc constrinfo {type action number "value {}"} {
2322    global expmap
2323    if {[lindex $expmap(phasetype) 0] == 4} {
2324        set mm 1
2325    } else {
2326        set mm 0
2327    }
2328    switch -glob ${type}-$action {
2329        atom-get {
2330            # does this constraint exist?
2331            set key [format "LNCN%4d%4d" $number 1]
2332            if {![existsexp $key]} {return -1}
2333            set clist {}
2334            for {set i 1} {$i < 999} {incr i} {
2335                set key [format "LNCN%4d%4d" $number $i]
2336                if {![existsexp $key]} break
2337                set line [readexp $key]
2338                set j1 2
2339                set j2 17
2340                set seg [string range $line $j1 $j2]
2341                while {[string trim $seg] != ""} {
2342                    set p [string range $seg 0 0]
2343                    if {$p == 1 && $mm} {
2344                        set atom [string trim [string range $seg 1 4]]
2345                        set var [string trim [string range $seg 5 7]]
2346                        if {$atom == "ALL"} {
2347                            set var UIS
2348                        } else {
2349                            scan $atom %x atom
2350                        }
2351                        lappend clist [list $p $atom $var \
2352                                [string trim [string range $seg 8 end]]]
2353                    } else {
2354                        lappend clist [list $p \
2355                                [string trim [string range $seg 1 3]] \
2356                                [string trim [string range $seg 4 7]] \
2357                                [string trim [string range $seg 8 end]]]
2358                    }
2359                    incr j1 16
2360                    incr j2 16
2361                    set seg [string range $line $j1 $j2]
2362                }
2363            }
2364            return $clist
2365        }
2366        atom-set {
2367            # delete records for current constraint
2368            for {set i 1} {$i < 999} {incr i} {
2369                set key [format "LNCN%4d%4d" $number $i]
2370                if {![existsexp $key]} break
2371                delexp $key
2372            }
2373            set line {}
2374            set i 1
2375            foreach tuple $value {
2376                set p [lindex $tuple 0]
2377                if {$p == 1 && $mm && \
2378                        [string toupper [lindex $tuple 1]] == "ALL"} {
2379                    set seg [format %1dALL UIS%8.4f \
2380                            [lindex $tuple 0] \
2381                            [lindex $tuple 3]]
2382                } elseif {$p == 1 && $mm} {
2383                    set seg [eval format %1d%.4X%-3s%8.4f $tuple]
2384                } elseif {[string toupper [lindex $tuple 1]] == "ALL"} {
2385                    set seg [format %1dALL%-4s%8.4f \
2386                            [lindex $tuple 0] \
2387                            [lindex $tuple 2] \
2388                            [lindex $tuple 3]]
2389                } else {
2390                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
2391                }
2392                append line $seg
2393                if {[string length $line] > 50} {
2394                    set key  [format "LNCN%4d%4d" $number $i]
2395                    makeexprec $key
2396                    setexp $key $line 3 68
2397                    set line {}
2398                    incr i
2399                }
2400            }
2401            if {$line != ""} {
2402                set key  [format "LNCN%4d%4d" $number $i]
2403                makeexprec $key
2404                setexp $key $line 3 68
2405            }
2406            return
2407        }
2408        atom-add {
2409            # loop over defined constraints
2410            for {set j 1} {$j < 9999} {incr j} {
2411                set key [format "LNCN%4d%4d" $j 1]
2412                if {![existsexp $key]} break
2413            }
2414            set number $j
2415            # save the constraint
2416            set line {}
2417            set i 1
2418            foreach tuple $value {
2419                set p [lindex $tuple 0]
2420                if {$p == 1 && $mm && \
2421                        [string toupper [lindex $tuple 1]] == "ALL"} {
2422                    set seg [format %1dALL UIS%8.4f \
2423                            [lindex $tuple 0] \
2424                            [lindex $tuple 3]]
2425                } elseif {$p == 1 && $mm} {
2426                    set seg [eval format %1d%.4X%-3s%8.4f $tuple]
2427                } elseif {[string toupper [lindex $tuple 1]] == "ALL"} {
2428                    set seg [format %1dALL%-4s%8.4f \
2429                            [lindex $tuple 0] \
2430                            [lindex $tuple 2] \
2431                            [lindex $tuple 3]]
2432                } else {
2433                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
2434                }
2435                append line $seg
2436                if {[string length $line] > 50} {
2437                    set key  [format "LNCN%4d%4d" $number $i]
2438                    makeexprec $key
2439                    setexp $key $line 3 68
2440                    set line {}
2441                    incr i
2442                }
2443            }
2444            if {$line != ""} {
2445                set key  [format "LNCN%4d%4d" $number $i]
2446                makeexprec $key
2447                setexp $key $line 3 68
2448            }
2449            return
2450        }
2451        atom-delete {
2452            for {set j $number} {$j < 9999} {incr j} {
2453                # delete records for current constraint
2454                for {set i 1} {$i < 999} {incr i} {
2455                    set key [format "LNCN%4d%4d" $j $i]
2456                    if {![existsexp $key]} break
2457                    delexp $key
2458                }
2459                # now copy records, from the next entry, if any
2460                set j1 $j
2461                incr j1
2462                set key1 [format "LNCN%4d%4d" $j1 1]
2463                # if there is no record, there is nothing to copy -- done
2464                if {![existsexp $key1]} return
2465                for {set i 1} {$i < 999} {incr i} {
2466                    set key1 [format "LNCN%4d%4d" $j1 $i]
2467                    if {![existsexp $key1]} break
2468                    set key  [format "LNCN%4d%4d" $j  $i]
2469                    makeexprec $key
2470                    setexp $key [readexp $key1] 1 68
2471                }
2472            }
2473        }
2474        profile*-delete {
2475            regsub profile $type {} term
2476            if {$term < 10} {
2477                set term " $term"
2478            }
2479            set key "LEQV PF$term   "
2480            # return nothing if no term exists
2481            if {![existsexp $key]} {return 0}
2482
2483            # number of constraint terms
2484            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2485            # don't delete a non-existing entry
2486            if {$number > $nterms} {return 0}
2487            set val [expr {$nterms - 1}]
2488            validint val 5
2489            setexp $key $val 1 5
2490            for {set i1 $number} {$i1 < $nterms} {incr i1} {
2491                set i2 [expr {1 + $i1}]
2492                # move the contents of constraint #i2 -> i1
2493                if {$i1 > 9} {
2494                    set k1 [expr {($i1+1)/10}]
2495                    set l1 $i1
2496                } else {
2497                    set k1 " "
2498                    set l1 " $i1"
2499                }
2500                set key1 "LEQV PF$term  $k1"
2501                # number of constraint lines for #i1
2502                set n1 [string trim [string range [readexp ${key1}] \
2503                        [expr {($i1%10)*5}] [expr {4+(($i1%10)*5)}]] ]
2504                if {$i2 > 9} {
2505                    set k2 [expr {($i2+1)/10}]
2506                    set l2 $i2
2507                } else {
2508                    set k2 " "
2509                    set l2 " $i2"
2510                }
2511                set key2 "LEQV PF$term  $k2"
2512                # number of constraint lines for #i2
2513                set n2 [string trim [string range [readexp ${key2}] \
2514                        [expr {($i2%10)*5}] [expr {4+(($i2%10)*5)}]] ]
2515                set val $n2
2516                validint val 5
2517                # move the # of terms
2518                setexp $key1 $val [expr {1+(($i1%10)*5)}] 5
2519                # move the terms
2520                for {set j 1} {$j <= $n2} {incr j 1} {
2521                    set key "LEQV PF${term}${l1}$j"
2522                    makeexprec $key
2523                    setexp $key [readexp "LEQV PF${term}${l2}$j"] 1 68
2524                }
2525                # delete any remaining lines
2526                for {set j [expr {$n2+1}]} {$j <= $n1} {incr j 1} {
2527                    delexp "LEQV PF${term}${l1}$j"
2528                }
2529            }
2530
2531            # clear the last term
2532            if {$nterms > 9} {
2533                set i [expr {($nterms+1)/10}]
2534            } else {
2535                set i " "
2536            }
2537            set key "LEQV PF$term  $i"
2538            set cb [expr {($nterms%10)*5}]
2539            set ce [expr {4+(($nterms%10)*5)}]
2540            set n2 [string trim [string range [readexp ${key}] $cb $ce] ]
2541            incr cb
2542            setexp $key "     " $cb 5
2543            # delete any remaining lines
2544            for {set j 1} {$j <= $n2} {incr j 1} {
2545                delexp "LEQV PF${term}${nterms}$j"
2546            }
2547        }
2548        profile*-set {
2549            regsub profile $type {} term
2550            if {$term < 10} {
2551                set term " $term"
2552            }
2553            set key "LEQV PF$term   "
2554            # get number of constraint terms
2555            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2556            # don't change a non-existing entry
2557            if {$number > $nterms} {return 0}
2558            if {$number > 9} {
2559                set k1 [expr {($number+1)/10}]
2560                set l1 $number
2561            } else {
2562                set k1 " "
2563                set l1 " $number"
2564            }
2565            set key1 "LEQV PF$term  $k1"
2566            # old number of constraint lines
2567            set n1 [string trim [string range [readexp ${key1}] \
2568                    [expr {($number%10)*5}] [expr {4+(($number%10)*5)}]] ]
2569            # number of new constraints
2570            set j2 [llength $value]
2571            # number of new constraint lines
2572            set val [set n2 [expr {($j2 + 2)/3}]]
2573            # store the new # of lines
2574            validint val 5
2575            setexp $key1 $val [expr {1+(($number%10)*5)}] 5
2576
2577            # loop over the # of lines in the old or new, whichever is greater
2578            set v0 0
2579            for {set j 1} {$j <= [expr {($n1 > $n2) ? $n1 : $n2}]} {incr j 1} {
2580                set key "LEQV PF${term}${l1}$j"
2581                # were there more lines in the old?
2582                if {$j > $n2} {
2583                    # this line is not needed
2584                    if {$j % 3 == 1} {
2585                        delexp %key
2586                    }
2587                    continue
2588                }
2589                # are we adding new lines?
2590                if {$j > $n1} {
2591                    makeexprec $key
2592                }
2593                # add the three constraints to the line
2594                foreach s {3 23 43} \
2595                        item [lrange $value $v0 [expr {2+$v0}]] {
2596                    if {$item != ""} {
2597                        set val [format %-10s%9.3f \
2598                                [lindex $item 0],[lindex $item 1] \
2599                                [lindex $item 2]]
2600                        setexp $key $val $s 19
2601                    } else {
2602                        setexp $key " " $s 19
2603                    }
2604                }
2605                incr v0 3
2606            }
2607        }
2608        profile*-add {
2609            regsub profile $type {} term
2610            if {$term < 10} {
2611                set term " $term"
2612            }
2613            set key "LEQV PF$term   "
2614            if {![existsexp $key]} {makeexprec $key}
2615            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2616            if {$nterms == ""} {
2617                set nterms 1
2618            } elseif {$nterms >= 99} {
2619                return 0
2620            } else {
2621                incr nterms
2622            }
2623            # store the new # of constraints
2624            set val $nterms
2625            validint val 5
2626            setexp $key $val 1 5
2627
2628            if {$nterms > 9} {
2629                set k1 [expr {($nterms+1)/10}]
2630                set l1 $nterms
2631            } else {
2632                set k1 " "
2633                set l1 " $nterms"
2634            }
2635            set key1 "LEQV PF$term  $k1"
2636
2637            # number of new constraints
2638            set j2 [llength $value]
2639            # number of new constraint lines
2640            set val [set n2 [expr {($j2 + 2)/3}]]
2641            # store the new # of lines
2642            validint val 5
2643            setexp $key1 $val [expr {1+(($nterms%10)*5)}] 5
2644
2645            # loop over the # of lines to be added
2646            set v0 0
2647            for {set j 1} {$j <= $n2} {incr j 1} {
2648                set key "LEQV PF${term}${l1}$j"
2649                makeexprec $key
2650                # add the three constraints to the line
2651                foreach s {3 23 43} \
2652                        item [lrange $value $v0 [expr {2+$v0}]] {
2653                    if {$item != ""} {
2654                        set val [format %-10s%9.3f \
2655                                [lindex $item 0],[lindex $item 1] \
2656                                [lindex $item 2]]
2657                        setexp $key $val $s 19
2658                    } else {
2659                        setexp $key " " $s 19
2660                    }
2661                }
2662                incr v0 3
2663            }
2664        }
2665        profile*-get {
2666            regsub profile $type {} term
2667            if {$term < 10} {
2668                set term " $term"
2669            }
2670            if {$number > 9} {
2671                set i [expr {($number+1)/10}]
2672            } else {
2673                set i " "
2674            }
2675            set key "LEQV PF$term  $i"
2676            # return nothing if no term exists
2677            if {![existsexp $key]} {return 0}
2678            # number of constraint lines
2679           
2680            set numline [string trim [string range [readexp ${key}] \
2681                    [expr {($number%10)*5}] [expr {4+(($number%10)*5)}]] ]
2682            if {$number == 0} {return $numline}
2683            set clist {}
2684            if {$number < 10} {
2685                set number " $number"
2686            }
2687            for {set i 1} {$i <= $numline} {incr i} {
2688                set key "LEQV PF${term}${number}$i"
2689                set line [readexp ${key}]
2690                foreach s {1 21 41} e {20 40 60} {
2691                    set seg [string range $line $s $e]
2692                    if {[string trim $seg] == ""} continue
2693                    # parse the string segment
2694                    set parse [regexp { *([0-9AL]+),([0-9AL]+) +([0-9.]+)} \
2695                            $seg junk phase hist mult]
2696                    # was parse successful
2697                    if {!$parse} {continue}
2698                    lappend clist [list $phase $hist $mult]
2699                }
2700            }
2701            return $clist
2702        }
2703        default {
2704            set msg "Unsupported constrinfo access: type=$type action=$action"
2705            tk_dialog .badexp "Error in readexp access" $msg error 0 OK
2706        }
2707
2708    }
2709}
2710
2711# read the default profile information for a histogram
2712# use: profdefinfo hist set# parm action
2713
2714#     proftype -- profile function number
2715#     profterms -- number of profile terms
2716#     pdamp -- damping value for the profile (*)
2717#     pcut -- cutoff value for the profile (*)
2718#     pterm$n -- profile term #n
2719#     pref$n -- refinement flag value for profile term #n (*)
2720
2721proc profdefinfo {hist set parm "action get"} {
2722    global expgui
2723    if {$hist < 10} {
2724        set key "HST  $hist"
2725    } else {
2726        set key "HST $hist"
2727    }
2728    switch -glob ${parm}-$action {
2729        proftype-get {
2730            set val [string range [readexp "${key}PRCF$set"] 0 4]
2731            if {$val == " "} {return 0}
2732            return $val
2733        }
2734        profterms-get {
2735            set val [string range [readexp "${key}PRCF$set"] 5 9]
2736            if {$val == " "} {return 0}
2737            return $val
2738        }
2739        pcut-get {
2740            return [string trim [string range [readexp "${key}PRCF$set"] 10 19]]
2741        }
2742        pdamp-get {
2743                set val [string range [readexp "${key}PRCF$set"] 24 24]
2744            if {$val == " "} {return 0}
2745            return $val
2746        }
2747        pterm*-get {
2748            regsub pterm $parm {} num
2749            set f1 [expr {15*(($num - 1) % 4)}]
2750            set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
2751            set line  [expr {1 + ($num - 1) / 4}]
2752            return [string trim [string range [\
2753                        readexp "${key}PRCF${set}$line"] $f1 $f2] ]
2754        }
2755        pref*-get {
2756            regsub pref $parm {} num
2757            set f [expr {24+$num}]
2758            if {[string toupper [string range [readexp "${key}PRCF$set"] $f $f]] == "Y"} {
2759                return 1
2760            }
2761            return 0
2762        }
2763        default {
2764            set msg "Unsupported profdefinfo access: parm=$parm action=$action"
2765            tk_dialog .badexp "Code Error" $msg error 0 Exit
2766        }
2767    }
2768}
2769
2770# get March-Dollase preferred orientation information
2771# use MDprefinfo hist phase axis-number parm action value
2772#    ratio    -- ratio of xtallites in PO direction vs random (>1 for more)
2773#    fraction -- fraction in this direction, when more than one axis is used
2774#    h k & l  -- indices of P.O. axis
2775#    ratioref -- flag to vary ratio
2776#    fracref  -- flag to vary fraction
2777#    damp     -- damping value
2778#    type     -- model type (0 = P.O. _|_ to beam, 1 = || to beam)
2779#    new      -- creates a new record with default values (set only)
2780proc MDprefinfo {histlist phaselist axislist parm "action get" "value {}"} {
2781    foreach phase $phaselist hist $histlist axis $axislist {
2782        if {$phase == ""} {set phase [lindex $phaselist end]}
2783        if {$hist == ""} {set hist [lindex $histlist end]}
2784        if {$axis == ""} {set axis [lindex $axislist end]}
2785        if {$hist < 10} {
2786            set hist " $hist"
2787        }
2788        if {$axis > 9} {
2789            set axis "0"
2790        }
2791        set key "HAP${phase}${hist}PREFO${axis}"
2792        switch -glob ${parm}-$action {
2793            ratio-get {
2794                return [string trim [string range [readexp $key] 0 9]]
2795            }
2796            ratio-set {
2797                if ![validreal value 10 6] {return 0}
2798                setexp $key $value 1 10
2799            }
2800            fraction-get {
2801                return [string trim [string range [readexp $key] 10 19]]
2802            }
2803            fraction-set {
2804                if ![validreal value 10 6] {return 0}
2805                setexp $key $value 11 10
2806            }
2807            h-get {
2808                set h [string trim [string range [readexp $key] 20 29]]
2809                # why not allow negative h values?
2810                #               if {$h < 1} {return 0}
2811                return $h
2812            }
2813            h-set {
2814                if ![validreal value 10 2] {return 0}
2815                setexp $key $value 21 10
2816            }
2817            k-get {
2818                set k [string trim [string range [readexp $key] 30 39]]
2819                #               if {$k < 1} {return 0}
2820                return $k
2821            }
2822            k-set {
2823                if ![validreal value 10 2] {return 0}
2824                setexp $key $value 31 10
2825            }
2826            l-get {
2827                set l [string trim [string range [readexp $key] 40 49]]
2828                #if {$l < 1} {return 0}
2829                return $l
2830            }
2831            l-set {
2832                if ![validreal value 10 2] {return 0}
2833                setexp $key $value 41 10
2834            }
2835            ratioref-get {
2836                if {[string toupper \
2837                        [string range [readexp $key] 53 53]] == "Y"} {
2838                    return 1
2839                }
2840                return 0
2841            }
2842            ratioref-set {
2843                if $value {
2844                    setexp $key "Y" 54 1
2845                } else {
2846                    setexp $key "N" 54 1
2847                }
2848            }
2849            fracref-get {
2850                if {[string toupper \
2851                        [string range [readexp $key] 54 54]] == "Y"} {
2852                    return 1
2853                }
2854                return 0
2855            }
2856            fracref-set {
2857                if $value {
2858                    setexp $key "Y" 55 1
2859                } else {
2860                    setexp $key "N" 55 1
2861              }
2862            }
2863            damp-get {
2864                set val [string trim [string range [readexp $key] 59 59]]
2865                if {$val == " "} {return 0}
2866                return $val
2867            }
2868            damp-set {
2869                setexp $key $value 60 1
2870            }
2871            type-get {
2872                set val [string trim [string range [readexp $key] 64 64]]
2873                if {$val == " "} {return 0}
2874                return $val
2875            }
2876            type-set {
2877                # only valid settings are 0 & 1
2878                if {$value != "0" && $value != "1"} {set value "0"}
2879                setexp $key $value 65 1
2880            }
2881            new-set {
2882                makeexprec $key
2883                setexp $key \
2884                        {  1.000000  1.000000  0.000000  0.000000  1.000000   NN    0    0} \
2885                        1 68
2886            }
2887            default {
2888                set msg "Unsupported MDprefinfo access: parm=$parm action=$action"
2889                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2890            }
2891
2892        }
2893
2894    }
2895}
2896
2897# get list of defined atom types
2898proc AtmTypList {} {
2899    set natypes [readexp " EXPR  NATYP"]
2900    if {$natypes == ""} return
2901    set j 0
2902    set typelist {}
2903    for {set i 1} {$i <= $natypes} {incr i} {
2904        set key {this should never be matched}
2905        while {![existsexp $key]} {
2906            incr j
2907            if {$j > 99} {
2908                return $typelist
2909            } elseif {$j <10} {
2910                set key " EXPR ATYP $j"
2911            } else {
2912                set key " EXPR ATYP$j"
2913            }
2914        }
2915        lappend typelist [string trim [string range $::exparray($key) 2 9]]
2916    }
2917    return $typelist
2918}
2919
2920# read information about atom types
2921#     distrad    atomic distance search radius (get/set)
2922#     angrad     atomic angle search radius (get/set)
2923proc AtmTypInfo {parm atmtype "action get" "value {}"} {
2924    # first, search through the records to find the record matching the type
2925    set natypes [readexp " EXPR  NATYP"]
2926    if {$natypes == ""} return
2927    set j 0
2928    set typelist {}
2929    for {set i 1} {$i <= $natypes} {incr i} {
2930        set key {this should never be matched}
2931        while {![existsexp $key]} {
2932            incr j
2933            if {$j > 99} {
2934                return $typelist
2935            } elseif {$j <10} {
2936                set key " EXPR ATYP $j"
2937            } else {
2938                set key " EXPR ATYP$j"
2939            }
2940        }
2941        if {[string toupper $atmtype] == \
2942                [string toupper [string trim [string range $::exparray($key) 2 9]]]} break
2943        set key {}
2944    }
2945    if {$key == ""} {
2946        # atom type not found
2947        return {}
2948    }
2949    switch -glob ${parm}-$action {
2950        distrad-get {
2951            return [string trim [string range [readexp $key] 15 24]]
2952        }
2953        distrad-set {
2954            if ![validreal value 10 2] {return 0}
2955            setexp $key $value 16 10
2956        }
2957        angrad-get {
2958            return [string trim [string range [readexp $key] 25 34]]
2959        }
2960        angrad-set {
2961            if ![validreal value 10 2] {return 0}
2962            setexp $key $value 26 10
2963        }
2964        default {
2965            set msg "Unsupported AtmTypInfo access: parm=$parm action=$action"
2966            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2967        }
2968    }
2969}
2970# read default information about atom types (records copied to the .EXP file
2971# from the gsas/data/atomdata.dat file as AFAC ...
2972#     distrad returns a list of atom types (one or two letters) and
2973#                the corresponding distance
2974# note that these values are read only (no set option)
2975proc DefAtmTypInfo {parm} {
2976    set keys [array names ::exparray " AFAC *_SIZ"]
2977    set elmlist {}
2978    if {[llength $keys] <= 0} {return ""}
2979    foreach key $keys {
2980        lappend elmlist [string trim [string range $key 6 7]]
2981    }
2982    switch -glob ${parm} {
2983        distrad {
2984            set out {}
2985            foreach key $keys elm $elmlist {
2986                set val [string range $::exparray($key) 0 9]
2987                lappend out "$elm [string trim $val]"
2988            }
2989            return $out
2990        }
2991        angrad {
2992            set out {}
2993            foreach key $keys elm $elmlist {
2994                set val [string range $::exparray($key) 10 19]
2995                lappend out "$elm [string trim $val]"
2996            }
2997            return $out
2998        }
2999        default {
3000            set msg "Unsupported DefAtmTypInfo access: parm=$parm"
3001            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3002        }
3003    }
3004}
3005# write the .EXP file
3006proc expwrite {expfile} {
3007    global exparray
3008    set blankline \
3009     "                                                                        "
3010    set fp [open ${expfile} w]
3011    fconfigure $fp -translation crlf -encoding ascii
3012    set keylist [lsort [array names exparray]]
3013    # reorder the keys so that VERSION comes 1st
3014    set pos [lsearch -exact $keylist {     VERSION}]
3015    set keylist "{     VERSION} [lreplace $keylist $pos $pos]"
3016    foreach key $keylist {
3017        puts $fp [string range \
3018                "$key$exparray($key)$blankline" 0 79]
3019    }
3020    close $fp
3021}
3022
3023# history commands -- delete all but last $keep history records,
3024# renumber if $renumber is true
3025proc DeleteHistory {keep renumber} {
3026    global exparray
3027    foreach y [lrange [lsort -decreasing \
3028            [array names exparray {    HSTRY*}]] $keep end] {
3029        unset exparray($y)
3030    }
3031    if !$renumber return
3032    # renumber
3033    set i 0
3034    foreach y [lsort -increasing \
3035            [array names exparray {    HSTRY*}]] {
3036        set key [format "    HSTRY%3d" [incr i]]
3037        set exparray($key) $exparray($y)
3038        unset exparray($y)
3039    }
3040    # list all history
3041    #    foreach y [lsort -decreasing [array names exparray {    HSTRY*}]] {puts "$y $exparray($y)"}
3042}
3043
3044proc CountHistory {} {
3045    global exparray
3046    return [llength [array names exparray {    HSTRY*}]]
3047}
3048
3049# set the phase flags for histogram $hist to $plist
3050proc SetPhaseFlag {hist plist} {
3051    # make a 2 digit key -- hh
3052    if {$hist < 10} {
3053        set hh " $hist"
3054    } else {
3055        set hh $hist
3056    }
3057    set key "HST $hh NPHAS"
3058    set str {}
3059    foreach iph {1 2 3 4 5 6 7 8 9} {
3060        if {[lsearch $plist $iph] != -1} {
3061            append str {    1}
3062        } else {
3063            append str {    0}     
3064        }
3065    }
3066    setexp $key $str 1 68
3067}
3068
3069# erase atom $atom from phase $phase
3070# update the list of atom types, erasing the record if not needed.
3071proc EraseAtom {atom phase} {
3072    set type [atominfo $phase $atom type]
3073    if {$type == ""} return
3074    if {$atom < 10} {
3075        set key "CRS$phase  AT  $atom"
3076    } elseif {$atom < 100} {
3077        set key "CRS$phase  AT $atom"
3078    } else {
3079        set key "CRS$phase  AT$atom"
3080    }
3081    # delete the records for the atom
3082    global exparray
3083    foreach k [array names exparray ${key}*] {
3084        delexp $k
3085    }
3086    # change the number of atoms in the phase
3087    phaseinfo $phase natoms set [expr {[phaseinfo $phase natoms] -1}]
3088
3089    # now adjust numbers in "EXPR ATYP" records and delete, if needed.
3090    set natypes [readexp " EXPR  NATYP"]
3091    if {$natypes == ""} return
3092    set j 0
3093    for {set i 1} {$i <= $natypes} {incr i} {
3094        incr j
3095        if {$j <10} {
3096            set key " EXPR ATYP $j"
3097        } else {
3098            set key " EXPR ATYP$j"
3099        }
3100        while {![existsexp $key]} {
3101            incr j
3102            if {$j > 99} {
3103                return
3104            } elseif {$j <10} {
3105                set key " EXPR ATYP $j"
3106            } else {
3107                set key " EXPR ATYP$j"
3108            }
3109        }
3110        set keytype [string trim [string range $exparray($key) 2 9]]
3111        if {$type == $keytype} {
3112            # found the type record
3113            set val [string trim [string range $exparray($key) 10 14]]
3114            incr val -1
3115            # if this is the last reference, remove the record,
3116            # otherwise, decrement the counter
3117            if {$val <= 0} {
3118                incr natypes -1 
3119                validint natypes 5
3120                setexp " EXPR  NATYP" $natypes 1 5
3121                delexp $key
3122            } else {
3123                validint val 5
3124                setexp $key $val 11 5
3125            }
3126            return
3127        }
3128    }
3129}
3130
3131# compute equivalent anisotropic temperature factor for Uequiv
3132proc CalcAniso {phase Uequiv} {
3133    foreach var {a b c alpha beta gamma} {
3134        set $var [phaseinfo $phase $var]
3135    }
3136
3137    set G(1,1) [expr {$a * $a}]
3138    set G(2,2) [expr {$b * $b}]
3139    set G(3,3) [expr {$c * $c}]
3140    set G(1,2) [expr {$a * $b * cos($gamma*0.017453292519943)}]
3141    set G(2,1) $G(1,2)
3142    set G(1,3) [expr {$a * $c * cos($beta *0.017453292519943)}]
3143    set G(3,1) $G(1,3)
3144    set G(2,3) [expr {$b * $c * cos($alpha*0.017453292519943)}]
3145    set G(3,2) $G(2,3)
3146
3147    # Calculate the volume**2
3148    set v2 0.0
3149    foreach i {1 2 3} {
3150        set J [expr {($i%3) + 1}]
3151        set K [expr {(($i+1)%3) + 1}]
3152        set v2 [expr {$v2+ $G(1,$i)*($G(2,$J)*$G(3,$K)-$G(3,$J)*$G(2,$K))}]
3153    }
3154    if {$v2 > 0} {
3155        set v [expr {sqrt($v2)}]
3156        foreach i {1 2 3} {
3157            set i1 [expr {($i%3) + 1}]
3158            set i2 [expr {(($i+1)%3) + 1}]
3159            foreach j {1 2 3} {
3160                set j1 [expr {($j%3) + 1}]
3161                set j2 [expr {(($j+1)%3) + 1}]
3162                set C($j,$i) [expr {(\
3163                        $G($i1,$j1) * $G($i2,$j2) - \
3164                        $G($i1,$j2)  * $G($i2,$j1)\
3165                        )/ $v}]
3166            }
3167        }
3168        set A(1,2) [expr {0.5 * ($C(1,2)+$C(2,1)) / sqrt( $C(1,1)* $C(2,2) )}]
3169        set A(1,3) [expr {0.5 * ($C(1,3)+$C(3,1)) / sqrt( $C(1,1)* $C(3,3) )}]
3170        set A(2,3) [expr {0.5 * ($C(2,3)+$C(3,2)) / sqrt( $C(2,2)* $C(3,3) )}]
3171        foreach i {1 1 2} j {2 3 3} {
3172            set A($i,$j) [expr {0.5 * ($C($i,$j) + $C($j,$i)) / \
3173                    sqrt( $C($i,$i)* $C($j,$j) )}]
3174            # clean up roundoff
3175            if {abs($A($i,$j)) < 1e-5} {set A($i,$j) 0.0}
3176        }
3177    } else {
3178        set A(1,2) 0.0
3179        set A(1,3) 0.0
3180        set A(2,3) 0.0
3181    }
3182    return "$Uequiv $Uequiv $Uequiv \
3183            [expr {$Uequiv * $A(1,2)}] \
3184            [expr {$Uequiv * $A(1,3)}] \
3185            [expr {$Uequiv * $A(2,3)}]"
3186}
3187
3188# read/edit soft (distance) restraint info
3189#  parm:
3190#    weight -- histogram weight (factr) *
3191#    restraintlist -- list of restraints *
3192#  action: get (default) or set
3193#  value: used only with set
3194#  * =>  read+write supported
3195proc SoftConst {parm "action get" "value {}"} {
3196    set HST {}
3197    # look for RSN record
3198    set n 0
3199    for {set i 0} {$i < $::expmap(nhst)} {incr i} {
3200        set ihist [expr {$i + 1}]
3201        if {[expr {$i % 12}] == 0} {
3202            incr n
3203            set line [readexp " EXPR  HTYP$n"]
3204            if {$line == ""} {
3205                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
3206                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3207            }
3208            set j 0
3209        } else {
3210            incr j
3211        }
3212        if {[string range $line [expr 2+5*$j] [expr 5*($j+1)]] == "RSN "} {
3213            set HST $ihist
3214        }
3215    }
3216    if {$HST <=9} {
3217        set key "HST  $HST"
3218    } else {
3219        set key "HST $HST"
3220    }
3221    if {$HST == "" && $action == "set"} {
3222        # no RSN found need to add the soft constr. histogram
3223        # increment number of histograms
3224        set hst [string trim [string range [readexp { EXPR  NHST }] 0 4]]
3225        incr hst
3226        set HST $hst
3227        if ![validint hst 5] {return 0}
3228        setexp  { EXPR  NHST } $hst 1 5
3229        # add to EXPR HTYPx rec, creating if needed
3230        set n [expr { 1+ (($HST - 1) / 12) }]
3231        set key " EXPR  HTYP$n"
3232        if {[array names ::exparray $key] == ""} {
3233            makeexprec $key
3234        }
3235        setexp $key "RSN " [expr 3 + 5*(($HST-1) % 12)] 5
3236        # create other HST  xx recs
3237        if {$HST <=9} {
3238            set key "HST  $HST"
3239        } else {
3240            set key "HST $HST"
3241        }
3242        makeexprec "$key  HNAM"
3243        setexp "$key  HNAM" "Bond distance restraints" 3 24
3244        makeexprec "$key FACTR"
3245        makeexprec "$key NBNDS"
3246        mapexp
3247    } elseif {$HST == ""} {
3248        if $::expgui(debug) {puts "no restraints"}
3249        return "1"
3250    }
3251
3252    switch -glob ${parm}-$action {
3253        weight-get {
3254            return [string trim [string range [readexp "$key FACTR"] 0 14]]
3255        }
3256        weight-set {
3257            # update FACTR
3258            if ![validreal value 15 6] {return 0}
3259            setexp "$key FACTR" $value 1 15
3260        }
3261        restraintlist-get {
3262            set ncons [string trim [string range [readexp "$key NBNDS"] 0 4]]
3263            set conslist {}
3264            for {set i 1} {$i <= $ncons} {incr i} {
3265                set fi [string toupper [format %.4x $i]]
3266                set line [readexp "${key}BD$fi"]
3267                set const {}
3268                foreach len {3 5 5 3 3 3 3 3 6 6} {
3269                  set lenm1 [expr {$len - 1}]
3270                  lappend const [string trim [string range $line 0 $lenm1]]
3271                  set line [string range $line $len end]
3272                }
3273                lappend conslist $const
3274            }
3275            return $conslist
3276        }
3277        restraintlist-set {
3278            set num [llength $value]
3279            if ![validint num 5] {return 0}
3280            setexp "$key NBNDS" $num 1 5
3281            # delete all old records
3282            foreach i [array names ::exparray "${key}BD*"] {unset ::exparray($i)}
3283            set i 0
3284            foreach cons $value {
3285                incr i
3286                set fi [string toupper [format %.4x $i]]
3287                makeexprec "${key}BD$fi"
3288                set pos 1
3289                foreach num $cons len {3 5 5 3 3 3 3 3 -6 -6} {
3290                    if {$len > 0} {
3291                        validint num $len
3292                        setexp "${key}BD$fi" $num $pos $len
3293                    } else {
3294                        set len [expr abs($len)]
3295                        validreal num $len 3
3296                        setexp "${key}BD$fi" $num $pos $len
3297                    }
3298                    incr pos $len
3299                }
3300            }
3301        }
3302        default {
3303            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
3304            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3305        }
3306    return 1
3307    }
3308}
3309
3310#======================================================================
3311# conversion routines
3312#======================================================================
3313
3314# convert x values to d-space
3315proc tod {xlist hst} {
3316    global expmap
3317    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3318        return [toftod $xlist $hst]
3319    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3320        return [tttod $xlist $hst]
3321    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3322        return [engtod $xlist $hst]
3323    } else {
3324        return {}
3325    }
3326}
3327
3328# convert tof to d-space
3329proc toftod {toflist hst} {
3330    set difc [expr {[histinfo $hst difc]/1000.}]
3331    set difc2 [expr {$difc*$difc}]
3332    set difa [expr {[histinfo $hst difa]/1000.}]
3333    set zero [expr {[histinfo $hst zero]/1000.}]
3334    set ans {}
3335    foreach tof $toflist {
3336        if {$tof == 0.} {
3337            lappend ans 0.
3338        } elseif {$tof == 1000.} {
3339            lappend ans 1000.
3340        } else {
3341            set td [expr {$tof-$zero}]
3342            lappend ans [expr {$td*($difc2+$difa*$td)/ \
3343                    ($difc2*$difc+2.0*$difa*$td)}]
3344        }
3345    }
3346    return $ans
3347}
3348
3349# convert two-theta to d-space
3350proc tttod {twotheta hst} {
3351    set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3352    set zero [expr [histinfo $hst zero]/100.]
3353    set ans {}
3354    set cnv [expr {acos(0.)/180.}]
3355    foreach tt $twotheta {
3356        if {$tt == 0.} {
3357            lappend ans 99999.
3358        } elseif {$tt == 1000.} {
3359            lappend ans 0.
3360        } else {
3361            lappend ans [expr {$lamo2 / sin($cnv*($tt-$zero))}]
3362        }
3363    }
3364    return $ans
3365}
3366
3367# convert energy (edx-ray) to d-space
3368# (note that this ignores the zero correction)
3369proc engtod {eng hst} {
3370    set lam [histinfo $hst lam1]
3371    set zero [histinfo $hst zero]
3372    set ans {}
3373    set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3374    foreach e $eng {
3375        if {$e == 0.} {
3376            lappend ans 1000.
3377        } elseif {$e == 1000.} {
3378            lappend ans 0.
3379        } else {
3380            lappend ans [expr {$v/$e}]
3381        }
3382    }
3383    return $ans
3384}
3385
3386# convert x values to Q
3387proc toQ {xlist hst} {
3388    global expmap
3389    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3390        return [toftoQ $xlist $hst]
3391    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3392        return [tttoQ $xlist $hst]
3393    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3394        return [engtoQ $xlist $hst]
3395    } else {
3396        return {}
3397    }
3398}
3399# convert tof to Q
3400proc toftoQ {toflist hst} {
3401    set difc [expr {[histinfo $hst difc]/1000.}]
3402    set difc2 [expr {$difc*$difc}]
3403    set difa [expr {[histinfo $hst difa]/1000.}]
3404    set zero [expr {[histinfo $hst zero]/1000.}]
3405    set 2pi [expr {4.*acos(0.)}]
3406    set ans {}
3407    foreach tof $toflist {
3408        if {$tof == 0.} {
3409            lappend ans 99999.
3410        } elseif {$tof == 1000.} {
3411            lappend ans 0.
3412        } else {
3413            set td [expr {$tof-$zero}]
3414            lappend ans [expr {$2pi * \
3415                    ($difc2*$difc+2.0*$difa*$td)/($td*($difc2+$difa*$td))}]
3416        }
3417    }
3418    return $ans
3419}
3420
3421# convert two-theta to Q
3422proc tttoQ {twotheta hst} {
3423    set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3424    set zero [expr [histinfo $hst zero]/100.]
3425    set ans {}
3426    set cnv [expr {acos(0.)/180.}]
3427    set 2pi [expr {4.*acos(0.)}]
3428    foreach tt $twotheta {
3429        if {$tt == 0.} {
3430            lappend ans 0.
3431        } elseif {$tt == 1000.} {
3432            lappend ans 1000.
3433        } else {
3434            lappend ans [expr {$2pi * sin($cnv*($tt-$zero)) / $lamo2}]
3435        }
3436    }
3437    return $ans
3438}
3439# convert energy (edx-ray) to Q
3440# (note that this ignores the zero correction)
3441proc engtoQ {eng hst} {
3442    set lam [histinfo $hst lam1]
3443    set zero [histinfo $hst zero]
3444    set ans {}
3445    set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3446    set 2pi [expr {4.*acos(0.)}]
3447    foreach e $eng {
3448        if {$e == 0.} {
3449            lappend ans 0.
3450        } elseif {$e == 1000.} {
3451            lappend ans 1000.
3452        } else {
3453            lappend ans [expr {$2pi * $e / $v}]
3454        }
3455    }
3456    return $ans
3457}
3458proc sind {angle} {
3459    return [expr {sin($angle*acos(0.)/90.)}]
3460}
3461
3462# convert d-space values to 2theta, TOF or KeV
3463proc fromd {dlist hst} {
3464    global expmap
3465    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3466        set difc [expr {[histinfo $hst difc]/1000.}]
3467        set difa [expr {[histinfo $hst difa]/1000.}]
3468        set zero [expr {[histinfo $hst zero]/1000.}]
3469        set ans {}
3470        foreach d $dlist {
3471            if {$d == 0.} {
3472                lappend ans 0.
3473            } elseif {$d == 1000.} {
3474                lappend ans 1000.
3475            } else {
3476                lappend ans [expr {$difc*$d + $difa*$d*$d + $zero}]
3477            }
3478        }
3479        return $ans
3480    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3481        set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3482        set zero [expr [histinfo $hst zero]/100.]
3483        set ans {}
3484        set cnv [expr {180./acos(0.)}]
3485        foreach d $dlist {
3486            if {$d == 99999.} {
3487                lappend ans 0
3488            } elseif {$d == 0.} {
3489                lappend ans 1000.
3490            } else {
3491                lappend ans [expr {$cnv*asin($lamo2/$d) + $zero}]
3492            }
3493        }
3494        return $ans
3495    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3496        set lam [histinfo $hst lam1]
3497        set zero [histinfo $hst zero]
3498        set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3499        set ans {}
3500        foreach d $dlist {
3501            if {$d == 1000.} {
3502                lappend ans 0
3503            } elseif {$d == 0.} {
3504                lappend ans 1000.
3505            } else {
3506                lappend ans [expr {$v/$d}]
3507            }
3508        }
3509        return $ans
3510    } else {
3511        return {}
3512    }
3513}
3514
3515# convert Q values to 2theta, TOF or KeV
3516proc fromQ {Qlist hst} {
3517    global expmap
3518    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3519        set difc [expr {[histinfo $hst difc]/1000.}]
3520        set difa [expr {[histinfo $hst difa]/1000.}]
3521        set zero [expr {[histinfo $hst zero]/1000.}]
3522        set ans {}
3523        foreach Q $Qlist {
3524            if {$Q == 0.} {
3525                lappend ans 1000.
3526            } elseif {$Q == 99999.} {
3527                lappend ans 1000.
3528            } else {
3529                set d [expr {4.*acos(0.)/$Q}]
3530                lappend ans [expr {$difc*$d + $difa*$d*$d + $zero}]
3531            }
3532        }
3533        return $ans
3534    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3535        set lamo4pi [expr {[histinfo $hst lam1]/(8.*acos(0.))}]
3536        set zero [expr [histinfo $hst zero]/100.]
3537        set ans {}
3538        set cnv [expr {180./acos(0.)}]
3539        foreach Q $Qlist {
3540            if {$Q == 0.} {
3541                lappend ans 0
3542            } elseif {$Q == 1000.} {
3543                lappend ans 1000.
3544            } else {
3545                lappend ans [expr {$cnv*asin($Q*$lamo4pi) + $zero}]
3546            }
3547        }
3548        return $ans
3549    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3550        set lam [histinfo $hst lam1]
3551        set zero [histinfo $hst zero]
3552        set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3553        set ans {}
3554        set 2pi [expr {4.*acos(0.)}]
3555        foreach Q $Qlist {
3556            if {$Q == 1000.} {
3557                lappend ans 0
3558            } elseif {$Q == 0.} {
3559                lappend ans 1000.
3560            } else {
3561                lappend ans [expr {$Q * $v/$2pi}]
3562            }
3563        }
3564        return $ans
3565    } else {
3566        return {}
3567    }
3568}
3569#============================================================================
3570# rigid body EXP editing routines (to move into readexp.tcl)
3571# RigidBodyList -- returns a list of the defined rigid body types
3572# SetRigidBodyVar -- set variables and damping for rigid body type multipliers
3573# ReadRigidBody  -- # of times a body is mapped, scaling factors, var #s & coordinates
3574# RigidBodyMappingList - return a list instances where a RB is mapped in phase
3575# RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping
3576# RigidBodySetTLS  -- change the TLS values for a rigid body mapping
3577# RigidBodySetDamp -- change the damping values for a rigid body mapping
3578# RigidBodyVary    -- set refinement variable numbers for a rigid body mapping
3579# RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping
3580# AddRigidBody -- defines a new rigid body type
3581# DeleteRigidBody -- remove a rigid body definition
3582# ReplaceRigidBody -- replaces a previous rigid body type
3583# ReadRigidBodyMapping  -- get parameters for a rigid body mapping
3584# MapRigidBody -- map a rigid body type into a phase
3585# EditRigidBodyMapping -- change the parameters in a rigid body mapping
3586# UnMapRigidBody --remove a rigid body constraint by removing a RB "instance"
3587#----- note that these older routines should not be used ------
3588# RigidBodyCount -- returns the number of defined rigid bodies (body types)
3589#    use RigidBodyList instead
3590# RigidBodyMappingCount -- # of times a rigid body is mapped in phase
3591#    use RigidBodyMappingList instead
3592#============================================================================
3593# returns the number of defined rigid bodies
3594proc RigidBodyCount {} {
3595    set n [string trim [readexp "RGBD  NRBDS"]]
3596    if {$n == ""} {
3597        set n 0
3598    }
3599    return $n
3600}
3601
3602# returns a list of the defined rigid body types
3603proc RigidBodyList {} {
3604    set n [string trim [readexp "RGBD  NRBDS"]]
3605    if {$n == ""} {
3606        set n 0
3607    }
3608    set rblist {}
3609    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
3610        set value $rbnum
3611        validint value 2
3612        set key "RGBD${value}"
3613        if {[existsexp "$key NATR "]} {
3614            lappend rblist $rbnum
3615        }
3616        if {[llength $rblist] == $n} break
3617    }
3618    return $rblist
3619}
3620
3621# ReadRigidBody provides all information associated with a rigid body type
3622#  rbnum is the rigid body type number
3623# it returns two items:
3624#   the number of times the rigid body is mapped
3625#   a list containing an element for each scaling factor in rigid body #rbnum.
3626# in each element there are four items:
3627#    a multiplier value for the rigid body coordinates
3628#    a damping value (0-9) for the refinement of the multiplier
3629#    a variable number if the multiplier will be refined
3630#    a list of cartesian coordinates coordinates
3631# each cartesian coordinate contains 4 items: x,y,z and a label
3632#  note that the label is present only when the RB is created in EXPGUI and is
3633#  not used in GSAS.
3634proc ReadRigidBody {rbnum} {
3635    if {[lsearch [RigidBodyList] $rbnum] == -1} {
3636        return ""
3637    }
3638    set value $rbnum
3639    validint value 2
3640    set key "RGBD${value}"
3641    set n [string trim [string range [readexp "$key NATR"] 0 4]]
3642    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
3643    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
3644    set out {}
3645    for {set i 1} {$i <= $nmult} {incr i} {
3646        set line [readexp "${key}${i}PARM"]
3647        set mult [string trim [string range $line 0 9]]
3648        set var [string trim [string range $line 10 14]]
3649        set damp [string trim [string range $line 15 19]]
3650        set coordlist {}
3651        for {set j 1} {$j <= $n} {incr j} {
3652            set value $j
3653            validint value 3
3654            set line [readexp "${key}${i}SC$value"]
3655            set x [string trim [string range $line 0 9]]
3656            set y [string trim [string range $line 10 19]]
3657            set z [string trim [string range $line 20 29]]
3658            set lbl [string trim [string range $line 30 39]]
3659            lappend coordlist [list $x $y $z $lbl]
3660        }
3661        lappend out [list $mult $damp $var $coordlist]
3662    }
3663    return [list $used $out]
3664}
3665
3666# SetRigidBodyVar
3667#   rbnum is the rigid body type number
3668#   varnumlist is a list of variable numbers
3669#      note that if this list is shorter than the number of actual multipliers
3670#      for the body, the unspecified variable will not be changed
3671#   damplist   is a list of damping values (0-9)
3672#      note that if the damplist is shorter than the number of actual multipliers
3673#      the unspecified values are not changed
3674#  SetRigidBodVar 2 {1 2 3} {}
3675#       will vary the (first 3) translations in body #3 and will not change the
3676#       damping values
3677#  SetRigidBodVar 3 {} {0 0 0}
3678#       will not change variable settings but will change the (first 3) damping values
3679#  SetRigidBodVar 4 {11 11} {8 8}
3680#      changes both variable numbers and damping at the same time
3681# Nothing is returned
3682proc SetRigidBodyVar {rbnum varnumlist damplist} {
3683    if {[lsearch [RigidBodyList] $rbnum] == -1} {
3684        return ""
3685    }
3686    set value $rbnum
3687    validint value 2
3688    set key "RGBD${value}"
3689    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
3690    for {set i 1} {$i <= $nmult} {incr i} {
3691        set j $i
3692        incr j -1
3693        set var [lindex $varnumlist $j]
3694        if {$var != ""} {
3695            validint var 5
3696            setexp "${key}${i}PARM" $var 11 15
3697        }
3698        set damp [lindex $damplist $j]
3699        if {$damp != ""} {
3700            if {$damp > 9} {set damp 9}
3701            if {$damp < 0} {set damp 0}
3702            validint damp 5
3703        }
3704        setexp "${key}${i}PARM" $damp 16 20
3705    }
3706}
3707
3708
3709# return the number of times rigid body $bodytyp is mapped in phase $phase
3710proc RigidBodyMappingCount {phase bodytyp} {
3711    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
3712    if {! [existsexp "$key  NBDS"]} {return 0}
3713    set n [string trim [readexp "$key  NBDS"]]
3714    if {$n == ""} {
3715        set n 0
3716    }
3717    return $n
3718}
3719# return a list of the instances where rigid body $bodytyp is mapped in phase $phase
3720proc RigidBodyMappingList {phase bodytyp} {
3721    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
3722    if {! [existsexp "$key  NBDS"]} {return {}}
3723    set n [string trim [readexp "$key  NBDS"]]
3724    if {$n == ""} {
3725        set n 0
3726    }
3727    set rblist {}
3728    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
3729        set value $rbnum
3730        validint value 2
3731        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
3732        if {[existsexp "$key  NDA"]} {
3733            lappend rblist $rbnum
3734        }
3735        if {[llength $rblist] == $n} break
3736    }
3737    return $rblist
3738}
3739
3740
3741
3742# reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num)
3743# returns a list of items (most lists) as follows:
3744#   1) sequence # of first atom in body
3745#   2) origin of body in fractional coordinates (3 elements)
3746#   3) Euler angles as 6 pairs of numbers (see below)
3747#   4) variable numbers for the 9 position variables (origin followed by rotations)
3748#   5) damping vals for the 9 position variables (origin followed by rotations)
3749#   6) the TLS values, in order below (empty list if TLS is not in use)
3750#   7) the variable numbers for each TLS values, in order below (or empty)
3751#   8) three damping values for the T, L and S terms.
3752# returns an empty list if no such body exists.
3753#
3754# Euler angles are a list of axes and angles to rotate:
3755#   { {axis1 angle1} {axis2 angle2} ...}
3756# where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
3757#
3758# The 20 TLS terms are ordered:
3759#    T11, T22, T33, T12, T13, T23
3760#    L11, L22, L33, L12, L13, L23
3761#    S12, S13, S21, S23, S31, S32, SAA, SBB
3762#
3763proc ReadRigidBodyMapping {phase bodytyp num} {
3764    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3765        return ""
3766    }
3767    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3768    set first [string trim [string range [readexp "$key  NDA"] 0 4]]
3769    set line [readexp "$key BDFL"]
3770    set varlist {}
3771    set damplist {}
3772    foreach i {0 1 2 3 4 5 6 7 8} {
3773        lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]]
3774        lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]]
3775    }
3776    set TLSdamplist {}
3777    foreach i {54 55 56} {
3778        lappend TLSdamplist [string trim [string range $line $i $i ]]
3779    }
3780    set line [readexp "${key} BDLC"]
3781    set x [string trim [string range $line 0 9]]
3782    set y [string trim [string range $line 10 19]]
3783    set z [string trim [string range $line 20 29]]
3784    set origin [list $x $y $z]
3785    set line [readexp "${key} BDOR"]
3786    set rotations {}
3787    foreach i {0 10 20 30 40 50} {
3788        set angle [string trim [string range $line $i [expr {$i+7}]]]
3789        set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]]
3790        lappend rotations [list $angle $axis]
3791    }
3792    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
3793    set tlsvars {}
3794    set tlsvals {}
3795    if {$TLS != 0} {
3796        set line [readexp "${key}TLSF1"]
3797        for {set j 0} {$j < 20} {incr j} {
3798            set var [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]]
3799            if {$var == ""} {set var 0}
3800            lappend tlsvars $var
3801        }
3802        for {set j 0} {$j < 20} {incr j} {
3803            set i 0
3804            if {$j == 0} {
3805                set i 1
3806            } elseif {$j == 8} {
3807                set i 2
3808            } elseif {$j == 16} {
3809                set i 3
3810            }
3811            if {$i != 0} {
3812                set line [readexp "${key}TLSP$i"]
3813                set i 0
3814                set j1 0
3815                set j2 7
3816            } else {
3817                incr j1 8
3818                incr j2 8
3819            }
3820            set val [string trim [string range $line $j1 $j2]]
3821            if {$val == ""} {set val 0}
3822            lappend tlsvals $val
3823        }
3824    }
3825    return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist]
3826}
3827
3828# Control TLS representation for phase, body # and instance number of a Rigid body mapping
3829#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
3830# Enable TLS use if TLS is non-zero (true). Disable if zero
3831proc RigidBodyEnableTLS {phase bodytyp num TLS} {
3832    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3833        return ""
3834    }
3835    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3836    if {$TLS} {
3837        setexp "${key} LSTF" [format "%5d" 1] 1 5
3838        if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
3839        if {![existsexp "${key}TLSP1"]} {
3840            makeexprec "${key}TLSP1"
3841            set str {}
3842            foreach v {.01 .01 .01 0 0 0 0 0} d {4 4 4 4 4 4 2 2} {
3843                validreal v 8 $d
3844                append str $v
3845            }
3846            setexp "${key}TLSP1" $str 1 64
3847        }
3848        if {![existsexp "${key}TLSP2"]} {
3849            makeexprec "${key}TLSP2"
3850            set str {}
3851            set v 0
3852            foreach d {2 2 2 2 4 4 4 4} {
3853                validreal v 8 $d
3854                append str $v
3855            }
3856            setexp "${key}TLSP2" $str 1 64
3857        }
3858        if {![existsexp "${key}TLSP3"]} {
3859            makeexprec "${key}TLSP3"
3860            set str {}
3861            set v 0
3862            foreach d {4 4 4 4} {
3863                validreal v 8 $d
3864                append str $v
3865            }
3866            setexp "${key}TLSP3" $str 1 64
3867        }
3868    } else {
3869        setexp "${key} LSTF" [format "%5d" 0] 1 5
3870    }
3871    return 1
3872}
3873
3874# Control the TLS values for Rigid body mapping for mapping with
3875#    phase ($phase), body type # ($bodytyp) and instance # ($num)
3876# set the 20 TLS values to the values in TLSvals
3877# There must be exactly 20 TLS terms, which are ordered:
3878#    T11, T22, T33, T12, T13, T23
3879#    L11, L22, L33, L12, L13, L23
3880#    S12, S13, S21, S23, S31, S32, SAA, SBB
3881proc RigidBodySetTLS {phase bodytyp num TLSvals} {
3882    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3883        return ""
3884    }
3885    if {[llength $TLSvals] != 20} {return ""}
3886    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3887    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
3888    if {$TLS == 0} {return ""}
3889    if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
3890    foreach n {1 2 3} {
3891        if {![existsexp "${key}TLSP$n"]} {makeexprec "${key}TLSP$n"}
3892    }
3893    set str {}
3894    set n 1
3895    set i 0
3896    foreach v $TLSvals d {4 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4} {
3897        incr i
3898        validreal v 8 $d
3899        append str $v
3900        if {$i == 8} {
3901            set i 0
3902            setexp "${key}TLSP$n" $str 1 64
3903            incr n
3904            set str {}
3905        }
3906    }
3907    setexp "${key}TLSP$n" $str 1 64
3908    return 1
3909}
3910
3911# set damping values for a Rigid body mapping
3912#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
3913# there must be 9 damping values in RBdamp for the 9 position variables (origin followed by rotations)
3914# Use of TLSdamp is optional, but to be used, TLS representation must be enabled and there must be
3915# three damping terms (for all T terms; for all L terms and for all S terms)
3916proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} {
3917    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3918        return ""
3919    }
3920    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3921    if {[llength $RBdamp] != 9} {return ""}
3922    set str {}
3923    foreach v $RBdamp {
3924        if {[validint v 1] != 1} {set v " "}
3925        append str $v
3926    }
3927    setexp "$key BDFL" $str 46 9
3928    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
3929    if {$TLS != 0 &&  [llength $TLSdamp] == 3} {
3930        set str {}
3931        foreach v $TLSdamp {
3932        if {[validint v 1] != 1} {set v " "}
3933            append str $v
3934        }
3935        setexp "$key BDFL" $str 55 3
3936    }
3937    return 1
3938}
3939
3940# set refinement variable numbers for a Rigid body mapping
3941#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
3942# there must be 9 variable values in RBvar for the 9 position variables (origin followed by rotations)
3943# note that the variable values should be unique integers
3944proc RigidBodyVary {phase bodytyp num RBvar} {
3945    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3946        return ""
3947    }
3948    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3949    if {[llength $RBvar] != 9} {return ""}
3950    set str {}
3951    foreach v $RBvar {
3952        if {[validint v 5] != 1} {set v " "}
3953        append str $v
3954    }
3955    setexp "$key BDFL" $str 1 45   
3956}
3957
3958# set TLS refinement variable numbers for a Rigid body mapping
3959#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
3960# there must be 20 variable values in TLSvar for the 20 parameters:
3961#    T11, T22, T33, T12, T13, T23
3962#    L11, L22, L33, L12, L13, L23
3963#    S12, S13, S21, S23, S31, S32, SAA, SBB
3964# note that the variable values should be unique integers
3965proc RigidBodyTLSVary {phase bodytyp num TLSvar} {
3966    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
3967        return ""
3968    }
3969    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
3970    if {[llength $TLSvar] != 20} {return ""}
3971    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
3972    if {$TLS == 0} {return ""}
3973    set str {}
3974    foreach v $TLSvar {
3975        if {[validint v 3] != 1} {set v " "}
3976        append str $v
3977    }
3978    setexp "${key}TLSF1" $str 1 60
3979
3980# AddRigidBody: add a new rigid body definition into the .EXP file
3981# arguments are:
3982#   multlist: defines a list of multipliers for each set of coordinates. In the
3983#             simplest case this will be {1}
3984#   coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } }
3985# note that when the length of multlist > 1 then coordlist must have the same length.
3986# for input where
3987#     multlist = {s1 s2} and
3988#     coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...}
3989#                     {0 0 0} {1 1 0} {2 1 2} ...}
3990#                 }
3991# the cartesian coordinates are defined from the input as
3992#    atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)]
3993#    atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)]
3994#    atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)]
3995#    ...
3996# Returns the number of the rigid body that has been created
3997proc AddRigidBody {multlist coordlist} {
3998    # find the first unused body #
3999    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
4000        set value $rbnum
4001        validint value 2
4002        set key "RGBD${value}"
4003        if {! [existsexp "$key NATR "]} {break}
4004    }
4005    # did we go too far?
4006    if {$rbnum == 16} {return ""}
4007    # increment the RB counter
4008    set n [string trim [readexp "RGBD  NRBDS"]]
4009    if {$n == ""} {
4010        makeexprec "RGBD  NRBDS"
4011        set n 0
4012    }
4013    incr n
4014    validint n 5
4015    setexp "RGBD  NRBDS" $n 1 5
4016    SetRigidBody $rbnum $multlist $coordlist
4017    return $rbnum
4018}
4019
4020# DeleteRigidBody: remove a rigid body definition from the .EXP file
4021# The body may not be mapped. I am not sure if GSAS allows more than 9 bodies,
4022# but if it does, the simplifed approach used here will fail, so this
4023# is not allowed.
4024# Input:
4025#   Rigid body number
4026# Returns:
4027#   1 on success
4028#   -1 if the body number is 11 or greater
4029#   -2 if the body is mapped
4030#   -3 if the body is not defined
4031proc DeleteRigidBody {rbnum} {
4032    # can't delete bodies with numbers higher than 10, since the key prefix
4033    # (RGBD11... will overlap with rigid body instance records, which would be
4034    # deleted below
4035    if {$rbnum > 10} {
4036        return -1
4037    }
4038    set value $rbnum
4039    validint value 2
4040    set key "RGBD${value}"
4041    if {![existsexp "$key NATR "]} {
4042        return -2
4043    }
4044    # make sure the body is not mapped
4045    if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} {
4046        return -3
4047    }
4048    # delete the records starting with "RGBD x" or "RGBD10"
4049    foreach key [array names ::exparray "${key}*"] {
4050        #puts $key
4051        delexp $key
4052    }
4053    # decrement the RB counter
4054    set n [string trim [readexp "RGBD  NRBDS"]]
4055    if {$n == ""} {
4056        set n 0
4057    }
4058    incr n -1
4059    validint n 5
4060    if {$n > 0} {
4061        setexp "RGBD  NRBDS" $n 1 5
4062    } else {
4063        delexp "RGBD  NRBDS"
4064    }
4065    return 1
4066}
4067
4068# ReplaceRigidBody: replace all the information for rigid body #rbnum
4069# Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather
4070# than added.
4071# Note that count of the # of times the body is used is preserved
4072proc ReplaceRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
4073    set value $rbnum
4074    validint value 2
4075    set key "RGBD${value}"
4076    set line [readexp "$key NBDS"]
4077    foreach key [array names ::exparray "${key}*"] {
4078        #puts $key
4079        delexp $key
4080    }
4081    SetRigidBody $rbnum $multlist $coordlist $varlist $damplist
4082    setexp "$key NBDS" $line 1 68
4083}
4084
4085# Edit the parameters for rigid body #rbnum
4086# (normally called from ReplaceRigidBody or AddRigidBody)
4087proc SetRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
4088    set value $rbnum
4089    validint value 2
4090    set key "RGBD${value}"
4091    # number of atoms
4092    set value [llength [lindex $coordlist 0]]
4093    validint value 5
4094    makeexprec "$key NATR"
4095    setexp "$key NATR" $value 1 5
4096    # number of times used
4097    set value 0
4098    validint value 5
4099    makeexprec "$key NBDS"
4100    setexp "$key NBDS" $value 1 5
4101    # number of coordinate matrices
4102    set value [llength $multlist]
4103    validint value 5
4104    makeexprec "$key NSMP"
4105    setexp "$key NSMP" $value 1 5
4106    set i 0
4107    foreach mult $multlist coords $coordlist {
4108        set var [lindex $varlist $i]
4109        if {$var == ""} {set var 0}
4110        set damp [lindex $damplist $i]
4111        if {$damp == ""} {set damp 0}
4112        incr i
4113        makeexprec "${key}${i}PARM"
4114        setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult $var $damp] 1 20
4115        set j 0
4116        foreach item $coords {
4117            #puts $item
4118            incr j
4119            set value $j
4120            validint value 3
4121            makeexprec "${key}${i}SC$value"
4122            if {[llength $item] == 4} {
4123                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
4124            } elseif {[llength $item] == 3} {
4125                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
4126            } else {
4127                return -code 3 "Invalid number of coordinates"
4128            }
4129        }
4130    }
4131}
4132
4133# convert a decimal to the GSAS hex encoding with a field $digits long.
4134proc ToHex {num digits} {
4135    return [string toupper [format "%${digits}x" $num]]
4136}
4137
4138# convert a GSAS hex encoding to a decimal integer
4139proc FromHex {hex} {
4140    return [scan $hex "%x"]
4141}
4142
4143# MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates
4144# (and optionally U values) for a set of atoms will be generated from the rigid body
4145# arguments:
4146#   phase: phase number (1-9)
4147#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4148#   firstatom: sequence number of the first atom in phase (note that atoms may
4149#              not be numbered sequentially)
4150#   position: list of three fractional coordinates for the origin of the rigid body coordinates
4151#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
4152#           cartesian system before the body is translated to position.
4153# returns the instance # (number of times body $bodytyp has been used in phase $phase)
4154proc MapRigidBody {phase bodytyp firstatom position angles} {
4155    # find the first unused body # for this phase & type
4156    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
4157        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
4158        if {! [existsexp "$key  NDA"]} {break}
4159    }
4160    # did we go too far?
4161    if {$rbnum == 16} {return ""}
4162    # increment number of mapped bodies of this type overall
4163    set value $bodytyp
4164    validint value 2
4165    set key "RGBD${value}"
4166    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
4167    incr used
4168    set value $used
4169    validint value 5
4170    setexp "$key NBDS" $value 1 5
4171    # increment number of mapped bodies of this type in this phase
4172    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4173    if {[existsexp "$key  NBDS"]} {
4174        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
4175    } else {
4176        makeexprec "$key  NBDS"
4177        set used 0
4178    }
4179    incr used
4180    set value $used
4181    validint value 5
4182    setexp "$key  NBDS" $value 1 5
4183    # now write the mapping parameters
4184    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
4185    set value $firstatom
4186    validint value 5
4187    makeexprec "$key  NDA"
4188    setexp "$key  NDA" $value 1 5
4189    set l1 {}
4190    set l2 {}
4191    for {set i 0} {$i < 9} {incr i} {
4192        append l1 [format %5d 0]
4193        append l2 [format %1d 0]
4194    }
4195    makeexprec "$key BDFL"
4196    setexp "$key BDFL" $l1$l2 1 54
4197    makeexprec "${key} BDLC"
4198    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
4199    makeexprec "${key} BDOR"
4200    set l1 {}
4201    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
4202        append l1 [format "%8.2f%2d" $val $dir]
4203    }
4204    setexp "${key} BDOR" $l1 1 60
4205    makeexprec "${key} LSTF"
4206    setexp "${key} LSTF" [format "%5d" 0] 1 5
4207    # turn off the X refinement flags for the new body
4208    set st [lsearch $::expmap(atomlist_$phase) $firstatom]
4209    set natoms [llength [lindex [lindex [lindex [ReadRigidBody $bodytyp] 1] 0] 3]]
4210    set en [expr {$st+$natoms-1}]
4211    set atomlist [lrange $::expmap(atomlist_$phase) $st $en]
4212    atominfo $phase $atomlist xref set 0
4213    # redo the mapping to capture the newly mapped atoms
4214    mapexp
4215    return $rbnum
4216}
4217
4218# EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody)
4219# arguments:
4220#   phase: phase number (1-9)
4221#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4222#   bodynum: instance number, as returned by MapRigidBody
4223#   position: list of three fractional coordinates for the origin of the rigid body coordinates
4224#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
4225#           cartesian system before the body is translated to position.
4226#
4227proc EditRigidBodyMapping {phase bodytyp bodynum position angles} {
4228    # number of bodies of this type in this phase
4229    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
4230    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
4231    set l1 {}
4232    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
4233        append l1 [format "%8.2f%2d" $val $dir]
4234    }
4235    setexp "${key} BDOR" $l1 1 60
4236}
4237
4238# UnMapRigidBody: remove a rigid body constraint by removing a RB "instance"
4239# (undoes MapRigidBody)
4240# arguments:
4241#   phase: phase number (1-9)
4242#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4243#   bodynum: instance number, as returned by MapRigidBody
4244proc UnMapRigidBody {phase bodytyp mapnum} {
4245    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} {
4246        return ""
4247    }
4248    # decrement number of mapped bodies of this type overall
4249    set value $bodytyp
4250    validint value 2
4251    set key "RGBD${value}"
4252    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
4253    incr used -1
4254    set value $used
4255    validint value 5
4256    setexp "$key NBDS" $value 1 5
4257    # decrement number of mapped bodies of this type in this phase
4258    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4259    if {[existsexp "$key  NBDS"]} {
4260        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
4261    } else {
4262        set used 0
4263    }
4264    incr used -1
4265    set value $used
4266    validint value 5
4267    if {$used > 0} {
4268        setexp "$key  NBDS" $value 1 5
4269    } else {
4270        delexp "$key  NBDS"
4271    }
4272    # now delete the mapping parameter records
4273    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]"
4274    foreach key [array names ::exparray "${key}*"] {
4275        delexp $key
4276    }
4277    return $used
4278}
4279
Note: See TracBrowser for help on using the repository browser.