#!/bin/bash
# AzimuthalEquidistant reference transformations generated by cs2cs from PROJ.4
#
# Output lines should be included AzimuthalEquidistant.txt
#
# Requires cs2cs from PROJ.4, bash, perl, and head

set -e

# convert cs2cs output to source and target lines
reformat () {
   perl -p -e 's/^(\S+)\s+(\S+)\s+(\S+)\s+(\S+).*$/source pt = (\1, \2)\ntarget pt = (\3, \4)/' 
}

echo "// AzimuthalEquidistant reference transformations generated by cs2cs from PROJ.4 with \"AzimuthalEquidistant-cs2cs.sh\""
echo "// cs2cs version: `cs2cs 2>&1 | head -1`"
echo
echo "// Repeat the Snyder International ellipsoid examples in higher precision to compare with cs2cs"
echo "source crs = _INTL_"
echo "target crs = _INTL_NORTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=intl +to +proj=aeqd +ellps=intl +lon_0=0 +lat_0=90 +no_defs"
echo "// $CS2CS"
$CS2CS <<EOF | reformat
0 90
0 75
180 75
90 75
45 75
-45 75
135 75
-135 75
EOF
echo "source crs = _INTL_"
echo "target crs = _INTL_SOUTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=intl +to +proj=aeqd +ellps=intl +lon_0=0 +lat_0=-90 +no_defs"
echo "// $CS2CS"
$CS2CS <<EOF | reformat
0 -90
0 -75
180 -75
90 -75
45 -75
-45 -75
135 -75
-135 -75
EOF

echo
echo "// Test the same 16 points for each of the eight WGS84 and spherical projections"

# 16 random points generated in Python with:
#
# import random
# for i in range(16):
#     print random.randint(-180,180), random.randint(-90,90)
#
POINTS=`cat <<EOF
-146 -46
176 4
-123 37
20 -69
143 44
-3 4
-104 21
66 6
-101 -85
-8 -36
111 79
80 37
-31 -18
110 34
-139 -71
-171 -48
EOF`

echo
echo "source crs = _WGS84_"
echo "target crs = _WGS84_NORTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=WGS84 +to +proj=aeqd +ellps=WGS84 +lon_0=15 +lat_0=90 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _WGS84_"
echo "target crs = _WGS84_SOUTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=WGS84 +to +proj=aeqd +ellps=WGS84 +lon_0=-35 +lat_0=-90 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _WGS84_"
echo "target crs = _WGS84_EQT_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=WGS84 +to +proj=aeqd +ellps=WGS84 +lon_0=105 +lat_0=0 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _WGS84_"
echo "target crs = _WGS84_OBL_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=WGS84 +to +proj=aeqd +ellps=WGS84 +lon_0=-55 +lat_0=40 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _SPHERE_"
echo "target crs = _SPHERE_NORTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=sphere +to +proj=aeqd +ellps=sphere +lon_0=15 +lat_0=90 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _SPHERE_"
echo "target crs = _SPHERE_SOUTH_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=sphere +to +proj=aeqd +ellps=sphere +lon_0=-35 +lat_0=-90 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _SPHERE_"
echo "target crs = _SPHERE_EQT_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=sphere +to +proj=aeqd +ellps=sphere +lon_0=105 +lat_0=0 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat

echo
echo "source crs = _SPHERE_"
echo "target crs = _SPHERE_OBL_"
CS2CS="cs2cs -E -f %.6f +proj=longlat +ellps=sphere +to +proj=aeqd +ellps=sphere +lon_0=-55 +lat_0=40 +no_defs"
echo "// $CS2CS"
echo -n "$POINTS" | $CS2CS | reformat
