' Name:  USPED
' Author: Bill Brown wmbrown@uiuc.edu
' Date: Nov 2001
' Revision Date: Tue Dec 4 2001
' --------------------------------------------------------------------
' Description: Computes Soiloss grid from 3 input grids representing 
' elevation, cfactor, and kfactor
' --------------------------------------------------------------------
' ------------------------------------
' VARIABLES : EDIT THESE FOR YOUR USE!
'     RESOLUTION MUST BE IN METERS

resolution = 3.0
rain = -280
prevention = 1.0

' RILL EXPONENTS - ALTER IF NEEDED:

mparam = 1.6
nparam = 1.3

' ------------------------------------


' GENERAL DATA GATHERING
theView = av.GetActiveDoc
theThemes = theView.GetThemes
thePrj = theView.GetProjection

' Get the Analysis Environment
ae = theView.GetExtension(AnalysisEnvironment)

' Create list for available grid themes. Make sure at least one
' grid is available.
theGridList = List.Make
for each aTheme in theThemes
  if (aTheme.Is(GTheme)) then
      theGridList.Add(aTheme)
  end
end
if (theGridList.Count < 1) then
  MsgBox.Error( "Not enough grid themes available", "Oops")
  Return NIL
end


' Select the elevation grid theme 
theGTheme_dem = MsgBox.List(theGridList,
            "Select ELEVATION grid: ",
            "Grid Theme(s)")
if (theGTheme_dem = NIL) then Return NIL end
demGrid = theGTheme_dem.GetGrid

' Select the c-factor grid theme
theGTheme_cfac = MsgBox.List(theGridList,
            "Select C-FACTOR grid: ",
            "Grid Theme(s)")
if (theGTheme_cfac = NIL) then Return NIL end
cfacGrid = theGTheme_cfac.GetGrid

' Select the k-factor grid theme 
theGTheme_kfac = MsgBox.List(theGridList,
            "Select K-FACTOR grid: ",
            "Grid Theme(s)")
if (theGTheme_kfac = NIL) then Return NIL end
kfacGrid = theGTheme_kfac.GetGrid

outroot = MsgBox.Input ("Enter root filename to use for outputs", 
          "Output", theGTheme_dem.GetName)

slpGrid = demGrid.Slope (1.0, TRUE)
aspGrid = demGrid.Aspect
accGrid =(demGrid.FlowDirection(FALSE)).FlowAccumulation(NIL)

' Calculate for prevailing rill erosion
sflow1 = 
    ((accGrid * resolution.AsGrid).Pow(mparam)) * 
    ((slpGrid * (0.01745).AsGrid).Sin).Pow(nparam)

qsx1 = ((((aspGrid * (-1).AsGrid)+450) * (0.01745).AsGrid).Cos) * sflow1 * 
	kfacGrid * cfacGrid * rain.AsGrid * prevention.AsGrid

qsy1 = ((((aspGrid * (-1).AsGrid)+450) * (0.01745).AsGrid).Sin) * sflow1 * 
	kfacGrid * cfacGrid * rain.AsGrid * prevention.AsGrid


qsx_slope1 = qsx1.Slope (1.0, FALSE)
qsx_aspect1 = qsx1.Aspect
qsy_slope1 = qsy1.Slope (1.0, FALSE)
qsy_aspect1 = qsy1.Aspect

qsx_dx1 = ((((450).AsGrid - qsx_aspect1) * (0.01745).AsGrid).Cos) * 
    (qsx_slope1 * (0.01745).AsGrid).Tan

qsy_dy1 = ((((450).AsGrid - qsy_aspect1) * (0.01745).AsGrid).Sin) * 
    (qsy_slope1 * (0.01745).AsGrid).Tan

erdprill = qsx_dx1 + qsy_dy1       

' Calculate for prevailing sheet erosion
sflow2 = (accGrid * resolution.AsGrid)* ((slpGrid * (0.01745).AsGrid).Sin)

qsx2 = ((((aspGrid * (-1).AsGrid)+450) * (0.01745).AsGrid).Cos) * 
	sflow2 * kfacGrid * cfacGrid * rain.AsGrid * prevention.AsGrid

qsy2 = ((((aspGrid * (-1).AsGrid)+450) * (0.01745).AsGrid).Sin) * 
	sflow2 * kfacGrid * cfacGrid * rain.AsGrid * prevention.AsGrid


qsx_slope2 = qsx2.Slope (1.0, FALSE)
qsx_aspect2 = qsx2.Aspect
qsy_slope2 = qsy2.Slope (1.0, FALSE)
qsy_aspect2 = qsy2.Aspect

qsx_dx2 = ((((450).AsGrid - qsx_aspect2) * (0.01745).AsGrid).Cos) * 
	(qsx_slope2 * (0.01745).AsGrid).Tan

qsy_dy2 = ((((450).AsGrid - qsy_aspect2) * (0.01745).AsGrid).Sin) * 
	(qsy_slope2 * (0.01745).AsGrid).Tan

erdpsheet = (qsx_dx2 + qsy_dy2) * (resolution).AsGrid
        


' Write output grid (prevailing rill erosion) and add to view
rillname = av.GetProject.GetWorkDir.AsString + "\" + outroot + "_usp_rill"
if ( TRUE = MsgBox.YesNo
         ("Save (prevailing rill) result as: " + rillname + "?", 
	 "Confirm Output", TRUE) ) 
         then
             erdprill.SaveDataSet((rillname).AsFileName)
             outTheme1 = Gtheme.make(erdprill)
             outTheme1.SetName( outroot + "_usp_rill" )
             theView.AddTheme(outTheme1)
end

' Write output grid (prevailing sheet erosion) and add to view
sheetname = av.GetProject.GetWorkDir.AsString + "\" + outroot + "_uspedsheet"
if ( TRUE = MsgBox.YesNo
         ("Save (prevailing sheet) result as: " + sheetname + "?", 
	 "Confirm Output", TRUE) ) 
         then
             erdpsheet.SaveDataSet((sheetname).AsFileName)
             outTheme2 = Gtheme.make(erdpsheet)
             outTheme2.SetName( outroot + "_usp_sheet" )
             theView.AddTheme(outTheme2)
end



' Delete temporary files
' 
' 
' File.Delete(xxx)
                                        
Msgbox.Info("Calculations completed successfully","Goodbye")