#! /bin/sh
#
#helpini
# _________________________________________________________________________________________________
#
#     Script annual_wave.sh
#
#     Target: adjusting input time serie(s) to a sinusoid
#    
#     Comments: Missing data in the input data are allowed. Missing data must be -9999.9004
#
#               Data input and output in plain txt format, sorted by columns. See format of
#               of sample data used in García-Pereira et al., 2023.
#
#
# Program created by Félix García-Pereira
# Last version: April 19, 2022
# Contact info: felgar03@ucm.es
# 
# If this program is used, please cite:
#
# García-Pereira, F., González-Rouco, J. F., Schmid, T., Melo-Aguilar, C, Vegas-Cañas, C., 
# Steinert, N. J., Roldán-Gómez, P. J., Cuesta-Valero, F. J., García-García, A., Beltrami, H., 
# and de Vrese, H.: "Thermodynamic and hydrological drivers of the subsurface thermal regime
# in Central Spain". Earth Surf. Dynam., submitted, 2023. 
#__________________________________________________________________________________________________
#helpfin        
#
# This script is run as follows:
# ./annual_wave.sh ${file} ${ncol}
#
#

dirework="<INSERT PATH WHERE SOIL TEMPERATURE SERIES WILL BE PROCESSED>"

cd ${dirework}

file=$1
ncol=$2

nfields=$(($ncol+2))
filein_name=$(echo "${file}" | awk -F "." '{print $1}' )

tau=365   # Daily data
# tau=8760  # Hourly data

for col in $(seq 3 ${nfields}) ; do

awk -v col=${col} '{print $1,$2,$col}' ${file} > temp${col}.dat

cat>annual_wave.r<<EOF
#
library(propagate)
#
ser<-read.table("temp${col}.dat", header=FALSE, dec = ".")
ser$V4[ser$V4 < -9000]<-NA
#
x<-ser[,1]
y<-ser[,3]

# Data adjustment to the function:
reslm <- lm(y ~ sin((2*pi/${tau})*x)+cos((2*pi/${tau})*x))

# To extract the adjusted timeseries:
pred <- predict(reslm, newdata=data.frame(x=x))

m <- summary(reslm)
coeff <- m\$coefficients
#coeff<-coef(reslm)

# From coefficients obtained from amplitude and phase:

b0 <- coeff[1]
alpha <- coeff[2]
beta <- coeff[3]
# da (db) is the error of alpha (beta)
da <- coeff[5]
db <- coeff[6]

# Alpha and beta vectors creation, with their respective errors.
al <- c(alpha,da)
bel <- c(beta,db)

# Matriz creation with cbind
DAT <- cbind(al, bel)

# Amplitud and phase determination
EXPR <- expression(sqrt(al^2 + bel^2))
EXPR2 <- expression(atan2(bel,al))
#amp <- sqrt(alpha^2 + beta^2)
#pha <- atan2(beta,alpha)
res <- propagate(EXPR, DAT)
rr <- res\$prop
amp <- c(rr[1], rr[4])
#
res <- propagate(EXPR2, DAT)
rr <- res\$prop
pha <- c(rr[1], rr[4])
#
ampln <- log(amp[1])

# Output creation with the coefficients of the adjustment
sali = data.frame(round((b0),3), round((amp[1]),3), round((amp[2]),3) , round((ampln),3), round((pha[1]),3), round((pha[2]),3))
#
write.table(sali, file="${filein_name}_sinfit_col${col}_coeffs.dat", col.name = FALSE, row.name = FALSE)
write.table(pred, file="${filein_name}_sinfit_col${col}_data.dat", append = FALSE, col.name = FALSE, eol = "\t", row.name = FALSE)

EOF
#
R --vanilla --slave < annual_wave.r

done # Column of data for sinusoidal adjustment

# The output file is ${filein_name}_sinfit_coeffs.dat. Columns inside are: b0 A error(A) lnA fi error(fi)
cat ${filein_name}_sinfit_col?_coeffs.dat ${filein_name}_sinfit_col??_coeffs.dat > ${filein_name}_sinfit_coeffs.dat

rm ${filein_name}_sinfit_col?_coeffs.dat ${filein_name}_sinfit_col??_coeffs.dat

