# ir.R: generate plots of IR spectra from dipole data
# use: on the linux command line:
# $ Rscript ir.R

# dip.dat: dipole data at each time step
# three numbers per line (x,y,z)
# read data into a data frame
df<-read.table("dip.dat")

# frequency of an oscillation of period 1 au_t in cm-1
f_au_t = 1378999.4
# dt: simulation time step in a.u.
dt = 100
# frequency scaling factor
fac = f_au_t/dt
# fmax: upper bound of frequency in plot in cm-1
fmax = 3500

# spectrum using Fourier transform
spx<-spectrum(df[[1]],plot=FALSE)
spy<-spectrum(df[[2]],plot=FALSE)
spz<-spectrum(df[[3]],plot=FALSE)

f<-spx$freq
s<-(spx$freq)^2 * (spx$spec+spy$spec+spz$spec)

# plot the spectrum on a cm-1 frequency scale in the range [0,fmax]
# plot frequency times spectrum vs frequency
# plot(fac*f,s,type="l",xlab=expression("f (cm"^'-1'*')'),ylab="",xlim=c(0,fmax))

# filtered Fourier spectrum

# normalized gaussian kernel for filter
gausswidth <- 25
kernelwidth <- 80
xg<- -kernelwidth:kernelwidth
g<-exp(-(xg/gausswidth)^2)
g<-g/sum(g)

# plot of filtered spectrum
x<-fac*f
y<-filter(s,g)

plot(x,y,type="l",col="red",lwd=2,
  xlab=expression("f (cm"^'-1'*')'),ylab="IR intensity (arb. units)",
  xlim=c(0,fmax),yaxt="n")
grid()

# write data on file ir_spectrum.dat for use with gnuplot
file<-"ir_spectrum.dat"
write(paste("# ir_spectrum",date()),file,1)
# the vector y contains "NA" values at beginning and end due to filter
# the plot range is from kernelwidth+1 to length(x)-kernelwidth
data_range<-(kernelwidth+1):(length(x)-kernelwidth)
xy<-c(x[data_range],y[data_range])
dim(xy)<-c(length(data_range),2)
write(t(xy),file,2,append=TRUE)