# raman.R: generate plots of Raman spectra from polarizability data
# use: on the linux command line:
# $ Rscript raman.R

# a_iso.dat: isotropic polarizability data at each time step
# one number per line
a_iso<-read.table("a_iso.dat")[[1]]

# remove mean
a_iso<-a_iso-mean(a_iso)

# 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
sp<-spectrum(a_iso,plot=FALSE)

f<-sp$freq
s<-sp$spec

# plot the spectrum on a cm-1 frequency scale in the range [0,fmax]
# plot frequency times spectrum vs frequency
# plot(fac*f,f*s,type="l",xlab="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<-f*filter(s,g)

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

# write data on file a_iso_spectrum.dat for use with gnuplot
file<-"a_iso_spectrum.dat"
write(paste("# a_iso_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)