#!/bin/sh
# the next line restarts using wish \
exec wish "$0" "$@"

# thermal.tcl stripped down version
#                                   David J C MacKay (1998)

set ownwindow 0
set verbose 0
if {$verbose>=1} {
    puts "Classical, Boson and Fermi version 2.1                David J C MacKay          1998"
    puts "                                             written under tcl 8.0 on linux"
}
# changes
#########################################################################
#
######################################################################
# contents:
#
######################################################################

frame .th
pack  .th
set w .th
set top .
if {$ownwindow>=1} {
    wm title . "Boltzmann distribution"
    wm iconname . "Boltzmann"
    wm geometry . +10+10
}

#######################################################################
#
#  variables
#
#######################################################################

set types "classical bosons fermions" ;# classical, boson, fermi
set color(classical) "red"
set color(bosons) "green"
set color(fermions) "yellow"
set palecolor(classical) "pink"
set palecolor(bosons) "palegreen"
set palecolor(fermions) "lightgoldenrod1"
global L ; set L 5 ;# number of levels
global N ; set N 3 ;# number of particles
global M ; set M 10 ;# number of microstates
foreach t [concat $types] {
    set Active($t) 1
}
global Ecutoff ;
set Ecutoff 15 ;# highest microstate energy worth considering when 
                #  building the display
global Emmax ; set Emmax 15 ;# highest energy of a microstate to display
global EcutoffF ; set EcutoffF 40 ;# highest energy of a fermi microstate to create
set bfmin 0.0 ;# boltzmann factors / probabilities
set bfmax 1.0 ;
set pmin 0.0 ;# boltzmann factors / probabilities
set pmax 1.0 ;
set betamin 0.05 ;#minimum temperature energy

#######################################################################
#
#                      procedures  for  packing
#
#######################################################################

proc makeControls { w controls } {
    global types palecolor color
    pack [frame $controls.rA] -side top -anchor w
    pack [frame $controls.rB] -side top -anchor w
    pack [button $w.restart -text "Restart" -borderwidth 1 \
	    -command {restart}] -in $controls.rA  -side left \
	    -pady 0 -padx 2 -anchor w
    # what is active
    foreach t [concat $types] {
        checkbutton $w.$t -text $t -variable Active($t) \
		-command ""  -background $palecolor($t)
	pack $w.$t -in $controls.rA  -side left -pady 2 -anchor w
    }
    pack [button $w.restart2 -text "Restart" -borderwidth 1 \
	    -command {restart}] -in $controls.rA  -side left \
	    -pady 0 -padx 2 -anchor w

    # how many of this and that 

    adjustableInteger $w $controls.rB "l" "L" "L"
    adjustableInteger $w $controls.rB "n" "N" "N"
    adjustableInteger $w $controls.rB "m" "M" "M"
    adjustableInteger $w $controls.rB "max" "Ecutoff" "Cut"
    adjustableInteger $w $controls.rB "emmax" "Emmax" "MaxDisplay"
    adjustableInteger $w $controls.rB "maxf" "EcutoffF" "CutF"
}
# make a frame called l within frame controls, and associate these buttons 
# with an integer called L
proc adjustableInteger { w controls l Lname Lstring } {
    frame $w.$l 
    pack $w.$l -in $controls -side left  -pady 2  -padx 2 -anchor w
    button $w.$l.l -text "$Lstring:" -padx 0 -pady 0 -borderwidth 1 -command {}
    button $w.$l.up -text ">" -padx 0 -pady 0 -borderwidth 1 -command "incr $Lname"
    button $w.$l.dn -text "<" -padx 0 -pady 0 -borderwidth 1 -command "incr $Lname -1"
    bind  $w.$l.up  <3> "$w.$l.dn invoke"
    bind  $w.$l.dn  <3> "$w.$l.up invoke"
    entry $w.$l.n -textvariable $Lname -width 3 -borderwidth 1
    pack $w.$l.l $w.$l.dn $w.$l.up $w.$l.n -in $w.$l -side left

}

