#!/usr/bin/gnuplot
# These mechanisms have been extensively studied, empirically quantified, and codified into an
# international standard for calculation: ANSI Standard S1-26:1995, or ISO 9613-1:1996.
T0=293.15;
set decimalsign ',';
#Directivitx index in dB (0=Freiraum, 3=über perfekt reflektierender Oberfläche)
#DI=2
# Sauerstoff relaxationsfrequenz abhängig von der prozentualen molaren Luftfeuchte h
FRelO(h) = 24.0 + 40400.0*h*(h+0.02)/(h+0.391)
# molare Feuchte zu rel Feuchte: h_rel = h_mol *42,7966. (2,33663 %mol = 100 %rel  )
# Stickstoffrelaxationsfrquenz
FRelN(h,T)=(T/T0)**(-0.5) * (9 + 280*h*exp(-4.17* ((T/T0)**(-1.0/3.0)-1)))
alpha(f,h,T)=869 * f**2.0 * \
  ( 1.84*10**(-11) * (T/T0)**(0.5) + \
      (T/T0)**(-2.5) * (\
         0.01275*exp(-2239.1/T)  /  (FRelO(h)+f**2.0/(FRelO(h))) +\
         0.10680*exp(-3352.0/T)  /  (FRelN(h,T)+f**2.0/(FRelN(h,T)))\
                       ) \
  )
Aabs(s,f,h,T)=alpha(f,h,T)*s/100        # absolute Dämpfung in DB/km  abhängig von der Entfernung s

#
set terminal pdfcairo size 21cm,29.7cm
set output 'AtmosphericSoundDamping.pdf'
#set terminal cairolatex standalone pdf size 21cm,29.7cm
#set output 'AtmosphericSoundDamping.tex'
set size 0.95,0.95
set origin 0.03,0.03
set title "Atmosphärische Schallabsorption abhängig von Temperatur und Luftfeuchte"
set style line 51 lt 1 dt 3 lc rgb 'blue'
set logscale x
set logscale y
xmin=100
xmax=1e6
set xrange [xmin:xmax]
set xlabel "Frequenz in Hz"
set ylabel "Dämpfung in dB / 100m"
ymin=0.01
ymax=10000
set yrange [ymin:ymax]
set for [j=0:5] for [i=2:8] ytics add (ymin*i*10**j ymin*i*10**j)       # Setzen von tics in logarithmischen Plots ist etwas umständlich da nicht einfach faktoren angegebenwerden können
set for [j=0:3] for [i in "2 5"] xtics add (xmin*i*10**j xmin*i*10**j)
set format x "%.0s %cHz"                      # Zahlendarstellung mit Si Prefixen, sonst werden die Zahlen sehr breit
set grid xtics ytics mxtics mytics ls 51, ls 0     # andere Farben für die Gitterlinien
set key top left                                  # Legende auf die linke Seite
plot for [kh=0:10] h=kh*10.0/42.7996 T=T0+20.0 alpha(x,h,T) title sprintf("h_{rel}=%d\%",kh*10.0)
# ----------------------
set title "Atmosphärische Schalldämpfung abhängig von Temperatur und Luftfeuchte"
s=1.0
set ytics 1
set ylabel sprintf("Dämpfung in dB / %d m",s)
unset logscale y
set mytics 10
set key top left
set multiplot
set xrange [100:100000]
set yrange [-5:0]
plot for [kh=0:10] h=kh*10.0/42.7996 T=T0+20.0 (-1.0)*Aabs(s,x,h,T) title sprintf("h_{rel}=%d\%",kh*10.0)
#--------- gestapelte Plots -----
set size 0.7,0.5
set origin 0.1,0.1
set xrange [100000:1000000]
set yrange [-140:0]
set ytics 10
unset ylabel
unset xlabel
unset title
unset key
#    Der gestapelte Plot soll eine Hintergrundfarbe bekommen, sonst ist er transparent
set object 1 rectangle from screen 0.1,0.1 to screen 0.8,0.6 behind fillcolor rgb "#22ffffaa"  fs noborder
set object 2 rectangle from graph 0, graph 0 to graph 1, graph 1 behind fillcolor rgb "#22ffffff" fs noborder
plot for [kh=0:10] h=kh*10.0/42.7996 T=T0+20.0 (-1.0)*Aabs(s,x,h,T) title sprintf("h_{rel}=%d\%",kh*10.0)
unset multiplot
system 
#h=0.7
#T=300.0
#s=100 
#plot Aabs(s,x,h,T) title sprintf("s=%d m, h=%d\%",s,h)
#set print
#system "pdflatex AtmosphericSoundDamping.tex "