#!/bin/bash

#################################################################
# Script che calcola la distanza ortodromica su file path di	#
# google earth esportato in formato kml (ellissoide WGS84)	#
# Necessita:							#
#		bc basic calculator				#
#								#
# Autore:	Pierpaolo Garofalo				#
# Data:		29.08.2010					#
# Script originale java:					#
# http://www.movable-type.co.uk/scripts/latlong-vincenty.html	#
# Documentazione:						#
#	http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf		#
#################################################################


if [ "$#" != "1" ]; then

    echo
    echo "Error: Provide a kml google path file"
    echo
    exit 1
    
fi

export KMLFILE="$1"


if [ -f "$KMLFILE" ]; then

    :
    
else

    echo
    echo "Error: file $KMLFILE not found"
    echo
    exit 1
    
fi

#################################################################
# Verifica esistenza di bc					#
#################################################################


if which bc &> /dev/null; then

    :

else
    
    echo
    echo "Error: bc not found"
    echo
    exit 1
    
fi

#################################################################
# Dati geometrici ellissoide WGS-84				#
#################################################################

a="6378137"		# m (±2 m)
b="6356752.314245"	# m
rf="298.257223563"	# 1/f=rf

pi=`echo "4*a(1)" | bc -l`

MAXERR=0.00000000001

export TMPFILE="/tmp/coordinate.txt"

cat $KMLFILE | grep -A 1 "<coordinates>" | sed -e 's/ /\n/g' | grep -v "</coordinates>" | grep -v "<coordinates>"> $TMPFILE

export NPOINTS=$(wc -l ${TMPFILE} | awk '{print $1}')

echo "Numero punti " $NPOINTS

DISTANZA=0

exec 3< $TMPFILE

    read <&3
    
    P0_lon=$(echo ${REPLY} | cut -d ',' -f1)
    P0_lat=$(echo ${REPLY} | cut -d ',' -f2)
    
    while read <&3
    do

	P1_lon=$(echo ${REPLY} | cut -d ',' -f1)
	P1_lat=$(echo ${REPLY} | cut -d ',' -f2)
	
	L=$(echo "(${P1_lon} - ${P0_lon}) / 180 * ${pi}" | bc -l)
	lambda=$L
	
	U1=$(echo "a((1-1/${rf})*s(${P0_lat}/180.0*${pi})/c(${P0_lat}/180.0*${pi}))" | bc -l)
	U2=$(echo "a((1-1/${rf})*s(${P1_lat}/180.0*${pi})/c(${P1_lat}/180.0*${pi}))" | bc -l)
	
	sinU1=$(echo "s(${U1})" | bc -l)
	sinU2=$(echo "s(${U2})" | bc -l)
	cosU1=$(echo "c(${U1})" | bc -l)
	cosU2=$(echo "c(${U2})" | bc -l)
	
	MAXITER=100
	COUNTER=1
	r=0
	DD=0
	DDNULL=0
	
	until [[ ($COUNTER -gt $MAXITER) || ("$r" == "1") ]]; do
	
	
	    sinLambda="$(echo "s(${lambda})" | bc -l)"
	    cosLambda=$(echo "c(${lambda})" | bc -l)
	
	    sinSigma=$(echo "sqrt((${cosU2}*${sinLambda})^2+(${cosU1}*${sinU2} - ${sinU1}*${cosU2}*${cosLambda})^2)" | bc -l)
	
	    if [ "$sinSigma" == "0" ]; then
	    
		DD=0
		DDNULL=1
		break
		
	    fi
	
	    cosSigma=$(echo "${sinU1}*${sinU2} + ${cosU1}*${cosU2}*${cosLambda}" | bc -l)
	    
	    quad=$(echo "if(${cosSigma}>0) if(${sinSigma}>0) 0 else 360; if(${cosSigma}<0) 180" | bc -l)
	    
	    sigma=$(echo "a(${sinSigma} / ${cosSigma}) + ${quad}" | bc -l)
	    
	    sinAlfa=$(echo "${cosU1} * ${cosU2} * ${sinLambda} / ${sinSigma}" | bc -l)
	
	    cosSqAlfa=$(echo "1 - ${sinAlfa}^2" | bc -l)
	
	    if [ "$cosSqAlfa" == "0" ]; then
	    
		cos2SigmaM=0
		
	    else
	    
		cos2SigmaM=$(echo "${cosSigma} - 2 * ${sinU1} * ${sinU2} / ${cosSqAlfa}" | bc -l)
		
	    fi
	
	    C=$(echo "1/(${rf} * 16) * ${cosSqAlfa} * (4 + 1/${rf}) * (4 - 3 * ${cosSqAlfa})" | bc -l)
	
	    lambdaP=$lambda
	    
	    lambda=$(echo "${L} + (1 - ${C}) * 1 / ${rf} * ${sinAlfa} * (${sigma} + ${C} * ${sinSigma} * (${cos2SigmaM} + ${C} * ${cosSigma} * (-1 + 2 * ${cos2SigmaM}^2)))" | bc -l)
	
	    r=$(echo "e=sqrt((${lambda} - ${lambdaP})^2);if(e < ${MAXERR}) 1 else 0" | bc -l)
	    
	    let COUNTER+=1
	    
	done
	
	if [ $COUNTER -gt $MAXITER ]; then
	
	    echo
	    echo "Error: lambda calculation did not converge"
	    echo
	    exit 1
	    
	fi
	
	if [ $DDNULL == "0" ]; then
	
	    uSq=$(echo "${cosSqAlfa}*(${a}^2 - ${b}^2)/${b}^2" | bc -l)
	
	    A=$(echo "1+${uSq}/16384*(4096+${uSq}*(-768+${uSq}*(320-175*${uSq})))" | bc -l)
	
	    B=$(echo "${uSq}/1024*(256+${uSq}*(-128+${uSq}*(74-47*${uSq})))" | bc -l)
	
	    deltaSigma=$(echo "${B}*${sinSigma}*(${cos2SigmaM}+${B}/4*(${cosSigma}*(-1+2*${cos2SigmaM}^2)-${B}/6*${cos2SigmaM}*(-3+4*${sinSigma}^2)*(-3+4*${cos2SigmaM}^2)))" | bc -l)
	
	    DD=$(echo "${b}*${A}*(${sigma} - ${deltaSigma})" | bc -l)
	
	    DISTANZA=$(echo "$DISTANZA+$DD" | bc -l)
	    
	fi
	
	P0_lon="$P1_lon"
	P0_lat="$P1_lat"

    done
    
    echo "TOTAL DISTANCE: " $DISTANZA "m"

exec 3>&-