Re: adding illumination to hemisphere (-JG)

[email protected]
Am 26.03.2014 16:16 schrieb "keithpickering ." <[email protected]>:

> 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]

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