#Searcarb training

#----------------------------------------------------------------------
# Simple calculations

#----------------------------------------------------------------------
# Using up-arrow to repeat or edit previous commands (try)

#----------------------------------------------------------------------
# Create vectors of data
c(1,2,3,4,5)
1:5
seq(1,5,by=1)

#----------------------------------------------------------------------
# Variable assignment
test1 <- 1:5
test1

#----------------------------------------------------------------------
# computation with vectors
test1*10
test2 <- 1:5
test1-test2

#----------------------------------------------------------------------
# How to get help?
help(c)
help(log)
?log

#----------------------------------------------------------------------
# Other functions useful for beginners
example(log)
demo()
demo(graphics)

#----------------------------------------------------------------------
# Additional help in the help menu of R

#----------------------------------------------------------------------
# Go to CRAN to show packages
# get seacarb.pdf

# CRITICAL step to do: load seacarb!
library(seacarb)

# List of functions
help(package="seacarb")

# example of a function
?buffer

# look through the functions:
example(Kspc)

# look through the functions
example(bjerrum)

# look through the functions
example(rho)

# explore carb
?carb

# note syntax at the bottom of the R Console
# note units

#----------------------------------------------------------------------
# AT and DIC are known:
carb(flag=15, var1=2400e-6, var2=2000e-6, S=35, T=25, P=0, Pt=0, Sit=0, k1k2="l", kf="pf", pHscale="T")

# AT and DIC are known but one uses the Roy et al. constants and the free scale:
carb(flag=15, var1=2400e-6, var2=2000e-6, S=35, T=25, P=0, Pt=0, Sit=0, k1k2="r", kf="dg", pHscale="F")

# pH and AT are known:
carb(flag=8, var1=8.2, var2=0.00234, S=35, T=25, P=0, Pt=0, Sit=0,
	k1k2="l", kf="pf", pHscale="T")

# pCO2 and AT are known:
carb(flag=24, var1=370, var2=2400e-6, S=35, T=25, P=0, Pt=0, Sit=0,
	k1k2="l", kf="pf", pHscale="T")

# AT is known and one wants to know pCO2 at various pHs:
seq(7.9, 8.2, 0.05)
pH <- seq(7.9, 8.2, 0.05)
carb(flag=8, var1=pH, var2=2400e-6, S=35, T=25, P=0, Pt=0, Sit=0,
	k1k2="l", kf="pf", pHscale="T")

#----------------------------------------------------------------------
# Other examples using multiple input lines
data(seacarb_test_P0)
data <- seacarb_test_P0
data
carb(flag=data$flag, var1=data$var1, var2=data$var2, S=data$S, T=data$T, P=data$P, Sit=data$Sit, Pt=data$Pt)

#----------------------------------------------------------------------
# Store, extract and re-use data (pH and at known)
carb <- carb(flag=8, var1=pH, var2=2400e-6, S=35, T=25, P=0, Pt=0, Sit=0, k1k2="l", kf="pf", pHscale="T")
carb$pCO2
carb$OmegaCalcite

#----------------------------------------------------------------------
# Simple xy plot
plot(carb$pH,carb$pCO2)

#----------------------------------------------------------------------
# Customize simple xy plot
plot(carb$pH,carb$pCO2,col="red")
plot(carb$pH,carb$pCO2,col="red", xlab="pH (total scale)", 
ylab="pCO2 (uatm)")
	
#----------------------------------------------------------------------
# One wants to know pCO2 at various pHs and ATs:
pH <- seq(7.9, 8.2, length.out=10)
AT <- seq(2200e-6, 2600e-6, length.out=10)
dat <- expand.grid(pH, AT)
dat
carb <- carb(flag=8, var1=dat$Var1, var2=dat$Var2, S=35, T=25, P=0, Pt=0, Sit=0, k1k2="l", kf="pf", pHscale="T")
carb

# Do a simple contour plot
pCO2 <- carb$pCO2
dim(pCO2) <- c(length(pH), length(AT)) # 
contour(pH, AT, pCO2)

# Do a more elaborated contour plot
contour(pH, AT*1e6, pCO2, 
	xlab="pH (total scale)",
	ylab="Total alkalinity (umol/kg)",
	levels=seq(200, 700, by=100),  
	labcex=1.5,
	method="edge",
	col="red",
	lwd=2,
	lty="solid",
	)

#----------------------------------------------------------------------
# Convert pH to a different scale
# SWS to total scale
pHconv(flag=1, pH=8.10, S=35, T=25, P=0)

# One can convert many values in one go
pHsws <- seq(7.7, 8.3, by=0.01)	
pHsws
pHtotal <- pHconv(flag=1, pH=pHsws, S=35, T=25, P=0)
pHtotal

#----------------------------------------------------------------------
# Read data from a file
# IMPORTANT: path is needed. In R: menu "Misc:Change directory"
# Or one can use the "setwd" function: setwd("/Users/gattuso/Documents/...")

# look at the file
cc <- c(rep("numeric", 8), rep("character",3))
inp <- read.table("comparison_of_software_inp.csv", header=T, sep=",", colClasses=cc)
inp
out <- carb(flag=inp$flag, var1=inp$pCO2, var2=inp$AT, S=inp$S, T=inp$T, P=inp$Pressure, Sit=inp$Sil, Pt=inp$PO4, k1k2=inp$k1k2, kf=inp$kf, pHscale=inp$scale)
out

#----------------------------------------------------------------------
# Save data to a file
write.table(out, file="comparison_of_software_out.csv", sep=",")