source: trunk/readexp.tcl @ 1227

Last change on this file since 1227 was 1227, checked in by toby, 8 years ago

fixes for anom and fourier GUIs

  • Property svn:keywords set to Author Date Revision Id
File size: 144.5 KB
Line 
1# $Id: readexp.tcl 1227 2012-12-07 22:44:51Z toby $
2# Routines to deal with the .EXP "data structure"
3set expmap(Revision) {$Revision: 1227 $ $Date: 2012-12-07 22:44:51 +0000 (Fri, 07 Dec 2012) $}
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                if {![existsexp ${key}HSCALE]} {
1464                    # fix missing scale factor record
1465                    makeexprec ${key}HSCALE
1466                    set value 1.0
1467                    validreal value 15 6
1468                    setexp ${key}HSCALE $value 1 15
1469                    catch {incr ::expgui(changed)}
1470                    setexp ${key}HSCALE "N" 20 1
1471                }
1472                return [string trim [string range [readexp ${key}HSCALE] 0 14]]
1473            }
1474            scale-set {
1475                if ![validreal value 15 6] {return 0}
1476                setexp ${key}HSCALE $value 1 15
1477            }
1478            sref-get {
1479                if {[string toupper [string range [readexp ${key}HSCALE] 19 19]] == "Y"} {
1480                    return 1
1481                }
1482                return 0
1483            }
1484            sref-set {
1485                if $value {
1486                    setexp ${key}HSCALE "Y" 20 1
1487                } else {
1488                    setexp ${key}HSCALE "N" 20 1
1489                }           
1490            }
1491            sdamp-get {
1492                set val [string range [readexp ${key}HSCALE] 24 24]
1493                if {$val == " "} {return 0}
1494                return $val
1495            }
1496            sdamp-set {
1497                setexp ${key}HSCALE $value 25 1
1498            }
1499
1500            difc-get -
1501            lam1-get {
1502                return [string trim [string range [readexp "${key} ICONS"] 0 9]]
1503            }
1504            difc-set -
1505            lam1-set {
1506                if ![validreal value 10 7] {return 0}
1507                setexp "${key} ICONS" $value 1 10
1508                # set the powpref warning (1 = suggested)
1509                catch {
1510                    global expgui
1511                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1512                    set msg "Diffractometer constants" 
1513                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1514                        append expgui(needpowpref_why) "\t$msg were changed\n"
1515                    }
1516                }
1517            }
1518            difa-get -
1519            lam2-get {
1520                return [string trim [string range [readexp "${key} ICONS"] 10 19]]
1521            }
1522            difa-set -
1523            lam2-set {
1524                if ![validreal value 10 7] {return 0}
1525                setexp "${key} ICONS" $value 11 10
1526                # set the powpref warning (1 = suggested)
1527                catch {
1528                    global expgui
1529                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1530                    set msg "Diffractometer constants" 
1531                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1532                        append expgui(needpowpref_why) "\t$msg were changed\n"
1533                    }
1534                }
1535            }
1536            zero-get {
1537                return [string trim [string range [readexp "${key} ICONS"] 20 29]]
1538            }
1539            zero-set {
1540                if ![validreal value 10 5] {return 0}
1541                setexp "${key} ICONS" $value 21 10
1542                # set the powpref warning (1 = suggested)
1543                catch {
1544                    global expgui
1545                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1546                    set msg "Diffractometer constants" 
1547                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1548                        append expgui(needpowpref_why) "\t$msg were changed\n"
1549                    }
1550                }
1551            }
1552            ipola-get {
1553                return [string trim [string range [readexp "${key} ICONS"] 54 54]]
1554            }
1555            ipola-set {
1556                if ![validint value 1] {return 0}
1557                setexp "${key} ICONS" $value 55 1
1558            }
1559            pola-get {
1560                return [string trim [string range [readexp "${key} ICONS"] 40 49]]
1561            }
1562            pola-set {
1563                if ![validreal value 10 5] {return 0}
1564                setexp "${key} ICONS" $value 41 10
1565            }
1566            kratio-get {
1567                set val [string trim [string range [readexp "${key} ICONS"] 55 64]]
1568                if {$val == ""} {set val 0}
1569                # N.B. this code is used w/CW, where Kratio may not be 0.0
1570                set lam2 [string trim [string range [readexp "${key} ICONS"] 10 19]]
1571                if {$lam2 == ""} {set lam2 0}
1572                # Change kratio & flag the change (this is rather kludged)
1573                if {$val == 0 && $lam2 != 0} {
1574                    set val 0.5
1575                    validreal val 10 5
1576                    setexp "${key} ICONS" $val 56 10
1577                    catch {incr ::expgui(changed)}
1578                }
1579                return $val
1580            }
1581            kratio-set {
1582                if ![validreal value 10 5] {return 0}
1583                setexp "${key} ICONS" $value 56 10
1584            }
1585
1586            wref-get {
1587            #------------------------------------------------------
1588            # col 33: refine flag for lambda, difc, ratio and theta
1589            #------------------------------------------------------
1590                if {[string toupper [string range \
1591                        [readexp "${key} ICONS"] 32 32]] == "L"} {
1592                    return 1
1593                }
1594                return 0
1595            }
1596            wref-set {
1597                if $value {
1598                    setexp "${key} ICONS" "L" 33 1
1599                } else {
1600                    setexp "${key} ICONS" " " 33 1
1601                }           
1602            }
1603            ratref-get {
1604                if {[string toupper [string range \
1605                        [readexp "${key} ICONS"] 32 32]] == "R"} {
1606                    return 1
1607                }
1608                return 0
1609            }
1610            ratref-set {
1611                if $value {
1612                    setexp "${key} ICONS" "R" 33 1
1613                } else {
1614                    setexp "${key} ICONS" " " 33 1
1615                }           
1616            }
1617            dcref-get {
1618                if {[string toupper [string range \
1619                        [readexp "${key} ICONS"] 32 32]] == "C"} {
1620                    return 1
1621                }
1622                return 0
1623            }
1624            dcref-set {
1625                if $value {
1626                    setexp "${key} ICONS" "C" 33 1
1627                } else {
1628                    setexp "${key} ICONS" " " 33 1
1629                }           
1630            }
1631            ttref-get {
1632                if {[string toupper [string range \
1633                        [readexp "${key} ICONS"] 32 32]] == "T"} {
1634                    return 1
1635                }
1636                return 0
1637            }
1638            ttref-set {
1639                if $value {
1640                    setexp "${key} ICONS" "T" 33 1
1641                } else {
1642                    setexp "${key} ICONS" " " 33 1
1643                }           
1644            }
1645
1646
1647            pref-get {
1648            #------------------------------------------------------
1649            # col 34: refine flag for POLA & DIFA
1650            #------------------------------------------------------
1651                if {[string toupper [string range \
1652                        [readexp "${key} ICONS"] 33 33]] == "P"} {
1653                    return 1
1654                }
1655                return 0
1656            }
1657            pref-set {
1658                if $value {
1659                    setexp "${key} ICONS" "P" 34 1
1660                } else {
1661                    setexp "${key} ICONS" " " 34 1
1662                }           
1663            }
1664            daref-get {
1665                if {[string toupper [string range \
1666                        [readexp "${key} ICONS"] 33 33]] == "A"} {
1667                    return 1
1668                }
1669                return 0
1670            }
1671            daref-set {
1672                if $value {
1673                    setexp "${key} ICONS" "A" 34 1
1674                } else {
1675                    setexp "${key} ICONS" " " 34 1
1676                }           
1677            }
1678
1679            zref-get {
1680            #------------------------------------------------------
1681            # col 34: refine flag for zero correction
1682            #------------------------------------------------------
1683                if {[string toupper [string range [readexp "${key} ICONS"] 34 34]] == "Z"} {
1684                    return 1
1685                }
1686                return 0
1687            }
1688            zref-set {
1689                if $value {
1690                    setexp "${key} ICONS" "Z" 35 1
1691                } else {
1692                    setexp "${key} ICONS" " " 35 1
1693                }           
1694            }
1695
1696            ddamp-get {
1697                set val [string range [readexp "${key} ICONS"] 39 39]
1698                if {$val == " "} {return 0}
1699                return $val
1700            }
1701            ddamp-set {
1702                setexp "${key} ICONS" $value 40 1
1703            }
1704
1705            backtype-get {
1706                set val [string trim [string range [readexp "${key}BAKGD "] 0 4]]
1707                if {$val == " "} {return 0}
1708                return $val
1709            }
1710            backtype-set {
1711                if ![validint value 5] {return 0}
1712                setexp "${key}BAKGD " $value 1 5
1713            }
1714            backterms-get {
1715                set val [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1716                if {$val == " "} {return 0}
1717                return $val
1718            }
1719            backterms-set {
1720                # this takes a bit of work -- if terms are added, add lines as needed to the .EXP
1721                set oldval [string trim [string range [readexp "${key}BAKGD "] 5 9]]
1722                if ![validint value 5] {return 0}
1723                if {$oldval < $value} {
1724                    set line1  [expr {2 + ($oldval - 1) / 4}]
1725                    set line2  [expr {1 + ($value - 1) / 4}]
1726                    for {set i $line1} {$i <= $line2} {incr i} {
1727                        # create a blank entry if needed
1728                        makeexprec ${key}BAKGD$i
1729                    }
1730                    incr oldval
1731                    for {set num $oldval} {$num <= $value} {incr num} {
1732                        set f1 [expr {15*(($num - 1) % 4)}]
1733                        set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
1734                        set line  [expr {1 + ($num - 1) / 4}]
1735                        if {[string trim [string range [readexp ${key}BAKGD$line] $f1 $f2]] == ""} {
1736                            set f1 [expr {15*(($num - 1) % 4)+1}]
1737                            setexp ${key}BAKGD$line 0.0 $f1 15                 
1738                        }
1739                    }
1740                }
1741                setexp "${key}BAKGD " $value 6 5
1742
1743            }
1744            bref-get {
1745                if {[string toupper [string range [readexp "${key}BAKGD"] 14 14]] == "Y"} {
1746                    return 1
1747                }
1748                return 0
1749            }
1750            bref-set {
1751                if $value {
1752                    setexp "${key}BAKGD "  "Y" 15 1
1753                } else {
1754                    setexp "${key}BAKGD "  "N" 15 1
1755                }
1756            }
1757            bdamp-get {
1758                set val [string range [readexp "${key}BAKGD "] 19 19]
1759                if {$val == " "} {return 0}
1760                return $val
1761            }
1762            bdamp-set {
1763                setexp "${key}BAKGD " $value 20 1
1764            }
1765            bterm*-get {
1766                regsub bterm $parm {} num
1767                set f1 [expr {15*(($num - 1) % 4)}]
1768                set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
1769                set line  [expr {1 + ($num - 1) / 4}]
1770                return [string trim [string range [readexp ${key}BAKGD$line] $f1 $f2] ]
1771            }
1772            bterm*-set {
1773                regsub bterm $parm {} num
1774                if ![validreal value 15 6] {return 0}
1775                set f1 [expr {15*(($num - 1) % 4)+1}]
1776                set line  [expr {1 + ($num - 1) / 4}]
1777                setexp ${key}BAKGD$line $value $f1 15
1778            }
1779            bank-get {
1780                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1781            }
1782            tofangle-get {
1783                return [string trim [string range [readexp "${key}BNKPAR"] 10 19]]
1784            }
1785            tmin-get {
1786                return [string trim [string range [readexp "${key} TRNGE"] 0 9]]
1787            }
1788            tmax-get {
1789                return [string trim [string range [readexp "${key} TRNGE"] 10 19]]
1790            }
1791            excl-get {
1792                set n [string trim [string range [readexp "${key} NEXC"] 0 4]]
1793                set exlist {}
1794                for {set i 1} {$i <= $n} {incr i} {
1795                    set line [readexp [format "${key}EXC%3d" $i]]
1796                    lappend exlist [list \
1797                            [string trim [string range $line  0  9]] \
1798                            [string trim [string range $line 10 19]]]
1799                }
1800                return $exlist
1801            }
1802            excl-set {
1803                set n [llength $value]
1804                if ![validint n 5] {return 0}
1805                setexp "${key} NEXC" $n 1 5
1806                set i 0
1807                foreach p $value {
1808                    incr i
1809                    foreach {r1 r2} $p {}
1810                    validreal r1 10 3
1811                    validreal r2 10 3
1812                    set k [format "${key}EXC%3d" $i]
1813                    if {![existsexp $k]} {
1814                        makeexprec $k
1815                    }
1816                    setexp $k ${r1}${r2} 1 20
1817                }
1818                # set the powpref warning (2 = required)
1819                catch {
1820                    global expgui
1821                    set expgui(needpowpref) 2
1822                    set msg "Excluded regions" 
1823                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1824                        append expgui(needpowpref_why) "\t$msg were changed\n"
1825                    }
1826                }
1827            }
1828            file-get {
1829                return [string trim [readexp "${key}  HFIL"] ]
1830            }
1831            file-set {
1832                setexp "${key}  HFIL" $value 3 65
1833            }
1834            bank-get {
1835                return [string trim [string range [readexp "${key} BANK"] 0 4]]
1836            }
1837            dmin-get {
1838                return [string trim [string range [readexp "${key} NREF"] 5 14]]
1839            }
1840            dmin-set {
1841                if ![validreal value 10 4] {return 0}
1842                setexp "${key} NREF" $value 6 10
1843                # set the powpref warning (2 = required)
1844                catch {
1845                    global expgui
1846                    set expgui(needpowpref) 2
1847                    set msg "Dmin (reflection range)" 
1848                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1849                        append expgui(needpowpref_why) "\t$msg was changed\n"
1850                    }
1851                }
1852            }
1853            use-get {
1854                set k [expr {($hist+11)/12}]
1855                set line [readexp " EXPR  HTYP$k"]
1856                set j [expr {((($hist-1) % 12)+1)*5}]
1857                if {[string range $line $j $j] == "*"} {return 0}
1858                return 1
1859            }
1860            use-set {
1861                set k [expr {($hist+11)/12}]
1862                set line [readexp " EXPR  HTYP$k"]
1863                set j [expr {((($hist-1) % 12)+1)*5+1}]
1864                if {$value} {
1865                    setexp " EXPR  HTYP$k" " " $j 1
1866                } else {
1867                    setexp " EXPR  HTYP$k" "*" $j 1
1868                }
1869                # set the powpref warning (2 = required)
1870                catch {
1871                    global expgui
1872                    set expgui(needpowpref) 2
1873                    set msg "Histogram use flags" 
1874                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1875                        append expgui(needpowpref_why) "\t$msg were changed\n"
1876                    }
1877                }
1878            }
1879            dstart-get {
1880                return [string trim [string range [readexp "${key} DUMMY"] 20 29]]
1881            }
1882            dstart-set {
1883                if ![validreal value 10 3] {return 0}
1884                setexp "${key} DUMMY" $value 21 10
1885                # set the powpref warning (1 = suggested)
1886                catch {
1887                    global expgui
1888                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1889                    set msg "Dummy histogram parameters" 
1890                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1891                        append expgui(needpowpref_why) "\t$msg were changed\n"
1892                    }
1893                }
1894            }
1895            dstep-get {
1896                return [string trim [string range [readexp "${key} DUMMY"] 30 39]]
1897            }
1898            dstep-set {
1899                if ![validreal value 10 3] {return 0}
1900                setexp "${key} DUMMY" $value 31 10
1901                catch {
1902                    global expgui
1903                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1904                    set msg "Dummy histogram parameters" 
1905                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1906                        append expgui(needpowpref_why) "\t$msg were changed\n"
1907                    }
1908                }
1909            }
1910            dpoints-get {
1911                return [string trim [string range [readexp "${key} DUMMY"] 0 9]]
1912            }
1913            dpoints-set {
1914                if ![validint value 10] {return 0}
1915                setexp "${key} DUMMY" $value 1 10
1916                catch {
1917                    global expgui
1918                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
1919                    set msg "Dummy histogram parameters" 
1920                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
1921                        append expgui(needpowpref_why) "\t$msg were changed\n"
1922                    }
1923                }
1924            }
1925            dtype-get {
1926                return [string trim [string range [readexp "${key} DUMMY"] 10 19]]
1927            }
1928            abscor1-get {
1929                return [string trim [string range [readexp "${key}ABSCOR"] 0 14]]
1930            }
1931            abscor1-set {
1932                if ![validreal value 15 7] {return 0}
1933                setexp "${key}ABSCOR" $value 1 15
1934            }
1935            abscor2-get {
1936                return [string trim [string range [readexp "${key}ABSCOR"] 15 29]]
1937            }
1938            abscor2-set {
1939                # this must have a decimal as the 5th character, so that we end up with a
1940                # decimal point in column 20.
1941                set tmp $value
1942                if ![validreal tmp 12 7] {return 0}
1943                set pos [string first "." $tmp]
1944                while {$pos < 4} {
1945                    set tmp " $tmp"
1946                    set pos [string first "." $tmp]
1947                }
1948                if {$pos == 4} {
1949                    setexp "${key}ABSCOR" $tmp 16 15
1950                    return
1951                }
1952                catch {
1953                    set tmp [format "%12.6E" $value]
1954                    set pos [string first "." $tmp]
1955                    while {$pos < 4} {
1956                        set tmp " $tmp"
1957                        set pos [string first "." $tmp]
1958                    }
1959                    if {$pos == 4} {
1960                        setexp "${key}ABSCOR" $tmp 16 15
1961                        return
1962                    }
1963                }
1964                return 0
1965            }
1966            abstype-get {
1967                set val [string trim [string range [readexp "${key}ABSCOR"] 40 44]]
1968                if {$val == ""} {set val 0}
1969                return $val
1970            }
1971            abstype-set {
1972                if ![validint value 5] {return 0}
1973                setexp "${key}ABSCOR" $value 41 5
1974            }
1975            absdamp-get {
1976                set val [string range [readexp "${key}ABSCOR"] 39 39]
1977                if {$val == " "} {return 0}
1978                return $val
1979            }
1980            absdamp-set {
1981                if ![validint value 5] {return 0}
1982                setexp "${key}ABSCOR" $value 36 5
1983            }
1984            absref-get {
1985                if {[string toupper \
1986                        [string range [readexp "${key}ABSCOR"] 34 34]] == "Y"} {
1987                    return 1
1988                }
1989                return 0
1990            }
1991            absref-set {
1992                if $value {
1993                    setexp "${key}ABSCOR" "    Y" 31 5
1994                } else {
1995                    setexp "${key}ABSCOR" "    N" 31 5
1996                }
1997            }
1998            ITYP-get {
1999                return [string trim [readexp "${key}I ITYP"]]
2000            }
2001            proftbl-get {
2002                set line [readexp "${key}PAB3"]
2003                if {$line == ""} {return 0}
2004                set val [string trim [string range $line 0 4]]
2005                if {$val == ""} {return 0}
2006                return $val
2007            }
2008            anomff-get {
2009                set l {}
2010                foreach i {1 2 3 4 5 6 7 8 9} {
2011                    if {![existsexp "${key}FFANS$i"]} continue
2012                    set line [readexp "${key}FFANS$i"]
2013                    set elem [string trim [string range $line 2 9]]
2014                    set fp [string trim [string range $line 10 19]]
2015                    set fpp [string trim [string range $line 20 29]]
2016                    lappend l [list $elem $fp $fpp]
2017                }
2018                return $l
2019            }
2020            anomff-set {
2021                # match up input against elements in list.
2022                # change elements included, return any elements that are
2023                # not found.
2024                set errorlist {}
2025                foreach triplet $value {
2026                    foreach {e fp fpp} $triplet {}               
2027                    foreach i {1 2 3 4 5 6 7 8 9} {
2028                        if {![existsexp "${key}FFANS$i"]} continue
2029                        # note that the element name is not used or validated
2030                        set elem [string trim [string range \
2031                                                   [readexp "${key}FFANS$i"] 2 9]]
2032                        if {[string match -nocase $e $elem]} { 
2033                            if ![validreal fp 10 3] {return 0}
2034                            setexp "${key}FFANS$i" $fp 11 10
2035                            if ![validreal fpp 10 3] {return 0}
2036                            setexp "${key}FFANS$i" $fpp 21 10
2037                            set e {}
2038                            break
2039                        }
2040                    }
2041                    if {$e != ""} {
2042                        # oops, no match
2043                        lappend errorlist $e
2044                    }
2045                }
2046                if {$errorlist != ""} {return [list 0 $errorlist]}
2047            }
2048            default {
2049                set msg "Unsupported histinfo access: parm=$parm action=$action"
2050                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2051            }
2052        }
2053    }
2054    return 1
2055}
2056
2057proc add_anomff {histlist type {fp 0} {fpp 0}} {
2058    global expgui
2059    foreach hist $histlist {
2060        if {$hist < 10} {
2061            set key "HST  $hist"
2062        } else {
2063            set key "HST $hist"
2064        }
2065        if ![validreal fp 10 3] {return 0}
2066        if ![validreal fpp 10 3] {return 0}
2067        foreach i {1 2 3 4 5 6 7 8 9} {
2068            if {![existsexp "${key}FFANS$i"]} {
2069                makeexprec "${key}FFANS$i"
2070                setexp "${key}FFANS$i" [string trim $type] 3 8
2071                setexp "${key}FFANS$i" $fp 11 10
2072                setexp "${key}FFANS$i" $fpp 21 10
2073                setexp "${key}FFANS$i" "NN" 34 2
2074                return 1
2075            }
2076        }
2077        return 0
2078    }
2079}
2080
2081# read the information that differs by both histogram and phase (profile & phase fraction)
2082# use: hapinfo hist phase parm action value
2083
2084#     frac -- phase fraction (*)
2085#     frref/frdamp -- refinement flag/damping value for the phase fraction (*)
2086#     proftype -- profile function number (*)
2087#     profterms -- number of profile terms (*)
2088#     pdamp -- damping value for the profile (*)
2089#     pcut -- cutoff value for the profile (*)
2090#     pterm$n -- profile term #n (*)
2091#     pref$n -- refinement flag value for profile term #n (*)
2092#     extmeth -- Fobs extraction method (*)
2093#     POnaxis -- number of defined M-D preferred axes
2094proc hapinfo {histlist phaselist parm "action get" "value {}"} {
2095    foreach phase $phaselist hist $histlist {
2096        if {$phase == ""} {set phase [lindex $phaselist end]}
2097        if {$hist == ""} {set hist [lindex $histlist end]}
2098        if {$hist < 10} {
2099            set hist " $hist"
2100        }
2101        set key "HAP${phase}${hist}"
2102        switch -glob ${parm}-$action {
2103            extmeth-get {
2104                set i1 [expr {($phase - 1)*5}]
2105                set i2 [expr {$i1 + 4}]
2106                return [string trim [string range [readexp "HST $hist EPHAS"] $i1 $i2]]
2107            }
2108            extmeth-set {
2109                set i1 [expr {($phase - 1)*5 + 1}]
2110                if ![validint value 5] {return 0}
2111                setexp "HST $hist EPHAS" $value $i1 5
2112            }
2113            frac-get {
2114                return [string trim [string range [readexp ${key}PHSFR] 0 14]]
2115            }
2116            frac-set {
2117                if ![validreal value 15 6] {return 0}
2118                setexp ${key}PHSFR $value 1 15
2119            }
2120            frref-get {
2121                if {[string toupper [string range [readexp ${key}PHSFR] 19 19]] == "Y"} {
2122                    return 1
2123                }
2124                return 0
2125            }
2126            frref-set {
2127                if $value {
2128                    setexp ${key}PHSFR "Y" 20 1
2129                } else {
2130                    setexp ${key}PHSFR "N" 20 1
2131                }           
2132            }
2133            frdamp-get {
2134                set val [string range [readexp ${key}PHSFR] 24 24]
2135                if {$val == " "} {return 0}
2136                return $val
2137            }
2138            frdamp-set {
2139                setexp ${key}PHSFR $value 25 1
2140            }
2141            proftype-get {
2142                set val [string range [readexp "${key}PRCF "] 0 4]
2143                if {$val == " "} {return 0}
2144                return $val
2145            }
2146            proftype-set {
2147                if ![validint value 5] {return 0}
2148                setexp "${key}PRCF " $value 1 5
2149                # set the powpref warning (1 = suggested)
2150                catch {
2151                    global expgui
2152                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2153                    set msg "Profile parameters" 
2154                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2155                        append expgui(needpowpref_why) "\t$msg were changed\n"
2156                    }
2157                }
2158            }
2159            profterms-get {
2160                set val [string range [readexp "${key}PRCF "] 5 9]
2161                if {$val == " "} {return 0}
2162                return $val
2163            }
2164            profterms-set {
2165                if ![validint value 5] {return 0}
2166                setexp "${key}PRCF " $value 6 5
2167                # now check that all needed entries exist
2168                set lines [expr {1 + ($value - 1) / 4}]
2169                for {set i 1} {$i <= $lines} {incr i} {
2170                    makeexprec "${key}PRCF $i"
2171                }
2172                # set the powpref warning (1 = suggested)
2173                catch {
2174                    global expgui
2175                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2176                    set msg "Profile parameters" 
2177                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2178                        append expgui(needpowpref_why) "\t$msg were changed\n"
2179                    }
2180                }
2181            }
2182            pcut-get {
2183                return [string trim [string range [readexp "${key}PRCF "] 10 19]]
2184            }
2185            pcut-set {
2186                if ![validreal value 10 5] {return 0}
2187                setexp "${key}PRCF " $value 11 10
2188                # set the powpref warning (1 = suggested)
2189                catch {
2190                    global expgui
2191                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2192                    set msg "Profile parameters" 
2193                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2194                        append expgui(needpowpref_why) "\t$msg were changed\n"
2195                    }
2196                }
2197            }
2198            pdamp-get {
2199                set val [string range [readexp "${key}PRCF "] 24 24]
2200                if {$val == " "} {return 0}
2201                return $val
2202            }
2203            pdamp-set {
2204                setexp "${key}PRCF   " $value 25 1
2205            }
2206            pterm*-get {
2207                regsub pterm $parm {} num
2208                set f1 [expr {15*(($num - 1) % 4)}]
2209                set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
2210                set line  [expr {1 + ($num - 1) / 4}]
2211                return [string trim [string range [readexp "${key}PRCF $line"] $f1 $f2] ]
2212            }
2213            pterm*-set {
2214                if ![validreal value 15 6] {return 0}
2215                regsub pterm $parm {} num
2216                set f1 [expr {1+ 15*(($num - 1) % 4)}]
2217                set line  [expr {1 + ($num - 1) / 4}]
2218                setexp "${key}PRCF $line" $value $f1 15
2219                # set the powpref warning (1 = suggested)
2220                catch {
2221                    global expgui
2222                    if {$expgui(needpowpref) == 0} {set expgui(needpowpref) 1}
2223                    set msg "Profile parameters" 
2224                    if {[string first $msg $expgui(needpowpref_why)] == -1} {
2225                        append expgui(needpowpref_why) "\t$msg were changed\n"
2226                    }
2227                }
2228            }
2229            pref*-get {
2230                regsub pref $parm {} num
2231                set f [expr {24+$num}]
2232                if {[string toupper [string range [readexp "${key}PRCF  "] $f $f]] == "Y"} {
2233                    return 1
2234                }
2235                return 0
2236            }
2237            pref*-set {
2238                regsub pref $parm {} num
2239                set f [expr {25+$num}]
2240                if $value {
2241                    setexp ${key}PRCF "Y" $f 1
2242                } else {
2243                    setexp ${key}PRCF "N" $f 1
2244                }           
2245            }
2246            POnaxis-get {
2247                set val [string trim \
2248                        [string range [readexp "${key}NAXIS"] 0 4]]
2249                if {$val == ""} {return 0}
2250                return $val
2251            }
2252            POnaxis-set {
2253                if ![validint value 5] {return 0}
2254                # there should be a NAXIS record, but if not make one
2255                if {![existsexp "${key}NAXIS"]} {
2256                    makeexprec "${key}NAXIS"
2257                }
2258                setexp "${key}NAXIS  " $value 1 5
2259            }
2260            default {
2261                set msg "Unsupported hapinfo access: parm=$parm action=$action"
2262                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
2263            }
2264        }
2265    }
2266    return 1
2267}
2268
2269#  read fixed constraints for a phase
2270proc atom_constraint_read {phase} {
2271    set fixlist ""
2272    foreach k {1 2 3 4 5 6 7 8 9} {
2273        set key [format "LEQV HOLD%1d%2d" $phase $k]
2274        set line [readexp $key]
2275        foreach j {2 10 18 26 34 42 50 58} {
2276            set fix_param [string range $line $j [expr $j+7]]
2277            if {[string trim $fix_param] == ""} {return $fixlist}
2278            lappend fixlist $fix_param
2279        }
2280    }
2281    return $fixlist
2282}
2283
2284# load all atom constraints into global array fix_param
2285proc atom_constraint_load { } {
2286    catch {unset ::fix_param}
2287    foreach i $::expmap(phaselist) {
2288        set temp [atom_constraint_read $i]
2289        foreach j $temp {
2290            set atomnum [string trim [string range $j 2 3]]
2291            set param [string trim [string range $j 4 6]]
2292            set ::fix_param($i,$atomnum,$param) 1   
2293        }
2294    }
2295}
2296
2297# returns 1 if the specified variable is fixed
2298proc atom_constraint_get {phase atom type} {
2299    if {[array names ::fix_param "$phase,$atom,$type"] == ""} {
2300        return 0
2301    }
2302    return 1
2303}
2304
2305proc atom_constraint_set {phase atomlist type mode} {
2306    foreach atom $atomlist {
2307        set key "$phase,$atom,$type"
2308        if {$mode} {
2309            set ::fix_param($key) 1
2310        } else {
2311            array unset ::fix_param $key
2312        }
2313    } 
2314    set fixlist {}
2315    foreach key [array names ::fix_param "$phase,*"] {
2316        foreach {j atom parm} [split $key ","] {}
2317        lappend fixlist \
2318            [format "%1s %+2s%-4s" $phase $atom $parm]
2319    }
2320    foreach key [array names ::exparray "LEQV HOLD$phase*"] {
2321        delexp $key
2322    }
2323    set k 0
2324    set j 1
2325    set line ""
2326    foreach fix $fixlist {
2327        incr k
2328        append line $fix
2329        if {$k == 8} {
2330            set key [format "LEQV HOLD%1d%2d" $phase $j]
2331            makeexprec $key
2332            setexp $key $line 3 [expr ($k * 8) + 2]
2333            set k 0
2334            incr j
2335            set line ""
2336        }
2337    }
2338    if {$line != ""} {
2339        set key [format "LEQV HOLD%1d%2d" $phase $j]
2340        makeexprec $key
2341        setexp $key $line 3 [expr ($k * 8) + 2]
2342    }   
2343}
2344
2345
2346#  get a logical constraint
2347#
2348#  type action
2349#  -----------
2350#  atom get  number        returns a list of constraints.
2351#   "   set  number value  replaces a list of constraints
2352#                          (value is a list of constraints)
2353#   "   add  number value  inserts a new list of constraints
2354#                          (number is ignored)
2355#   "   delete number      deletes a set of constraint entries
2356# Each item in the list of constraints is composed of 4 items:
2357#              phase, atom, variable, multiplier
2358# If variable=UISO atom can be ALL, otherwise atom is a number
2359# legal variable names: FRAC, X, Y, Z, UISO, U11, U22, U33, U12, U23, U13,
2360#                       MX, MY, MZ
2361#
2362#  type action
2363#  -----------
2364#  profileXX get number         returns a list of constraints for term XX=1-36
2365#                               use number=0 to get # of defined
2366#                                  constraints for term XX
2367#   "        set number value   replaces a list of constraints
2368#                               (value is a list of constraints)
2369#   "        add number value   inserts a new list of constraints
2370#                               (number is ignored)
2371#   "        delete number      deletes a set of constraint entries
2372# Each item in the list of constraints is composed of 3 items:
2373#              phase-list, histogram-list, multiplier
2374# Note that phase-list and/or histogram-list can be ALL
2375#
2376#  type action
2377#  -----------
2378#  absorbX get number         returns a list of constraints for term X=1 or 2
2379#   returns a the number of constraints for number = 0
2380#   returns a list of lists {{hist mult} {hist mult} ...}
2381
2382#  absorbX set number value   replaces a list of constraints
2383#      number corresponds to a specific constraint see "absorbX get 0"
2384#      value is a list of lists {{hist mult} {hist mult} ...}
2385#  absorbX add number value   inserts a new list of constraints
2386#                               (number is ignored)
2387#  absorbX  delete number      deletes a set of constraint entries and renumbers
2388# note that hist can be:
2389#      a histogram number (such as 2) or
2390#      range of histograms (such as 1:10 or 11:99, etc.) or
2391#      the string "ALL"
2392
2393proc constrinfo {type action number "value {}"} {
2394    global expmap
2395    if {[lindex $expmap(phasetype) 0] == 4} {
2396        set mm 1
2397    } else {
2398        set mm 0
2399    }
2400    switch -glob ${type}-$action {
2401        atom-get {
2402            # does this constraint exist?
2403            set key [format "LNCN%4d%4d" $number 1]
2404            if {![existsexp $key]} {return -1}
2405            set clist {}
2406            for {set i 1} {$i < 999} {incr i} {
2407                set key [format "LNCN%4d%4d" $number $i]
2408                if {![existsexp $key]} break
2409                set line [readexp $key]
2410                set j1 2
2411                set j2 17
2412                set seg [string range $line $j1 $j2]
2413                while {[string trim $seg] != ""} {
2414                    set p [string range $seg 0 0]
2415                    if {$p == 1 && $mm} {
2416                        set atom [string trim [string range $seg 1 4]]
2417                        set var [string trim [string range $seg 5 7]]
2418                        if {$atom == "ALL"} {
2419                            set var UIS
2420                        } else {
2421                            scan $atom %x atom
2422                        }
2423                        lappend clist [list $p $atom $var \
2424                                [string trim [string range $seg 8 end]]]
2425                    } else {
2426                        lappend clist [list $p \
2427                                [string trim [string range $seg 1 3]] \
2428                                [string trim [string range $seg 4 7]] \
2429                                [string trim [string range $seg 8 end]]]
2430                    }
2431                    incr j1 16
2432                    incr j2 16
2433                    set seg [string range $line $j1 $j2]
2434                }
2435            }
2436            return $clist
2437        }
2438        atom-set {
2439            # delete records for current constraint
2440            for {set i 1} {$i < 999} {incr i} {
2441                set key [format "LNCN%4d%4d" $number $i]
2442                if {![existsexp $key]} break
2443                delexp $key
2444            }
2445            set line {}
2446            set i 1
2447            foreach tuple $value {
2448                set p [lindex $tuple 0]
2449                if {$p == 1 && $mm && \
2450                        [string toupper [lindex $tuple 1]] == "ALL"} {
2451                    set seg [format %1dALL UIS%8.4f \
2452                            [lindex $tuple 0] \
2453                            [lindex $tuple 3]]
2454                } elseif {$p == 1 && $mm} {
2455                    set seg [eval format %1d%.4X%-3s%8.4f $tuple]
2456                } elseif {[string toupper [lindex $tuple 1]] == "ALL"} {
2457                    set seg [format %1dALL%-4s%8.4f \
2458                            [lindex $tuple 0] \
2459                            [lindex $tuple 2] \
2460                            [lindex $tuple 3]]
2461                } else {
2462                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
2463                }
2464                append line $seg
2465                if {[string length $line] > 50} {
2466                    set key  [format "LNCN%4d%4d" $number $i]
2467                    makeexprec $key
2468                    setexp $key $line 3 68
2469                    set line {}
2470                    incr i
2471                }
2472            }
2473            if {$line != ""} {
2474                set key  [format "LNCN%4d%4d" $number $i]
2475                makeexprec $key
2476                setexp $key $line 3 68
2477            }
2478            return
2479        }
2480        atom-add {
2481            # loop over defined constraints
2482            for {set j 1} {$j < 9999} {incr j} {
2483                set key [format "LNCN%4d%4d" $j 1]
2484                if {![existsexp $key]} break
2485            }
2486            set number $j
2487            # save the constraint
2488            set line {}
2489            set i 1
2490            foreach tuple $value {
2491                set p [lindex $tuple 0]
2492                if {$p == 1 && $mm && \
2493                        [string toupper [lindex $tuple 1]] == "ALL"} {
2494                    set seg [format %1dALL UIS%8.4f \
2495                            [lindex $tuple 0] \
2496                            [lindex $tuple 3]]
2497                } elseif {$p == 1 && $mm} {
2498                    set seg [eval format %1d%.4X%-3s%8.4f $tuple]
2499                } elseif {[string toupper [lindex $tuple 1]] == "ALL"} {
2500                    set seg [format %1dALL%-4s%8.4f \
2501                            [lindex $tuple 0] \
2502                            [lindex $tuple 2] \
2503                            [lindex $tuple 3]]
2504                } else {
2505                    set seg [eval format %1d%3d%-4s%8.4f $tuple]
2506                }
2507                append line $seg
2508                if {[string length $line] > 50} {
2509                    set key  [format "LNCN%4d%4d" $number $i]
2510                    makeexprec $key
2511                    setexp $key $line 3 68
2512                    set line {}
2513                    incr i
2514                }
2515            }
2516            if {$line != ""} {
2517                set key  [format "LNCN%4d%4d" $number $i]
2518                makeexprec $key
2519                setexp $key $line 3 68
2520            }
2521            return
2522        }
2523        atom-delete {
2524            for {set j $number} {$j < 9999} {incr j} {
2525                # delete records for current constraint
2526                for {set i 1} {$i < 999} {incr i} {
2527                    set key [format "LNCN%4d%4d" $j $i]
2528                    if {![existsexp $key]} break
2529                    delexp $key
2530                }
2531                # now copy records, from the next entry, if any
2532                set j1 $j
2533                incr j1
2534                set key1 [format "LNCN%4d%4d" $j1 1]
2535                # if there is no record, there is nothing to copy -- done
2536                if {![existsexp $key1]} return
2537                for {set i 1} {$i < 999} {incr i} {
2538                    set key1 [format "LNCN%4d%4d" $j1 $i]
2539                    if {![existsexp $key1]} break
2540                    set key  [format "LNCN%4d%4d" $j  $i]
2541                    makeexprec $key
2542                    setexp $key [readexp $key1] 1 68
2543                }
2544            }
2545        }
2546        profile*-delete {
2547            regsub profile $type {} term
2548            if {$term < 10} {
2549                set term " $term"
2550            }
2551            set key "LEQV PF$term   "
2552            # return nothing if no term exists
2553            if {![existsexp $key]} {return 0}
2554
2555            # number of constraint terms
2556            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2557            # don't delete a non-existing entry
2558            if {$number > $nterms} {return 0}
2559            set val [expr {$nterms - 1}]
2560            validint val 5
2561            setexp $key $val 1 5
2562            for {set i1 $number} {$i1 < $nterms} {incr i1} {
2563                set i2 [expr {1 + $i1}]
2564                # move the contents of constraint #i2 -> i1
2565                if {$i1 > 9} {
2566                    set k1 [expr {($i1+1)/10}]
2567                    set l1 $i1
2568                } else {
2569                    set k1 " "
2570                    set l1 " $i1"
2571                }
2572                set key1 "LEQV PF$term  $k1"
2573                # number of constraint lines for #i1
2574                set n1 [string trim [string range [readexp ${key1}] \
2575                        [expr {($i1%10)*5}] [expr {4+(($i1%10)*5)}]] ]
2576                if {$i2 > 9} {
2577                    set k2 [expr {($i2+1)/10}]
2578                    set l2 $i2
2579                } else {
2580                    set k2 " "
2581                    set l2 " $i2"
2582                }
2583                set key2 "LEQV PF$term  $k2"
2584                # number of constraint lines for #i2
2585                set n2 [string trim [string range [readexp ${key2}] \
2586                        [expr {($i2%10)*5}] [expr {4+(($i2%10)*5)}]] ]
2587                set val $n2
2588                validint val 5
2589                # move the # of terms
2590                setexp $key1 $val [expr {1+(($i1%10)*5)}] 5
2591                # move the terms
2592                for {set j 1} {$j <= $n2} {incr j 1} {
2593                    set key "LEQV PF${term}${l1}$j"
2594                    makeexprec $key
2595                    setexp $key [readexp "LEQV PF${term}${l2}$j"] 1 68
2596                }
2597                # delete any remaining lines
2598                for {set j [expr {$n2+1}]} {$j <= $n1} {incr j 1} {
2599                    delexp "LEQV PF${term}${l1}$j"
2600                }
2601            }
2602
2603            # clear the last term
2604            if {$nterms > 9} {
2605                set i [expr {($nterms+1)/10}]
2606            } else {
2607                set i " "
2608            }
2609            set key "LEQV PF$term  $i"
2610            set cb [expr {($nterms%10)*5}]
2611            set ce [expr {4+(($nterms%10)*5)}]
2612            set n2 [string trim [string range [readexp ${key}] $cb $ce] ]
2613            incr cb
2614            setexp $key "     " $cb 5
2615            # delete any remaining lines
2616            for {set j 1} {$j <= $n2} {incr j 1} {
2617                delexp "LEQV PF${term}${nterms}$j"
2618            }
2619        }
2620        profile*-set {
2621            regsub profile $type {} term
2622            if {$term < 10} {
2623                set term " $term"
2624            }
2625            set key "LEQV PF$term   "
2626            # get number of constraint terms
2627            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2628            # don't change a non-existing entry
2629            if {$number > $nterms} {return 0}
2630            if {$number > 9} {
2631                set k1 [expr {($number+1)/10}]
2632                set l1 $number
2633            } else {
2634                set k1 " "
2635                set l1 " $number"
2636            }
2637            set key1 "LEQV PF$term  $k1"
2638            # old number of constraint lines
2639            set n1 [string trim [string range [readexp ${key1}] \
2640                    [expr {($number%10)*5}] [expr {4+(($number%10)*5)}]] ]
2641            # number of new constraints
2642            set j2 [llength $value]
2643            # number of new constraint lines
2644            set val [set n2 [expr {($j2 + 2)/3}]]
2645            # store the new # of lines
2646            validint val 5
2647            setexp $key1 $val [expr {1+(($number%10)*5)}] 5
2648
2649            # loop over the # of lines in the old or new, whichever is greater
2650            set v0 0
2651            for {set j 1} {$j <= [expr {($n1 > $n2) ? $n1 : $n2}]} {incr j 1} {
2652                set key "LEQV PF${term}${l1}$j"
2653                # were there more lines in the old?
2654                if {$j > $n2} {
2655                    # this line is not needed
2656                    if {$j % 3 == 1} {
2657                        delexp $key
2658                    }
2659                    continue
2660                }
2661                # are we adding new lines?
2662                if {$j > $n1} {
2663                    makeexprec $key
2664                }
2665                # add the three constraints to the line
2666                foreach s {3 23 43} \
2667                        item [lrange $value $v0 [expr {2+$v0}]] {
2668                    if {$item != ""} {
2669                        set val [format %-10s%9.3f \
2670                                [lindex $item 0],[lindex $item 1] \
2671                                [lindex $item 2]]
2672                        setexp $key $val $s 19
2673                    } else {
2674                        setexp $key " " $s 19
2675                    }
2676                }
2677                incr v0 3
2678            }
2679        }
2680        profile*-add {
2681            regsub profile $type {} term
2682            if {$term < 10} {
2683                set term " $term"
2684            }
2685            set key "LEQV PF$term   "
2686            if {![existsexp $key]} {makeexprec $key}
2687            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2688            if {$nterms == ""} {
2689                set nterms 1
2690            } elseif {$nterms >= 99} {
2691                return 0
2692            } else {
2693                incr nterms
2694            }
2695            # store the new # of constraints
2696            set val $nterms
2697            validint val 5
2698            setexp $key $val 1 5
2699
2700            if {$nterms > 9} {
2701                set k1 [expr {($nterms+1)/10}]
2702                set l1 $nterms
2703            } else {
2704                set k1 " "
2705                set l1 " $nterms"
2706            }
2707            set key1 "LEQV PF$term  $k1"
2708
2709            # number of new constraints
2710            set j2 [llength $value]
2711            # number of new constraint lines
2712            set val [set n2 [expr {($j2 + 2)/3}]]
2713            # store the new # of lines
2714            validint val 5
2715            setexp $key1 $val [expr {1+(($nterms%10)*5)}] 5
2716
2717            # loop over the # of lines to be added
2718            set v0 0
2719            for {set j 1} {$j <= $n2} {incr j 1} {
2720                set key "LEQV PF${term}${l1}$j"
2721                makeexprec $key
2722                # add the three constraints to the line
2723                foreach s {3 23 43} \
2724                        item [lrange $value $v0 [expr {2+$v0}]] {
2725                    if {$item != ""} {
2726                        set val [format %-10s%9.3f \
2727                                [lindex $item 0],[lindex $item 1] \
2728                                [lindex $item 2]]
2729                        setexp $key $val $s 19
2730                    } else {
2731                        setexp $key " " $s 19
2732                    }
2733                }
2734                incr v0 3
2735            }
2736        }
2737        profile*-get {
2738            regsub profile $type {} term
2739            if {$term < 10} {
2740                set term " $term"
2741            }
2742            if {$number > 9} {
2743                set i [expr {($number+1)/10}]
2744            } else {
2745                set i " "
2746            }
2747            set key "LEQV PF$term  $i"
2748            # return nothing if no term exists
2749            if {![existsexp $key]} {return 0}
2750            # number of constraint lines
2751           
2752            set numline [string trim [string range [readexp ${key}] \
2753                    [expr {($number%10)*5}] [expr {4+(($number%10)*5)}]] ]
2754            if {$number == 0} {return $numline}
2755            set clist {}
2756            if {$number < 10} {
2757                set number " $number"
2758            }
2759            for {set i 1} {$i <= $numline} {incr i} {
2760                set key "LEQV PF${term}${number}$i"
2761                set line [readexp ${key}]
2762                foreach s {1 21 41} e {20 40 60} {
2763                    set seg [string range $line $s $e]
2764                    if {[string trim $seg] == ""} continue
2765                    # parse the string segment
2766                    set parse [regexp { *([0-9AL]+),([0-9AL]+) +([0-9.]+)} \
2767                            $seg junk phase hist mult]
2768                    # was parse successful
2769                    if {!$parse} {continue}
2770                    lappend clist [list $phase $hist $mult]
2771                }
2772            }
2773            return $clist
2774        }
2775        absorb*-delete {
2776            regsub absorb $type {} term
2777            set key "LEQV ABS$term   "
2778            if {! [existsexp $key]} {return 0}
2779            # current number of constraints
2780            set nterm [string trim [string range [readexp $key] 0 5]]
2781            if {$nterm == ""} {return 0}
2782            # does the entry exist?
2783            if {$nterm < $number} {
2784                puts "deleted!"
2785                return $nterm
2786            }
2787            for {set target $number} {$target < $nterm} {incr target} {
2788                set source [expr {$target + 1}]
2789                set recs [GetAbsCount $term $source]
2790                SetAbsCount $term $target [expr {3*$recs}]
2791                validint source 2
2792                validint target 2
2793                for {set i 1} {$i <= $recs} {incr i} {
2794                    set keyin "LEQV ABS${term}${source}$i"
2795                    set keyout "LEQV ABS${term}${target}$i"
2796                    set ::exparray($keyout) $::exparray($keyin)
2797                }
2798            }
2799            SetAbsCount $term $nterm 0
2800            # delete the last entry
2801            validint nterm 2
2802            foreach i {1 2 3 4 5 6 7 8 9} {
2803                set key "LEQV ABS${term}${nterm}$i"
2804                delexp $key
2805            }
2806            # decrease the count by one
2807            set nterm [expr {[string trim $nterm] - 1}]
2808            if {$nterm == 0} {
2809                delexp "LEQV ABS$term   "
2810            } else {
2811                validint nterm 5
2812                setexp "LEQV ABS$term   " $nterm 1 5                   
2813            }
2814            return [string trim $nterm]
2815
2816            if {$term < 10} {
2817                set term " $term"
2818            }
2819            set key "LEQV PF$term   "
2820            # return nothing if no term exists
2821            if {![existsexp $key]} {return 0}
2822
2823            # number of constraint terms
2824            set nterms [string trim [string range [readexp ${key}] 0 4] ]
2825            # don't delete a non-existing entry
2826            if {$number > $nterms} {return 0}
2827            set val [expr {$nterms - 1}]
2828            validint val 5
2829            setexp $key $val 1 5
2830            for {set i1 $number} {$i1 < $nterms} {incr i1} {
2831                set i2 [expr {1 + $i1}]
2832                # move the contents of constraint #i2 -> i1
2833                if {$i1 > 9} {
2834                    set k1 [expr {($i1+1)/10}]
2835                    set l1 $i1
2836                } else {
2837                    set k1 " "
2838                    set l1 " $i1"
2839                }
2840                set key1 "LEQV PF$term  $k1"
2841                # number of constraint lines for #i1
2842                set n1 [string trim [string range [readexp ${key1}] \
2843                        [expr {($i1%10)*5}] [expr {4+(($i1%10)*5)}]] ]
2844                if {$i2 > 9} {
2845                    set k2 [expr {($i2+1)/10}]
2846                    set l2 $i2
2847                } else {
2848                    set k2 " "
2849                    set l2 " $i2"
2850                }
2851                set key2 "LEQV PF$term  $k2"
2852                # number of constraint lines for #i2
2853                set n2 [string trim [string range [readexp ${key2}] \
2854                        [expr {($i2%10)*5}] [expr {4+(($i2%10)*5)}]] ]
2855                set val $n2
2856                validint val 5
2857                # move the # of terms
2858                setexp $key1 $val [expr {1+(($i1%10)*5)}] 5
2859                # move the terms
2860                for {set j 1} {$j <= $n2} {incr j 1} {
2861                    set key "LEQV PF${term}${l1}$j"
2862                    makeexprec $key
2863                    setexp $key [readexp "LEQV PF${term}${l2}$j"] 1 68
2864                }
2865                # delete any remaining lines
2866                for {set j [expr {$n2+1}]} {$j <= $n1} {incr j 1} {
2867                    delexp "LEQV PF${term}${l1}$j"
2868                }
2869            }
2870
2871            # clear the last term
2872            if {$nterms > 9} {
2873                set i [expr {($nterms+1)/10}]
2874            } else {
2875                set i " "
2876            }
2877            set key "LEQV PF$term  $i"
2878            set cb [expr {($nterms%10)*5}]
2879            set ce [expr {4+(($nterms%10)*5)}]
2880            set n2 [string trim [string range [readexp ${key}] $cb $ce] ]
2881            incr cb
2882            setexp $key "     " $cb 5
2883            # delete any remaining lines
2884            for {set j 1} {$j <= $n2} {incr j 1} {
2885                delexp "LEQV PF${term}${nterms}$j"
2886            }
2887        }
2888        absorb*-set {
2889            regsub absorb $type {} term
2890            if {$number < 1} return   
2891            # delete old records
2892            set l [GetAbsCount $term $number]
2893            set num $number
2894            validint num 2
2895            for {set i 1} {$i <= $l} {incr i} {
2896                delexp "LEQV ABS${term}${num}$i"
2897            }
2898            # record the new number of records
2899            SetAbsCount $term $number [llength $value]
2900            # save the new records
2901            set i 1
2902            set offh 2
2903            set offm 14
2904            foreach set $value {
2905                set hist [string trim [lindex $set 0]]
2906                set mult [string trim [lindex $set 1]]
2907                validreal mult 8 4
2908                set key "LEQV ABS${term}${num}$i"
2909                if {$offh == 2} {
2910                    makeexprec $key
2911                }
2912                setexp $key $hist [expr {$offh+1}] 11
2913                setexp $key $mult [expr {$offm+1}] 8
2914                incr offh 21
2915                incr offm 21
2916                if {$offm > 67} {
2917                    incr i
2918                    set offh 2
2919                    set offm 14
2920                }
2921            }
2922            return
2923        }
2924        absorb*-add {
2925            regsub absorb $type {} term
2926            set key "LEQV ABS$term   "
2927            if {! [existsexp $key]} {makeexprec $key}
2928            # current number of constraints
2929            set nterm [string trim [string range [readexp $key] 0 5]]
2930            if {$nterm == ""} {set nterm 0}
2931            if {$nterm >= 99} {
2932                return $nterm
2933            }
2934            incr nterm
2935            validint nterm 5
2936            setexp $key $nterm 1 5
2937            constrinfo $type set [string trim $nterm] $value
2938            return [string trim $nterm]
2939        }
2940        absorb*-get {
2941            regsub absorb $type {} term
2942            # no constraints, return blank
2943            set key "LEQV ABS$term   "
2944            if {! [existsexp $key]} {return ""}
2945            # requesting number of constraints
2946            if {$number == 0} {
2947                set l [string trim [string range [readexp ${key}] 0 5]]
2948                if {$l == ""} {return 0}
2949                return $l
2950            }
2951            #
2952            if {$number > 9} {
2953                set num $number
2954                set i [expr {($number+1)/10}]
2955                set off [expr {5*($number % 10)}]
2956                set key "LEQV ABS$term  $i"
2957            } else {
2958                set num " $number"
2959                set i " "
2960                set off [expr {5*($number % 10)}]
2961            }
2962            set off1 [expr {$off + 5}]
2963            set l [string trim [string range [readexp ${key}] $off $off1]]
2964            if {$l == ""} {return {}}
2965            # now look up those records
2966            set res {}
2967            for {set i 1} {$i <= $l} {incr i} {
2968                set key "LEQV ABS${term}${num}$i"
2969                for {set j 0} {$j < 3} {incr j} {
2970                    set off [expr {2 + 21*$j}]
2971                    set off1 [expr {$off + 11}]
2972                    set hist [string trim [string range [readexp ${key}] $off $off1]]
2973                    set off [expr {14 + 21*$j}]
2974                    set off1 [expr {$off + 7}]
2975                    set mult [string trim [string range [readexp ${key}] $off $off1]]
2976                    if {$mult == ""} break
2977                    lappend res [list $hist $mult]
2978                }
2979            }
2980            return $res
2981        }
2982        default {
2983            set msg "Unsupported constrinfo access: type=$type action=$action"
2984            tk_dialog .badexp "Error in readexp access" $msg error 0 OK
2985        }
2986
2987    }
2988}
2989proc GetAbsCount {term number} {
2990    if {$number > 9} {
2991        set num $number
2992        set off [expr {5*($number % 10)}]
2993        set i [expr {($number+1)/10}]
2994        set key "LEQV ABS$term  $i"
2995    } else {
2996        set num " $number"
2997        set off [expr {5*($number % 10)}]
2998        set key "LEQV ABS$term   "
2999    }
3000    set off1 [expr {$off + 5}]
3001    set l [string trim [string range [readexp ${key}] $off $off1]]
3002    if {$l == ""} {set l 0}
3003    return $l
3004}
3005proc SetAbsCount {term number len} {
3006    if {$number > 9} {
3007        set num $number
3008        set off [expr {1 + 5*($number % 10)}]
3009        set i [expr {($number+1)/10}]
3010        set key "LEQV ABS$term  $i"
3011    } else {
3012        set num " $number"
3013        set off [expr {1 + 5*($number % 10)}]
3014        set key "LEQV ABS$term   "
3015    }
3016    set l [expr {($len + 2)/3}]
3017    set val $l
3018    validint val 5
3019    setexp $key $val $off 5
3020}
3021
3022# read the default profile information for a histogram
3023# use: profdefinfo hist set# parm action
3024
3025#     proftype -- profile function number
3026#     profterms -- number of profile terms
3027#     pdamp -- damping value for the profile (*)
3028#     pcut -- cutoff value for the profile (*)
3029#     pterm$n -- profile term #n
3030#     pref$n -- refinement flag value for profile term #n (*)
3031
3032proc profdefinfo {hist set parm "action get"} {
3033    global expgui
3034    if {$hist < 10} {
3035        set key "HST  $hist"
3036    } else {
3037        set key "HST $hist"
3038    }
3039    switch -glob ${parm}-$action {
3040        proftype-get {
3041            set val [string range [readexp "${key}PRCF$set"] 0 4]
3042            if {$val == " "} {return 0}
3043            return $val
3044        }
3045        profterms-get {
3046            set val [string range [readexp "${key}PRCF$set"] 5 9]
3047            if {$val == " "} {return 0}
3048            return $val
3049        }
3050        pcut-get {
3051            return [string trim [string range [readexp "${key}PRCF$set"] 10 19]]
3052        }
3053        pdamp-get {
3054                set val [string range [readexp "${key}PRCF$set"] 24 24]
3055            if {$val == " "} {return 0}
3056            return $val
3057        }
3058        pterm*-get {
3059            regsub pterm $parm {} num
3060            set f1 [expr {15*(($num - 1) % 4)}]
3061            set f2 [expr {15*(1 + ($num - 1) % 4)-1}]
3062            set line  [expr {1 + ($num - 1) / 4}]
3063            return [string trim [string range [\
3064                        readexp "${key}PRCF${set}$line"] $f1 $f2] ]
3065        }
3066        pref*-get {
3067            regsub pref $parm {} num
3068            set f [expr {24+$num}]
3069            if {[string toupper [string range [readexp "${key}PRCF$set"] $f $f]] == "Y"} {
3070                return 1
3071            }
3072            return 0
3073        }
3074        default {
3075            set msg "Unsupported profdefinfo access: parm=$parm action=$action"
3076            tk_dialog .badexp "Code Error" $msg error 0 Exit
3077        }
3078    }
3079}
3080
3081# get March-Dollase preferred orientation information
3082# use MDprefinfo hist phase axis-number parm action value
3083#    ratio    -- ratio of xtallites in PO direction vs random (>1 for more)
3084#    fraction -- fraction in this direction, when more than one axis is used
3085#    h k & l  -- indices of P.O. axis
3086#    ratioref -- flag to vary ratio
3087#    fracref  -- flag to vary fraction
3088#    damp     -- damping value
3089#    type     -- model type (0 = P.O. _|_ to beam, 1 = || to beam)
3090#    new      -- creates a new record with default values (set only)
3091proc MDprefinfo {histlist phaselist axislist parm "action get" "value {}"} {
3092    foreach phase $phaselist hist $histlist axis $axislist {
3093        if {$phase == ""} {set phase [lindex $phaselist end]}
3094        if {$hist == ""} {set hist [lindex $histlist end]}
3095        if {$axis == ""} {set axis [lindex $axislist end]}
3096        if {$hist < 10} {
3097            set hist " $hist"
3098        }
3099        if {$axis > 9} {
3100            set axis "0"
3101        }
3102        set key "HAP${phase}${hist}PREFO${axis}"
3103        switch -glob ${parm}-$action {
3104            ratio-get {
3105                return [string trim [string range [readexp $key] 0 9]]
3106            }
3107            ratio-set {
3108                if ![validreal value 10 6] {return 0}
3109                setexp $key $value 1 10
3110            }
3111            fraction-get {
3112                return [string trim [string range [readexp $key] 10 19]]
3113            }
3114            fraction-set {
3115                if ![validreal value 10 6] {return 0}
3116                setexp $key $value 11 10
3117            }
3118            h-get {
3119                set h [string trim [string range [readexp $key] 20 29]]
3120                # why not allow negative h values?
3121                #               if {$h < 1} {return 0}
3122                return $h
3123            }
3124            h-set {
3125                if ![validreal value 10 2] {return 0}
3126                setexp $key $value 21 10
3127            }
3128            k-get {
3129                set k [string trim [string range [readexp $key] 30 39]]
3130                #               if {$k < 1} {return 0}
3131                return $k
3132            }
3133            k-set {
3134                if ![validreal value 10 2] {return 0}
3135                setexp $key $value 31 10
3136            }
3137            l-get {
3138                set l [string trim [string range [readexp $key] 40 49]]
3139                #if {$l < 1} {return 0}
3140                return $l
3141            }
3142            l-set {
3143                if ![validreal value 10 2] {return 0}
3144                setexp $key $value 41 10
3145            }
3146            ratioref-get {
3147                if {[string toupper \
3148                        [string range [readexp $key] 53 53]] == "Y"} {
3149                    return 1
3150                }
3151                return 0
3152            }
3153            ratioref-set {
3154                if $value {
3155                    setexp $key "Y" 54 1
3156                } else {
3157                    setexp $key "N" 54 1
3158                }
3159            }
3160            fracref-get {
3161                if {[string toupper \
3162                        [string range [readexp $key] 54 54]] == "Y"} {
3163                    return 1
3164                }
3165                return 0
3166            }
3167            fracref-set {
3168                if $value {
3169                    setexp $key "Y" 55 1
3170                } else {
3171                    setexp $key "N" 55 1
3172              }
3173            }
3174            damp-get {
3175                set val [string trim [string range [readexp $key] 59 59]]
3176                if {$val == " "} {return 0}
3177                return $val
3178            }
3179            damp-set {
3180                setexp $key $value 60 1
3181            }
3182            type-get {
3183                set val [string trim [string range [readexp $key] 64 64]]
3184                if {$val == " "} {return 0}
3185                return $val
3186            }
3187            type-set {
3188                # only valid settings are 0 & 1
3189                if {$value != "0" && $value != "1"} {set value "0"}
3190                setexp $key $value 65 1
3191            }
3192            new-set {
3193                makeexprec $key
3194                setexp $key \
3195                        {  1.000000  1.000000  0.000000  0.000000  1.000000   NN    0    0} \
3196                        1 68
3197            }
3198            default {
3199                set msg "Unsupported MDprefinfo access: parm=$parm action=$action"
3200                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3201            }
3202
3203        }
3204
3205    }
3206}
3207
3208# get list of defined atom types
3209proc AtmTypList {} {
3210    set natypes [readexp " EXPR  NATYP"]
3211    if {$natypes == ""} return
3212    set j 0
3213    set typelist {}
3214    for {set i 1} {$i <= $natypes} {incr i} {
3215        set key {this should never be matched}
3216        while {![existsexp $key]} {
3217            incr j
3218            if {$j > 99} {
3219                return $typelist
3220            } elseif {$j <10} {
3221                set key " EXPR ATYP $j"
3222            } else {
3223                set key " EXPR ATYP$j"
3224            }
3225        }
3226        lappend typelist [string trim [string range $::exparray($key) 2 9]]
3227    }
3228    return $typelist
3229}
3230
3231# read information about atom types
3232#     distrad    atomic distance search radius (get/set)
3233#     angrad     atomic angle search radius (get/set)
3234proc AtmTypInfo {parm atmtype "action get" "value {}"} {
3235    # first, search through the records to find the record matching the type
3236    set natypes [readexp " EXPR  NATYP"]
3237    if {$natypes == ""} return
3238    set j 0
3239    set typelist {}
3240    for {set i 1} {$i <= $natypes} {incr i} {
3241        set key {this should never be matched}
3242        while {![existsexp $key]} {
3243            incr j
3244            if {$j > 99} {
3245                return $typelist
3246            } elseif {$j <10} {
3247                set key " EXPR ATYP $j"
3248            } else {
3249                set key " EXPR ATYP$j"
3250            }
3251        }
3252        if {[string toupper $atmtype] == \
3253                [string toupper [string trim [string range $::exparray($key) 2 9]]]} break
3254        set key {}
3255    }
3256    if {$key == ""} {
3257        # atom type not found
3258        return {}
3259    }
3260    switch -glob ${parm}-$action {
3261        distrad-get {
3262            return [string trim [string range [readexp $key] 15 24]]
3263        }
3264        distrad-set {
3265            if ![validreal value 10 2] {return 0}
3266            setexp $key $value 16 10
3267        }
3268        angrad-get {
3269            return [string trim [string range [readexp $key] 25 34]]
3270        }
3271        angrad-set {
3272            if ![validreal value 10 2] {return 0}
3273            setexp $key $value 26 10
3274        }
3275        default {
3276            set msg "Unsupported AtmTypInfo access: parm=$parm action=$action"
3277            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3278        }
3279    }
3280}
3281# read default information about atom types (records copied to the .EXP file
3282# from the gsas/data/atomdata.dat file as AFAC ...
3283#     distrad returns a list of atom types (one or two letters) and
3284#                the corresponding distance
3285# note that these values are read only (no set option)
3286proc DefAtmTypInfo {parm} {
3287    set keys [array names ::exparray " AFAC *_SIZ"]
3288    set elmlist {}
3289    if {[llength $keys] <= 0} {return ""}
3290    foreach key $keys {
3291        lappend elmlist [string trim [string range $key 6 7]]
3292    }
3293    switch -glob ${parm} {
3294        distrad {
3295            set out {}
3296            foreach key $keys elm $elmlist {
3297                set val [string range $::exparray($key) 0 9]
3298                lappend out "$elm [string trim $val]"
3299            }
3300            return $out
3301        }
3302        angrad {
3303            set out {}
3304            foreach key $keys elm $elmlist {
3305                set val [string range $::exparray($key) 10 19]
3306                lappend out "$elm [string trim $val]"
3307            }
3308            return $out
3309        }
3310        default {
3311            set msg "Unsupported DefAtmTypInfo access: parm=$parm"
3312            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3313        }
3314    }
3315}
3316# write the .EXP file
3317proc expwrite {expfile} {
3318    global exparray
3319    set blankline \
3320     "                                                                        "
3321    set fp [open ${expfile} w]
3322    fconfigure $fp -translation crlf -encoding ascii
3323    set keylist [lsort [array names exparray]]
3324    # reorder the keys so that VERSION comes 1st
3325    set pos [lsearch -exact $keylist {     VERSION}]
3326    set keylist "{     VERSION} [lreplace $keylist $pos $pos]"
3327    foreach key $keylist {
3328        puts $fp [string range \
3329                "$key$exparray($key)$blankline" 0 79]
3330    }
3331    close $fp
3332}
3333
3334# history commands -- delete all but last $keep history records,
3335# renumber if $renumber is true
3336proc DeleteHistory {keep renumber} {
3337    global exparray
3338    foreach y [lrange [lsort -decreasing \
3339            [array names exparray {    HSTRY*}]] $keep end] {
3340        unset exparray($y)
3341    }
3342    if !$renumber return
3343    # renumber
3344    set i 0
3345    foreach y [lsort -increasing \
3346            [array names exparray {    HSTRY*}]] {
3347        set key [format "    HSTRY%3d" [incr i]]
3348        set exparray($key) $exparray($y)
3349        unset exparray($y)
3350    }
3351    # list all history
3352    #    foreach y [lsort -decreasing [array names exparray {    HSTRY*}]] {puts "$y $exparray($y)"}
3353}
3354
3355proc CountHistory {} {
3356    global exparray
3357    return [llength [array names exparray {    HSTRY*}]]
3358}
3359
3360# set the phase flags for histogram $hist to $plist
3361proc SetPhaseFlag {hist plist} {
3362    # make a 2 digit key -- hh
3363    if {$hist < 10} {
3364        set hh " $hist"
3365    } else {
3366        set hh $hist
3367    }
3368    set key "HST $hh NPHAS"
3369    set str {}
3370    foreach iph {1 2 3 4 5 6 7 8 9} {
3371        if {[lsearch $plist $iph] != -1} {
3372            append str {    1}
3373        } else {
3374            append str {    0}     
3375        }
3376    }
3377    setexp $key $str 1 68
3378}
3379
3380# erase atom $atom from phase $phase
3381# update the list of atom types, erasing the record if not needed.
3382proc EraseAtom {atom phase} {
3383    set type [atominfo $phase $atom type]
3384    if {$type == ""} return
3385    if {$atom < 10} {
3386        set key "CRS$phase  AT  $atom"
3387    } elseif {$atom < 100} {
3388        set key "CRS$phase  AT $atom"
3389    } else {
3390        set key "CRS$phase  AT$atom"
3391    }
3392    # delete the records for the atom
3393    global exparray
3394    foreach k [array names exparray ${key}*] {
3395        delexp $k
3396    }
3397    # change the number of atoms in the phase
3398    phaseinfo $phase natoms set [expr {[phaseinfo $phase natoms] -1}]
3399
3400    # now adjust numbers in "EXPR ATYP" records and delete, if needed.
3401    set natypes [readexp " EXPR  NATYP"]
3402    if {$natypes == ""} return
3403    set j 0
3404    for {set i 1} {$i <= $natypes} {incr i} {
3405        incr j
3406        if {$j <10} {
3407            set key " EXPR ATYP $j"
3408        } else {
3409            set key " EXPR ATYP$j"
3410        }
3411        while {![existsexp $key]} {
3412            incr j
3413            if {$j > 99} {
3414                return
3415            } elseif {$j <10} {
3416                set key " EXPR ATYP $j"
3417            } else {
3418                set key " EXPR ATYP$j"
3419            }
3420        }
3421        set keytype [string trim [string range $exparray($key) 2 9]]
3422        if {$type == $keytype} {
3423            # found the type record
3424            set val [string trim [string range $exparray($key) 10 14]]
3425            incr val -1
3426            # if this is the last reference, remove the record,
3427            # otherwise, decrement the counter
3428            if {$val <= 0} {
3429                incr natypes -1 
3430                validint natypes 5
3431                setexp " EXPR  NATYP" $natypes 1 5
3432                delexp $key
3433            } else {
3434                validint val 5
3435                setexp $key $val 11 5
3436            }
3437            return
3438        }
3439    }
3440}
3441
3442# compute equivalent anisotropic temperature factor for Uequiv
3443proc CalcAniso {phase Uequiv} {
3444    foreach var {a b c alpha beta gamma} {
3445        set $var [phaseinfo $phase $var]
3446    }
3447
3448    set G(1,1) [expr {$a * $a}]
3449    set G(2,2) [expr {$b * $b}]
3450    set G(3,3) [expr {$c * $c}]
3451    set G(1,2) [expr {$a * $b * cos($gamma*0.017453292519943)}]
3452    set G(2,1) $G(1,2)
3453    set G(1,3) [expr {$a * $c * cos($beta *0.017453292519943)}]
3454    set G(3,1) $G(1,3)
3455    set G(2,3) [expr {$b * $c * cos($alpha*0.017453292519943)}]
3456    set G(3,2) $G(2,3)
3457
3458    # Calculate the volume**2
3459    set v2 0.0
3460    foreach i {1 2 3} {
3461        set J [expr {($i%3) + 1}]
3462        set K [expr {(($i+1)%3) + 1}]
3463        set v2 [expr {$v2+ $G(1,$i)*($G(2,$J)*$G(3,$K)-$G(3,$J)*$G(2,$K))}]
3464    }
3465    if {$v2 > 0} {
3466        set v [expr {sqrt($v2)}]
3467        foreach i {1 2 3} {
3468            set i1 [expr {($i%3) + 1}]
3469            set i2 [expr {(($i+1)%3) + 1}]
3470            foreach j {1 2 3} {
3471                set j1 [expr {($j%3) + 1}]
3472                set j2 [expr {(($j+1)%3) + 1}]
3473                set C($j,$i) [expr {(\
3474                        $G($i1,$j1) * $G($i2,$j2) - \
3475                        $G($i1,$j2)  * $G($i2,$j1)\
3476                        )/ $v}]
3477            }
3478        }
3479        set A(1,2) [expr {0.5 * ($C(1,2)+$C(2,1)) / sqrt( $C(1,1)* $C(2,2) )}]
3480        set A(1,3) [expr {0.5 * ($C(1,3)+$C(3,1)) / sqrt( $C(1,1)* $C(3,3) )}]
3481        set A(2,3) [expr {0.5 * ($C(2,3)+$C(3,2)) / sqrt( $C(2,2)* $C(3,3) )}]
3482        foreach i {1 1 2} j {2 3 3} {
3483            set A($i,$j) [expr {0.5 * ($C($i,$j) + $C($j,$i)) / \
3484                    sqrt( $C($i,$i)* $C($j,$j) )}]
3485            # clean up roundoff
3486            if {abs($A($i,$j)) < 1e-5} {set A($i,$j) 0.0}
3487        }
3488    } else {
3489        set A(1,2) 0.0
3490        set A(1,3) 0.0
3491        set A(2,3) 0.0
3492    }
3493    return "$Uequiv $Uequiv $Uequiv \
3494            [expr {$Uequiv * $A(1,2)}] \
3495            [expr {$Uequiv * $A(1,3)}] \
3496            [expr {$Uequiv * $A(2,3)}]"
3497}
3498
3499# read/edit soft (distance) restraint info
3500#  parm:
3501#    weight -- histogram weight (factr) *
3502#    restraintlist -- list of restraints *
3503#  action: get (default) or set
3504#  value: used only with set
3505#  * =>  read+write supported
3506proc SoftConst {parm "action get" "value {}"} {
3507    set HST {}
3508    # look for RSN record
3509    set n 0
3510    for {set i 0} {$i < $::expmap(nhst)} {incr i} {
3511        set ihist [expr {$i + 1}]
3512        if {[expr {$i % 12}] == 0} {
3513            incr n
3514            set line [readexp " EXPR  HTYP$n"]
3515            if {$line == ""} {
3516                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
3517                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3518            }
3519            set j 0
3520        } else {
3521            incr j
3522        }
3523        if {[string range $line [expr 2+5*$j] [expr 5*($j+1)]] == "RSN "} {
3524            set HST $ihist
3525        }
3526    }
3527    if {$HST <=9} {
3528        set key "HST  $HST"
3529    } else {
3530        set key "HST $HST"
3531    }
3532    if {$HST == "" && $action == "set"} {
3533        # no RSN found need to add the soft constr. histogram
3534        # increment number of histograms
3535        set hst [string trim [string range [readexp { EXPR  NHST }] 0 4]]
3536        incr hst
3537        set HST $hst
3538        if ![validint hst 5] {return 0}
3539        setexp  { EXPR  NHST } $hst 1 5
3540        # add to EXPR HTYPx rec, creating if needed
3541        set n [expr { 1+ (($HST - 1) / 12) }]
3542        set key " EXPR  HTYP$n"
3543        if {[array names ::exparray $key] == ""} {
3544            makeexprec $key
3545        }
3546        setexp $key "RSN " [expr 3 + 5*(($HST-1) % 12)] 5
3547        # create other HST  xx recs
3548        if {$HST <=9} {
3549            set key "HST  $HST"
3550        } else {
3551            set key "HST $HST"
3552        }
3553        makeexprec "$key  HNAM"
3554        setexp "$key  HNAM" "Bond distance restraints" 3 24
3555        makeexprec "$key FACTR"
3556        makeexprec "$key NBNDS"
3557        mapexp
3558    } elseif {$HST == ""} {
3559        if $::expgui(debug) {puts "no restraints"}
3560        return "1"
3561    }
3562
3563    switch -glob ${parm}-$action {
3564        weight-get {
3565            return [string trim [string range [readexp "$key FACTR"] 0 14]]
3566        }
3567        weight-set {
3568            # update FACTR
3569            if ![validreal value 15 6] {return 0}
3570            setexp "$key FACTR" $value 1 15
3571        }
3572        restraintlist-get {
3573            set ncons [string trim [string range [readexp "$key NBNDS"] 0 4]]
3574            set conslist {}
3575            for {set i 1} {$i <= $ncons} {incr i} {
3576                set fi [string toupper [format %.4x $i]]
3577                set line [readexp "${key}BD$fi"]
3578                set const {}
3579                foreach len {3 5 5 3 3 3 3 3 6 6} {
3580                  set lenm1 [expr {$len - 1}]
3581                  lappend const [string trim [string range $line 0 $lenm1]]
3582                  set line [string range $line $len end]
3583                }
3584                lappend conslist $const
3585            }
3586            return $conslist
3587        }
3588        restraintlist-set {
3589            set num [llength $value]
3590            if ![validint num 5] {return 0}
3591            setexp "$key NBNDS" $num 1 5
3592            # delete all old records
3593            foreach i [array names ::exparray "${key}BD*"] {unset ::exparray($i)}
3594            set i 0
3595            foreach cons $value {
3596                incr i
3597                set fi [string toupper [format %.4x $i]]
3598                makeexprec "${key}BD$fi"
3599                set pos 1
3600                foreach num $cons len {3 5 5 3 3 3 3 3 -6 -6} {
3601                    if {$len > 0} {
3602                        validint num $len
3603                        setexp "${key}BD$fi" $num $pos $len
3604                    } else {
3605                        set len [expr abs($len)]
3606                        validreal num $len 3
3607                        setexp "${key}BD$fi" $num $pos $len
3608                    }
3609                    incr pos $len
3610                }
3611            }
3612        }
3613        default {
3614            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
3615            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3616        }
3617    return 1
3618    }
3619}
3620
3621# read/edit chemical restraint info
3622#  parm:
3623#    weight -- histogram weight (factr) *
3624#    restraintlist -- list of restraints *
3625#  action: get (default) or set
3626#  value: used only with set
3627#      value is a list of constraints
3628#      each constrain contains {sum esd cons1 cons2...}
3629#      each consN contains {phase atomnum multiplier}
3630#  * =>  read+write supported
3631# Examples:
3632#
3633#ChemConst restraintlist set { {10 1.1 {1 1 2} {2 2 3}} {0 1 {1 1 1} {1 2 -2}} }
3634#
3635#ChemConst restraintlist get
3636#{10.00000 1.10000 {1 1 2.00000} {2 2 3.00000}} {0.00000 1.00000 {1 1 1.00000} {1 2 -2.00000}}
3637# constraint one 2*(1:1) + 3*(2:2) = 10(1.1)
3638# constraint two 1*(1:1) - 2*(1:2) = 0(1)
3639#   where (1:2) is the total number of atoms (multiplicity*occupancy) for atom 2 in phase 1
3640
3641proc ChemConst {parm "action get" "value {}"} {
3642    set HST {}
3643    # look for CMP record
3644    set n 0
3645    for {set i 0} {$i < $::expmap(nhst)} {incr i} {
3646        set ihist [expr {$i + 1}]
3647        if {[expr {$i % 12}] == 0} {
3648            incr n
3649            set line [readexp " EXPR  HTYP$n"]
3650            if {$line == ""} {
3651                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
3652                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3653            }
3654            set j 0
3655        } else {
3656            incr j
3657        }
3658        if {[string range $line [expr 2+5*$j] [expr 5*($j+1)]] == "CMP "} {
3659            set HST $ihist
3660        }
3661    }
3662    if {$HST <=9} {
3663        set key "HST  $HST"
3664    } else {
3665        set key "HST $HST"
3666    }
3667    if {$HST == "" && $action == "set"} {
3668        # no CMP found need to add the soft constr. histogram
3669        # increment number of histograms
3670        set hst [string trim [string range [readexp { EXPR  NHST }] 0 4]]
3671        incr hst
3672        set HST $hst
3673        if ![validint hst 5] {return 0}
3674        setexp  { EXPR  NHST } $hst 1 5
3675        # add to EXPR HTYPx rec, creating if needed
3676        set n [expr { 1+ (($HST - 1) / 12) }]
3677        set key " EXPR  HTYP$n"
3678        if {[array names ::exparray $key] == ""} {
3679            makeexprec $key
3680        }
3681        setexp $key "CMP " [expr 3 + 5*(($HST-1) % 12)] 5
3682        # create other HST  xx recs
3683        if {$HST <=9} {
3684            set key "HST  $HST"
3685        } else {
3686            set key "HST $HST"
3687        }
3688        makeexprec "$key  HNAM"
3689        setexp "$key  HNAM" "Chemical composition restraints" 3 31
3690        makeexprec "$key FACTR"
3691#       makeexprec "$key NBNDS"
3692        makeexprec "$key NCMPS"
3693        mapexp
3694    } elseif {$HST == ""} {
3695        if $::expgui(debug) {puts "no restraints"}
3696        return "1"
3697    }
3698
3699    switch -glob ${parm}-$action {
3700        weight-get {
3701            return [string trim [string range [readexp "$key FACTR"] 0 14]]
3702        }
3703        weight-set {
3704            # update FACTR
3705            if ![validreal value 15 6] {return 0}
3706            setexp "$key FACTR" $value 1 15
3707        }
3708        restraintlist-get {
3709            set ncons [string trim [string range [readexp "$key NCMPS"] 0 4]]
3710            set conslist {}
3711            for {set i 1} {$i <= $ncons} {incr i} {
3712                set const {}
3713                set line [readexp "${key} CM$i  "]
3714                # number of terms
3715                set nterm [string trim [string range $line 0 4]]
3716                if {$nterm == ""} {set nterm 0}
3717                # chemical sum and esd
3718                lappend const [string trim [string range $line 5 14]]
3719                lappend const [string trim [string range $line 15 24]]
3720                for {set j 1} {$j <= $nterm} {incr j} {
3721                    set n [expr {($j + 2)/3}]
3722                    set o1 [expr {20*(($j-1)%3)}]
3723                    set o2 [expr {19 + 20*(($j-1)%3)}]
3724                    validint n 2
3725                    if {$o1 == 0} {
3726                        set line [readexp "${key} CM${i}${n}"]
3727                    }
3728                    set frag [string range $line $o1 $o2]                   
3729                    lappend const [list \
3730                                       [string trim [string range $frag 0 4]] \
3731                                       [string trim [string range $frag 5 9]] \
3732                                       [string trim [string range $frag 10 19]] \
3733                                   ]
3734                }
3735                lappend conslist $const
3736            }
3737            return $conslist
3738        }
3739        restraintlist-set {
3740            set num [llength $value]
3741            if ![validint num 5] {return 0}
3742            setexp "$key NCMPS" $num 1 5
3743            # delete all old records
3744            foreach i [array names ::exparray "${key} CM*"] {
3745                unset ::exparray($i)
3746            }
3747            set i 0
3748            foreach cons $value {
3749                incr i
3750                set sum [lindex $cons 0]
3751                set esd [lindex $cons 1]
3752                set terms [lrange $cons 2 end]
3753                set nterms [llength $terms]
3754                validint nterms 5
3755                validreal sum 10 5
3756                validreal esd 10 5
3757                makeexprec "${key} CM$i  " 
3758                setexp "${key} CM$i  " "${nterms}${sum}${esd}" 1 25
3759                set j 0
3760                set str {}
3761                foreach term $terms {
3762                    incr j
3763                    set n [expr {($j + 2)/3}]
3764                    if {$n > 99} break
3765                    validint n 2
3766                    foreach {phase atom mult} $term {}
3767                    validint phase 5
3768                    validint atom 5
3769                    validreal mult 10 5
3770                    append str "${phase}${atom}${mult}"
3771                    if {[expr {$j%3}] == 0} {
3772                        #puts [readexp "${key} CM${i}${n}"]
3773                        makeexprec "${key} CM${i}${n}"
3774                        setexp "${key} CM${i}${n}" $str 1 60
3775                        set str {}
3776                    }
3777                }
3778                if {[string length $str] > 0} {
3779                    makeexprec "${key} CM${i}${n}"
3780                    setexp "${key} CM${i}${n}" $str 1 60
3781                }
3782            }
3783        }
3784        default {
3785            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
3786            puts $msg
3787            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
3788        }
3789    return 1
3790    }
3791}
3792
3793#======================================================================
3794# conversion routines
3795#======================================================================
3796
3797# convert x values to d-space
3798proc tod {xlist hst} {
3799    global expmap
3800    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3801        return [toftod $xlist $hst]
3802    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3803        return [tttod $xlist $hst]
3804    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3805        return [engtod $xlist $hst]
3806    } else {
3807        return {}
3808    }
3809}
3810
3811# convert tof to d-space
3812proc toftod {toflist hst} {
3813    set difc [expr {[histinfo $hst difc]/1000.}]
3814    set difc2 [expr {$difc*$difc}]
3815    set difa [expr {[histinfo $hst difa]/1000.}]
3816    set zero [expr {[histinfo $hst zero]/1000.}]
3817    set ans {}
3818    foreach tof $toflist {
3819        if {$tof == 0.} {
3820            lappend ans 0.
3821        } elseif {$tof == 1000.} {
3822            lappend ans 1000.
3823        } else {
3824            set td [expr {$tof-$zero}]
3825            lappend ans [expr {$td*($difc2+$difa*$td)/ \
3826                    ($difc2*$difc+2.0*$difa*$td)}]
3827        }
3828    }
3829    return $ans
3830}
3831
3832# convert two-theta to d-space
3833proc tttod {twotheta hst} {
3834    set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3835    set zero [expr [histinfo $hst zero]/100.]
3836    set ans {}
3837    set cnv [expr {acos(0.)/180.}]
3838    foreach tt $twotheta {
3839        if {$tt == 0.} {
3840            lappend ans 99999.
3841        } elseif {$tt == 1000.} {
3842            lappend ans 0.
3843        } else {
3844            lappend ans [expr {$lamo2 / sin($cnv*($tt-$zero))}]
3845        }
3846    }
3847    return $ans
3848}
3849
3850# convert energy (edx-ray) to d-space
3851# (note that this ignores the zero correction)
3852proc engtod {eng hst} {
3853    set lam [histinfo $hst lam1]
3854    set zero [histinfo $hst zero]
3855    set ans {}
3856    set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3857    foreach e $eng {
3858        if {$e == 0.} {
3859            lappend ans 1000.
3860        } elseif {$e == 1000.} {
3861            lappend ans 0.
3862        } else {
3863            lappend ans [expr {$v/$e}]
3864        }
3865    }
3866    return $ans
3867}
3868
3869# convert x values to Q
3870proc toQ {xlist hst} {
3871    global expmap
3872    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3873        return [toftoQ $xlist $hst]
3874    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3875        return [tttoQ $xlist $hst]
3876    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3877        return [engtoQ $xlist $hst]
3878    } else {
3879        return {}
3880    }
3881}
3882# convert tof to Q
3883proc toftoQ {toflist hst} {
3884    set difc [expr {[histinfo $hst difc]/1000.}]
3885    set difc2 [expr {$difc*$difc}]
3886    set difa [expr {[histinfo $hst difa]/1000.}]
3887    set zero [expr {[histinfo $hst zero]/1000.}]
3888    set 2pi [expr {4.*acos(0.)}]
3889    set ans {}
3890    foreach tof $toflist {
3891        if {$tof == 0.} {
3892            lappend ans 99999.
3893        } elseif {$tof == 1000.} {
3894            lappend ans 0.
3895        } else {
3896            set td [expr {$tof-$zero}]
3897            lappend ans [expr {$2pi * \
3898                    ($difc2*$difc+2.0*$difa*$td)/($td*($difc2+$difa*$td))}]
3899        }
3900    }
3901    return $ans
3902}
3903
3904# convert two-theta to Q
3905proc tttoQ {twotheta hst} {
3906    set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3907    set zero [expr [histinfo $hst zero]/100.]
3908    set ans {}
3909    set cnv [expr {acos(0.)/180.}]
3910    set 2pi [expr {4.*acos(0.)}]
3911    foreach tt $twotheta {
3912        if {$tt == 0.} {
3913            lappend ans 0.
3914        } elseif {$tt == 1000.} {
3915            lappend ans 1000.
3916        } else {
3917            lappend ans [expr {$2pi * sin($cnv*($tt-$zero)) / $lamo2}]
3918        }
3919    }
3920    return $ans
3921}
3922# convert energy (edx-ray) to Q
3923# (note that this ignores the zero correction)
3924proc engtoQ {eng hst} {
3925    set lam [histinfo $hst lam1]
3926    set zero [histinfo $hst zero]
3927    set ans {}
3928    set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3929    set 2pi [expr {4.*acos(0.)}]
3930    foreach e $eng {
3931        if {$e == 0.} {
3932            lappend ans 0.
3933        } elseif {$e == 1000.} {
3934            lappend ans 1000.
3935        } else {
3936            lappend ans [expr {$2pi * $e / $v}]
3937        }
3938    }
3939    return $ans
3940}
3941proc sind {angle} {
3942    return [expr {sin($angle*acos(0.)/90.)}]
3943}
3944
3945# convert d-space values to 2theta, TOF or KeV
3946proc fromd {dlist hst} {
3947    global expmap
3948    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
3949        set difc [expr {[histinfo $hst difc]/1000.}]
3950        set difa [expr {[histinfo $hst difa]/1000.}]
3951        set zero [expr {[histinfo $hst zero]/1000.}]
3952        set ans {}
3953        foreach d $dlist {
3954            if {$d == 0.} {
3955                lappend ans 0.
3956            } elseif {$d == 1000.} {
3957                lappend ans 1000.
3958            } else {
3959                lappend ans [expr {$difc*$d + $difa*$d*$d + $zero}]
3960            }
3961        }
3962        return $ans
3963    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
3964        set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
3965        set zero [expr [histinfo $hst zero]/100.]
3966        set ans {}
3967        set cnv [expr {180./acos(0.)}]
3968        foreach d $dlist {
3969            if {$d == 99999.} {
3970                lappend ans 0
3971            } elseif {$d == 0.} {
3972                lappend ans 1000.
3973            } else {
3974                lappend ans [expr {$cnv*asin($lamo2/$d) + $zero}]
3975            }
3976        }
3977        return $ans
3978    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
3979        set lam [histinfo $hst lam1]
3980        set zero [histinfo $hst zero]
3981        set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
3982        set ans {}
3983        foreach d $dlist {
3984            if {$d == 1000.} {
3985                lappend ans 0
3986            } elseif {$d == 0.} {
3987                lappend ans 1000.
3988            } else {
3989                lappend ans [expr {$v/$d}]
3990            }
3991        }
3992        return $ans
3993    } else {
3994        return {}
3995    }
3996}
3997
3998# convert Q values to 2theta, TOF or KeV
3999proc fromQ {Qlist hst} {
4000    global expmap
4001    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
4002        set difc [expr {[histinfo $hst difc]/1000.}]
4003        set difa [expr {[histinfo $hst difa]/1000.}]
4004        set zero [expr {[histinfo $hst zero]/1000.}]
4005        set ans {}
4006        foreach Q $Qlist {
4007            if {$Q == 0.} {
4008                lappend ans 1000.
4009            } elseif {$Q == 99999.} {
4010                lappend ans 1000.
4011            } else {
4012                set d [expr {4.*acos(0.)/$Q}]
4013                lappend ans [expr {$difc*$d + $difa*$d*$d + $zero}]
4014            }
4015        }
4016        return $ans
4017    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
4018        set lamo4pi [expr {[histinfo $hst lam1]/(8.*acos(0.))}]
4019        set zero [expr [histinfo $hst zero]/100.]
4020        set ans {}
4021        set cnv [expr {180./acos(0.)}]
4022        foreach Q $Qlist {
4023            if {$Q == 0.} {
4024                lappend ans 0
4025            } elseif {$Q == 1000.} {
4026                lappend ans 1000.
4027            } else {
4028                lappend ans [expr {$cnv*asin($Q*$lamo4pi) + $zero}]
4029            }
4030        }
4031        return $ans
4032    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
4033        set lam [histinfo $hst lam1]
4034        set zero [histinfo $hst zero]
4035        set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
4036        set ans {}
4037        set 2pi [expr {4.*acos(0.)}]
4038        foreach Q $Qlist {
4039            if {$Q == 1000.} {
4040                lappend ans 0
4041            } elseif {$Q == 0.} {
4042                lappend ans 1000.
4043            } else {
4044                lappend ans [expr {$Q * $v/$2pi}]
4045            }
4046        }
4047        return $ans
4048    } else {
4049        return {}
4050    }
4051}
4052#============================================================================
4053# rigid body EXP editing routines (to move into readexp.tcl)
4054# RigidBodyList -- returns a list of the defined rigid body types
4055# SetRigidBodyVar -- set variables and damping for rigid body type multipliers
4056# ReadRigidBody  -- # of times a body is mapped, scaling factors, var #s & coordinates
4057# RigidBodyMappingList - return a list instances where a RB is mapped in phase
4058# RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping
4059# RigidBodySetTLS  -- change the TLS values for a rigid body mapping
4060# RigidBodySetDamp -- change the damping values for a rigid body mapping
4061# RigidBodyVary    -- set refinement variable numbers for a rigid body mapping
4062# RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping
4063# AddRigidBody -- defines a new rigid body type
4064# DeleteRigidBody -- remove a rigid body definition
4065# ReplaceRigidBody -- replaces a previous rigid body type
4066# ReadRigidBodyMapping  -- get parameters for a rigid body mapping
4067# MapRigidBody -- map a rigid body type into a phase
4068# EditRigidBodyMapping -- change the parameters in a rigid body mapping
4069# UnMapRigidBody --remove a rigid body constraint by removing a RB "instance"
4070#----- note that these older routines should not be used ------
4071# RigidBodyCount -- returns the number of defined rigid bodies (body types)
4072#    use RigidBodyList instead
4073# RigidBodyMappingCount -- # of times a rigid body is mapped in phase
4074#    use RigidBodyMappingList instead
4075#============================================================================
4076# returns the number of defined rigid bodies
4077proc RigidBodyCount {} {
4078    set n [string trim [readexp "RGBD  NRBDS"]]
4079    if {$n == ""} {
4080        set n 0
4081    }
4082    return $n
4083}
4084
4085# returns a list of the defined rigid body types
4086proc RigidBodyList {} {
4087    set n [string trim [readexp "RGBD  NRBDS"]]
4088    if {$n == ""} {
4089        set n 0
4090    }
4091    set rblist {}
4092    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
4093        set value $rbnum
4094        validint value 2
4095        set key "RGBD${value}"
4096        if {[existsexp "$key NATR "]} {
4097            lappend rblist $rbnum
4098        }
4099        if {[llength $rblist] == $n} break
4100    }
4101    return $rblist
4102}
4103
4104# ReadRigidBody provides all information associated with a rigid body type
4105#  rbnum is the rigid body type number
4106# it returns two items:
4107#   the number of times the rigid body is mapped
4108#   a list containing an element for each scaling factor in rigid body #rbnum.
4109# in each element there are four items:
4110#    a multiplier value for the rigid body coordinates
4111#    a damping value (0-9) for the refinement of the multiplier
4112#    a variable number if the multiplier will be refined
4113#    a list of cartesian coordinates coordinates
4114# each cartesian coordinate contains 4 items: x,y,z and a label
4115#  note that the label is present only when the RB is created in EXPGUI and is
4116#  not used in GSAS.
4117proc ReadRigidBody {rbnum} {
4118    if {[lsearch [RigidBodyList] $rbnum] == -1} {
4119        return ""
4120    }
4121    set value $rbnum
4122    validint value 2
4123    set key "RGBD${value}"
4124    set n [string trim [string range [readexp "$key NATR"] 0 4]]
4125    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
4126    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
4127    set out {}
4128    for {set i 1} {$i <= $nmult} {incr i} {
4129        set line [readexp "${key}${i}PARM"]
4130        set mult [string trim [string range $line 0 9]]
4131        set var [string trim [string range $line 10 14]]
4132        set damp [string trim [string range $line 15 19]]
4133        set coordlist {}
4134        for {set j 1} {$j <= $n} {incr j} {
4135            set value $j
4136            validint value 3
4137            set line [readexp "${key}${i}SC$value"]
4138            set x [string trim [string range $line 0 9]]
4139            set y [string trim [string range $line 10 19]]
4140            set z [string trim [string range $line 20 29]]
4141            set lbl [string trim [string range $line 30 39]]
4142            lappend coordlist [list $x $y $z $lbl]
4143        }
4144        lappend out [list $mult $damp $var $coordlist]
4145    }
4146    return [list $used $out]
4147}
4148
4149# SetRigidBodyVar
4150#   rbnum is the rigid body type number
4151#   varnumlist is a list of variable numbers
4152#      note that if this list is shorter than the number of actual multipliers
4153#      for the body, the unspecified variable will not be changed
4154#   damplist   is a list of damping values (0-9)
4155#      note that if the damplist is shorter than the number of actual multipliers
4156#      the unspecified values are not changed
4157#  SetRigidBodVar 2 {1 2 3} {}
4158#       will vary the (first 3) translations in body #3 and will not change the
4159#       damping values
4160#  SetRigidBodVar 3 {} {0 0 0}
4161#       will not change variable settings but will change the (first 3) damping values
4162#  SetRigidBodVar 4 {11 11} {8 8}
4163#      changes both variable numbers and damping at the same time
4164# Nothing is returned
4165proc SetRigidBodyVar {rbnum varnumlist damplist} {
4166    if {[lsearch [RigidBodyList] $rbnum] == -1} {
4167        return ""
4168    }
4169    set value $rbnum
4170    validint value 2
4171    set key "RGBD${value}"
4172    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
4173    for {set i 1} {$i <= $nmult} {incr i} {
4174        set j $i
4175        incr j -1
4176        set var [lindex $varnumlist $j]
4177        if {$var != ""} {
4178            validint var 5
4179            setexp "${key}${i}PARM" $var 11 15
4180        }
4181        set damp [lindex $damplist $j]
4182        if {$damp != ""} {
4183            if {$damp > 9} {set damp 9}
4184            if {$damp < 0} {set damp 0}
4185            validint damp 5
4186        }
4187        setexp "${key}${i}PARM" $damp 16 20
4188    }
4189}
4190
4191
4192# return the number of times rigid body $bodytyp is mapped in phase $phase
4193proc RigidBodyMappingCount {phase bodytyp} {
4194    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4195    if {! [existsexp "$key  NBDS"]} {return 0}
4196    set n [string trim [readexp "$key  NBDS"]]
4197    if {$n == ""} {
4198        set n 0
4199    }
4200    return $n
4201}
4202# return a list of the instances where rigid body $bodytyp is mapped in phase $phase
4203proc RigidBodyMappingList {phase bodytyp} {
4204    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4205    if {! [existsexp "$key  NBDS"]} {return {}}
4206    set n [string trim [readexp "$key  NBDS"]]
4207    if {$n == ""} {
4208        set n 0
4209    }
4210    set rblist {}
4211    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
4212        set value $rbnum
4213        validint value 2
4214        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
4215        if {[existsexp "$key  NDA"]} {
4216            lappend rblist $rbnum
4217        }
4218        if {[llength $rblist] == $n} break
4219    }
4220    return $rblist
4221}
4222
4223
4224
4225# reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num)
4226# returns a list of items (most lists) as follows:
4227#   1) sequence # of first atom in body
4228#   2) origin of body in fractional coordinates (3 elements)
4229#   3) Euler angles as 6 pairs of numbers (see below)
4230#   4) variable numbers for the 9 position variables (origin followed by rotations)
4231#   5) damping vals for the 9 position variables (origin followed by rotations)
4232#   6) the TLS values, in order below (empty list if TLS is not in use)
4233#   7) the variable numbers for each TLS values, in order below (or empty)
4234#   8) three damping values for the T, L and S terms.
4235# returns an empty list if no such body exists.
4236#
4237# Euler angles are a list of axes and angles to rotate:
4238#   { {axis1 angle1} {axis2 angle2} ...}
4239# where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
4240#
4241# The 20 TLS terms are ordered:
4242#    T11, T22, T33, T12, T13, T23
4243#    L11, L22, L33, L12, L13, L23
4244#    S12, S13, S21, S23, S31, S32, SAA, SBB
4245#
4246proc ReadRigidBodyMapping {phase bodytyp num} {
4247    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4248        return ""
4249    }
4250    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4251    set first [string trim [string range [readexp "$key  NDA"] 0 4]]
4252    set line [readexp "$key BDFL"]
4253    set varlist {}
4254    set damplist {}
4255    foreach i {0 1 2 3 4 5 6 7 8} {
4256        lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]]
4257        lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]]
4258    }
4259    set TLSdamplist {}
4260    foreach i {54 55 56} {
4261        lappend TLSdamplist [string trim [string range $line $i $i ]]
4262    }
4263    set line [readexp "${key} BDLC"]
4264    set x [string trim [string range $line 0 9]]
4265    set y [string trim [string range $line 10 19]]
4266    set z [string trim [string range $line 20 29]]
4267    set origin [list $x $y $z]
4268    set line [readexp "${key} BDOR"]
4269    set rotations {}
4270    foreach i {0 10 20 30 40 50} {
4271        set angle [string trim [string range $line $i [expr {$i+7}]]]
4272        set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]]
4273        lappend rotations [list $angle $axis]
4274    }
4275    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
4276    set tlsvars {}
4277    set tlsvals {}
4278    if {$TLS != 0} {
4279        set line [readexp "${key}TLSF1"]
4280        for {set j 0} {$j < 20} {incr j} {
4281            set var [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]]
4282            if {$var == ""} {set var 0}
4283            lappend tlsvars $var
4284        }
4285        for {set j 0} {$j < 20} {incr j} {
4286            set i 0
4287            if {$j == 0} {
4288                set i 1
4289            } elseif {$j == 8} {
4290                set i 2
4291            } elseif {$j == 16} {
4292                set i 3
4293            }
4294            if {$i != 0} {
4295                set line [readexp "${key}TLSP$i"]
4296                set i 0
4297                set j1 0
4298                set j2 7
4299            } else {
4300                incr j1 8
4301                incr j2 8
4302            }
4303            set val [string trim [string range $line $j1 $j2]]
4304            if {$val == ""} {set val 0}
4305            lappend tlsvals $val
4306        }
4307    }
4308    return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist]
4309}
4310
4311# Control TLS representation for phase, body # and instance number of a Rigid body mapping
4312#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
4313# Enable TLS use if TLS is non-zero (true). Disable if zero
4314proc RigidBodyEnableTLS {phase bodytyp num TLS} {
4315    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4316        return ""
4317    }
4318    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4319    if {$TLS} {
4320        setexp "${key} LSTF" [format "%5d" 1] 1 5
4321        if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
4322        if {![existsexp "${key}TLSP1"]} {
4323            makeexprec "${key}TLSP1"
4324            set str {}
4325            foreach v {.01 .01 .01 0 0 0 0 0} d {4 4 4 4 4 4 2 2} {
4326                validreal v 8 $d
4327                append str $v
4328            }
4329            setexp "${key}TLSP1" $str 1 64
4330        }
4331        if {![existsexp "${key}TLSP2"]} {
4332            makeexprec "${key}TLSP2"
4333            set str {}
4334            set v 0
4335            foreach d {2 2 2 2 4 4 4 4} {
4336                validreal v 8 $d
4337                append str $v
4338            }
4339            setexp "${key}TLSP2" $str 1 64
4340        }
4341        if {![existsexp "${key}TLSP3"]} {
4342            makeexprec "${key}TLSP3"
4343            set str {}
4344            set v 0
4345            foreach d {4 4 4 4} {
4346                validreal v 8 $d
4347                append str $v
4348            }
4349            setexp "${key}TLSP3" $str 1 64
4350        }
4351    } else {
4352        setexp "${key} LSTF" [format "%5d" 0] 1 5
4353    }
4354    return 1
4355}
4356
4357# Control the TLS values for Rigid body mapping for mapping with
4358#    phase ($phase), body type # ($bodytyp) and instance # ($num)
4359# set the 20 TLS values to the values in TLSvals
4360# There must be exactly 20 TLS terms, which are ordered:
4361#    T11, T22, T33, T12, T13, T23
4362#    L11, L22, L33, L12, L13, L23
4363#    S12, S13, S21, S23, S31, S32, SAA, SBB
4364proc RigidBodySetTLS {phase bodytyp num TLSvals} {
4365    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4366        return ""
4367    }
4368    if {[llength $TLSvals] != 20} {return ""}
4369    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4370    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
4371    if {$TLS == 0} {return ""}
4372    if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
4373    foreach n {1 2 3} {
4374        if {![existsexp "${key}TLSP$n"]} {makeexprec "${key}TLSP$n"}
4375    }
4376    set str {}
4377    set n 1
4378    set i 0
4379    foreach v $TLSvals d {4 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4} {
4380        incr i
4381        validreal v 8 $d
4382        append str $v
4383        if {$i == 8} {
4384            set i 0
4385            setexp "${key}TLSP$n" $str 1 64
4386            incr n
4387            set str {}
4388        }
4389    }
4390    setexp "${key}TLSP$n" $str 1 64
4391    return 1
4392}
4393
4394# set damping values for a Rigid body mapping
4395#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
4396# there must be 9 damping values in RBdamp for the 9 position variables (origin followed by rotations)
4397# Use of TLSdamp is optional, but to be used, TLS representation must be enabled and there must be
4398# three damping terms (for all T terms; for all L terms and for all S terms)
4399proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} {
4400    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4401        return ""
4402    }
4403    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4404    if {[llength $RBdamp] != 9} {return ""}
4405    set str {}
4406    foreach v $RBdamp {
4407        if {[validint v 1] != 1} {set v " "}
4408        append str $v
4409    }
4410    setexp "$key BDFL" $str 46 9
4411    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
4412    if {$TLS != 0 &&  [llength $TLSdamp] == 3} {
4413        set str {}
4414        foreach v $TLSdamp {
4415        if {[validint v 1] != 1} {set v " "}
4416            append str $v
4417        }
4418        setexp "$key BDFL" $str 55 3
4419    }
4420    return 1
4421}
4422
4423# set refinement variable numbers for a Rigid body mapping
4424#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
4425# there must be 9 variable values in RBvar for the 9 position variables (origin followed by rotations)
4426# note that the variable values should be unique integers
4427proc RigidBodyVary {phase bodytyp num RBvar} {
4428    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4429        return ""
4430    }
4431    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4432    if {[llength $RBvar] != 9} {return ""}
4433    set str {}
4434    foreach v $RBvar {
4435        if {[validint v 5] != 1} {set v " "}
4436        append str $v
4437    }
4438    setexp "$key BDFL" $str 1 45   
4439}
4440
4441# set TLS refinement variable numbers for a Rigid body mapping
4442#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
4443# there must be 20 variable values in TLSvar for the 20 parameters:
4444#    T11, T22, T33, T12, T13, T23
4445#    L11, L22, L33, L12, L13, L23
4446#    S12, S13, S21, S23, S31, S32, SAA, SBB
4447# note that the variable values should be unique integers
4448proc RigidBodyTLSVary {phase bodytyp num TLSvar} {
4449    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
4450        return ""
4451    }
4452    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
4453    if {[llength $TLSvar] != 20} {return ""}
4454    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
4455    if {$TLS == 0} {return ""}
4456    set str {}
4457    foreach v $TLSvar {
4458        if {[validint v 3] != 1} {set v " "}
4459        append str $v
4460    }
4461    setexp "${key}TLSF1" $str 1 60
4462
4463# AddRigidBody: add a new rigid body definition into the .EXP file
4464# arguments are:
4465#   multlist: defines a list of multipliers for each set of coordinates. In the
4466#             simplest case this will be {1}
4467#   coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } }
4468# note that when the length of multlist > 1 then coordlist must have the same length.
4469# for input where
4470#     multlist = {s1 s2} and
4471#     coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...}
4472#                     {0 0 0} {1 1 0} {2 1 2} ...}
4473#                 }
4474# the cartesian coordinates are defined from the input as
4475#    atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)]
4476#    atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)]
4477#    atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)]
4478#    ...
4479# Returns the number of the rigid body that has been created
4480proc AddRigidBody {multlist coordlist} {
4481    # find the first unused body #
4482    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
4483        set value $rbnum
4484        validint value 2
4485        set key "RGBD${value}"
4486        if {! [existsexp "$key NATR "]} {break}
4487    }
4488    # did we go too far?
4489    if {$rbnum == 16} {return ""}
4490    # increment the RB counter
4491    set n [string trim [readexp "RGBD  NRBDS"]]
4492    if {$n == ""} {
4493        makeexprec "RGBD  NRBDS"
4494        set n 0
4495    }
4496    incr n
4497    validint n 5
4498    setexp "RGBD  NRBDS" $n 1 5
4499    SetRigidBody $rbnum $multlist $coordlist
4500    return $rbnum
4501}
4502
4503# DeleteRigidBody: remove a rigid body definition from the .EXP file
4504# The body may not be mapped. I am not sure if GSAS allows more than 9 bodies,
4505# but if it does, the simplifed approach used here will fail, so this
4506# is not allowed.
4507# Input:
4508#   Rigid body number
4509# Returns:
4510#   1 on success
4511#   -1 if the body number is 11 or greater
4512#   -2 if the body is mapped
4513#   -3 if the body is not defined
4514proc DeleteRigidBody {rbnum} {
4515    # can't delete bodies with numbers higher than 10, since the key prefix
4516    # (RGBD11... will overlap with rigid body instance records, which would be
4517    # deleted below
4518    if {$rbnum > 10} {
4519        return -1
4520    }
4521    set value $rbnum
4522    validint value 2
4523    set key "RGBD${value}"
4524    if {![existsexp "$key NATR "]} {
4525        return -2
4526    }
4527    # make sure the body is not mapped
4528    if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} {
4529        return -3
4530    }
4531    # delete the records starting with "RGBD x" or "RGBD10"
4532    foreach key [array names ::exparray "${key}*"] {
4533        #puts $key
4534        delexp $key
4535    }
4536    # decrement the RB counter
4537    set n [string trim [readexp "RGBD  NRBDS"]]
4538    if {$n == ""} {
4539        set n 0
4540    }
4541    incr n -1
4542    validint n 5
4543    if {$n > 0} {
4544        setexp "RGBD  NRBDS" $n 1 5
4545    } else {
4546        delexp "RGBD  NRBDS"
4547    }
4548    return 1
4549}
4550
4551# ReplaceRigidBody: replace all the information for rigid body #rbnum
4552# Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather
4553# than added.
4554# Note that count of the # of times the body is used is preserved
4555proc ReplaceRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
4556    set value $rbnum
4557    validint value 2
4558    set key "RGBD${value}"
4559    set line [readexp "$key NBDS"]
4560    foreach key [array names ::exparray "${key}*"] {
4561        #puts $key
4562        delexp $key
4563    }
4564    SetRigidBody $rbnum $multlist $coordlist $varlist $damplist
4565    setexp "$key NBDS" $line 1 68
4566}
4567
4568# Edit the parameters for rigid body #rbnum
4569# (normally called from ReplaceRigidBody or AddRigidBody)
4570proc SetRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
4571    set value $rbnum
4572    validint value 2
4573    set key "RGBD${value}"
4574    # number of atoms
4575    set value [llength [lindex $coordlist 0]]
4576    validint value 5
4577    makeexprec "$key NATR"
4578    setexp "$key NATR" $value 1 5
4579    # number of times used
4580    set value 0
4581    validint value 5
4582    makeexprec "$key NBDS"
4583    setexp "$key NBDS" $value 1 5
4584    # number of coordinate matrices
4585    set value [llength $multlist]
4586    validint value 5
4587    makeexprec "$key NSMP"
4588    setexp "$key NSMP" $value 1 5
4589    set i 0
4590    foreach mult $multlist coords $coordlist {
4591        set var [lindex $varlist $i]
4592        if {$var == ""} {set var 0}
4593        set damp [lindex $damplist $i]
4594        if {$damp == ""} {set damp 0}
4595        incr i
4596        makeexprec "${key}${i}PARM"
4597        setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult $var $damp] 1 20
4598        set j 0
4599        foreach item $coords {
4600            #puts $item
4601            incr j
4602            set value $j
4603            validint value 3
4604            makeexprec "${key}${i}SC$value"
4605            if {[llength $item] == 4} {
4606                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
4607            } elseif {[llength $item] == 3} {
4608                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
4609            } else {
4610                return -code 3 "Invalid number of coordinates"
4611            }
4612        }
4613    }
4614}
4615
4616# convert a decimal to the GSAS hex encoding with a field $digits long.
4617proc ToHex {num digits} {
4618    return [string toupper [format "%${digits}x" $num]]
4619}
4620
4621# convert a GSAS hex encoding to a decimal integer
4622proc FromHex {hex} {
4623    return [scan $hex "%x"]
4624}
4625
4626# MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates
4627# (and optionally U values) for a set of atoms will be generated from the rigid body
4628# arguments:
4629#   phase: phase number (1-9)
4630#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4631#   firstatom: sequence number of the first atom in phase (note that atoms may
4632#              not be numbered sequentially)
4633#   position: list of three fractional coordinates for the origin of the rigid body coordinates
4634#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
4635#           cartesian system before the body is translated to position.
4636# returns the instance # (number of times body $bodytyp has been used in phase $phase)
4637proc MapRigidBody {phase bodytyp firstatom position angles} {
4638    # find the first unused body # for this phase & type
4639    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
4640        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
4641        if {! [existsexp "$key  NDA"]} {break}
4642    }
4643    # did we go too far?
4644    if {$rbnum == 16} {return ""}
4645    # increment number of mapped bodies of this type overall
4646    set value $bodytyp
4647    validint value 2
4648    set key "RGBD${value}"
4649    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
4650    incr used
4651    set value $used
4652    validint value 5
4653    setexp "$key NBDS" $value 1 5
4654    # increment number of mapped bodies of this type in this phase
4655    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4656    if {[existsexp "$key  NBDS"]} {
4657        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
4658    } else {
4659        makeexprec "$key  NBDS"
4660        set used 0
4661    }
4662    incr used
4663    set value $used
4664    validint value 5
4665    setexp "$key  NBDS" $value 1 5
4666    # now write the mapping parameters
4667    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
4668    set value $firstatom
4669    validint value 5
4670    makeexprec "$key  NDA"
4671    setexp "$key  NDA" $value 1 5
4672    set l1 {}
4673    set l2 {}
4674    for {set i 0} {$i < 9} {incr i} {
4675        append l1 [format %5d 0]
4676        append l2 [format %1d 0]
4677    }
4678    makeexprec "$key BDFL"
4679    setexp "$key BDFL" $l1$l2 1 54
4680    makeexprec "${key} BDLC"
4681    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
4682    makeexprec "${key} BDOR"
4683    set l1 {}
4684    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
4685        append l1 [format "%8.2f%2d" $val $dir]
4686    }
4687    setexp "${key} BDOR" $l1 1 60
4688    makeexprec "${key} LSTF"
4689    setexp "${key} LSTF" [format "%5d" 0] 1 5
4690    # turn off the X refinement flags for the new body
4691    set st [lsearch $::expmap(atomlist_$phase) $firstatom]
4692    set natoms [llength [lindex [lindex [lindex [ReadRigidBody $bodytyp] 1] 0] 3]]
4693    set en [expr {$st+$natoms-1}]
4694    set atomlist [lrange $::expmap(atomlist_$phase) $st $en]
4695    atominfo $phase $atomlist xref set 0
4696    # redo the mapping to capture the newly mapped atoms
4697    mapexp
4698    return $rbnum
4699}
4700
4701# EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody)
4702# arguments:
4703#   phase: phase number (1-9)
4704#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4705#   bodynum: instance number, as returned by MapRigidBody
4706#   position: list of three fractional coordinates for the origin of the rigid body coordinates
4707#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
4708#           cartesian system before the body is translated to position.
4709#
4710proc EditRigidBodyMapping {phase bodytyp bodynum position angles} {
4711    # number of bodies of this type in this phase
4712    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
4713    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
4714    set l1 {}
4715    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
4716        append l1 [format "%8.2f%2d" $val $dir]
4717    }
4718    setexp "${key} BDOR" $l1 1 60
4719}
4720
4721# UnMapRigidBody: remove a rigid body constraint by removing a RB "instance"
4722# (undoes MapRigidBody)
4723# arguments:
4724#   phase: phase number (1-9)
4725#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
4726#   bodynum: instance number, as returned by MapRigidBody
4727proc UnMapRigidBody {phase bodytyp mapnum} {
4728    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} {
4729        return ""
4730    }
4731    # decrement number of mapped bodies of this type overall
4732    set value $bodytyp
4733    validint value 2
4734    set key "RGBD${value}"
4735    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
4736    incr used -1
4737    set value $used
4738    validint value 5
4739    setexp "$key NBDS" $value 1 5
4740    # decrement number of mapped bodies of this type in this phase
4741    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
4742    if {[existsexp "$key  NBDS"]} {
4743        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
4744    } else {
4745        set used 0
4746    }
4747    incr used -1
4748    set value $used
4749    validint value 5
4750    if {$used > 0} {
4751        setexp "$key  NBDS" $value 1 5
4752    } else {
4753        delexp "$key  NBDS"
4754    }
4755    # now delete the mapping parameter records
4756    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]"
4757    foreach key [array names ::exparray "${key}*"] {
4758        delexp $key
4759    }
4760    return $used
4761}
4762
4763# return a list of defined Fourier maps
4764proc listFourier {} {
4765    set l {}
4766    foreach i {1 2 3 4 5 6 7 8 9} {
4767        if {[existsexp "  FOUR CDAT$i"]} {
4768            lappend l $i
4769        }
4770    }
4771    return $l
4772}
4773
4774# read a Fourier map entry
4775# returns five values:
4776#   0: type of map (DELF,FCLC,FOBS,*FDF,PTSN,DPTS)
4777#   1: section (X,Y or Z)
4778#   2: phase (1-9)
4779#   3: DMIN (usually 0.0)
4780#   4: DMAX (usually 999.99)
4781proc readFourier {num} {
4782    set key "  FOUR CDAT$num"
4783    if {![existsexp $key]} {
4784        return {}
4785    }
4786    set vals {}
4787    # 0: type of map (DELF,FCLC,FOBS,[2-9]FDF,PTSN,DPTS)
4788    lappend vals [string trim [string range [readexp $key] 2 6]]
4789    # 1: section (X,Y or Z)
4790    lappend vals [string trim [string range [readexp $key] 7 8]]
4791    # 2: phase (1-9)
4792    lappend vals [string trim [string range [readexp $key] 8 13]]
4793    # 3: DMIN (usually 0.0)
4794    lappend vals [string trim [string range [readexp $key] 18 25]]
4795    # 4: DMAX (usually 999.99)
4796    lappend vals [string trim [string range [readexp $key] 30 37]]
4797    return $vals
4798}
4799
4800# add a new Fourier map computation type
4801#   arguments:
4802#      phase: (1-9)
4803#      type: type of map (DELF,FCLC,FOBS,*FDF,PTSN,DPTS) - default DELF
4804#      section: (X,Y or Z) - default Z
4805#   returns the number of the map that is added
4806proc addFourier {phase {type "DELF"} {section "Z"}} {
4807    set num {}
4808    foreach i {1 2 3 4 5 6 7 8 9} {
4809        set key "  FOUR CDAT$i"
4810        if {! [existsexp "  FOUR CDAT$i"]} {
4811            set num $i
4812            break
4813        }
4814    }
4815    if {$num == ""} {return {}}
4816    set key "  FOUR CDAT$num"
4817    makeexprec $key
4818    setexp $key $type 3 4
4819    setexp $key $section 8 1
4820    validint phase 5
4821    setexp $key $phase 9 5
4822    setexp $key "NOPR   0.00      999.99" 15 23
4823    return $num
4824}
4825
4826# delete all Fourier map computations
4827proc delFourier {} {
4828    foreach i {1 2 3 4 5 6 7 8 9} {
4829        set key "  FOUR CDAT$i"
4830        delexp $key
4831    }
4832}
4833
4834# read/set a Fourier computation value
4835# use: Fourierinfo num parm
4836#  or: Fourierinfo num parm set value
4837#
4838#  num is the Fourier entry
4839#  parm is one of the following
4840#     type    -- type of map (DELF,FCLC,FOBS,*FDF,PTSN,DPTS)
4841#     section -- last running map direction (X,Y or Z)
4842#     phase   -- phase (1-9)
4843#     dmin    -- d-space for highest order reflection to use (usually 0.0)
4844#     dmax    -- d-space for lowest order reflection to use (usually 999.99)
4845# all parameters may be read or set
4846proc Fourierinfo {num parm "action get" "value {}"} {
4847    set key "  FOUR CDAT$num"
4848    if {![existsexp $key]} {
4849        return {}
4850    }
4851    switch -glob ${parm}-$action {
4852        type-get {
4853            # type of map (DELF,FCLC,FOBS,*FDF,PTSN,DPTS)
4854            return [string trim [string range [readexp $key] 2 6]]
4855        }
4856        type-set {
4857            set found 0
4858            foreach val {DELF FCLC FOBS 2FDF 3FDF 4FDF 5FDF 6FDF 7FDF 8FDF 9FDF PTSN DPTS} {
4859                if {$val == $value} {
4860                    set found 1
4861                    break
4862                }
4863            }
4864            if $found {
4865                setexp $key $value 3 4
4866            }
4867        }
4868        section-get {
4869            # section (X,Y or Z)
4870            return [string range [readexp $key] 7 8]
4871        }
4872        section-set {
4873            set found 0
4874            foreach val {X Y Z} {
4875                if {$val == $value} {
4876                    set found 1
4877                    break
4878                }
4879            }
4880            if $found {
4881                setexp $key $value 8 1
4882            }
4883        }
4884        phase-get {
4885            # phase (1-9)
4886            return [string trim [string range [readexp $key] 8 13]]
4887        }
4888        phase-set {
4889            validint value 5
4890            setexp $key $value 9 5
4891        }
4892        dmin-get {
4893            # DMIN (usually 0.0)
4894            return [string trim [string range [readexp $key] 18 25]]
4895        }
4896        dmin-set {
4897            validreal value 7 2
4898            setexp $key $value 19 7
4899        }
4900        dmax-get {
4901            # DMAX (usually 999.99)
4902            return [string trim [string range [readexp $key] 30 37]]
4903        }
4904        dmax-set {
4905            validreal value 7 2
4906            setexp $key $value 31 7
4907        }
4908        default {
4909            set msg "Unsupported Fourierinfo access: parm=$parm action=$action"
4910            puts $msg
4911            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
4912        }
4913    }
4914}
4915
4916# set histograms used in Fourier computation
4917#  use:
4918#     FourierHists $phase
4919#     FourierHists $phase set {4 3 2 1}
4920# returns a list of histograms to be used to compute that phase's Fourier map
4921# or sets a list of histograms to be used to compute that phase's Fourier map
4922#
4923# Note that the histograms are loaded in the order specified with reflections in
4924# the last histogram overwriting those in earlier ones, where a reflection
4925# occurs in more than one place
4926proc FourierHists {phase "action get" "value {}"} {
4927    # note that in theory one can have more than one CRSm  FMHSTn record
4928    # if more than 22 histograms are used but we will ignore this
4929    set key "CRS$phase  FMHST1"
4930    if {![existsexp $key]} {
4931        makeexprec $key
4932    }
4933    if {$action == "get"} {
4934        return [string trim [readexp $key]]
4935    } else {
4936        set hlist {}
4937        foreach hist $value {
4938            validint hist 3
4939            append hlist $hist
4940        }
4941        setexp $key $hlist 0 67
4942    }
4943}
4944# get the Fourier map computation step and limits
4945# returns 4 lists:
4946#   {stepx stepy stepz} : step size in Angstroms
4947#   {xmin xmax} : min and max x in fractional coordinates
4948#   {ymin ymax} : min and max y in fractional coordinates
4949#   {zmin zmax} : min and max z in fractional coordinates
4950proc getFourierLimits {phase} {
4951    set key "CRS$phase  FMPCTL"
4952    if {![existsexp $key]} {
4953        setFourierLimits $phase
4954    }
4955    set i 0
4956    set line [readexp $key]
4957    foreach v {x y z} cell {a b c} {
4958        set cell_$v [phaseinfo $phase $cell]
4959    }
4960    foreach typ {step min max} {
4961        foreach v {x y z} {
4962            set val [string trim [string range $line $i [expr $i+5]]]
4963            if {$val == ""} {set val 0}
4964            set ${typ}_${v} $val
4965            incr i 5
4966        }           
4967    }
4968    set steps {}
4969    if {[catch {
4970        foreach v {x y z} {
4971            set range_$v {}
4972            lappend steps [expr {[set cell_$v] / [set step_$v]}]
4973            lappend range_$v [expr {[set min_$v] * 1. / [set step_$v] }]
4974            lappend range_$v [expr {[set max_$v] * 1. / [set step_$v] }]
4975        }
4976    }]} {
4977        return [list {.2 .2 .2} {0 1} {0 1} {0 1}]
4978    } else {
4979        return [list $steps $range_x $range_y $range_z]
4980    }
4981}
4982
4983# set the Fourier map computation step and limits
4984#   Asteps contains {stepx stepy stepz} : step size in Angstroms
4985#   range_x contains {xmin xmax} : min and max x in fractional coordinates
4986#   range_y contains {ymin ymax} : min and max y in fractional coordinates
4987#   range_z contains {zmin zmax} : min and max z in fractional coordinates
4988proc setFourierLimits {phase \
4989                           {Asteps {.2 .2 .2}} \
4990                           {range_x {0 1}} \
4991                           {range_y {0 1}} \
4992                           {range_z {0 1}} } {
4993    set key "CRS$phase  FMPCTL"
4994    if {![existsexp $key]} {
4995        makeexprec $key
4996    }
4997    set i 1
4998    # steps across map
4999    foreach v {x y z} cell {a b c} As $Asteps {
5000        set s [expr {1 + int([phaseinfo $phase $cell] / $As)}]
5001        set s [expr {$s + ($s % 2)}]
5002        set step_$v $s
5003        lappend steps [set step_$v]
5004        validint s 5
5005        setexp $key $s $i 5
5006        incr i 5
5007    }
5008    # x,y,z min in steps
5009    foreach v {x y z} {
5010        foreach {min max} [set range_$v] {}
5011        set s [expr {int($min * [set step_$v]-.5)}]
5012        validint s 5
5013        setexp $key $s $i 5
5014        incr i 5
5015    }
5016    # x,y,z max in steps
5017    foreach v {x y z} {
5018        foreach {min max} [set range_$v] {}
5019        set s [expr {int($max * [set step_$v]+.5)}]
5020        validint s 5
5021        setexp $key $s $i 5
5022        incr i 5
5023    }
5024}
Note: See TracBrowser for help on using the repository browser.