;$Id. doshade.pro v 5.5  15/3/2002
;+
; NAME:
;   DOSHADE
;
; PURPOSE:
;   Computes terrain shading, both cast shadows and illumination intensity
;
; CATEGORY:
;   terrain parameters
;   solar radiation
;
; CALLING SEQUENCE:
;   Result = DOSHADE(Dem, Sun, Dl,[SHADING=shading])
;
;
; INPUTS:
;   Dem:  2D array with elevation data
;          (digital elevation model with regular grid spacing)
;
;   Sun:  Position of the sun:
;      3 dimensional unit vector [SunX,SunY,SunZ]
;      or 2 dimensional vector [azimuth,elevation] in degrees
;
;   Dl:     Cell size of the DEM
;
; OPTIONAL INPUTS:
;   none
;
; KEYWORD PARAMETERS:
;   SHADING:  Set this keyword to a named array variable to return
;      array of illumination intensities
;
; OUTPUTS:
;   This function returns an array of the same dimensions of DEM containing
;   binary values: 0 for shaded cell, 1 for in-sun cells
;
; OPTIONAL OUTPUTS:
;   SHADING = shading:  contains an array of the same dimensions of DEM
;     with illumination intensities, depending on the sun-terrain
;     geometry
;
; SIDE EFFECTS:
;   On some very regular DEMs or with low relief it produces some
;   spatial aliasing artifacts at certain shallow elevation
;   angles and oblique azimuthal angles.
;   Antialiasing procedure under research
;
; REQUIREMENTS
;   function CELL_GRADIENT if keyword SHADING set
;
; REMARKS
;   doshade without the SHADING keyword computes cast shadows from other cells.
;     To include self-shading (slope steeper than the local solar elevation
;     angle) set shading to a named variable and set to zero the zero values
;     of shading.
;     Example:  sh=doshade(DEM, dl,SHADING=hillshading)
;            sh = sh*(hillshading GT 0)
;
;   It follows the general convention for Solar position:
;      The sun vector is positive
;      * eastwards
;      * southwards
;      * upwards
;
; PROCEDURE:
;   See reference
;
; EXAMPLE:
;   shad = DOSHADE(demtopo, [145,35], 50, SHADING=hillshading)
;   shad = DOSHADE(demtopo, [0.46984631,0.67101007,0.57357644], 50)
;
;
; EXAMPLE with synthetic surface
;   ;creates an empty array
;   z = FLTARR(360,360)
;    ;;arrays of x and y values ranging from [-4*!PI,4*!PI]
;   x=(replicate(1,360)#FINDGEN(360))*2-360
;   Y=(FINDGEN(360)#replicate(1,360))*2-360
;   x=2*x*!dtor
;   y=2*y*!dtor
;     ;;synthetic surface
;   z=cos(x)*cos(y)*exp( (-1)/4.*SQRT(x^2+y^2) )
;
;     ;vertical exageration:
;   z=z*5
;     ; shading from 135 degrees at 30 degrees elevation
;   sh=doshade(z,[135,30],2*!dtor,shading=shading)
;     ;display the results
;   shade_surf,z,shades=bytscl(shading*sh),ax=55,az=260,xstyle=4,ystyle=4,zstyle=4
;
;     ;if borders are rugged or shaggy on very regular surfaces,
;     ;due to spatial aliasing,pass a filter:
;   k=fltarr(3,3)+.5
;   k[1,1]=1
;   shcon=convol(sh*shading,k,total(k))
;     ; display the results
;   shade_surf,z,shades=bytscl(shcon),ax=55,az=260,xstyle=4,ystyle=4,zstyle=4
;
;
; REFERENCE
;   Corripio, J.G. 2003: "Vectorial algebra algorithms for calculating terrain
;      parameters from DEMs and the position of the sun for solar radiation
;      modelling in mountainous terrain" International Journal of Geographical
;      Information Science. 17(1), 1-23.
;
; MODIFICATION HISTORY:
;   Written by:  Javier G. Corripio.
;   March, 2002
;-
;  Copyright (c) 2002 Javier G Corripio.
;  jgc@geo.ed.ac.uk  http://www.geo.ed.ac.uk/~jgc/


FUNCTION DOSHADE, dem, Sun, dl, SHADING=shading
  On_error,2
nparams=N_PARAMS()
    ;check that the argument exists
IF (nparams EQ 0) THEN BEGIN
    ;message=[ ['USAGE: Result = DOSHADE(Dem, Sun, dl, [SHADING=hillshading])'] ]
    ;print,message
     message, 'USAGE: Result = DOSHADE(Dem, Sun, dl, [SHADING=hillshading])   -------     '  + $
     			' Sun = [x, y, z] or Sun = [azimuth, elevation] in degrees', /INFO

        RETURN, !VALUES.F_NAN
ENDIF

SunVector = double(sun) ;*** for the 3D case
SunVector = SunVector / SQRT(TOTAL(SunVector^2)) ; in case is not unit vector

IF Size(SunVector,/dimensions) EQ 2 THEN BEGIN  ; for the 2D case
    azimuth = double(sun[0] * (!dPI/180.0d) )
    elevation = double(sun[1] * (!dPI/180.0d) )
    SunVector = dblarr(3)
    SunVector[0] = SIN(azimuth)*COS(elevation)     ;positive Eastwards
    SunVector[1] = -COS(azimuth)*COS(elevation) ;positive Southwards
    SunVector[2] = SIN(elevation)
ENDIF

; If operated with a different reference system, eg: Y is positive northwards,
;   simply make SunVector[1] = -SunVector[1]


; solvector is opposite to sunvector and used for integer steps scanning
;   along the solar direction
; this makes the largest coordinate of solvector =  1

 SolVector = -SunVector/MAX(ABS(SunVector(0:1)))

SunVectorPXY = SunVector    ; horizontal pojection sunvector, step1
SunVectorPXY[2] = 0.0d   ; horizontal pojection sunvector, step2

; vector normal to sun in the horizontal plane
SunVectorH = CROSSP(SunVectorPXY,SunVector)
SunVectorH = SunVectorH/SQRT(TOTAL(SunVectorH^2)) ;** Unit SunVectorH

; unit vector normal to sun upwards
SunVectorU = CROSSP(SunVectorH,SunVector)

; direction of view of sun
SunVectorN = -SunVector


sz=size(dem,/dimensions)
cols = FIX(sz[0])   ; number of columns
rows = FIX(sz[1])   ; number of rows

;Vector to Origin of coordinates, for projection onto plane perpendicular to sun
 ; It doesn't need to be (0,0,0) any arbitrary point is fine (0,0,0) is simpler

xvalues = INDGEN(cols) # REPLICATE(1,rows)
yvalues = REPLICATE(1,cols) # INDGEN(rows)
xvalues = xvalues * dl
yvalues = yvalues * dl
zvalues = dem

;VectortoOrigin[i] is [xvalues[i],yvalues[i],zvalues[i]

;dot product [x,y,z] * SunVectorU
zprojection =   xvalues[*,*] * SunVectorU[0] +$
       yvalues[*,*] * SunVectorU[1] +$
       zvalues[*,*] * SunVectorU[2]


; Determine origin of scanning lines in the direction of the sun
CASE 1 OF
    (SunVector[0] LE 0.0): Xend=cols-1 ; sun is on the West
    (SunVector[0] GT 0.0): Xend=0   ; sun is on the East
    ELSE:
ENDCASE

CASE 1 OF
    (SunVector[1] LE 0.0): Yend = rows-1   ; sun is on the North
    (SunVector[1] GT 0.0): Yend = 0         ; sun is on the South
    ELSE:
ENDCASE

; sombra is the array to store binary shadow values
sombra = replicate(1,cols,rows)

;------------- SCANNING CASE SUN at  0, 90, 180, 270 degrees -----------

; 0 or 180 degrees
IF (Sunvector[0] EQ 0) THEN BEGIN
    lastloop = cols-1
    lengthSC = rows-1
    Firsti = 0
    FirstJ = rows-Yend-1
    scanindex=[0,1]
    ScanVector = INTARR(2,lengthSC)
    ScanVector[1,*] = Firstj + indgen(lengthSC) * SolVector[1]
    GOTO, JumpNoZero
ENDIF

; 90 or 270 degrees
IF (Sunvector[1] EQ 0) THEN BEGIN
    lastloop = rows-1
    lengthSC = cols-1
    Firsti = Xend-cols-1
    FirstJ = 0
    scanindex=[1,0]
    ScanVector = INTARR(2,lengthSC)
    ScanVector[0,*] = Firsti + indgen(lengthSC) * SolVector[0]
    GOTO, JumpNoZero
ENDIF

;-------------------------  SCANNING BLOCK I ----------------------------------

;***************** scanning along X-AXIS  ****************

FOR i = 0, cols-1 DO BEGIN
; length of scanning line (vector)
    lengthX1 = ABS(Xend - i) / SolVector[0]    ; solve for: l*Solvector=Xmax
    lengthX2 = (rows-1) / SolVector[1]     ; solve for: l*Sv=Ymax
    lengthSC = FIX(MIN(ABS([lengthX1,lengthX2]) ))

    IF (lengthSC EQ 0 ) THEN CONTINUE ; this means a corner = sunny

    vectorlength = INDGEN(lengthSC)

; X-origin of scan line
    firsti = i
; Y-origin of scan line
    firstj = rows-Yend-1

; ScanVectorX: arrays of xvalues along scan line  ScanVectorY: Yvalues
    ScanVectorX = FIX(i + vectorlength * SolVector[0])

; ScanVectorY: Yvalues
    ScanVectorY = FIX(Firstj + vectorlength * SolVector[1])

; Array of Z_Projections along scan line
    ZPdxdy = zprojection[ScanVectorX,ScanVectorY]

; Some antialiasing device is required.  A temporary improvement is
    ;achieved by averaging the cell projections with the previous value
    ;along the scan line

    IF (lengthsc GT 2) THEN $
    ZPdxdy[1:lengthsc-1] = 0.5*(ZPdxdy[0:lengthsc-2] + ZPdxdy[1:lengthsc-1] )

; initial zpcompare smaller than any possible zprojection
    zpcompare =-1e30

    FOR n=0,lengthSC-1 DO BEGIN
       IF (ZPdxdy[n] LT zpcompare) THEN $
         sombra[ScanvectorX[n],ScanvectorY[n]] = 0 ELSE $
         zpcompare = ZPdxdy[n]
    ENDFOR
ENDFOR

xend = ABS(XEND-i)-1

;==============================================================================
;-------------------------  SCANNING BLOCK J ----------------------------------

;*** scanning   Along Y-AXIS
lenarray=fltarr(rows)
FOR j = 0, rows-1 DO BEGIN
;length of scanning line (vector)
    lengthY1 = ABS(Yend - j) / SolVector[1]    ; solve for: l*Sv=Ymax
    lengthY2 = (cols-1) / SolVector[0] ; solve for: l*Sv=Xmax
    lengthSC = FIX(MIN(ABS([lengthY1,lengthY2]) ))

    IF (lengthSC EQ 0 ) THEN CONTINUE ; this means a corner = sunny

    vectorlength = INDGEN(lengthSC)
lenarray[j]=lengthSC
; X-origin of scan line
    firsti = Xend
; Y-origin of scan line
    firstj =  j

 ; ScanVectorX: arrays of xvalues along scan line
    ScanVectorX = FIX(Xend + vectorlength * SolVector[0])

; ScanVectorY: Yvalues
    ScanVectorY = FIX(Firstj + vectorlength * SolVector[1])

;Array of Z_Projections along scan line
    ZPdxdy = zprojection[ScanVectorX,ScanVectorY]

; Some antialiasing device is required.  A temporary improvement is
    ;achieved by averaging the cell projections with the previous value
    ;along the scan line

    IF (lengthsc GT 2) THEN $
    ZPdxdy[1:lengthsc-1] = 0.5*(ZPdxdy[0:lengthsc-2] + ZPdxdy[1:lengthsc-1] )

; initial zpcompare smaller than any possible zprojection
    zpcompare =-1e30

    FOR n=0,lengthSC-1 DO BEGIN
       IF (ZPdxdy[n] LT zpcompare) THEN $
         sombra[ScanvectorX[n],ScanvectorY[n]] = 0 ELSE $
         zpcompare = ZPdxdy[n]
    ENDFOR

ENDFOR

;==============================================================================

GOTO, JumpZero

;-------------------------  SCANNING BLOCK ZERO -----------------------------

JUMPNOZERO:

FOR i = 0,lastloop DO BEGIN

    ScanVector[scanindex[0],*] = i
; Array of Z_Projections along scan line
    ZPdxdy = zprojection[ScanVector[scanindex[0],*],ScanVector[scanindex[1],*]]

; Some antialiasing device is required.  A temporary improvement is
    ;achieved by averaging the cell projections with the previous value
    ;along the scan line

    ZPdxdy[1:lengthsc-1] = 0.5*(ZPdxdy[0:lengthsc-2] + ZPdxdy[1:lengthsc-1] )

; initial zpcompare smaller than any possible zprojection
    zpcompare =-1e30

    FOR n=0,lengthSC-1 DO BEGIN
      IF (ZPdxdy[n] LT zpcompare) THEN $
       sombra[Scanvector[scanindex[0],n],Scanvector[scanindex[1],n]] = 0 $
       ELSE zpcompare = ZPdxdy[n]
    ENDFOR
ENDFOR

;==============================================================================

JUMPZERO:


; border columns and rows are set to 1st (or 2nd) inner column and row
;sombra[1,*] = sombra[2,*]
sombra[0,*] = sombra[1,*]
;sombra[cols-2,*] = sombra[cols-3,*]
sombra[cols-1,*] = sombra[cols-2,*]
;sombra[*,1] = sombra[*,2]
sombra[*,0] = sombra[*,1]
;sombra[*,rows-2] = sombra[*,rows-3]
sombra[*,rows-1] = sombra[*,rows-2]

IF  ARG_PRESENT(shading) THEN BEGIN
    gr = cell_gradient(dem,Dl)
    shading=fltarr(cols,rows)
    shading[*,*] = gr[0,*,*]* SunVector[0] + $
         gr[1,*,*]* SunVector[1] + $
         gr[2,*,*]* SunVector[2]
    shading =  shading > 0.
ENDIF


RETURN, sombra

END