function Hillshds=hillshademe(CurrentAzimuth, CurrentZenith, CurrentZ_fact, dElev_dy, dElev_dx)
%Function to compute a hillshade based on arguements of illumination
%azimuth, zenith, and vertical exaggeration (CurrentZ_fact). Needs the
%elevation gradients in x and y as computed from function loadfile.
%
%Azimuth is 0 towards south and increases counterclockwise so 135 is
%NE and 225 is NW
%
%Originally written by Olaf Zielke and modifies slightly by Ramon
%Arrowsmith in 2010
CurrentAzimuth = 360.0-CurrentAzimuth+90; %convert to mathematic unit
CurrentAzimuth(CurrentAzimuth>=360)=CurrentAzimuth-360;
CurrentAzimuth = CurrentAzimuth * (pi/180); % convert to radians
%lighting altitude
CurrentZenith = (90-CurrentZenith) * (pi/180); % convert to zenith angle in radians
[asp,grad] = cart2pol(dElev_dy,dElev_dx); % convert to carthesian coordinates
grad = atan(CurrentZ_fact*grad); %steepest slope
% convert asp
asp(asp