In [17]:
options(repr.plot.width=7, repr.plot.height=5)
par(mai=c(0.4,0.4,0.2,0.4),bty="l",mgp=c(1.5,0.5,0))
In [18]:
library(chron)
#********************#
# CO2 concentrations #
#********************#
dat1=read.table(file="datCO2_1.txt", header=FALSE)
dat2=read.table(file="datCO2_2.txt", header=FALSE)
# Conversion to DOY
dayno1=julian(dat1$V2,dat1$V1,dat1$V3,origin=c(1,1,2017))
dayno2=julian(dat2$V2,dat2$V1,dat2$V3,origin=c(1,1,2017))
# fractional day
fday1=dayno1+(dat1$V4/24)+(dat1$V5/(24*60))+(dat1$V6/(24*60*60))
fday2=dayno2+(dat2$V4/24)+(dat2$V5/(24*60))+(dat2$V6/(24*60*60))
# check for continuity of record [0 --> continuous record]; result=number of gaps
# 30-sec measurements; 1 day--> 24x60x2=2880 values
ngaps1=length(which(c(1,round(diff(fday1)*2880))!=1))
ngaps2=length(which(c(1,round(diff(fday2)*2880))!=1))
# Insert NA for continuous record
x=fday1
xx=cbind(x,dat1$V7)
y=c(1,round(diff(x)*2880))
z=rep(x,times=(y))
ind=which(diff(z)==0)
z[which(diff(z)==0)]=NA
zz=rep(xx[,2],times=(y))
zz[ind]=NA
time.1=seq(from=fday1[1],to=(fday1[length(fday1)]+0.0001),by=1/(2880))
time.2=fday2
hour.1=24*(time.1-floor(time.1))
hour.2=24*(time.2-floor(time.2))
co2.1=zz
co2.2=dat2$V7
#********************#
# Alpha counts #
#********************#
datA1=read.table(file="datAlpha_1.txt",header=FALSE)
datA2=read.table(file="datAlpha_2.txt",header=FALSE)
timeA.1=datA1$V4
timeA.2=datA2$V4
# conversion from counts to Bq/m3
rad.1=(datA1$V7)*50.2
#
# Plots
#
plot(hour.1,co2.1,type="b",pch=20,main="14/06/2017",xlab="hour",ylab="CO2 concentration",xlim=c(10,24))
par(new=TRUE)
plot(timeA.1,rad.1,col=2,pch=4,xaxt="n",yaxt="n",xlim=c(10,24),xlab="",ylab="",ylim=c(0,250))
axis(side=4,col=2,col.axis=2,col.lab=2)
mtext("Radon (Bq/m3)",side=4,col=2,line=-2)
# conversion from counts to Bq/m3 (1-h and 2-h integration time)
rad.2=c(datA2$V7[1:6]*50.2, datA2$V7[7:14]*25.1)
plot(hour.2,co2.2,type="b",pch=20,main="15/06/2017",xlab="hour",ylab="CO2 concentration",xlim=c(10,24))
par(new=TRUE)
plot(timeA.2,rad.2,col=2,pch=4,xaxt="n",yaxt="n",xlim=c(10,24),xlab="",ylab="",ylim=c(0,250))
axis(side=4,col=2,col.axis=2,col.lab=2)
mtext("Radon (Bq/m3)",side=4,col=2,line=-2)
In [ ]: