; $Id: Cell_Gradient.pro,v 5.4 2001/04/4
;
; Copyright (c) 2001 Javier G Corripio.
;
;+
; NAME:
;   Cell_Gradient
;
; PURPOSE:
;   Returns the gradient vector of every grid cell in a 2D variable
;   Optionally returns the surface area of every grid cell
;
; CALLING SEQUENCE:
;   Result = CellGradient(DEM, Dl [,AREA=area])
;
; INPUTS:
;   DEM:  input variable m cols x n rows
;   DL:   cell size
;
; KEYWORD PARAMETERS:
;   AREA:  returns an array containing the surface area of every grid cell
;
; OUTPUTS:
;   Result: a 3D array [3,m,n] containing the X,Y,Z coordinates
;    of the vector normal to the grid cell surface
;
;   area:  a 2D array [m,n] containing the surface area of every
;    grid cell in the inputDEM
;
; EXAMPLE
;   z = read_tiff('demdata.tif')   ; z is a 2D array with elevation data
;   Result=Cell_Gradient(z,50,AREA=surface)
;
; MODIFICATION HISTORY:
;   Written by:  Javier Corripio, May 2001.
;
;-



FUNCTION unitMvector, matrixM, LENGTH=length
On_error,2  ;Return to caller if error

typearray = size(matrixM)

IF (typearray(4)=2) THEN matMf=FLOAT(matrixM) ELSE matMf=matrixM
    ;if matrixM is integer, then it is converted to float to operate with enough accuracy

length = SQRT(TOTAL(matMf^2,1))

FOR i=0,2 DO BEGIN
  matMf(i,*,*) = matMf(i,*,*)/length
ENDFOR

RETURN, matMf
END


FUNCTION Cell_Gradient,DEM, Dl, AREA=area
On_error,2  ;Return to caller if error
nparams=N_PARAMS()
IF (nparams EQ 0) THEN BEGIN
MESSAGE, 'USAGE: Result = CellGradient(DEM, Dl, [,AREA=area])',/INFO
    return,!VALUES.F_NaN
ENDIF

sz=size(DEM,/dimensions)
m=fix(sz(0))    ; number of columns
n=fix(sz(1))    ; number of rows

;******** COMPUTE NORMAL VECTOR   ************
cellgr = FLTARR(3,m,n)*0.0
cellarea = FLTARR(m,n)*0.0

Zr=FLTARR(m,n)
Zd=FLTARR(m,n)
Zrd=FLTARR(m,n)

Zr(0:m-2,*)=DEM(1:m-1,*)
Zd(*,0:n-2)=DEM(*,1:n-1)
Zrd(0:m-2,0:n-2)=DEM(1:m-1,1:n-1)

cellgr(0,*,*) = .5*dl*(DEM+Zd-Zr-Zrd)
cellgr(1,*,*) = .5*dl*(DEM-Zd+Zr-Zrd)
cellgr(2,*,*) = (dl^2)

cellgr(*,m-1,0:n-1) = cellgr(*,m-2,0:n-1)
cellgr(*,0:m-1,n-1) = cellgr(*,0:m-1,n-2)
cellgr(*,m-1,n-1) = cellgr(*,m-2,n-2)

; norm_matvector() returns the matrix of unit vectors and the length of
; every vector  (area, in this case)



IF ARG_PRESENT(area) THEN BEGIN
    area = SQRT(TOTAL(cellgr^2,1))
    ; If the DEM has null or no data values such as  -9999
    ;  at the boundary it may produce unrealistic large surface areas
    ; uncomment next lines to set them to zero:

    ;maxarea = 3600.* dl
       ; this is a cell containing a wall twice the Eiger Nordwand
    ;too_big = where(area GT maxarea, count)
    ;if (count NE 0) THEN area(toobig) = 0
ENDIF


CellGr = unitMvector(Cellgr)


RETURN, cellgr
END