# 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 grid
set key bottom right
set width 10
set logscale x1y1x2x3
set title 'Simulated infrared dust spectrum for an \
{\mbox{\normalsize H\thinspace\footnotesize II}\kern3pt} region'

set xlabel '$\lambda$/$\mu$m'
set x2range [phy_c/(min*1e-6):phy_c/(max*1e-6)]
set x2label '$\nu$/Hz'
set x3range [phy_h*phy_c/(min*1e-6)/phy_ev:phy_h*phy_c/ \
(max*1e-6)/phy_ev]
set x3label 'Photon Energy / eV'

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}$'
# 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] 'example1.dat' using 1:2 \
t 'Nikoli\v{c}-Ford Dust Code' with lines, \
'example2.dat' t 'IRAS Photometry' using \
($1):(($2)/3e8*((($1)*1e-6)**2)*1.375191e+13/3.668333e+17) \
with points, \
f(x)/f(60)*1.375191e+13/(3e8/(60e-6**2)) t '$\beta=2$ Greybody' with lines