yout184<-dget("yout184-dget.Data")
time<-seq(1,16384)*2/1000
postscript(file="scr.ps",horiz=FALSE)
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")
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")
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")
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")
aa<-Re(fft(Conj(a)*a,inverse=TRUE))
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)")
graphics.off()
sink("scr.txt")
cat("\n")
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")
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")
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")
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")
cat("freq. of Brownian model oscillation, from C & tau: w1 = ", w1, "ps^-1")
cat("\n")
sink()