# Gnuplot script plotting the Planck function for -20, 0 and 20 degC # and its integrals # It uses logscale x-axis, the y-axis being lambda/micrometer # times the spectral irradiance -- this is the recommended way in the # wavelength domain, see Marr & Wilkin 2012: # http://link.aip.org/link/AJPIAS/v80/i5/p399/s1&Agg=doi # The integrating part is taken from an example at # http://gnuplot.sourceforge.net/demo/bivariat.html # Jan Hollan, 2013. Titles in Czech language. set encoding iso_8859_2 set term postscript eps color enhanced font "Arial-Medium" set out "plan_sb_cz.eps" set size 0.7 set grid sigma = 5.67e-8 sfb(x) = sigma * (x+273.15)**4 Planck(x,t) = \ 1E-6*3.741E-16/(( exp(5*log((x*1E-6))) )*(exp(1.439E-2/(x*t*1E-6))-1)) abst(x)=x+273.15 set title \ "Záření černého tělesa \ - Planckovy funkce a jejich integrály \ \n (je vyznačeno pásmo LWIR od 7 {/Symbol m}m do 15 {/Symbol m}m)" set ylabel \ "({/Symbol l}/{/Symbol m}m) {/Symbol ×} \ spektrální intenzita vyzařování\n \ /(W{/Symbol ×}m^{-{/*0.8 2}}{/Symbol m}m^{-{/*0.8 1}})" set xlabel "ln({/Symbol l} / {/Symbol m}m)" set y2label \ "integrál spektrální intenzity vyzařování / (W{/Symbol ×}m^{-{/*0.8 2}})" set xrange [0.5:5.5] set yrange [0:450] set xtics 1 nomirror set x2label "{/Symbol l} / {/Symbol m}m" set x2tics ("2" log(2), "3" log(3), "5" log(5), "10" log(10), \ "20" log(20), "30" log(30), "50" log(50), "100" log(100), \ "200" log(200)) delta = 0.2 # delta can be set to 0.025 for non-MSDOS machines # # integral_f(x) takes one variable, the upper limit. 0 is the lower limit. # calculate the integral of function f(t) from 0 to x # choose a step size no larger than delta such that an integral number of # steps will cover the range of integration. integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta)) int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.) int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.) # # integral2_f(x,y) takes two variables; x is the lower limit, and y the upper. # calculate the integral of function f(t) from x to y integral2_f(x,y) = (xy-d*.5) ? 0 : (int2(x+d,y,d) + (f(x)+4*f(x+d*.5)+f(x+d))*d/6.) f(x)=exp(x)*Planck(exp(x),abst(0)) set samples 50 set key bottom right set multiplot set arrow from log(7), graph 0 to log(7), graph 1 nohead set arrow from log(15), graph 0 to log(15), graph 1 nohead plot \ exp(x)*Planck(exp(x),abst(20)) title '20 {/Symbol \260}C' lt 1 lw 5, \ exp(x)*Planck(exp(x),abst(0)) title '0 {/Symbol \260}C' lt 2 lw 5 , \ exp(x)*Planck(exp(x),abst(-20)) title '-20 {/Symbol \260}C' lt 3 lw 5 ,\ integral_f(x) title "" lt 2 lw 3 f(x)=exp(x)*Planck(exp(x),abst(-20)) plot integral_f(x) title "" lt 3 lw 3 f(x)=exp(x)*Planck(exp(x),abst(20)) plot integral_f(x) title "" lt 1 lw 3 unset multiplot quit