#uncomment to retrieve the data set, 16384 data at 2 fs intervals
yout184<-dget("yout184-dget.Data")
time<-seq(1,16384)*2/1000

#turn on display device, x11 or postscript
#	(if latter, may want to comment out the readline() waits)
x11()

# assume file yout184 = 16384 point trajectory (2 fs/pt) for chi2-35;
# signal from Brownian motion of Tyr-35 ring about chi-2
# 2 fs/frame, 16384 frames
plot(time,yout184,type="l",ylab="chi-2 (degrees)",xlab="time (ps)")
title("Trajectory of chi-2 of Tyr-35, from 32 ps MD simulation")

readline()

plot(time[1:500], yout184[1:500],type="l",ylab="chi-2 (degrees)",xlab="time (ps)")
title("Trajectory of chi-2 of Tyr-35, from 32 ps MD simulation")

readline()

# a = fft of signal
a<-fft(yout184)/16384
f<-c(seq(0,8192)*1/32.768,seq(8191,1)*1/32.768)
plot(f[1:8193],Mod(a)[1:8193],type="l",ylab="Mod(fft(chi-2))",xlab="frequency (1/ps)")
title("Fourier transform of trajectory of chi-2 of Tyr-35")

readline()

plot(f[1:750],Mod(a[1:750]),type="l",ylab="Mod(fft(chi-2))",xlab="frequency (1/ps)")
title("Fourier transform of trajectory of chi-2 of Tyr-35")

readline()

# construct correlation in freq. space and invert to obtain aa=C(t)
aa<-Re(fft(Conj(a)*a,inverse=TRUE))

# fit Debye relaxation modulated by periodic term to C(t)
# = from model: tethered Brownian oscillator
# derive by inversion of chi'', dissipative (Im) part of susceptibility, chi
aa.df<-data.frame(cbind(time=(1:500)*.002,Ct=aa[1:500]))

aa3.nls<-nls(Ct ~ C*exp(-time/tau)*(cos(w*time) + 1/(tau*w)*sin(w*time)),aa.df[1:500,],start=c(tau=.1, w=20, C=20))

aa2.nls<-nls(Ct ~ aa.df[1,"Ct"]*exp(-time/tau)*(cos(w*time) + 1/(tau*w)*sin(w*time)),aa.df[1:500,],start=c(tau=.1, w=20))

I=7.7*10^15
XX=1.98*100*(180/pi)^2*4.184*10^7/(7.7*10^-15)*10^-24
aa2c.nls<-nls(Ct ~ C*exp(-time/tau)*(cos(sqrt(XX/C - 1/tau^2)*time) + 1/(tau*sqrt(XX/C - 1/tau^2))*sin(sqrt(XX/C - 1/tau^2)*time)),aa.df[1:500,],start=c(tau=.15, C=20))

aafit.nls<-aa2c.nls
tau<-coef(aafit.nls)["tau"]
C<-coef(aafit.nls)["C"]
w1<-sqrt(XX/C - 1/tau^2)

plot(time[1:500], aa.df[1:500,2],type="l",ylab="C(t)",xlab="time (ps)")
lines(time[1:500], fitted.values(aafit.nls),type="l",lty=2)
title("Correlation func. C(t) for trajectory of chi-2 of Tyr-35\n(dashed line, fitted values from Brownian relaxation model)")

#readline()

cat("\n")

# zero mean, 19 deg var, 4+ degree stddev
cat("mean chi-2 =     ", mean(yout184))
cat("\n")
cat("variance chi-2 = ", var(yout184))
cat("\n")
cat("std.dev. chi-2 = ", sqrt(var(yout184)))
cat("\n")

cat("\n")


# Parvseval's theorem - size signal vector invariant to change basis
# var(yout184) = sum(yout184^2)/16384
cat("Parvseval's theorem - size signal vector invariant to change basis\n")
cat("var(yout184) = sum(yout184^2)/16384 = sum(Mod(a)^2)) = C(0) = ", sum(Mod(a)^2))
cat("\n")
cat("\n")

cat("first 20 values of C(t); note C(0) = variance\n")
print(aa[1:20])

cat("\n")
cat("\n")

cat("summary of fit to C(t)\n")
print(summary(aafit.nls))

cat("\n")

# torsional spring (force) constant
# K(kcal/mol) = var(rad^2) * R(kcal/mol) * T(K) / 1000
# K(kcal/mol) = C(0)(rad^2) * R(kcal/mol) * T(K) / 1000
# NB: temperature of run = T = 100K
K<-(1/var(yout184)*(180/pi)^2*1.98*100/1000)
K<-(1/C*(180/pi)^2*1.98*100/1000)
cat("torsional spring constant, from fit value C(0): K  = ", K, " Kcal/mol")

cat("\n")
cat("\n")

# frequency corresponding to torsional spring constant K, w0=sqrt(K/I)
w0<-sqrt(K*1000*4.184*10^7/(7.7*10^-15))/10^12
cat("freq. of torsional oscillation, from K: w0 = ", w0, "ps^-1")

cat("\n")
cat("\n")

# frequency of periodic term in tethered Brownian model, w1=sqrt(K/I - 1/tau^2)
cat("freq. of Brownian model oscillation, from C & tau: w1 = ", w1, "ps^-1")

cat("\n")