#! /bin/bash


f1='archivio_filtrato_01.dat'
f2='archivio_filtrato_2b.dat'
f3='archivio_filtrato_2c.dat'


stat=`cat $f1 | grep -v '#' | awk 'BEGIN{i=0; c1=0; c2=0; c11=0; c12=0} \
{i+=1; c1+=$1; c2+=$1*$1; c11+=$2; c12+=$2*$2} \
END {printf "%3d %10.4f %10.4f %10.4f %10.4f",i,c1/i,sqrt((c2/i)-((c1/i)*(c1/i))),c11/i,sqrt((c12/i)-((c11/i)*(c11/i)))}'`
npop01=`echo $stat | awk '{print $1}'`
avg_u01=`echo $stat | awk '{print $2}'`
sd_u01=`echo $stat | awk '{print $3}'`
avg_d01=`echo $stat | awk '{print $4}'`
sd_d01=`echo $stat | awk '{print $5}'`
echo " Meteo station: 01"
echo "  SIZE = "$npop01
echo "   AVG U (m/s) = "$avg_u01
echo "ST.DEV U (m/s) = "$sd_u01
echo "   AVG D (degN)= "$avg_d01
echo "ST.DEV D (degN)= "$sd_d01
echo "=============================="

stat=`cat $f2 | grep -v '#' | awk 'BEGIN{i=0; c1=0; c2=0; c11=0; c12=0} \
{i+=1; c1+=$1; c2+=$1*$1; c11+=$2; c12+=$2*$2} \
END {printf "%3d %10.4f %10.4f %10.4f %10.4f",i,c1/i,sqrt((c2/i)-((c1/i)*(c1/i))),c11/i,sqrt((c12/i)-((c11/i)*(c11/i)))}'`
npop2b=`echo $stat | awk '{print $1}'`
avg_u2b=`echo $stat | awk '{print $2}'`
sd_u2b=`echo $stat | awk '{print $3}'`
avg_d2b=`echo $stat | awk '{print $4}'`
sd_d2b=`echo $stat | awk '{print $5}'`
echo " Meteo station: 2b"
echo "  SIZE = "$npop2b
echo "   AVG U (m/s) = "$avg_u2b
echo "ST.DEV U (m/s) = "$sd_u2b
echo "   AVG D (degN)= "$avg_d2b
echo "ST.DEV D (degN)= "$sd_d2b
echo "=============================="

stat=`cat $f3 | grep -v '#' | awk 'BEGIN{i=0; c1=0; c2=0; c11=0; c12=0} \
{i+=1; c1+=$1; c2+=$1*$1; c11+=$2; c12+=$2*$2} \
END {printf "%3d %10.4f %10.4f %10.4f %10.4f",i,c1/i,sqrt((c2/i)-((c1/i)*(c1/i))),c11/i,sqrt((c12/i)-((c11/i)*(c11/i)))}'`
npop2c=`echo $stat | awk '{print $1}'`
avg_u2c=`echo $stat | awk '{print $2}'`
sd_u2c=`echo $stat | awk '{print $3}'`
avg_d2c=`echo $stat | awk '{print $4}'`
sd_d2c=`echo $stat | awk '{print $5}'`
echo " Meteo station: 2c"
echo "  SIZE = "$npop2c
echo "   AVG U (m/s) = "$avg_u2c
echo "ST.DEV U (m/s) = "$sd_u2c
echo "   AVG D (degN)= "$avg_d2c
echo "ST.DEV D (degN)= "$sd_d2c
echo "=============================="
of="histo_windspeed.png"
of2="histo_directions.png"


cat > plotta << EOF

set zero 1.E-30
set grid
set border 31 lw 2
set style data lines
set style line 1 lt 1 lw 4 lc rgb "blue" pt 1 ps 1
set style line 2 lt 1 lw 4 lc rgb "cyan" pt 1 ps 1
set style line 3 lt 1 lw 4 lc rgb "green" pt 1 ps 1
set style line 4 lt 1 lw 4 lc rgb "yellow" pt 1 ps 1
set style line 5 lt 1 lw 4 lc rgb "orange" pt 1 ps 1
set style line 6 lt 1 lw 4 lc rgb "red" pt 1 ps 1
set style line 8 lt 8 lw 4 pt 1 ps 1
set style line 10 lt 0 lw 3 lc rgb "black" pt 1 ps 1
set style line 11 lt -1 lw 3 lc rgb "gray" pt 1 ps 1