proc makeEnergylevels { w energylevels } {
    global L Elmin Elmax ocmax ocmin betamin ;
    global energy types color meanoc verbose

    # make our own local frame
    set e $energylevels.e
    catch { destroy $e }
    pack [frame $e]

    frame $e.labels 
    pack $e.labels -in $e -side top -pady 2 -expand 1 -fill x
    label $e.labels.e -text "Energy levels" -anchor w
    label $e.labels.o -text "Occupancies" -anchor w
    pack $e.labels.e -in $e.labels -side left -pady 0
    pack $e.labels.o -in $e.labels -side right -pady 0
    for { set l 1 } { $l <= $L } { incr l } {
	set energy($l) [expr $l-1]
	catch { destroy $e.ee$l  ; destroy $e.e$l }
	if {$verbose>=2} { puts "elmax = $Elmax" }
	set ee$l [frame $e.ee$l] ; pack [set ee$l] -in $e -side top -pady 1
	set e$l [scale $e.e$l -orient horizontal -length 280  \
		-from $Elmin -to $Elmax   -width 8 -sliderlength 8 -background blue \
		-borderwidth 0 -showvalue 0  \
		-variable energy($l) -tickinterval 0 -resolution 0.1  \
		-bigincrement 0     -command {recomputeMicrostateEnergies}]
	pack [set e$l] -in [set ee$l] -side left -pady 1 -padx 1 -expand 1 -fill x
	# show occupancies too.
	foreach  t [concat $types] {
	    catch { destroy $e.o$l$t  }
	    set meanoc($l,$t) 1
	    set o$l$t [scale $e.o$l$t -orient horizontal -length 280  \
		    -from $ocmin -to $ocmax   -width 6 -sliderlength 6 \
		    -borderwidth 0 -showvalue 0  -background \
		    $color($t) \
		    -variable meanoc($l,$t) -tickinterval 0 \
		    -resolution 0.001 ]
	    pack [set o$l$t] -in [set ee$l] -side top  -padx 1 -expand 1 -fill x
	}
    }
    raise    $e.labels
    catch { destroy $e.betalabel ; destroy $e.ebeta  }
    label $e.betalabel -text "Temperature" -anchor w
    pack $e.betalabel -in $e -side top -pady 2 -anchor w
    set ebeta [scale $e.ebeta -orient horizontal -length 284  \
	    -from $betamin -to $Elmax   -width 8 \
	    -borderwidth 0 -showvalue 0  -sliderlength 8 -background blue \
	    -variable temp -tickinterval 0 -resolution $betamin  \
	    -bigincrement 0     -command {set_temp}]
    pack [set ebeta] -in $e -side top -pady 1 -padx 1 -anchor w

    catch { destroy $e.zframe }
# partition functions (these don't quite belong here... 
    set zframe [frame $e.zframe]
    pack $zframe -in $e -side right
    global Z palecolor
    pack [label $zframe.zl -text "Partition functions"]  \
	    -side left -pady 0 -padx 10
    foreach t [concat $types] {
	set Z($t) 1.0 ;
	set z$t [entry $zframe.z$t  -width 10  -background \
			$palecolor($t) -borderwidth 0 -textvariable Z($t)]
	pack [set z$t]  -side left -pady 0 -padx 6
    }

}

#######################################################################
#
#                      procedures
#
#######################################################################

global beta; set beta 1
global temp; set temp 1
proc set_temp { t } {
    global beta
    set beta [expr 1.0/$t]
    computeMicrostateProbs
}

proc makeFactorials { N } {
    global factorial
    set factorial(0) 1 
    set now 1
    for { set m 1 } { $m <= $N } { incr m } {
	set now [expr $now*$m]
	set factorial($m) $now
    }
}

