/* AML to do topographic shielding for cosmogenic-isotope production
/* calculations at a point. 

/* You need a DEM of the area, preferably in UTM coordinates, x,y,and z
/* denominated in meters. The name is specified as an argument.

&args dem

/* run this at the arc prompt. First go to grid.

grid

/* the next command is our shortcut to open a display window. Fix as needed.

d9

/* first, enter a point at which to determine shielding.  

mape %dem%
image %dem%

&ty 'Input point at which to calculate shielding...'

&getpoint &map &current
&sv currentx %pnt$x%
&sv currenty %pnt$y%
&lv currentx
&lv currenty

/* find the elevation of this point

&sv currentz [show cellvalue %dem% %currentx% %currenty%]

/* start shielding calcs. First, develop x and y coordinate grids

&describe %dem%
setcell %dem%
setwindow %dem%

&ty 'Making coordinate grids...'
&if [exists ygrid -grid] &then kill ygrid all
ygrid = (%grd$ymax% - (%grd$dy% / 2)) - ( $$rowmap * %grd$dy% )
&if [exists xgrid -grid] &then kill xgrid all
xgrid = (%grd$xmin% + (%grd$dx% / 2)) + ( $$colmap * %grd$dx% )

/* create grid showing r, i.e. distance to current point from each other cell

&ty 'Making r and theta grids...'
&if [exists rgrid -grid] &then kill rgrid all
rgrid = sqrt(sqr(xgrid - %currentx%) + sqr(ygrid - %currenty%))

/* use this to create azimuth and angular-elevation-of-each-point grids

&if [exists azgrid -grid] &then kill azgrid all
azgrid = con((xgrid - %currentx%) GE 0, acos((ygrid - %currenty%) / rgrid), 2 * pi - acos((ygrid - %currenty%) / rgrid))

&ty 'Computing vertical angles...'
&if [exists elgrid -grid] &then kill elgrid all
elgrid = atan((%dem% - %currentz%) / rgrid)

/* bin azimuths every 2 degrees, which is massive overkill

&ty 'Binning by azimuth...'
&if [exists azbin -grid] &then kill azbin all
azbin = reclass(azgrid, azimuth.rc)

/* apply integration formula 

&ty 'Calculating shielding factors...'
&sv bind 2
&sv binr [calc ( %bind% * 2 * 3.1416 / 360 ) ]

&if [exists s_factor -grid] &then kill s_factor all
s_factor = con(elgrid > 0,(( %binr% / ( 2 * pi ) ) * pow(sin(elgrid),3.3)),0)

/* apply zonalstats to obtain maximum shielding factor in each azimuth bin

&ty 'Calculating total shielding factor...'
&if [exists horiz1.dat -info] &then kill horiz1 info

horiz1.dat = zonalstats(azbin,int(s_factor * 100000),max)

/* calculate total shielding factor in tables

q
tables

sel horiz1.dat
&sv total 0
&do i = 1 &to 180
      &s  a [show RECORD %i%  max]
      &s  total [calc ( %total% + %a% ) ] 
&end 

&s total [calc %total% / 100000 ]
&ty 'Total shielding factor is:'
&ty %total%