# updated 8/13/2010

# This S-PLUS script computes left-censored summary statistics for NAWQA Major Ion data
# It also computes summary statistics and quantiles for Physical Properties and Major Ions
# that do not have censored data

# Text output will be written to the file specified in the sink command.
sink("report.majors.summary2.out.txt")



# Read data file
# Data file updated 8/13/10 to add alkalinity and bicarbonate for 4 WMIC wells
# that were missing overlooked values because they were for unfiltered water samples
attach.data.frame(majors.nawqa.mdl)


#start of wilson function that we create as 
wilson<-function(target, target2, target3) {


cat("Processing target", target3, "\n")

#Parametric Summary Statistics for left censored data
print(mdlpar(lcens(target, target2 == "<")))

cat("Processing target", target3, "\n")

#Parametric Summary Statistics and Percentiles for left censored data
print(summary(mdlpar(lcens(target, target2 == "<"))))



# End of wilson function
}

# for every trace element generate a new line below
# following the example

wilson(bromide, bromide.rmk, "bromide")
wilson(chloride, chloride.rmk, "chloride")
wilson(fluoride, fluoride.rmk, "fluoride")
wilson(do, do.rmk, "do")




#start of wilson2 function that we create as 
wilson2<-function(target, target5) {


cat("Processing target", target5, "\n")

# do summary statistics for listed analytes
print(summary(target))

cat("Processing target", target5, "\n")

# print 10th and 90th percentiles for listed analytes
print(quantile(target,c(0.10,0.90), na.rm=T))



# End of wilson2 function
}



wilson2(bromide, "bromide")
wilson2(chloride, "chloride")
wilson2(fluoride, "fluoride")
wilson2(do, "do")
wilson2(calcium, "calcium")
wilson2(magnesium, "magnesium")
wilson2(potassium, "potassium")
wilson2(TDS, "TDS")
wilson2(silica, "silica")
wilson2(sodium, "sodium")
wilson2(sulfate, "sulfate")
wilson2(alkalinity, "alkalinity")
wilson2(HCO3, "HCO3")
wilson2(sc, "sc")
wilson2(temperature, "temperature")
wilson2(ph, "ph")
wilson2(hardness, "hardness")





#Detach the data frame
detach(majors.nawqa.mdl)

# close output file
sink()