#import modules arcpy and string
import arcpy, string

#import the environment settings and spatial analyst extension
from arcpy import env
from arcpy.sa import *

#Check out the Spatial Analyst extension (must be available to check out!)
arcpy.CheckOutExtension('spatial')

#Add your workspace here
env.workspace = r'C:WORKSPACE'
env.overwriteOutput = True

#Define inputs from TIFFs
Band2 = 'LANDSAT_B2.tif' #Band 2
Band4 = 'LC08_L1TP_019036_20130604_20170310_01_T1_B4.tif' #Band 4
Band5 = 'LC08_L1TP_019036_20130604_20170310_01_T1_B5.tif' #Band 5

#Band Ratios
R42 = arcpy.sa.Divide(Band4, Band2)
R52 = arcpy.sa.Divide(Band5, Band2)
R54 = arcpy.sa.Divide(Band5, Band4)
print 'Ratioed Bands'

#Slopes
a = 52
b = 11
c = 33.5
d = 120

#Algorithm processing
ALG = a - (b * R42) + (c * R52) - (d * R54)

#Saving final product
result = 'ALG.tif'
ALG.save(result)
print 'Success'

arcpy.env.workspace=r'C:/WORKSPACE/Band2'
B2rasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'C:/WORKSPACE/Band4'
B4rasters=arcpy.ListRasters('*.tif*')

for idx in range(0, len(B2rasters) - 1):
    b2 = B2rasters[idx]
    b4 = B4rasters[idx]
output = b4 / b2
result = output + idx
output.save(result)
print result