/* This AML program was used for work published as: Kumar, L., Skidmore, A.K., /* Knowles, E., (1997) Modelling Topographic Variation in Solar Radiation in a /* GIS Environment. International Journal of Geographical Information Science, /* 11(5): 475-497. /* This program calculates the shortwave radiation received at the /* surface of the earth over a period of time. /* For the given day(s), it calculates the sunset and sunrise times /* and intergrates solar radiation from sunrise to sunset each day. /* /* Written by: Lalit Kumar /* School of Geography /* University of New South Wales /* Sydney 2052 /* Australia /* /* Last Modified: 5/08/1997 /* /* Send comments/suggestions to p2114659@geog.unsw.edu.au /* &args cover outgrid latdeg daystart dayend timeint &if [show program] ne GRID &then grid &s cover = [response 'ENTER NAME of DEM OF SITE'] &if ^ [exists %cover% -grid] &then &return &error Grid %cover% does not exist. &s outgrid = [response 'ENTER OUTGRID NAME'] &if [exist %outgrid% -grid] &then &return &error Grid %outgrid% already exists. &s latdeg = [response 'ENTER LATITUDE IN DEGREES'] &if ^ [variable latdeg] &then &return &error need a valid latitude &s daystart = [response 'ENTER JULIAN DATE FOR DAY TO START CALCULATION'] &if ^ [variable daystart] or %daystart% < 1 or %daystart% > 365 &then &return &error please enter a start day between 1 and 365 &s dayend = [response 'ENTER JULIAN DATE FOR DAY TO END CALCULATION'] &if ^ [variable dayend] or %dayend% < 1 or %dayend% > 365 &then &return &error please enter an end day between 1 and 365 &s timeint = [response 'ENTER TIME INTERVAL IN MINUTES'] &if ^ [variable timeint] or %timeint% <= 0 &then &return &error please enter a time greater than zero &s time = %timeint% * 0.0043633 /* convert minutes into radians &s daynumber = %daystart% &s pi = 22 / 7 &s degtorad = 2 * %pi% / 360 &s latitude = %latdeg% * %degtorad% &if [exist slopegrid -grid] &then kill slopegrid all slopegrid = slope(%cover%) &if [exist slopegrid2 -grid] &then kill slopegrid2 all slopegrid2 = slopegrid * %degtorad% &if [exist aspectgrid -grid] &then kill aspectgrid all aspectgrid = aspect(%cover%) &if [exist aspectgrid1 -grid] &then kill aspectgrid1 all /* change so that surface azimuth is zero at south if (aspectgrid == -1) aspectgrid1 = 0 else if (aspectgrid le 180) aspectgrid1 = 180 - aspectgrid else aspectgrid1 = 540 - aspectgrid &if [exist aspectgrid2 -grid] &then kill aspectgrid2 all aspectgrid2 = aspectgrid1 * %degtorad% &if [exist outgrid0 -grid] &then kill outgrid0 all outgrid0 = hillshade (%cover%, 0, 0, shadow) &if [exist initialgrid1 -grid] &then kill initialgrid all initialgrid = outgrid0 * 0 kill outgrid0 all kill slopegrid all kill aspectgrid all kill aspectgrid1 all &if [exist sungrid -grid] &then kill sungrid all &if [exist wattsgrid -grid] &then kill wattsgrid all &if [exist wattsgrid2 -grid] &then kill wattsgrid2 all &if [exist shaded -grid] &then kill shaded all &if [exist cosi -grid] &then kill cosi all &do &until %daynumber% > %dayend% /* See Eq 7 &s Io = [calc 1.367 * ( 1 + 0.034 * [cos [calc 360 * %degtorad% * %daynumber% / 365]] )] /* See Eq 4 &s decl = [calc 23.45 * %degtorad% * [sin [calc %degtorad% * 360 * ( 284 + %daynumber% ) / 365]]] /* See Eq 6 &s sunrise = [acos [calc -1 * [tan %latitude%] * [tan %decl%]]] &s sunset = -1 * %sunrise% &s hourangle = %sunrise% - ( %time% / 2 ) /* To ensure calculations are at half the time interval &s pass = 1 &do &until %hourangle% < %sunset% &ty day : %daynumber% &ty pass: %pass% /* See Eq 2 &s solaralt = [asin [calc [sin %latitude%] * [sin %decl%] + [cos %latitude%] * [cos %decl%] * [cos %hourangle%]]] &s test = [tan %decl%] / [tan %latitude%] /* See Eq 3 &if ( [cos %hourangle%] > %test% ) &then &s solaraz = [asin [calc [cos %decl%] * [sin %hourangle%] / [cos %solaralt%]]] &else &if ( [cos %hourangle%] < %test% ) &then &s solaraz = %pi% - [asin [calc [cos %decl%] * [sin %hourangle%] / [cos %solaralt%]]] &else &if ( %test% = [cos %hourangle%] and %hourangle% >= 0 ) &then &s solaraz = %pi% / 2 &else &if ( %test% = [cos %hourangle%] and %hourangle% < 0 ) &then &s solaraz = -1 * %pi% / 2 &end /* don't ask :)) &if ( %solaraz% >= 0 ) &then &s solarazdeg = %solaraz% * 57.29578 &else &s solarazdeg = 360 - ( [abs %solaraz%] * 57.29578 ) /* See Eq 11 /* Below is the original code, which doesn't run under AI 8.0.2 (no raise to fractional powers) /*&s M = [calc ( 1229 + ( 614 * [sin %solaralt%] ) ** 2 ) ** 0.5 - 614 * [sin %solaralt%]] /* Calvin Tolkamp proposed the workaround below: &s M = [calc [sqrt [calc 1229 + ( 614 * [sin %solaralt%] ) ** 2 ]] - 614 * [sin %solaralt%]] &s e = 2.7182818 /* See Eq 15 /* Below is the original code, which doesn't run under AI 8.0.2 (no raise to fractional powers) /*&s Is = %Io% * 0.56 * ( %e% ** ( -0.65 * %M% ) + %e% ** ( -0.095 * %M% ) ) /* Calvin Tolkamp proposed the workaround below: &s Is = %Io% * 0.56 * ( [exp [calc -0.65 * %M% ]] + [exp [calc -0.095 * %M% ]] ) &s solaraltdeg = %solaralt% * 57.29578 sungrid = HILLSHADE ( %cover%, %solarazdeg%, %solaraltdeg%, shadow ) /* See Eq 16 cosi = sin(%decl%) * (sin(%latitude%) * cos(slopegrid2) - cos(%latitude%) * sin(slopegrid2) * cos(aspectgrid2)) + cos(%decl%) * cos(%hourangle%) * (cos(%latitude%) * cos(slopegrid2) + sin(%latitude%) * sin(slopegrid2) * cos(aspectgrid2) ) + cos(%decl%) * sin(slopegrid2) * sin(aspectgrid2) * sin(%hourangle%) /* Trust me on the following two lines (You always need to fudge a little). if (cosi lt 0) shaded = 0 else shaded = 1 wattsgrid = %Is% * cosi * sungrid * shaded * 60 * %timeint% wattsgrid2 = wattsgrid + initialgrid kill wattsgrid all kill initialgrid all kill sungrid all kill shaded all kill cosi all rename wattsgrid2 initialgrid &s hourangle = %hourangle% - %time% &s pass = %pass% + 1 &end &s daynumber = %daynumber% + 1 &end %outgrid% = int(initialgrid) /* NOTE : THE UNITS FOR RADIATION IN THE OUTGRID ARE kJ/m^2/timeperiod /* If these are to be converted to MJ/m^2/day, then they should be divided by 1000 * no. of days. kill initialgrid all kill aspectgrid2 all kill slopegrid2 all