# invoked when the energy level scales are touched
proc recomputeMicrostateEnergies { {junk 0} } {
    global energy L M e_m occupancy
    for { set m 1 } { $m <= $M } { incr m } {
	set e 0
	for { set l 1 } { $l <= $L } { incr l } {
	    set n $occupancy($m,$l)
	    set e [expr $e+$energy($l)*$n]
	}
	set e_m($m) $e
    }
    computeMicrostateProbs
}

# invoked when the energy level scales are touched
# finds probabilities and occupancies
proc computeMicrostateProbs { {junk 0} } {
    global energy L M e_m types beta degeneracy bf p meanoc occupancy Z
    global Active
    foreach  t [concat $types] {
	for { set l 1 } { $l <= $L } { incr l } {
	    set meanoc($l,$t) 0
	}
    }
    # compute boltz factors once only
    for { set m 1 } { $m <= $M } { incr m } {
	set bf($m) [expr exp(-1.0*$beta*$e_m($m))]
    }
    foreach  t [concat $types] {
	if {$Active($t)} {
	    set Z($t) 0.0
	    for { set m 1 } { $m <= $M } { incr m } {
		set d $degeneracy($m,$t)
		if {$d} {
		    set thisbf  [expr 1.0*$d*$bf($m)]
		    set Z($t) [expr $Z($t) + $thisbf]
		}
	    }
	    set IZ [expr 1.0/$Z($t)]
	    for { set m 1 } { $m <= $M } { incr m } {
		set d $degeneracy($m,$t)
		if {$d} {
		    set p($m,$t) [expr 1.0*$d*$bf($m)*$IZ]
		    for { set l 1 } { $l <= $L } { incr l } {
			set meanoc($l,$t) [expr \
				$meanoc($l,$t)+ $p($m,$t)*$occupancy($m,$l)]
		    }
		} 
	    }
	}
    }
}

# computes energies and degeneracies of all microstates
proc findMicrostates { } {
    global microstate ; set microstate 0
    global L N factorial  Fermions Active Ecutoff EcutoffF visited
    set mstring ""  ; 
    set l 1 ; set l_left $L ; set n_left $N ; set degen $factorial($N)
    set degenf 1 ; 
    set e 0 ;
    set Fermions 0 
    catch { unset visited }

    # initiate recursive function to make all states up to Ecutoff

    continueMicrostate $l $l_left $n_left $degen $mstring $e $degenf

    # repeat a recursive function to catch some more fermion states
    if {$Ecutoff&&($EcutoffF>$Ecutoff)} {
	set l 1 ; set l_left $L ; set n_left $N ; set degen $factorial($N)
	set degenf 1 ; 
	set e 0 ;
	continueMicrostate $l $l_left $n_left $degen $mstring $e $degenf 1
	# extra argument indicates fermion mode
    }
    # finished
    global M ; set M $microstate
    if {$Fermions<=0} {
	set Active(fermions) 0
    }
    return $M
}

