#read in the data, 1024 data points for chi-2(35) and the time in fs.
xxx<-read.table("exam1.data",header=T)

#fit the data to a model
aa<-nls(chi ~ A * sin((2 * pi)/128* t + B),data = xxx, start = c(A = max(xxx[,"chi"]), B = 0))

#print a summary
print(summary(aa))

#print the parameters from fft of signal
#normalize amplitude Mod(fft)[17] by N/2 = 512:
#   works for R or S fft algorithm
#adjust phase Arg(fft)[17] by +pi/2:
#   0 of expnl Fourier series is at argument pi/2 > 0 for sine series
cat("\nparameters from fft of signal\n")
cat("A = ", Mod(fft(xxx[,2]))[17]/512, " ;   B = ", (Arg(fft(xxx[,2]))[17] + pi/2), "\n")

#startup a plot display
x11()

#plot the data and the fitted values
plot(xxx[,"t"],xxx[,"chi"],type="l",ylab="chi-2 (degrees)",xlab="time (femtoseconds)")
lines(xxx[,"t"],fitted.values(aa),type="l",lty=2)
title(main="Motion resulting from torque applied at chi-2 of Tyr-35 of BPTI\n(dotted line, fit to data)")

#type "return" to continue
readline()

#plot the residuals vs time
plot(xxx[,"t"],residuals(aa),type="l",ylab="residuals",xlab="time (femtoseconds)")

#type "return" to continue
readline()

#analyze the residuals
qqnorm(residuals(aa))
qqline(residuals(aa))