loglog.py

#!/usr/bin/env python
'''
	loglog.py
	Eric Ayars
	10/4/16

	Graphing demonstration showing how to create log-log 
	plots and curve fits relevant to Advanced Lab.
'''

from pylab import *
from scipy.optimize import curve_fit

# Data (Made-up, with realistic error)
omega = array([ 1,  2,  3,  4,  5,  6,  7,  8,  9])
force = array([ 0.025, 0.233, 0.406, 0.776, 1.136, 1.579, 2.199, 2.932, 3.694])
Om = linspace(1,9,100)

# Linear plot, with curve fits
errorbar(omega, force, yerr=0.05, fmt='bd')
title('Centripetal Force')
xlabel(r'Angular Velocity $\omega$ (rad/s)')
ylabel('Force (N)')

# Equation for curve fit
def power(t, A, B):
    return A*t**B

guess = [1,2]   # Doesn't have to be a particulary good guess, usually.
params, cov = curve_fit(power, omega, force, guess, sigma=0.05)

# add the curve fit line
plot(Om, power(Om, params[0], params[1]), 'g-')

# add annotations
text(2,3, r'$F = A \omega^B$')
text(2,2.7, r'$A = %0.4f \pm %0.4f$' %(params[0], sqrt(cov[0,0])))
text(2,2.4, r'$B = %0.2f \pm %0.2f$' %(params[1], sqrt(cov[1,1])))

# limit x,y range on display if necessary
xlim(0,10)

# If you want to save the graph...
savefig('Force.pdf')

show()

# Here's a better plot: log-log plot!
loglog(omega,force,'bd')
errorbar(omega,force,yerr=0.05,fmt='bd')
title('Centripetal Force')
xlabel(r'Angular Velocity $\omega$ (rad/s)')
ylabel('Force (N)')

# add the curve fit line
plot(Om, power(Om, params[0], params[1]), 'g-')

# add annotations
text(1.05,7, r'$F = A \omega^B$')
text(1.05,5, r'$A = %0.4f \pm %0.4f$' %(params[0], sqrt(cov[0,0])))
text(1.05,3.5, r'$B = %0.2f \pm %0.2f$' %(params[1], sqrt(cov[1,1])))

savefig('log-log_force.pdf')
show()