Your equation is fine.
Here is a working implementation on bash.
# /bin/bash
DEG_PER_RAD=$(echo "scale=10; 180 / 3.141592653589" | bc)
function atan2 {
if [ $(echo "scale=10; $1 > 0" | bc -l) -eq 1 ]
then
if [ $(echo "scale=10; $2 > 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; a($1 / $2)" | bc -l)
fi
if [ $(echo "scale=10; $2 < 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; 3.141592653589 - a(-1 * $1 / $2)" | bc -l)
fi
if [ $(echo "scale=10; $2 == 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; 3.141592653589 / 2" | bc -l)
fi
fi
if [ $(echo "scale=10; $1 < 0" | bc -l) -eq 1 ]
then
if [ $(echo "scale=10; $2 > 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; -1 * a(-1 * $1 / $2)" | bc -l)
fi
if [ $(echo "scale=10; $2 < 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; a($1 / $2) - 3.141592653589" | bc -l)
fi
if [ $(echo "scale=10; $2 == 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; 3 * 3.141592653589 / 2" | bc -l)
fi
fi
if [ $(echo "scale=10; $1 == 0" | bc -l) -eq 1 ]
then
if [ $(echo "scale=10; $2 > 0" | bc -l) -eq 1 ]
then
a2=0
fi
if [ $(echo "scale=10; $2 < 0" | bc -l) -eq 1 ]
then
a2=$(echo "scale=10; 3.141592653589" | bc -l)
fi
if [ $(echo "scale=10; $2 == 0" | bc -l) -eq 1 ]
then
a2=0
fi
fi
}
function get_bearing {
rad_lat1=$(echo "scale=9; ("$1" / $DEG_PER_RAD)" | bc)
rad_lon1=$(echo "scale=9; ("$2" / $DEG_PER_RAD)" | bc)
rad_lat2=$(echo "scale=9; ("$3" / $DEG_PER_RAD)" | bc)
rad_lon2=$(echo "scale=9; ("$4" / $DEG_PER_RAD)" | bc)
dLon=$(echo "scale=10; $rad_lon2 - $rad_lon1" | bc -l)
y=$(echo "scale=10; s($dLon) * c($rad_lat2)" | bc -l)
x=$(echo "scale=10; (c($rad_lat1) * s($rad_lat2)) - (s($rad_lat1) * c($rad_lat2) * c($dLon))" | bc -l)
atan2 $y $x
bearing=$(echo "scale=10; $DEG_PER_RAD * $a2" | bc -l)
}
get_bearing 49.624196795 6.1256280094 49.6241653669 6.1242621755
echo "Bearing is "$bearing
Regarding PLC I have no experiance so I can't help you. Sorry