#! /usr/bin/env gnuplot

# This work is in the public domain.
#
# It was completed by Arun I in 2016.
#
# Though you are not legally obliged to do so, I would appreciate it
# if you credit me for this work by not removing this notice, and
# hopefully linking to my blog post at
# https://systemreboot.net/post/black-body-energy-as-a-sum-of-photon-oscillator-energies

# Base filename for output files
filename='blackbody-energies'

# Set terminal and output file
set term epslatex color colortext standalone
set output sprintf('%s.tex', filename)

# Axes' properties, title and key
set xrange [-0.01:5]
set yrange [0:1.05]
set mxtics 5
set ytics 0.2
set mytics 4
set xlabel '$ h\nu/kT $'
set ylabel 'Average oscillator energy (in units of $ kT $)'
set title "Planck's solution as a freeze-out of\nphoton oscillator-energies as $ h\\nu/kT \\to \\infty $"
unset key

# Line styles
set style line 1 lc 'blue' dt 4
set style line 2 lc 'red' dt 4
set style line 3 lc 'dark-yellow' dt 4
set style line 4 lc 'forest-green' dt 4

# Labels
set arrow from graph 0, first 1 to graph 1, first 1 nohead lc 'black' dt 6
set label '$ E_{\textrm{\tiny avg}}^{\textrm{\tiny Rayleigh}} = kT $' at graph 1, first 1 right offset graph -0.05, -0.05
set label '\footnotesize $ h\nu $' at 1.4, 0.15 tc ls 1
set label '\footnotesize $ 2h\nu $' at 1, 0.33 tc ls 2
set label '\footnotesize $ 3h\nu $' at 0.8, 0.45 tc ls 3
set label '\footnotesize $ 4h\nu $' at 0.65, 0.53 tc ls 4
set label '$ E_{\textrm{\tiny avg}}^{\textrm{\tiny Planck}} $' at 0.6, 0.8

# Expressions to plot
En(x,n) = n*x*exp(-n*x)*(1-exp(-x))
yN(x,N) = sum [n=1:N] En(x,n)

plot x/(exp(x) - 1) with lines lc 'black', \
     for [N=1:4] yN(x,N) with lines ls N

# Convert tex -> pdf -> svg
set output
system(sprintf('pdflatex %s.tex && pdf2svg %s.pdf %s.svg', filename, filename, filename))