set terminal png large size 800,800
set out '$of'

set xlabel 'Wind speed (m/s)'
set ylabel 'Relative frequency'
set xrange [-0.5:5.0]
set yrange [0:1]

#set xdata time
#set timefmt "%d-%b-%Y-%H:%M:%S"
#
# method 2:
#
# to make a histogram with *area* of a bin equal to the count
# in a bin, so the area under the curve is the number of data:
#
# Reference:
# http://www.inference.phy.cam.ac.uk/teaching/comput/C++/examples/gnuplot/index.shtml#four

# Creiamo l'istogramma delle velocita'
set title 'Wind speed histogram'
set label 1 "SIZE01 = $npop01"           at graph 0.6,0.80 
set label 2 "AVG u01 = $avg_u01"        at graph 0.6,0.77
set label 3 "ST.DEV u01 = $sd_u01"      at graph 0.6,0.74

set label 4 "SIZE2b = $npop2b"           at graph 0.6,0.71
set label 5 "AVG u2b = $avg_u2b"        at graph 0.6,0.68
set label 6 "ST.DEV u2b = $sd_u2b"      at graph 0.6,0.65

set label 7 "SIZE2c = $npop2c"           at graph 0.6,0.62
set label 8 "AVG u2c = $avg_u2c"        at graph 0.6,0.59
set label 9 "ST.DEV u2c = $sd_u2c"      at graph 0.6,0.56

bin_width = 0.5
bin_number(x) = floor(x/bin_width)
rounded(x) = bin_width * ( bin_number(x) + 0.5 )
UNITY = 1 #/ $npop
plot '$f1' u (rounded(\$1)):((UNITY)/$npop01.0) t 'Station 01' smooth frequency w histeps ls 1, \\
 '$f2' u (rounded(\$1)):((UNITY)/$npop2b.0) t 'Station 2b' smooth frequency w histeps ls 3, \\
 '$f3' u (rounded(\$1)):((UNITY)/$npop2c.0) t 'Station 2c' smooth frequency w histeps ls 6

#Creiamo l'istogramma delle direzioni

set out'$of2'
set title 'Wind direction histogram'
set xlabel 'Wind direction (degrees)'
set ylabel 'Relative frequency'
set xrange [0:360]
set yrange [0:0.5]

set label 1 "SIZE01 = $npop01"           at graph 0.6,0.80
set label 2 "AVG d01 = $avg_d01"        at graph 0.6,0.77
set label 3 "ST.DEV d01 = $sd_d01"      at graph 0.6,0.74

set label 4 "SIZE2b = $npop2b"           at graph 0.6,0.71
set label 5 "AVG d2b = $avg_d2b"        at graph 0.6,0.68
set label 6 "ST.DEV d2b = $sd_d2b"      at graph 0.6,0.65

set label 7 "SIZE2c = $npop2c"           at graph 0.6,0.62
set label 8 "AVG d2c = $avg_d2c"        at graph 0.6,0.59
set label 9 "ST.DEV d2c = $sd_d2c"      at graph 0.6,0.56

bin_width = 22.5
bin_number(x) = floor(x/bin_width)
rounded(x) = bin_width * ( bin_number(x) + 0.5 )
UNITY = 1 #/ $npop
plot '$f1' u (rounded(\$2)):((UNITY)/$npop01.0) t 'Station 01' smooth frequency w histeps ls 1, \\
 '$f2' u (rounded(\$2)):((UNITY)/$npop2b.0) t 'Station 2b' smooth frequency w histeps ls 3, \\
 '$f3' u (rounded(\$2)):((UNITY)/$npop2c.0) t 'Station 2c' smooth frequency w histeps ls 6



#pause -1 'ready 1'



EOF


gnuplot plotta
rm -f plotta 


exit 0









