#!/bin/sh
# hemispheres -- divide and dominate the world
# Copyright       : http://www.fsf.org/copyleft/gpl.html
# Original Author : Benjamin Horner-Johnson <ben@rice.edu>
# New Commander   : Dan Jacobson -- http://jidanni.org/
# Created On      : 2002-MAY-11
# Last Modified On: Sat Nov 21 08:17:32 2020
# Update Count    : 77

## sh script to find great circle dividing globe into hemispheres between two
## points.  2002-MAY-11   Benjamin Horner-Johnson <ben@rice.edu>

# > I live in Taiwan, my little brother lives in Sweden.  I want to
# > divide up the world into hemispheres of control.  OK, first draw world
# > map with two great circles centered with axes passing through each of
# > us.

#gmtset PAPER_MEDIA B5+ BASEMAP_TYPE PLAIN, well, cropped anyway in convert(1)

test $# -ne 4 && { cat 1>&2 <<EOF
usage: $0 longitude1 latitude1 longitude2 latitude2
e.g. Taiwan Chicago: $0 120.87 24.18 -88 42
EOF
exit 4
}
set -ue
lon1=$1 lat1=$2 lon2=$3 lat2=$4
c1=`echo ${lat1} | awk '{if ($1 < 0) {print $1+90} else {print $1-90}}'`
c2=`echo ${lat2} | awk '{if ($1 < 0) {print $1+90} else {print $1-90}}'`

# generate great circles ("equators") around points 1 and 2
project -T${lon1}/${lat1} -C${lon1}/${c1} -G111.195 -Q -L-20015/20015 > gc1
project -T${lon2}/${lat2} -C${lon2}/${c2} -G111.195 -Q -L-20015/20015 > gc2

# > OK, now a great circle is needed to pass through the two common points of
# > those two great circles.  This is the final boundary for our hemispheres
# > of control for world dominance, etc. muhahahaha, whatever.

# find the two intersections
gmtselect gc1 -C0.75/gc2 > gcInt
pole1=`head -1 gcInt`
pole2=`tail -1 gcInt`
intlon1=`echo ${pole1} | awk '{print $1}'`
intlat1=`echo ${pole1} | awk '{print $2}'`
intlon2=`echo ${pole2} | awk '{print $1}'`
intlat2=`echo ${pole2} | awk '{print $2}'`

# > OK, time for a meeting with little brother.  Need to pinpoint best
# > place [midpoint of great circle through us], and worst place [other
# > midpoint of said circle, antipode of first.]

# Use that midpoint and one of the intersections to create the great circle
# dividing your hemispheres of domination.  (Otherwise you could live in the
# same hemisphere because there are an infinite number of great circles that
# intersect two points).

# find the midpoint between point 1 and point 2
project -C${lon1}/${lat1} -E${lon2}/${lat2} -G5 -Q > gcMid12
mid=`tail -1 gcMid12 | awk '{print $3/2}'`
midpoint=`gawk -v m=$mid '((m-$3)<2.6 && ($3-m)<2.6)' gcMid12`
lonm=`echo $midpoint | awk '{print $1}'`
latm=`echo $midpoint | awk '{print $2}'`
antilonm=`echo $midpoint | awk '{ p=$1+180; if (p>360) {print p-180} else {print p}}'`
antilatm=`echo $midpoint | awk '{print -1*$2}'`

# generate the border line for hemispheres of control
project -C${intlon1}/${intlat1} -E${lonm}/${latm} -G111.195 -Q -L-20015/20015 > gcBorder

###
### mapping
###
## part 1 border
B=hemispheres shade=220
O=$B.ps
pscoast -JG${lonm}/${latm}/4.5i -R-180/180/-90/90 -Bg15 -K -G$shade -Dc -Y2.0i -X0.75i > $O

psxy -JG -R -K -O gc1 -W5/255/0/0 >> $O
psxy -JG -R -K -O gc2 -W5/0/255/0 >> $O
psxy -JG -R -K -O gcBorder -W10/0/0/0 >> $O
psxy -JG -R -K -O -Sc0.1i -M -W1 <<END>> $O
> -G255/0/0
${lon1} ${lat1}
> -G0/255/0
${lon2} ${lat2}
> -G0
${lonm} ${latm}
END
pstext -JG -R -K -O -D0.1i/0.0i <<END>> $O
${lon1} ${lat1} 12 0 0 6 1
${lon2} ${lat2} 12 0 0 6 2
${lonm} ${latm} 12 0 0 5 Best meeting
END

## part 2 - antipodes
pscoast -JG${antilonm}/${antilatm}/4.5i -R-180/180/-90/90 -Bg15 \
    -K -O -G$shade -Dc -X4.6i >> $O
#to change separation of the two globes, adjust -X above.

psxy -JG -R -K -O gc1 -W5/255/0/0 >> $O
psxy -JG -R -K -O gc2 -W5/0/255/0 >> $O
psxy -JG -R -K -O gcBorder -W10/0/0/0 >> $O
psxy -JG -R -K -O -Sc0.1i -W1 -M <<END>> $O
> -G255/0/0
${lon1} ${lat1}
> -G0/255/0
${lon2} ${lat2}
> -G255
${antilonm} ${antilatm}
END
echo "${antilonm} ${antilatm} 12 0 0 5 Worst meeting"|pstext -JG -R -K -O -D0.1i/0.0i >> $O
psxy -JG -R -O /dev/null >> $O #we have made the .ps
rm gc1 gc2 gcBorder gcMid12 gcInt

# > Wait, it's mother's day and forgot about Mom, who lives in Chicago.
# > She wants to come to the meeting.  Need to show best point for all
# > three parties to meet.  And worst too.  Wait, mom wants a cut of the
# > pie... so just in case anybody has spare brain cells, lets extend this
# > to a general case where say one has a list of long,lat pairs of all
# > players whom we should divide up the world into equal sized areas of
# > control... [or never mind, just do the two player case.]

# Worst point would probably be the intersection of your worst point and the
# great circle 90 degrees away from Chicago.  Best point should be the center
# of a triangle (spherical) with your three locations as the vertices.  Don't
# have enough time for the general case.

#OK, thanks Ben for the programming!