#!/bin/bash

# create a graph of the histograms of common time observations from two wind stations.

# create the distributions for each stationlet "totsta = 0"

let "totname = 0"
#for name in midland big-spring lubbock roswell carlsbad hobbs clovis
for name in harrisburg_1 harrisburg_2
#for name in fort_ord_1 fort_ord_2
do
	let "totname = totname +1"
	staname[totname]=${name}
	#echo "station number ${totname} is ${staname[totname]}"
done

let "totsta = 0"
#for stanum in 722650 722657 722670 722680 722687 722688 722689
for stanum in 725110 725118
#for stanum in 690070 724916
do
	./ncdc2weru < wnd2/${stanum}.wnd2 > ${stanum}.org
	let "totsta = totsta +1"
	let "station[totsta] = stanum"
	echo "station number ${totsta} is ${station[totsta]}"
done

if [[ totname -ne totsta ]]
then
	echo "Count of names and station numbers doesn't match"
	echo "Please Check and try again"
	exit
fi

# create gnuplot scripts to plot one station against the other
plotfile="sidebyside.plt"
outfile="sidebyside.ps"

# header for text comparison table
echo "Station_1 #_points Energy Common_energy #_common_points Correlation Common_energy energy #_points Station_2"

# create datafiles for pairs of stations
for nsta1 in $( seq 1 ${totsta} )
do
	# compare two stations only once, and not station with self
	let "base2 = nsta1 + 1"
	for nsta2 in $( seq ${base2} ${totsta} )
	do
		# create the files with records recorded at the same time
		./commonhours -f${station[nsta1]}.org -g${station[nsta2]}.org -p${station[nsta1]}_${station[nsta2]}.com -q${station[nsta2]}_${station[nsta1]}.com
		# create histograms for original records
		./wind_mdb -x -w ${station[nsta1]} -n ${staname[nsta1]} < ${station[nsta1]}.org > ${station[nsta1]}.his
		./wind_mdb -x -w ${station[nsta2]} -n ${staname[nsta2]} < ${station[nsta2]}.org > ${station[nsta2]}.his
		# create histograms reflecting the records from the common times
		./wind_mdb -x -w ${station[nsta1]} -n ${staname[nsta1]}_${staname[nsta2]} < ${station[nsta1]}_${station[nsta2]}.com > ${station[nsta1]}_${station[nsta2]}.his
		./wind_mdb -x -w ${station[nsta2]} -n ${staname[nsta2]}_${staname[nsta1]} < ${station[nsta2]}_${station[nsta1]}.com > ${station[nsta2]}_${station[nsta1]}.his

		# graph the histograms for comparison
		./windhisto_files ${station[nsta1]}.his ${station[nsta1]}_${station[nsta2]}.his ${station[nsta2]}.his ${station[nsta2]}_${station[nsta1]}.his
		mv windplot.ps ${staname[nsta1]}_${staname[nsta2]}.ps

		#echo "processing ${staname[nsta1]} vs. ${staname[nsta2]}"
		# create the file with records recorded at the same time
		./sidebyside -f${station[nsta1]}.org -g${station[nsta2]}.org > ${station[nsta1]}_${station[nsta2]}.sbs

		# find stats from histograms
		# full station 1, station 1 common with 2, station 2 common with 1, station 2
		energy_1=`./stats_from_histo -t8 < ${station[nsta1]}.his | grep ${station[nsta1]} | dm x9`
		energy_12=`./stats_from_histo -t8 < ${station[nsta1]}_${station[nsta2]}.his | grep ${station[nsta1]} | dm x9`
		energy_21=`./stats_from_histo -t8 < ${station[nsta2]}_${station[nsta1]}.his | grep ${station[nsta2]} | dm x9`
		energy_2=`./stats_from_histo -t8 < ${station[nsta2]}.his | grep ${station[nsta2]} | dm x9`

		# find number of data points in original files
		numpts1=`cat ${station[nsta1]}.org | colex 6 | stats n`
		numpts2=`cat ${station[nsta2]}.org | colex 6 | stats n`

		# repeat for wind velocities > 8 m/s in at least one of the pairs
		# find the number of common time data points
		cat ${station[nsta1]}_${station[nsta2]}.sbs | colex 7 9 | dm "if x1 >= 8 or x2 >= 8 then INPUT else SKIP"  > temp.sbs
		numcompts=`cat temp.sbs | pair | grep Analysis | dm x3`
		# find the correlation between wind velocities
		corr_coef=`cat temp.sbs | pair | grep -A 1 Correlation | grep -v Correlation | dm x1`
		# find intercept
		intercept=`cat temp.sbs | pair | grep -A 1 Intercept | grep -v Intercept | dm x1`
		# find slope
		slope=`cat temp.sbs | pair | grep -A 1 Slope | grep -v Slope | dm x2`

		echo "# GNUPLOT commands to plot wind station data records" > ${plotfile}
		echo "set terminal postscript color font 'Helvetica,5'" >> ${plotfile}
		echo "set output '${outfile}'" >> ${plotfile}
		echo "set title 'For ${numcompts} points, correlation = ${corr_coef}'" >> ${plotfile}
		echo "set xlabel '${staname[nsta1]} energy = ${energy_1}, Common energy = ${energy_12}'" >> ${plotfile}
		echo "set ylabel '${staname[nsta2]} energy = ${energy_2}, Common energy = ${energy_21}'" >> ${plotfile}
		echo "plot \\" >> ${plotfile}
		echo "'temp.sbs' using 1:2 with points notitle, \\" >> ${plotfile}
		echo "${intercept}+${slope}*x title 'trend', \\" >> ${plotfile}
		echo "x title '1:1'" >> ${plotfile}

		gnuplot ${plotfile}

		mv ${outfile} ${staname[nsta1]}_x8_${staname[nsta2]}.ps

		# text output for summary table
		echo "${staname[nsta1]} ${numpts1} ${energy_1} ${energy_12} ${numcompts} ${corr_coef} ${energy_21} ${energy_2} ${numpts2} ${staname[nsta2]}"

	done
done