# Example 3
# A rather more complicated plot to show off multiple axes,
# and LaTeXed labels on plots.
reset
# A few physical constants
min = 5
max = 200
phy_h = 6.626068e-34
phy_c = 3e8
phy_ev = 1.6e-19

# Set up plot basics...
set output 'examples/eps/example3.eps'
set terminal eps monochrome
set grid
set key bottom right
set width 10
set log x
set log y
set title 'Simulated infrared dust spectrum for an \
{\mbox{\normalsize H\thinspace\footnotesize II}\kern3pt} region'

# X-axis is wavelength, lambda
set xlabel '$\lambda$/$\mu$m'

# Y-axis is emitted flux, integrated over grainsize a
set ylabel '$\int F_{\nu}(a)\mathrm{d}a \cdot 4\pi r^2 / \
\mathrm{W} \, \mathrm{Hz}^{-1}\, \mathrm{m}^2 \, \
\mathrm{H}^{-1}$'

# Make a second X-axis, in units of frequency, nu
set x2range [phy_c/(min*1e-6):phy_c/(max*1e-6)]
set log x2
set x2label '$\nu$/Hz'

# And a third X-axis, in units of photon energy, in eV
set x3range [phy_h*phy_c/(min*1e-6)/phy_ev:phy_h*phy_c/ \
(max*1e-6)/phy_ev]
set log x3
set x3label 'Photon Energy / eV'

# Put an arrow and label on our plot, labelling one
# of the lines
set arrow 1 from 60, 2e-5 to 38, 1e-5
set label 1 "$F_\nu=\nu^{\beta}B_\nu(30\mathrm{K})$" \
at 62, 1.8e-5

# Make f(x) a 30K greybody
T=30.0
h=6.626e-34
k=1.38e-23
c=3e8
f(x)=((c/(x*1e-6))**(3+2))/(exp(h*c/(x*1e-6*k*T))-1.0)

# Finally, plot all of our data
plot [min:max][1e-7:1e-3] 'examples/example2a.dat' using 1:2 \
t 'Nikoli\v{c}-Ford Dust Code' with lines, \
'examples/example2b.dat' t 'IRAS Photometry' using \
($1):(($2)/3e8*((($1)*1e-6)**2)*1.375191e+13/3.668333e+17), \
f(x)/f(60)*1.375191e+13/(3e8/(60e-6**2)) t '$\beta=2$ Greybody'

# Produce a gif copy
set term gif
set dpi 168
set output 'examples/eps/example3.gif'
refresh