# recursive routine for finding all microstates.
# if mode=fermions then only fermion states are produced
proc continueMicrostate { l l_left n_left degen mstring e degenf {fermion_mode 0} } {
    global energy e_m microstring palecolor types degeneracy microstate 
    global current_ns occupancy L EcutoffF Ecutoff Fermions visited verbose
    if {$l_left==1} {
	# extend the string with n_left and make this a new microstate
	# if we are in normal mode 'all' or else the quanta remaining are 
	# valid for a fermion.
	if {!$fermion_mode||($n_left<=1)} {
	    extendMicrostate  $l $n_left $n_left $degen $mstring $e  $degenf 
	    # if  the energy is not too big, or else we are in fermion mode
	    # and the energy *is* too big, ...
	    if {$fermion_mode} {
		set energyok [expr (($EcutoffF==0)||($EcutoffF>=$eplus))]
	    } else {
		set energyok [expr (($Ecutoff==0)||($Ecutoff>=$eplus))]
	    }
	    if [info exists visited($mstringplus)] {
		if {$verbose>=2} {
		    puts "visited $mstringplus already" 
		}
		set beenhere 1
	    } else {
		set beenhere 0
	    }
	    if {(!$beenhere)&&($energyok)} {
		incr microstate ; set m $microstate
		set degeneracy($m,classical) $degenplus
		set degeneracy($m,bosons) 1
		set degeneracy($m,fermions) $degenfplus
		incr Fermions $degenfplus
		set e_m($m) $eplus
		set microstring($m) $mstringplus
		set visited($mstringplus) $m
		# store occupancies in a convenient array
		for { set l 1 } { $l <= $L } { incr l } {
		    set occupancy($m,$l) [string index $mstringplus [expr $l-1]]
		}
	    }
	}
    } else {
	# loop through all valid numbers that could be put here
	set n_validmax $n_left 
	if {($fermion_mode)&&($n_left>1)} {
	    set n_validmax 1
	}
	for { set n $n_validmax } { $n >= 0 } { incr n -1 } {
	    extendMicrostate  $l $n_left $n $degen $mstring $e $degenf 
	    if {$fermion_mode} {
		set energyok [expr (($EcutoffF==0)||($EcutoffF>=$eplus))]
	    } else {
		set energyok [expr (($Ecutoff==0)||($Ecutoff>=$eplus))]
	    }
	    if {$energyok} {
		continueMicrostate [expr $l+1] \
			[expr $l_left-1] $n_leftplus  $degenplus \
			$mstringplus  $eplus $degenfplus $fermion_mode
	    }

	}
    }
}

proc extendMicrostate { l n_left n degen mstring e degenf {fermion_mode 0}} {
    set spacer "" 
    global factorial energy  current_ns
    upvar mstringplus mstringplus
    upvar eplus eplus 
    upvar degenplus degenplus
    upvar n_leftplus n_leftplus
    upvar degenfplus degenfplus
    # extend the string with n
    set current_ns($l) $n
    set mstringplus "$mstring$spacer$n" 
    set eplus [expr $e+$energy($l)*$n]
    set degenfplus [expr (($n>1)?0:$degenf)]
    set n_leftplus [expr $n_left-$n]
    set f [set factorial($n)]
    set degenplus [expr $degen/$f]
}

# makes the widgets for showing the m's
proc makeMicrostates { w micro } {
    global M Emmin Emmax L verbose  bfmin bfmax pmin pmax ;
    global e_m microstring color palecolor types degeneracy ; 
    set tinyfont "Courier 8" ; set degenwidth 3

    set e $micro.e 
    catch { destroy $e }
    pack [frame $e]

    # top row
    label $e.mlabel -text "Microstates" -anchor w
    pack $e.mlabel -in $e -side top -pady 2
    #
    for { set m 1 } { $m <= $M } { incr m } {
	set thisrow [frame $e.e$m]
	pack $thisrow -in $e -side top -pady 0 -anchor w
	# show microstate energy
	set m$m [scale $e.m$m -orient horizontal -length 184  \
		-from $Emmin -to $Emmax   -width 6 \
		-borderwidth 0 -showvalue 0  -sliderlength 6 -background blue \
		-variable e_m($m) -tickinterval 0 -resolution 0.1  \
		-bigincrement 0     -command {}]
	set ml$m [label $e.ml$m -text $microstring($m) -font $tinyfont -width $L -anchor w -borderwidth 0]
	pack [set ml$m] [set m$m] -in $thisrow -side left -pady 0 -padx 1
	foreach t [concat $types] {
	    if {$verbose>=2} {
		puts " degeneracy($m,$t)"
		puts [set degeneracy($m,$t)]
	    }
	    set ml$m$t [entry $e.ml$m$t -textvariable degeneracy($m,$t) \
		    -font $tinyfont -width $degenwidth -background \
		    $palecolor($t) -borderwidth 0]
	    pack [set ml$m$t] -in $thisrow -side left -pady 0 -padx 1
	}
	# boltzmann factors:
	set mlbf$m$t [scale $e.mlbf$m$t -orient horizontal -length 86  \
		-from $bfmin -to $bfmax   -width 6  -background \
		blue	-borderwidth 0 -showvalue 0  -sliderlength 6 \
		-variable bf($m) -tickinterval 0 -resolution 0.001  \
		-bigincrement 0     -command {}]
	pack [set mlbf$m$t] -in $thisrow -side left -pady 0 -padx 5
	# probabilities
	foreach t [concat $types] {
	    if {$degeneracy($m,$t)} {
		set mlp$m$t [scale $e.mlp$m$t -orient horizontal -length 86  \
			-from $pmin -to $pmax   -width 6  -background \
			$color($t) -borderwidth 0 -showvalue 0 -sliderlength 6\
			-variable p($m,$t) -tickinterval 0 -resolution 0.001  \
			-bigincrement 0     -command {}]
		pack [set mlp$m$t] -in $thisrow -side left -pady 0 -padx 1
	    }
	}
    }
}

