Re: adding illumination to hemisphere (-JG)

Hi Pasquale,

I created a script to do that some time ago, which I include below. The
script depends on the "bc" utility to handle the math, but most Unix/Linux
installations should have that.

The script actually loops through various viewpoints and subsolar points to
create frames that form a "spinning earth" .gif, but you can easily modify
it to choose just one frame or view.

Best wishes,

Keith Pickering

*****************************

#!/bin/sh
##############################################################
#            Animated turning globe with terminator
#            Author: Keith Pickering
#
##############################################################
# Usage: set latitude and longitude of subsolar point, plus
#        latitude and longitude of viewpoint.
# Longitudes of the two must be marginally different, or an
# error is thrown. Also, direct 180 opposite longitudes
# throws the same error. Viewpoints over North and South poles
# give warnings, but produce good output. All other rotations
# should work with no errors or warnings.
###############################################################

# Create file name in tif format
height=1.0i
width=1.0i
dpi=100
. gmt_shell_functions.sh
gmtset DOTS_PR_INCH $dpi PAPER_MEDIA=Custom_${width}x${height}
name=globe3
# Set starting frame number, total number of frames, and number of degrees
of rotation between each frame
frame=1
nframes=90
framesd=4
mkdir $$
# Standard variables, non-looping:
# Define latitide and longitude of subsolar point:
slon=164.2152
slat=15.6865
# Define latitude and longitude of view center:
vlon=-83
vlat=41
# Mapping options and colors
D="-Dc -A600/0/1"
day="-G110/150/120 -S33/33/150"
nite="-G33/45/36 -S9/9/45"
dr=.017453293
#
# Main loop starts here:
#
while [ $frame -le $nframes ]; do
   ps=$$/$(eval "gmt_set_framename $name $frame").ps
   vlon=$(echo "-81.5-($frame*$framesd) "|bc)
   slon=$(echo "sl=$vlon+241.2; if(sl<-180) sl+=360;sl"|bc)
   echo "file: " $ps "vlon: " $vlon "slon: " $slon
OPT="-Rd -JG"$vlon/$vlat"/.9i -P"

## Standard terminator script starts here
tlat=$(echo "scale=4; $slat-90"|bc)
vlaty=$(echo "scale=4; $vlat-90"|bc)
ld=$(bc << EOF
scale=4
ld=$slon-($vlon)
if (ld>180) ld-=360
if (ld<-180) ld+=360
ld
EOF
)
# Gibbous/crescent switch: 1=crescent
# New routine: compute spherical angle between viewpoint and subsolar point.
ldi=$(bc -l << EOF
scale=4
ld=$slon-($vlon)
cd=s($dr*$vlat)*s($dr*$slat)+c($dr*$vlat)*c($dr*$slat)*c($dr*ld)
if (cd < 0) ldi=1 else ldi=0
ldi
EOF
)
# Angle for -L option in clip path.
l1=$(bc -l << EOF
scale=4
sb=c($slat*$dr)*c($vlat*$dr)*c($ld*$dr)+s($slat*$dr)*s($vlat*$dr)
cb=sqrt(1-sb*sb)
b=a(sb/cb)
y=s($slat*$dr)-sb*s($vlat*$dr)
x=c($slat*$dr)*s($ld*$dr)*c($vlat*$dr)
res=a(y/x)
# resolve arctan ambiguity:
if (x<0) if (y<0) res-=(180*$dr) else res+=(180*$dr)
res/$dr
EOF
)

# generate terminator clipping path
project -G5 -S -T$slon/$slat -C$slon/$tlat -L0/360 > clip1.txt

if [ $ldi -eq 1 ]; then
# cresent: central meridian is unlit
   l2=$(echo "scale=4; $l1-180"|bc -l)
   # GC2. At right angles to the view, covering the lit half of the globe.
   project -G5 -S -T$vlon/$vlat -C$vlon/$vlaty -L$l2/$l1 > clip2.txt
   # day hemisphere: whole globe, will be mostly overlain with night color
scheme.
   pscoast $OPT $day $D -K > $ps
   # Terminator ellipse
   psclip clip1.txt  $OPT -K -O >> $ps
   pscoast  $OPT $nite $D -K -O >> $ps
   psclip -C -O -K >> $ps
   # Unlit half.
   psclip clip2.txt $OPT -K -O >> $ps
   pscoast  $OPT $nite $D -K -O >> $ps
   psclip -C -O >> $ps
else
# gibbous: central meridian is lit
   l2=$(echo "scale=4; $l1+180"|bc -l)
   # GC2. At right angles to the view, covering the lit half of the globe.
   project -G5 -S -T$vlon/$vlat -C$vlon/$vlaty -L$l1/$l2 > clip2.txt
   # night hemisphere: whole globe, will be mostly overlain with day color
scheme.
   pscoast $OPT $nite $D -K > $ps
   # Terminator ellipse
   psclip clip1.txt  $OPT -K -O >> $ps
   pscoast  $OPT $day $D -K -O >> $ps
   psclip -C -O -K >> $ps
   # Lit half.
   psclip clip2.txt $OPT -K -O >> $ps
   pscoast  $OPT $day $D -K -O >> $ps
   psclip -C -O >>$ps
fi
   ps2raster -E$dpi -Tt $ps
   frame=$(eval "gmt_set_framenext $frame")
done
echo "converting to gif, this may take a while ..."
   convert -delay 4 -loop 0 $$/*.tif $name.gif
echo $name".gif conversion, done, cleaning up ..."
   rm $$/*.ps
cat << EOF > $name.html
<HTML>
<TITLE>Rotating Globe</TITLE>
<BODY>
<img src="$name.gif" border="1">
</BODY>
</HTML>
EOF


On Tue, Mar 25, 2014 at 8:58 PM, Pasquale Tricarico <[email protected]>wrote:

> Hi,
>
> I am using GMT 4.5.12, using grdimage and -JG projection to obtain an
> hemisphere view of a planet. I would like to add illumination, i.e. by
> specifying
> the lat/lon of the sun, so that the correct part of the hemisphere is
> illuminated and the other part is much darker, but still visible. It
> looks like grdgradient may not be able to do this as it would have a
> fixed azimuth for illumination, maybe by using grdmath? Any
> suggestions will be greatly appreciated.
>
> Regards,
> Pasquale
>
> --
> Pasquale Tricarico, Ph.D.
> Planetary Science Institute
> http://orbit.psi.edu/~tricaric
>
> To unsubscribe, send the message "signoff gmt-help" to
> [email protected]
>

To unsubscribe, send the message "signoff gmt-help" to [email protected]