#!/bin/sh
# $Id: shaded_globe,v 1.4 2005-02-16 09:46:15-07 braup Exp braup $
# Makes a shaded globe

originlon="-105.3"     # Boulder
originlat="40"         # Boulder
#originlon="-85.3"
#originlat="30"
#origin="-150 61.16"    # Anchorage
#origin="0 -90"          # South Pole
#origin="-75.867 0.12"    # somewhere near equator


sunposition="-130 45"
#sunposition="-170 35"
#sunposition=$origin

if [ $# -eq 2 ] ; then
  originlon="$1"
  originlat="$2"
  sunlon=`echo "$1 - 30" | bc`
  sunlat=`echo "$2 + 5"  | bc`
#  sunposition="$sunlon $sunlat"
  sunposition="$originlon $originlat"
fi

origin="$originlon $originlat"

proj=G$originlon/$originlat/5c
region=-180/180/-90/90
inc=1
grdfile=shade.grd
out=globe.ps
pngout=globe.png

echo "Origin = $origin" >&2
bndpen=0.5p/0/0/0
gmtset PAPER_MEDIA letter BASEMAP_TYPE plain FRAME_PEN 1p HEADER_FONT_SIZE 20p
gmtset D_FORMAT %lg DEGREE_FORMAT 1
grdmath -R$region -I$inc -F $sunposition GDIST COSD = $grdfile
grdmath $grdfile 0.2 ADD = land.grd
grd2cpt $grdfile -S-0.2/1.2/0.05 -Cgray -Z > shade.cpt
grdimage $grdfile -R$region -J$proj -Cshade.cpt -K -V > $out
pscoast -R$region -J$proj -A3000 -Dc -Gc -O -K -V >> $out
grdimage land.grd -R$region -J$proj -Cshade.cpt -O -K -V >> $out
pscoast -R$region -J$proj -Q -O -K -V >> $out
#pscoast -R$region -J$proj -A3000 -Dc -N1/$bndpen -W$bndpen -O -K -V >> $out
echo "$origin" | psxy -J$proj -R$region -Sc0.2c -G0/0/255 -O -V >> $out

convert -density 80 -rotate 90 -trim $out $pngout