#
# packing procedures
#
proc restart { } {
    global N L 
    upvar w w   energylevels energylevels  controls controls  
    upvar microstates microstates

    global Elmin ; set Elmin 0 ;
    global Emmin ; set Emmin 0 ;
    global Elmax ; set Elmax $L ;# highest energy for a level
    global ocmin ; set ocmin 0 ;
    global ocmax ; set ocmax $N;

    makeEnergylevels $w $energylevels
    makeFactorials $N 
    findMicrostates
    makeMicrostates $w $microstates
}

####################################################################
#
# set up windows 
# 
####################################################################

bind . <Control-x> "destroy ."
bind . <Control-q> "destroy ."
bind . <Control-c> "destroy ."
bind . <Control-z> "destroy ."

frame $w.controls
set controls $w.controls
pack $controls -side top
makeControls $w $controls

frame $w.middlerow
set middlerow $w.middlerow
pack $middlerow -side top

frame $w.energylevels
set energylevels $w.energylevels
pack $energylevels -in $middlerow -side top -expand y -fill x

frame $w.microstates
set microstates $w.microstates
pack $microstates -in $middlerow -side left

# buttons to do with overall control (quit, etc.)
frame $w.buttons
pack $w.buttons -side bottom -fill x -pady 2 -expand 1

#
# overall control buttons
#
button $w.dismiss -text Quit -command "destroy $top"
button $w.help -text Help -command "help"

pack $w.dismiss  $w.help \
	-in $w.buttons -side left -fill x -expand 1 -anchor w  -padx 3 -pady 1
############
#  end bottom row
############

# make it happen!
restart

#####################################################################

proc help { } {
    set w .help
    catch {destroy $w}
    toplevel $w
    wm geometry $w +10+10
    bind $w <Control-c> "destroy $w"
    frame $w.buttons
    pack $w.buttons -side bottom -fill x -pady 2 -expand 1 -padx 4
    button $w.buttons.dismiss -text Dismiss -command "destroy $w"
    pack $w.buttons.dismiss  -side left -fill x -expand 1 -padx 4

  text $w.t -background white -height 24 -wrap word\
            -xscrollcommand "$w.xscroll set" \
            -yscrollcommand "$w.yscroll set" \
            -setgrid 1 -highlightthickness 0 -pady 2 -padx 3
        scrollbar $w.xscroll -command "$w.t xview" \
            -highlightthickness 0 -orient horizontal
        scrollbar $w.yscroll -command "$w.t yview" \
            -highlightthickness 0 -orient vertical

pack $w.yscroll -side right -fill y
pack $w.xscroll -side bottom -fill x
pack $w.t -expand yes -fill both

    $w.t insert 0.0 \
{Boltzmann distributions          - author David J C MacKay  mackay@mrao.cam.ac.uk


 General layout:

 Shortcuts: 
     C-r     reset

}
}