' Name:  rusle3D
' 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

' EXPONENTS - ALTER IF NEEDED:

mparam = 0.4
nparam = 1.4

' ------------------------------------
' per Helena's note of Nov 30: 
' ...use the exponent 0.4 instead
' of 0.6 (and 1.4 instead of 1.6) - it works better with D8 
' (which overpredicts concentrated flow)


' 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 grid theme of interest
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 grid theme of interest
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 grid theme of interest
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, FALSE)
accGrid =(demGrid.FlowDirection(FALSE)).FlowAccumulation(NIL)
lsfacGrid = 
    ((accGrid * (resolution/22.13).AsGrid).Pow(mparam))*((((slpGrid * 
    (0.01745).AsGrid).Sin)/0.0896).Pow(nparam)) * (mparam + 1.0).AsGrid
lossGrid = rain.AsGrid * kfacGrid * cfacGrid * prevention.AsGrid * lsfacGrid

' Write output grid and add to view
lossname = av.GetProject.GetWorkDir.AsString + "\" + outroot + "_rusle"
if ( TRUE = MsgBox.YesNo
         ("Save result as: " + lossname + "?", "Confirm Output", TRUE) ) 
         then
             lossGrid.SaveDataSet((lossname).AsFileName)
             outTheme = Gtheme.make(lossGrid)
             outTheme.SetName( outroot + "_rusle" )
             theView.AddTheme(outTheme)
end

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