project 8 large amplitude pendulum

import numpy as np import matplotlib.pyplot as pl import pandas as pd from scipy.optimize import curve_fit # automated curve fitting

df_alex = pd.read_csv('alex1.csv') df_alex.head()

pl.plot(df_alex.time, df_alex.angle, 'b-')

good_time = (df_alex.time > 8.5) & (df_alex.time<16) time = df_alex[good_time].time angle = df_alex[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.title("Fit at around 50 degrees") pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_dimi1 = pd.read_csv('dimitris1.csv') df_dimi1.head()

pl.plot(df_dimi1.time, df_dimi1.angle, 'b-')

good_time = (df_dimi1.time > 8.5) & (df_dimi1.time<16) time = df_dimi1[good_time].time angle = df_dimi1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_dimi2 = pd.read_csv('dimitris2.csv') df_dimi2.head()

pl.plot(df_dimi2.time, df_dimi2.angle, 'b-')

good_time = (df_dimi2.time > 8.5) & (df_dimi2.time<16) time = df_dimi2[good_time].time angle = df_dimi2[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.title("Fit at around 50 degrees") pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_mitch1 = pd.read_csv('mitchell1.csv') df_mitch1.head()

pl.plot(df_mitch1.time, df_mitch1.angle, 'b-')

good_time = (df_mitch1.time > 8.5) & (df_mitch1.time<16) time = df_mitch1[good_time].time angle = df_mitch1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.title("Fit at around 50 degrees") pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_soren1 = pd.read_csv('soren1.csv') df_soren1.head()

pl.plot(df_soren1.time, df_soren1.angle, 'b-')

good_time = (df_soren1.time > 8.5) & (df_soren1.time<16) time = df_soren1[good_time].time angle = df_soren1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.title("Fit at around 50 degrees") pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_soren2 = pd.read_csv('soren2.csv') df_soren2.head()

pl.plot(df_soren2.time, df_soren2.angle, 'b-')

bad_time = (df_soren2.time > 8.5) & (df_soren2.time<16) time = df_soren2[bad_time].time angle = df_soren2[bad_time].angle pl.plot(time, angle ,'b.') pl.grid()

T0 = 2 A0 = 45 def model(t, A, omega, phi, off): return A*np.sin(omega*t + phi) + off par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0)) A, omega, phi, off = par dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov)) print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm)) pl.plot(time, angle ,'b.', label='data') pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit") pl.grid() pl.title("Fit at around 50 degrees") pl.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

import os import csv fname = 'rawdata.csv' dataFile = open(os.path.join('rawdata',fname)) reader = csv.reader(dataFile, delimiter=',') lineCnt=0 t1_list=[] t2_list=[] theta1_list=[] theta2_list=[] for line in reader: if lineCnt in [0]: pass else: try: t1_list.append(float(line[1])) theta1_list.append(float(line[2])*np.pi/180.0) # convert to radians t2_list.append(float(line[3])) theta2_list.append(float(line[4])*np.pi/180.0) # convert to radians except: print("Ack!", line) lineCnt+=1 N=len(t1_list) t1arr=np.array(t1_list) t2arr=np.array(t2_list) th1arr=np.array(theta1_list)+off # theta 1 corrected for offset th2arr=np.array(theta2_list)+off # theta 2 corrected for offset ratio = (t2arr-t1arr)/T0 # array of experimental periods devided by zero amplitude period thAvg = (th1arr+th2arr)/2.0 vsaa = 1.0*np.ones(N) saa = 1.0+thAvg*thAvg/16 def simpson(f, a, b, N, args=None): if args is None: args=[] x=np.linspace(a,b,2*N) # sample the function h=(b-a)/(2*N-1) y=f(*([x] +list(args))) evens = y[2:-2:2] odds = y[1:-1:2] return (y[0] + y[-1] + 2*odds.sum() + 4*evens.sum())*h/3.0 def K_int(phi, m): return 1.0/np.sqrt(1-m*np.sin(phi)**2) def K(m, N=1000): return simpson(K_int, 0, np.pi/2, N=N, args=(m,)) sraList=[] for theta in thAvg: sraList.append(2*K(np.sin(theta/2.0)**2)/np.pi) thDeg=thAvg*180/np.pi # array of average amplitudes (in radians) pl.title("Amplitude of a simple pendulum") pl.xlabel("$\\theta_{\\rm max}$ (degrees)") pl.ylabel("T/To") pl.axis([30,70,0.99,1.1]) pl.plot(thDeg, ratio,'b.', label="experiment") pl.plot(thDeg, saa, 'r-',label="small angle approx.") pl.plot(thDeg, vsaa, 'g-', label="very small angle approx.") pl.plot(thDeg, sraList, 'c-', label="exact theory") pl.legend(loc=2)