#!/bin/sh

# Markus Neteler 2007 + 2011, http://www.grassbook.org/
# Script to calculate global energy for a certain time period of a year

# from http://www.grassbook.org/examples_menu2nd.php
# updated to GRASS 6.4+

echo "Enter elevation map (e.g., elevation):"
read elev
g.copy rast=$elev,SOLelev

i=75
lastday=288

g.region rast=$elev

#generate an empty map for global radiation:
r.mapcalc "global.rad=0"

while [ $i -le $lastday ]
do
 # generate map names convenient for xganim and r.out.mpeg:
 DAY=`echo $i | awk '{printf "%03i", $1}'`

 echo "Computing radiation for day $DAY..."
 r.sun -s elevin=SOLelev day="$i"\
       beam_rad=b_rad.$DAY diff_rad=d_rad.$DAY\
                            refl_rad=r_rad.$DAY

 #add to (cell-wise) global energy:
 r.mapcalc "global.rad=global.rad + b_rad.$DAY +\
                         d_rad.$DAY + r_rad.$DAY"
 r.timestamp b_rad.$DAY date="$i days"
 r.colors b_rad.$DAY col=gyr
 r.timestamp d_rad.$DAY date="$i days"
 r.colors d_rad.$DAY col=gyr
 r.timestamp r_rad.$DAY date="$i days"
 r.colors r_rad.$DAY col=gyr

 i=`expr $i + 1`
done

#cleanup:
g.remove rast=SOLelev,SOLaspect,SOLslope
echo "Finished. Written <global.rad>"

exit 0