include("Path.jl")
include("Reading.jl")
include("Table.jl")
include("Plots.jl")
include("Option.jl")

using DataFrames
#=
	For New Zealand soil classification
	Webb abd Liburne Criteria for defining the soil family and soil sibling
=#
function START_PSD()
	println("=== START RUNNING PSD TEXTURE === \n")
		∑Psd, Id, N_∑Psd = reading.PSD(path.Input)

		Sand = Array{Float64}(undef, N_∑Psd)
		Silt = Array{Float64}(undef, N_∑Psd)
		Clay = Array{Float64}(undef, N_∑Psd)
		SoilClass = Array{String}(undef, N_∑Psd)

		for iSoil in 1:N_∑Psd
			Clay[iSoil] =  ∑Psd[iSoil, 1]
			Silt[iSoil] = ∑Psd[iSoil, 5] - Clay[iSoil]
			Sand[iSoil] = ∑Psd[iSoil,10] - Silt[iSoil]

			ClaySiltSand = Clay[iSoil] + Silt[iSoil] + Sand[iSoil]

			if 0.9999 < ClaySiltSand < 1.0001
				error("ID = $iSoil  Clay+Sand+Silt≠1")
			end  # if name
		end #  iSoil in 1:N_∑Psd

		if option.Classification == "NZsimple"
			N_Clayey = 0
			N_Silty = 0
			N_Loamy = 0
			N_Sandy = 0
			N_Error = 0

		elseif option.Classification == "NZ"
			N_Sand = 0
			N_LoamySand = 0
			N_SandyLoam = 0
			N_LoamySilt = 0
			N_Silt = 0
			N_SandyClayLoam = 0
			N_ClayLoam = 0
			N_SiltLoam = 0
			N_SiltyClay = 0
			N_LoamyClay = 0
			N_Clay = 0
			N_Error = 0
		end

		
		for iSoil in 1:N_∑Psd
			if option.Classification == "NZsimple"
				if Clay[iSoil] >= 0.35
					SoilClass[iSoil] = "Clayey"
					N_Clayey += 1
	
				elseif Clay[iSoil] < 0.35 && Silt[iSoil] >= 0.40
					SoilClass[iSoil] = "Silty"
					N_Silty += 1
	
				elseif Clay[iSoil] < 0.35 && Silt[iSoil] < 0.40  && Clay[iSoil] >= 0.08
					SoilClass[iSoil] = "Loamy"
					N_Loamy += 1
	
				elseif Clay[iSoil] < 0.35 && Silt[iSoil] < 0.40  && Clay[iSoil] < 0.08
					SoilClass[iSoil] = "Sandy"
					N_Sandy += 1
					
				else
					SoilClass[iSoil] = "Error"
					N_Error += 1
				end

			elseif option.Classification == "NZ" # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
				if Clay[iSoil] <= 0.08 && Sand[iSoil] >= 0.8
					SoilClass[iSoil] = "Sand"
					N_Sand += 1
	
				elseif Clay[iSoil] <= 0.08 && Silt[iSoil] < 0.40 &&  Sand[iSoil] < 0.8
					SoilClass[iSoil] = "LoamySand"
					N_LoamySand += 1
	
				elseif Clay[iSoil] > 0.08 && Silt[iSoil] < 0.40
					SoilClass[iSoil] = "SandyLoam"
					N_SandyLoam += 1
	
				elseif 0.4 < Silt[iSoil] < 0.82
					SoilClass[iSoil] = "LoamySilt"
					N_LoamySilt += 1
					
				elseif 0.82 < Silt[iSoil]
					SoilClass[iSoil] = "Silt"
					N_Silt += 1

				elseif Silt[iSoil] <= 0.15
					SoilClass[iSoil] = "SandyClayLoam"
					N_SandyClayLoam += 1

				elseif 0.15 < Silt[iSoil] < 0.4
					SoilClass[iSoil] = "ClayLoam"
					N_ClayLoam += 1

				elseif 0.4 <= Silt[iSoil]
					SoilClass[iSoil] = "SiltLoam"
					N_SiltLoam += 1

				elseif Clay[iSoil] < 0.60 && 0.30 < Silt[iSoil]
					SoilClass[iSoil] = "SiltyClay"
					N_SiltyClay += 1

				elseif Clay[iSoil] < 0.60 && Silt[iSoil] < 0.3
					SoilClass[iSoil] = "LoamyClay"
					N_LoamyClay += 1

				elseif Clay[iSoil] >= 0.60
					SoilClass[iSoil] = "Clay"
					N_Clay += 1

				else
					N_Error += 1
				end

			end 
		end # for iSoil in 1:N_∑Psd

		# OUTPUTS
			table.Output(path.Output, ∑Psd, Id, N_∑Psd, Sand, Silt, Clay, SoilClass)
			emptyDF = DataFrame(N_Clayey=Int64[], N_Silty=Int64[], N_Loamy=Int64[], N_Sandy=Int64[], N_Error=Int64[])

			if option.Classification == "NZsimple"
				Statistics_N = DataFrame(N_Clayey=N_Clayey, N_Silty=N_Silty, N_Loamy=N_Loamy, N_Sandy=N_Sandy, N_Error=N_Error)
				println(Statistics_N, "\n")
	
				emptyDF = DataFrame(Sand=Float64[], Silt=Float64[], Clay=Float64[])
				Statistics_SandSiltClay = DataFrame(Sand=Sand, Silt=Silt, Clay=Clay)
				println(describe(Statistics_SandSiltClay))
			
			elseif option.Classification == "NZ"
				Statistics_N = DataFrame(N_Sand = N_Sand, N_LoamySand = N_LoamySand, N_SandyLoam = N_SandyLoam, N_LoamySilt = N_LoamySilt, N_Silt = N_Silt, N_SandyClayLoam = N_SandyClayLoam, N_ClayLoam = N_ClayLoam,  N_SiltLoam = N_SiltLoam, N_SiltyClay = N_SiltyClay, N_LoamyClay = N_LoamyClay, N_Clay = N_Clay, N_Error=N_Error)
				println(Statistics_N, "\n")	
			end



			# Work in progress
			plots.TERNARY(path.Plots, Sand, Silt, Clay)

		println("\n === END RUNNING PSD TEXTURE ===")

end # function START_PSD()


START_PSD()