Changeset 537


Ignore:
Timestamp:
Dec 4, 2009 5:07:50 PM (13 years ago)
Author:
toby
Message:

# on 2002/01/22 22:16:51, toby did:
Major revision:

export to GRACE
BKGEDIT use linear least-squares & fits bkg functions 1-6

add label by phase name
edit phase name in "edit tickmarks"
make obs & calc pattern less obscuring in bkgedit

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/liveplot

    • Property rcs:date changed from 2001/12/16 19:38:52 to 2002/01/22 22:16:51
    • Property rcs:lines changed from +3 -1 to +616 -183
    • Property rcs:rev changed from 1.25 to 1.26
    r498 r537  
    4848set graph(color_chi2) magenta
    4949set graph(color_bkg) green
    50 set graph(color_calc) red
    5150set graph(color_obs) black
    5251set graph(color_input) magenta
    … …  
    6463set expgui(website) www.ncnr.nist.gov/xtal/software/expgui
    6564set peakinfo(obssym) scross
    66 set peakinfo(obssize) 1.0
     65if {$program == "bkgedit"}  {
     66    set peakinfo(obssize) 0.15
     67    set graph(color_calc) pink
     68} else {
     69    set peakinfo(obssize) 1.0
     70    set graph(color_calc) red
     71}
    6772set peakinfo(inpsym) triangle
    6873set peakinfo(inpsize) 1.0
    … …  
    7378    set peakinfo(min$i) -Inf
    7479    set peakinfo(dashes$i) 1
    75 }
     80    set graph(label$i) Phase$i
     81}
     82set expgui(RadiiList) {}
    7683
    7784proc waitmsg {message} {
    … …  
    593600            catch {$box element create phase$i}
    594601            catch {
    595                 $box element config phase$i -color $peakinfo(color$i)
     602                $box element config phase$i -color $peakinfo(color$i) \
     603                        -label $graph(label$i)
    596604            }
    597605        } else {
    … …  
    658666    pack [button $bx.c$i.1 -text "Color\nmenu" \
    659667            -command "setcolor $i"] -side left
     668
     669    pack [frame $bx.l$i -bd 2 -relief groove] -side top
     670   
     671    pack [label $bx.l$i.1 -text " Phase label:"] -side left
     672   
     673    pack [entry $bx.l$i.2 -textvariable graph(label$i) -width 20] \
     674            -side left
     675
    660676    pack [frame $bx.b] -side top
    661     #pack [button $bx.b.1 -command plotdata -text "Update Plot"] \
    662             #    -side left
    663     pack [button $bx.b.4 -command "destroy $bx" -text Close ] -side right
     677    pack [button $bx.b.4 -command "destroy $bx; plotdata" \
     678            -text Close ] -side right
    664679}
    665680
    … …  
    746761    }
    747762    setprintopt $box
     763}
     764
     765#-------------------------------------------------------------------------
     766# export current plot to Grace
     767#-------------------------------------------------------------------------
     768if {$tcl_platform(platform) == "unix"} {
     769    set graph(GraceFile) /tmp/grace_out.agr
     770} else {
     771    set graph(GraceFile) C:/graceout.agr
     772}
     773proc exportgrace {} {
     774    global graph box
     775    global tcl_platform graph
     776    catch {toplevel .export}
     777    raise .export
     778    eval destroy [grid slaves .export]
     779    set col 5
     780    grid [label .export.1a -text Title:] -column 1 -row 1
     781    set graph(title) [$box cget -title]
     782    grid [entry .export.1b -width 60 -textvariable graph(title)] \
     783            -column 2 -row 1 -columnspan 4
     784    grid [label .export.2a -text Subtitle:] -column 1 -row 2
     785    grid [entry .export.2b -width 60 -textvariable graph(subtitle)] \
     786            -column 2 -row 2 -columnspan 4
     787    grid [label .export.3a -text "File name:"] -column 1 -row 3
     788    grid [entry .export.3b -width 60 -textvariable graph(GraceFile)] \
     789            -column 2 -row 3 -columnspan 4
     790    grid [button .export.help -text Help -bg yellow \
     791            -command "MakeWWWHelp liveplot.html grace"] \
     792            -column [incr col -1] -row 4
     793    grid [button .export.c -text "Close" \
     794            -command "set graph(export) 0; destroy .export"] \
     795            -column [incr col -1] -row 4
     796    if {$tcl_platform(platform) == "unix" && [auto_execok xmgrace] != ""} {
     797        grid [button .export.d -text "Export & \nstart grace" \
     798            -command "set graph(export) 1; destroy .export"] \
     799                -column [incr col -1] -row 4
     800    }
     801    grid [button .export.e -text "Export" \
     802            -command "set graph(export) 2; destroy .export"] \
     803            -column [incr col -1] -row 4
     804    tkwait window .export
     805    if {$graph(export) == 0} return
     806    if {[catch {
     807        set fp [open $graph(GraceFile) w]
     808        puts $fp [output_grace $box $graph(title) $graph(subtitle)]
     809        close $fp
     810    } errmsg]} {
     811        MyMessageBox -parent . -title "Export Error" \
     812                -message "An error occured during the export: $errmsg" \
     813                -icon error -type Ignore -default ignore
     814        return
     815    }
     816
     817    if {$graph(export) == 1} {
     818        set err [catch {exec xmgrace $graph(GraceFile) &} errmsg]
     819        if $err {
     820        MyMessageBox -parent . -title "Grace Error" \
     821                -message "An error occured launching grace (xmgrace): $errmsg" \
     822                -icon error -type Ignore -default ignore
     823        }
     824    } else {
     825        MyMessageBox -parent . -title "OK" \
     826                -message "File $graph(GraceFile) created" \
     827                -type OK -default ok
     828    }
    748829}
    749830
    … …  
    9511032}
    9521033
    953 # Chebyschev fit, evaluate the LS vector, V_n = sum_i (T_n(X_i) * Y_i/sigma_i**2)
    954 proc ChebMakeV {xl yl o xmin xmax} {
    955     for {set i 0} {$i < $o} {incr i} {
    956         set sum($i) 0.
    957     }
    958     foreach y $yl x $xl {
    959         # rescale x
    960         set xs [expr {-1 + 2 * (1.*$x - $xmin) / (1.*$xmax - 1.*$xmin)}]
    961         # compute the Chebyschev term Tn(xs)
    962         set Tpp 0
    963         set Tp 0
    964         for {set i 0} {$i < $o} {incr i} {
    965             if {$Tpp == $Tp && $Tp == 0} {
    966                 set T 1
    967             } elseif {$Tpp == 0} {
    968                 set T $xs
    969             } else {
    970                 set T [expr {2. * $xs * $Tp - $Tpp}]
    971             }
    972             set sum($i) [expr {$sum($i) + $y * $T}]
    973 # weighted
    974             # set sum($i) [expr {$sum($i) + $y * $T / ($sigma*$sigma)}]
    975             set Tpp $Tp
    976             set Tp $T
    977         }
    978     }
    979     lappend vec 2 $o 0
    980     for {set i 0} {$i < $o} {incr i} {
    981         lappend vec $sum($i)
    982     }
    983     return $vec
    984 }
    985 
    986 # Chebyschev fit, evaluate the LS Hessian, A
    987 # where A_jk = sum_i {T_j(X_i) * T_k(X_i)/(sigma_i**2)}
    988 proc ChebMakeA {xl o xmin xmax} {
    989     for {set j 0} {$j < $o} {incr j} {
    990         for {set i 0} {$i <= $j} {incr i} {
    991             set sum(${i}_$j) 0.
    992         }
    993     }
    994     foreach x $xl {
    995         # rescale x
    996         set xs [expr {-1 + 2 * (1.*$x - $xmin) / (1.*$xmax - 1.*$xmin)}]
    997         # compute the Chebyschev term Tj(xs)
    998         set Tjpp 0
    999         set Tjp 0
    1000         for {set j 0} {$j < $o} {incr j} {
    1001             if {$Tjpp == $Tjp && $Tjp == 0} {
    1002                 set Tj 1
    1003             } elseif {$Tjpp == 0} {
    1004                 set Tj $xs
    1005             } else {
    1006                 set Tj [expr {2. * $xs * $Tjp - $Tjpp}]
    1007             }
    1008             set Tjpp $Tjp
    1009             set Tjp $Tj
    1010             # compute the Chebyschev term Ti(xs)
    1011             set Tipp 0
    1012             set Tip 0
    1013             for {set i 0} {$i <= $j} {incr i} {
    1014                 if {$Tipp == $Tip && $Tip == 0} {
    1015                     set Ti 1
    1016                 } elseif {$Tipp == 0} {
    1017                     set Ti $xs
    1018                 } else {
    1019                     set Ti [expr {2. * $xs * $Tip - $Tipp}]
    1020                 }
    1021                 set Tipp $Tip
    1022                 set Tip $Ti
    1023                 set sum(${i}_$j) [expr {$sum(${i}_$j) + $Ti * $Tj}]
    1024 # weighted
    1025                 # set sum(${i}_$j) [expr {$sum(${i}_$j) + $Ti * $Tj / ($sigma * $sigma)}]
    1026             }
    1027         }
    1028     }
    1029     lappend mat 2 $o $o
    1030     for {set i 0} {$i < $o} {incr i} {
    1031         for {set j 0} {$j < $o} {incr j} {
    1032             if {$j < $i} {
    1033                 lappend mat $sum(${j}_$i)
    1034             } else {
    1035                 lappend mat $sum(${i}_$j)
    1036             }
    1037         }
    1038     }
    1039     return $mat
    1040 }
    1041 
    10421034# change the binding of the mouse, based on the selected mode
    10431035proc bkgEditMode {b} {
    … …  
    10861078# plot the background points
    10871079proc bkgPointPlot {} {
    1088     global bkglist termmenu chebterms expnam hst tmin tmax
     1080    global bkglist termmenu expgui expnam hst tmin tmax
    10891081    set l {}
    10901082    set fp [open $expnam.bkg$hst w]
    … …  
    11071099        for {set i 2} {$i <= $l/1.5} {incr i 2} {
    11081100            $termmenu insert end radiobutton -label $i \
    1109                     -variable chebterms
     1101                    -variable expgui(FitOrder) -command "BkgFillTermBoxes nosave"
    11101102            set imax $i
    11111103        }
    1112         if {$imax < $chebterms} {set chebterms $imax}
     1104        if {$imax < $expgui(FitOrder)} {set expgui(FitOrder) $imax}
    11131105    } else {
    11141106        .bkg.f.fit1 config -state disabled
    11151107        .bkg.f.terms config -state disabled
    1116         set chebterms 2
     1108        set expgui(FitOrder) 2
    11171109    }
    11181110}
    … …  
    11641156# initialize the background plot
    11651157proc bkghstInit {} {
    1166     global bkglist tmin tmax hst expnam cheblist chebterms
     1158    global bkglist tmin tmax hst expnam termlist expgui
    11671159    set tmin [histinfo $hst tmin]
    11681160    set tmax [histinfo $hst tmax]
    … …  
    11921184    bkgPointPlot
    11931185    bkgFillPoints
    1194     set cheblist ""
    1195     BkgFillCheb
    1196     set chebterms 2
    1197 }
    1198 
    1199 # fit a Chebyshev polynomial to the selected background points
     1186    set termlist ""
     1187    set expgui(FitOrder) 2
     1188    BkgFillTermBoxes nosave
     1189}
     1190
    12001191proc bkgFit {button} {
    1201     global bkglist chebterms cheblist
     1192    global bkglist termlist expgui
    12021193    # keep the button down while working
    12031194    $button config -relief sunken
    … …  
    12121203
    12131204    # perform the Fit
    1214     global tmin tmax
    1215     set V [ChebMakeV $X $Y $chebterms $tmin $tmax]
    1216     #La::show $V
    1217     set A [ChebMakeA $X $chebterms $tmin $tmax]
    1218     #La::show $A
    1219     set cheblist [lrange [La::msolve $A $V] 3 end]
    1220     BkgFillCheb
    1221     bkgFillPoints
    1222     # compute the curve and display it
    1223     set calcb {}
    1224     foreach x [xvec range 0 end] {
    1225         lappend calcb [chebeval $cheblist $x $tmin $tmax]
    1226     }
    1227     .g element configure 11 -xdata xvec -ydata $calcb
    1228     update
     1205    set termlist [FitBkgFunc $X $Y $expgui(FitOrder) $expgui(FitFunction) \
     1206            $expgui(RadiiList)]
     1207    # set the bkg terms in the edit boxes & update the plot
     1208    BkgFillTermBoxes
    12291209    $button config -relief raised
    12301210}
    12311211
    1232 # put the Chebyshev coefficients into edit widgets
    1233 proc BkgFillCheb {} {
    1234     global cheblist
    1235     global chebedit
     1212# put the Background coefficients into edit widgets
     1213proc BkgFillTermBoxes {"save {}"} {
     1214    global termlist expgui
     1215    global bkgeditbox
    12361216    catch {destroy .bkg.canvas.fr}
    12371217    set top [frame .bkg.canvas.fr]
    12381218    .bkg.canvas create window 0 0 -anchor nw -window $top
    1239     # delete trace on chebedit
    1240     foreach v [ trace vinfo chebedit] {
    1241         eval trace vdelete chebedit $v
    1242     }
    1243     if {[llength $cheblist] == 0} {
    1244         grid [label $top.0 -text "(no terms defined)"] -col 1 -row 1
    1245         .bkg.cw config -state disabled
     1219    # delete trace on bkgeditbox
     1220    foreach v [ trace vinfo bkgeditbox] {
     1221        eval trace vdelete bkgeditbox $v
     1222    }
     1223
     1224    .bkg.cw config -state normal
     1225    set k 0
     1226    if {$expgui(FitFunction) == 3} {
     1227        # o is number of refinable terms
     1228        set o [expr {2 + ($expgui(FitOrder) - 2)/2}]
     1229        grid [label $top.lbl -text terms] -col $k -row 1
     1230        if {$expgui(FitOrder) >= 4} {
     1231            grid [label $top.rlbl -text radii] -col $k -row 2
     1232        }
     1233        incr k
     1234        set width 7
    12461235    } else {
    1247         set i -1
    1248         .bkg.cw config -state normal
    1249         foreach c $cheblist {
    1250             incr i
    1251             grid [frame $top.$i -relief groove -bd 3] -col $i -row 1
    1252             grid [label $top.$i.l -text "[expr 1+$i]"] -col 1 -row 1
    1253             grid [entry $top.$i.e -textvariable chebedit($i) -width 13] \
    1254                     -col 2 -row 1
    1255             set chebedit($i) $c
    1256         }
    1257         trace variable chebedit w "BkgRecalcCheb $top"
    1258     }
     1236        set o $expgui(FitOrder)
     1237        set width 10
     1238    }
     1239    for {set i 0} {$i < $o} {incr i} {
     1240        if {$i >= [llength $termlist]} {lappend termlist 0.}
     1241        set bkgeditbox($i) [lindex $termlist $i]
     1242        grid [frame $top.$i -relief groove -bd 3] -col $k -row 1
     1243        grid [label $top.$i.l -text "[expr 1+$i]"] -col 1 -row 1
     1244        grid [entry $top.$i.e -textvariable bkgeditbox($i) -width $width] \
     1245                -col 2 -row 1
     1246        if {$expgui(FitFunction) == 3 && $i > 1} {
     1247            set j [expr $i-2]
     1248            if {$j >= [llength $expgui(RadiiList)]} {lappend expgui(RadiiList) 0.}
     1249            set bkgeditbox(r$j) [lindex $expgui(RadiiList) $j]
     1250            if {$bkgeditbox(r$j) == 0} {
     1251                set bkgeditbox(r$j) ??
     1252            }
     1253            grid [frame $top.r$j -relief groove -bd 3] \
     1254                    -col [expr $k-2] -row 2
     1255            grid [label $top.r$j.l -text "[expr -1+$i]"] -col 1 -row 1
     1256            grid [entry $top.r$j.e -textvariable bkgeditbox(r$j) -width $width] \
     1257                    -col 2 -row 1       
     1258        }
     1259        incr k
     1260    }
     1261    trace variable bkgeditbox w "BkgRecalcPlot $top"
     1262    BkgRecalcPlot $top x x x
    12591263    update idletasks
    12601264    set sizes [grid bbox $top]
    12611265    .bkg.canvas config -scrollregion $sizes -height [lindex $sizes 3]
    1262 }
    1263 
    1264 # respond to edits made to Chebyshev terms
    1265 proc BkgRecalcCheb {top var i mode} {
    1266     global chebedit cheblist
    1267     if [catch {expr $chebedit($i)}] {
    1268         $top.$i.e config -fg red
     1266    # inhibit the save button, if requested
     1267    if {$save == "nosave"} {
     1268        .bkg.cw config -state disabled
     1269        .g element configure 11 -xdata {} -ydata {}
     1270        update
     1271    }
     1272}
     1273
     1274# respond to edits made to background terms
     1275proc BkgRecalcPlot {top var i mode} {
     1276    global bkgeditbox termlist expgui expgui(FitOrder)
     1277
     1278    set good 1
     1279
     1280    if {$expgui(FitFunction) == 3} {
     1281        set expgui(RadiiList) {}
     1282        for {set j 0} {$j < ($expgui(FitOrder) - 2)/2} {incr j} {
     1283            lappend expgui(RadiiList) $bkgeditbox(r$j)
     1284            if {[catch {expr $bkgeditbox(r$j)}]} {
     1285                $top.r$j.e config -fg red
     1286                set good 0
     1287            } elseif {$bkgeditbox(r$j) == 0} {
     1288                $top.r$j.e config -fg red
     1289                set good 0
     1290            } else {
     1291                $top.r$j.e config -fg black
     1292            }
     1293        }
     1294        set o [expr {2 + ($expgui(FitOrder) - 2)/2}]
    12691295    } else {
    1270         $top.$i.e config -fg black
    1271         set cheblist [lreplace $cheblist $i $i $chebedit($i)]
    1272         global tmin tmax
     1296        set o $expgui(FitOrder)
     1297    }
     1298    set termlist {}
     1299    for {set j 0} {$j < $o} {incr j} {
     1300        lappend termlist $bkgeditbox($j)
     1301        if {[catch {expr $bkgeditbox($j)}]} {
     1302            $top.$j.e config -fg red
     1303            set good 0
     1304        } else {
     1305            $top.$j.e config -fg black
     1306        }
     1307    }
     1308
     1309    # disable fit for invalid values
     1310    if {$good} {
     1311        .bkg.cw config -state normal
     1312        .bkg.f.fit1 config -state normal
    12731313        # plot it
    1274         set calcb {}
    1275         foreach x [xvec range 0 end] {
    1276             lappend calcb [chebeval $cheblist $x $tmin $tmax]
    1277         }
     1314        set calcb [BkgEval $termlist $expgui(FitFunction) \
     1315                [xvec range 0 end] $expgui(RadiiList)]
    12781316        .g element configure 11 -xdata xvec -ydata $calcb
     1317        update
     1318    } else {
     1319        .bkg.cw config -state disabled
     1320        .bkg.f.fit1 config -state disabled
     1321        .g element configure 11 -xdata {} -ydata {}
    12791322        update
    12801323    }
    … …  
    13371380}
    13381381
     1382# convert x values to Q
     1383proc toQ {xlist hst} {
     1384    global expmap
     1385    if {[string range $expmap(htype_$hst) 2 2] == "T"} {
     1386        return [toftoQ $xlist $hst]
     1387    } elseif {[string range $expmap(htype_$hst) 2 2] == "C"} {
     1388        return [tttoQ $xlist $hst]
     1389    } elseif {[string range $expmap(htype_$hst) 2 2] == "E"} {
     1390        return [engtoQ $xlist $hst]
     1391    } else {
     1392        return {}
     1393    }
     1394}
     1395# convert tof to Q
     1396proc toftoQ {toflist hst} {
     1397    set difc [expr {[histinfo $hst difc]/1000.}]
     1398    set difc2 [expr {$difc*$difc}]
     1399    set difa [expr {[histinfo $hst difa]/1000.}]
     1400    set zero [expr {[histinfo $hst zero]/1000.}]
     1401    set 2pi [expr {4.*acos(0.)}]
     1402    set ans {}
     1403    foreach tof $toflist {
     1404        if {$tof == 0.} {
     1405            lappend ans 99999.
     1406        } elseif {$tof == 1000.} {
     1407            lappend ans 0.
     1408        } else {
     1409            set td [expr {$tof-$zero}]
     1410            lappend ans [expr {$2pi * \
     1411                    ($difc2*$difc+2.0*$difa*$td)/($td*($difc2+$difa*$td))}]
     1412        }
     1413    }
     1414    return $ans
     1415}
     1416
     1417# convert two-theta to Q
     1418proc tttoQ {twotheta hst} {
     1419    set lamo2 [expr {0.5 * [histinfo $hst lam1]}]
     1420    set zero [expr [histinfo $hst zero]/100.]
     1421    set ans {}
     1422    set cnv [expr {acos(0.)/180.}]
     1423    set 2pi [expr {4.*acos(0.)}]
     1424    foreach tt $twotheta {
     1425        if {$tt == 0.} {
     1426            lappend ans 0.
     1427        } elseif {$tt == 1000.} {
     1428            lappend ans 1000.
     1429        } else {
     1430            lappend ans [expr {$2pi * sin($cnv*($tt-$zero)) / $lamo2}]
     1431        }
     1432    }
     1433    return $ans
     1434}
     1435# convert energy (edx-ray) to Q
     1436# (note that this ignores the zero correction)
     1437proc engtoQ {eng hst} {
     1438    set lam [histinfo $hst lam1]
     1439    set zero [histinfo $hst zero]
     1440    set ans {}
     1441    set v [expr {12.398/(2.0*[sind[expr ($lam/2.0)]])}]
     1442    set 2pi [expr {4.*acos(0.)}]
     1443    foreach e $eng {
     1444        if {$e == 0.} {
     1445            lappend ans 0.
     1446        } elseif {$e == 1000.} {
     1447            lappend ans 1000.
     1448        } else {
     1449            lappend ans [expr {$2pi * $e / $v}]
     1450        }
     1451    }
     1452    return $ans
     1453}
     1454
     1455proc BkgEval {terms num tlist "rlist {}"} {
     1456    global expmap hst
     1457    if {$num == 1} {
     1458        global tmin tmax
     1459        foreach x $tlist {
     1460            lappend blist [chebeval $terms $x $tmin $tmax]
     1461        }
     1462        return $blist
     1463    } elseif {$num == 2} {
     1464        set ts 1
     1465        if {[string range $expmap(htype_$hst) 2 2] == "T"} {
     1466            catch {
     1467                set line [histinfo $hst ITYP]
     1468                set ts [expr 180./ [lindex $line 2]]
     1469            }
     1470        }
     1471        foreach tof $tlist {
     1472            set tofm [expr {$tof * $ts}]
     1473            set bkg 0
     1474            set i -1
     1475            foreach t $terms {
     1476                incr i
     1477                set bkg [expr {$bkg + $t * cos($i * $tofm * 3.14159/180.)}]
     1478            }
     1479            lappend blist $bkg
     1480        }
     1481        return $blist
     1482    } elseif {$num == 3} {
     1483        set Qlist [toQ $tlist $hst]
     1484        foreach Q $Qlist tofm $tlist {
     1485            set i 0
     1486            set j -1
     1487            foreach t $terms {
     1488                incr i
     1489                if {$i == 1} {
     1490                    set bkg $t
     1491                } elseif {$i == 2} {
     1492                    set bkg [expr {$bkg + $tofm * $t}]
     1493                } else {
     1494                    incr j
     1495                    set r [lindex $rlist $j]
     1496                    set QR [expr {$Q * $r}]
     1497                    set bkg [expr {$bkg + $t * sin($QR)/$QR}]
     1498                }
     1499            }
     1500            lappend blist $bkg
     1501        }
     1502        return $blist
     1503    } elseif {$num == 4} {
     1504        set Qlist [toQ $tlist $hst]
     1505        foreach Q $Qlist {
     1506            set i -1
     1507            set QT 1
     1508            foreach t $terms {
     1509                incr i
     1510                if {$i == 0} {
     1511                    set bkg $t
     1512                } else {
     1513                    set QT [expr {$QT * $Q * $Q / $i}]
     1514                    set bkg [expr {$bkg + $t * $QT}]
     1515                }
     1516            }
     1517            lappend blist $bkg
     1518        }
     1519        return $blist
     1520    } elseif {$num == 5} {
     1521        set Qlist [toQ $tlist $hst]
     1522        foreach Q $Qlist {
     1523            set i -1
     1524            set QT 1
     1525            foreach t $terms {
     1526                incr i
     1527                if {$i == 0} {
     1528                    set bkg $t
     1529                } else {
     1530                    set QT [expr {$QT * $i /($Q * $Q)}]
     1531                    set bkg [expr {$bkg + $t * $QT}]
     1532                }
     1533            }
     1534            lappend blist $bkg
     1535        }
     1536        return $blist
     1537    } elseif {$num == 6} {
     1538        set Qlist [toQ $tlist $hst]
     1539        foreach Q $Qlist {
     1540            set i 0
     1541            set QT 1
     1542            foreach t $terms {
     1543                incr i
     1544                if {$i == 1} {
     1545                    set bkg $t
     1546                } elseif {$i % 2} {
     1547                    # odd
     1548                    set QT1 [expr {1./$QT}]
     1549                    set bkg [expr {$bkg + $t * $QT1}]
     1550                } else {
     1551                    # even
     1552                    set QT [expr {2*$QT*$Q*$Q/$i}]
     1553                    set QT1 $QT
     1554                    set bkg [expr {$bkg + $t * $QT1}]
     1555                }
     1556            }
     1557            lappend blist $bkg     
     1558        }
     1559        return $blist
     1560    }
     1561}
     1562
     1563proc backderivcal {nterms num tof "rlist {}"} {
     1564    global expmap hst
     1565    if {$num == 1} {
     1566        global tmin tmax
     1567        # rescale x
     1568        set xs [expr {-1 + 2 * (1.*$tof - $tmin) / (1.*$tmax - 1.*$tmin)}]
     1569        # compute the Chebyschev term Tn(xs)
     1570        set deriv {}
     1571        set Tpp 0
     1572        set Tp 0
     1573        for {set i 0} {$i < $nterms} {incr i} {
     1574            if {$Tpp == $Tp && $Tp == 0} {
     1575                set T 1
     1576            } elseif {$Tpp == 0} {
     1577                set T $xs
     1578            } else {
     1579                set T [expr {2. * $xs * $Tp - $Tpp}]
     1580            }
     1581            lappend deriv $T
     1582            set Tpp $Tp
     1583            set Tp $T
     1584        }
     1585        return $deriv
     1586    } elseif {$num == 2} {
     1587        set ts 1
     1588        if {[string range $expmap(htype_$hst) 2 2] == "T"} {
     1589            catch {
     1590                set line [histinfo $hst ITYP]
     1591                set ts [expr 180./ [lindex $line 2]]
     1592            }
     1593            set tofm [expr {$tof * $ts}]
     1594        } else {
     1595            set tofm $tof
     1596        }
     1597        set deriv {}
     1598        for {set i 0} {$i < $nterms} {incr i} {
     1599            lappend deriv [expr {cos($i * $tofm * 3.14159/180.)}]
     1600        }
     1601        return $deriv
     1602    } elseif {$num == 3} {
     1603        set Q [toQ $tof $hst]
     1604        set j -1
     1605        #set n [expr {2 + ($nterms - 2)/2}]
     1606        for {set i 1} {$i <= $nterms} {incr i} {
     1607            if {$i == 1} {
     1608                set deriv 1
     1609            } elseif {$i == 2} {
     1610                lappend deriv $tof
     1611            } else {
     1612                incr j
     1613                set r [lindex $rlist $j]
     1614                set QR [expr {$Q * $r}]
     1615                lappend deriv [expr {sin($QR)/$QR}]
     1616            }
     1617        }
     1618        return $deriv
     1619    } elseif {$num == 4} {
     1620        set Q [toQ $tof $hst]
     1621        set QT 1
     1622        for {set i 0} {$i < $nterms} {incr i} {
     1623            if {$i == 0} {
     1624                set deriv 1
     1625            } else {
     1626                lappend deriv [set QT [expr {$QT * $Q * $Q / $i}]]
     1627            }
     1628        }
     1629        return $deriv
     1630    } elseif {$num == 5} {
     1631        set Q [toQ $tof $hst]
     1632        set QT 1
     1633        for {set i 0} {$i < $nterms} {incr i} {
     1634            if {$i == 0} {
     1635                set deriv 1
     1636            } else {
     1637                lappend deriv [set QT [expr {$QT * $i /($Q * $Q)}]]
     1638            }
     1639        }
     1640        return $deriv
     1641    } elseif {$num == 6} {
     1642        set Q [toQ $tof $hst]
     1643        set QT 1
     1644        for {set i 1} {$i <= $nterms} {incr i} {
     1645            if {$i == 1} {
     1646                set deriv 1
     1647            } elseif {$i % 2} {
     1648                # odd
     1649                lappend deriv [set QT1 [expr {1./$QT}]]
     1650            } else {
     1651                # even
     1652                set QT [expr {2*$QT*$Q*$Q/$i}]
     1653                lappend deriv [set QT1 $QT]
     1654            }
     1655        }
     1656        return $deriv
     1657    }
     1658}
     1659
     1660# evaluate the best-fit background terms to fit GSAS background functions 1-6
     1661# to a set of X and Y values.
     1662# num is the function number,
     1663# order is the # of terms
     1664# rlist is used only for function type 3; there must be (order-2)/2 values
     1665proc FitBkgFunc {X Y order num "rlist {}"} {
     1666    if {$num == 3} {
     1667        set o [expr {2 + ($order - 2)/2}]
     1668    } else {
     1669        set o $order
     1670    }
     1671    # zero the matrix and vector
     1672    for {set j 0} {$j < $o} {incr j} {
     1673        set sum($j) 0.
     1674        for {set i 0} {$i <= $j} {incr i} {
     1675            set sum(${i}_$j) 0.
     1676        }
     1677    }
     1678global octave
     1679set octave {}
     1680append octave {des = [}
     1681    foreach y $Y x $X {
     1682        # compute derivatives at point x
     1683        set derivlist [backderivcal $o $num $x $rlist]
     1684append octave " $derivlist ;\n"
     1685        # compute matrix elements
     1686        for {set j 0} {$j < $o} {incr j} {
     1687            set Tj [lindex $derivlist $j]
     1688            # weighted
     1689            # set sum($j) [expr {$sum($j) + $y * $Tj / ($sigma*$sigma)}]
     1690            set sum($j) [expr {$sum($j) + $y * $Tj}]
     1691            for {set i 0} {$i <= $j} {incr i} {
     1692                set Ti [lindex $derivlist $i]
     1693                # weighted
     1694                # set sum(${i}_$j) [expr {$sum(${i}_$j) + $Ti * $Tj / ($sigma * $sigma)}]
     1695                set sum(${i}_$j) [expr {$sum(${i}_$j) + $Ti * $Tj}]
     1696            }
     1697        }
     1698    }
     1699    # populate the matrix & vector in La format
     1700    lappend V 2 $o 0
     1701    lappend A 2 $o $o
     1702    for {set i 0} {$i < $o} {incr i} {
     1703        lappend V $sum($i)
     1704        for {set j 0} {$j < $o} {incr j} {
     1705            if {$j < $i} {
     1706                lappend A $sum(${j}_$i)
     1707            } else {
     1708                lappend A $sum(${i}_$j)
     1709            }
     1710        }
     1711    }
     1712    set termlist {}
     1713    if {[catch {
     1714        set termlist [lrange [La::msolve $A $V] 3 end]
     1715    }]} {
     1716        tk_dialog .singlar "Singular Matrix" \
     1717            "Unable to fit function: singular matrix. Too many terms or something else is wrong." ""\
     1718            0 OK
     1719    }
     1720    return $termlist
     1721}
     1722
    13391723# save the Chebyshev terms in the .EXP file
    1340 proc bkgChebSave {} {
    1341     global hst cheblist expgui Revision expmap expnam
    1342     histinfo $hst backtype set 1
    1343     histinfo $hst backterms set [llength $cheblist]
     1724proc bkgSave {} {
     1725    global hst termlist expgui Revision expmap expnam
     1726    histinfo $hst backtype set $expgui(FitFunction)
     1727    # stick the r values into the list
     1728    if {$expgui(FitFunction) == 3} {
     1729        set t [lrange $termlist 0 1]
     1730        foreach a [lrange $termlist 2 end] b $expgui(RadiiList) {lappend t $a $b}
     1731    } else {
     1732        set t $termlist
     1733    }
     1734    histinfo $hst backterms set [llength $t]
    13441735    set num 0
    1345     foreach v $cheblist {
     1736    foreach v $t {
    13461737        set var "bterm[incr num]"
    13471738        histinfo $hst $var set $v
    … …  
    14261817if [file executable [file join $expgui(gsasexe) $expgui(tcldump)]] {
    14271818    set expgui(tcldump) [file join $expgui(gsasexe) $expgui(tcldump)]
    1428 #    puts "got tcldump"
    14291819} else {
    14301820    set expgui(tcldump) {}
    1431 #    puts "no tcldump"
    14321821}
    14331822
    … …  
    15051894            -line 0 -symbol $peakinfo(inpsym) -label "bkg pts"
    15061895    $box element configure 11 -color $graph(color_fit) \
    1507             -symbol none -label "Cheb fit" -dashes 5 -line 2
     1896            -symbol none -label "bkg fit" -dashes 5 -line 2
    15081897    $box element show "3 2 11 12"
    15091898}
    … …  
    15411930.a.file.menu add cascade -label Tickmarks -menu .a.file.menu.tick
    15421931menu .a.file.menu.tick
    1543 foreach num {1 2 3 4 5 6 7 8 9} {
    1544     .a.file.menu.tick add checkbutton -label "Phase $num" \
    1545             -variable  peakinfo(flag$num) \
    1546             -command plotdata
    1547 }
    15481932.a.file.menu add cascade -label Histogram -menu .a.file.menu.hist -state disabled
    15491933.a.file.menu add command -label "Update Plot" \
    15501934        -command {set cycle [getcycle];readdata .g}
    1551 .a.file.menu add command -label "Make PostScript" -command makepostscriptout
     1935.a.file.menu add cascade -label "Export plot" -menu .a.file.menu.export
     1936menu .a.file.menu.export
     1937.a.file.menu.export add command -label "to PostScript" \
     1938        -command makepostscriptout
     1939if {$blt_version > 2.3 && $blt_version != 8.0} {
     1940    source [file join $expgui(scriptdir) graceexport.tcl]
     1941    .a.file.menu.export add command -label "to Grace" -command exportgrace
     1942}
    15521943.a.file.menu add command -label Quit -command "destroy ."
    15531944
    … …  
    15621953        -value 1 -variable expgui(autotick) -command plotdata
    15631954.a.options.menu.tick add separator
    1564 foreach num {1 2 3 4 5 6 7 8 9} {
    1565     .a.options.menu.tick add command -label "Phase $num" \
    1566             -command "minioptionsbox $num"
    1567 }
     1955.a.options.menu.tick add command -label "Label by name" \
     1956        -command {
     1957    foreach p $expmap(phaselist) {
     1958        # 20 characters, max
     1959        set graph(label$p) [string range [phaseinfo $p name] 0 19]
     1960        plotdata
     1961    }
     1962}
     1963.a.options.menu.tick add separator
     1964
    15681965if {$program == "liveplot"} {
    15691966    .a.options.menu add command -label "Obs symbol" -command getsymopts
    … …  
    15811978if {$program != "liveplot"} {
    15821979    lappend l1 input fit
    1583     lappend l2 "Input points" "Cheb. fit"
     1980    lappend l2 "Input points" "bkg fit"
    15841981}
    15851982   
    … …  
    16382035} elseif {$program == "bkgedit"}  {
    16392036    catch {pack [frame .bkg -bd 3 -relief sunken] -side bottom -fill both}
    1640     grid [label .bkg.top -text "Background Point Editing"] \
    1641             -col 0 -row 0 -columnspan 4
    1642     grid [button .bkg.help -text Help -bg yellow \
    1643             -command "MakeWWWHelp liveplot.html bkgedit"] \
    1644             -column 5 -row 0 -rowspan 2 -sticky n
     2037#    grid [label .bkg.top -text "Background Point Editing"] \
     2038#           -col 0 -row 0 -columnspan 4
     2039#    grid [button .bkg.help -text Help -bg yellow \
     2040#           -command "MakeWWWHelp liveplot.html bkgedit"] \
     2041#           -column 5 -row 0 -rowspan 2 -sticky n
    16452042   
    16462043    grid [frame .bkg.l -bd 3 -relief groove] \
    … …  
    16512048                -col $c -row 0
    16522049    }
     2050    # leave a small blank space
     2051    grid columnconfigure .bkg 2 -pad 0 -min 10
    16532052    grid [frame .bkg.f -bd 3 -relief groove] \
    16542053            -col 3 -row 1 -columnspan 2 -sticky nsw
    … …  
    16562055            -col 1 -row 1
    16572056    grid [label .bkg.f.tl -text "with"] -col 3 -row 1
    1658     set termmenu [tk_optionMenu .bkg.f.terms chebterms 0]
     2057    set termmenu [tk_optionMenu .bkg.f.terms expgui(FitOrder) 0]
    16592058    grid .bkg.f.terms -col 4 -row 1
    16602059    grid [label .bkg.f.tl1 -text "terms"] -col 5 -row 1
    … …  
    16622061    grid [frame .bkg.c1 -bd 3 -relief groove] \
    16632062            -col 0 -row 5 -rowspan 2 -sticky nsew
    1664     grid [label .bkg.c1.1 -text "Chebyshev\nterms"] -col 0 -row 0
     2063    grid [label .bkg.c1.0 -text "Background\nfunction #"] -col 0 -row 0
     2064    set bkgmenu [tk_optionMenu .bkg.c1.1 expgui(FitFunction) stuff]
     2065    grid .bkg.c1.1 -col 0 -row 1
     2066    $bkgmenu delete 0 end
     2067    foreach item {
     2068        "1 - Shifted Chebyschev polynomial"
     2069        "2 - Cosine Fourier series"
     2070        "3 - Radial distribution peaks"
     2071        "4 - Power series in Q**2n/n!"
     2072        "5 - Power series in n!/Q**2n"
     2073        "6 - Power series in Q**2n/n! and n!/Q**2n"
     2074    } {
     2075        set val [lindex $item 0]
     2076        $bkgmenu insert end radiobutton -variable expgui(FitFunction) \
     2077                -label $item -value $val \
     2078                -command "set termlist {};BkgFillTermBoxes nosave"
     2079    }
     2080    set expgui(FitFunction) 1
     2081
    16652082    grid [canvas .bkg.canvas \
    16662083            -scrollregion {0 0 5000 500} -width 0 -height 0 \
    … …  
    16692086    grid [scrollbar .bkg.scroll -command ".bkg.canvas xview" \
    16702087            -orient horizontal] -column 1 -row 6 -columnspan 3 -sticky nsew
    1671     grid [button .bkg.cw -text "Save in EXP\nfile & Exit" \
    1672             -command "bkgChebSave;exit"] \
     2088    grid [button .bkg.cw -text "Save in\nEXP file\n& Exit" \
     2089            -command "bkgSave;exit"] \
    16732090            -col 4 -columnspan 2 -row 5 -rowspan 2 -sticky ns
    16742091
    … …  
    16932110pack [menubutton .a.help -text Help -underline 0 -menu .a.help.menu] -side right
    16942111menu .a.help.menu -tearoff 0
    1695 .a.help.menu add command -command "MakeWWWHelp liveplot.html" -label "Web page"
     2112if {$program == "bkgedit"}  {
     2113    .a.help.menu add command -command "MakeWWWHelp liveplot.html bkgedit" \
     2114            -label "Web page"
     2115} else {
     2116    .a.help.menu add command -command "MakeWWWHelp liveplot.html" \
     2117            -label "Web page"
     2118}
    16962119.a.help.menu add command -command aboutliveplot -label About
    16972120
    … …  
    17572180    set expgui(plotlist) [lindex $expmap(powderlist) 0]
    17582181}
     2182
     2183foreach num $expmap(phaselist) {
     2184    .a.file.menu.tick add checkbutton -label "Phase $num" \
     2185            -variable  peakinfo(flag$num) \
     2186            -command plotdata
     2187    bind . <Key-$num> ".a.file.menu.tick invoke $num"
     2188    .a.options.menu.tick add command -label "Phase $num opts" \
     2189            -command "minioptionsbox $num"
     2190}
     2191
    17592192# N = load next histogram
    17602193bind . <Key-n> {
Note: See TracChangeset for help on using the changeset viewer.