uncertainty.py

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt

class measurement:
    def __init__(self, val, err):
        self.val = val
        self.err = err

    def __repr__(self):
        return str(self.__dict__)

    def percent_err(self):
        if self.err == 0:
            return 0
        return self.err / self.val

    def to_list(self):
        return [self.val, self.err]

    def __add__(lhs, rhs):
        return measurement(lhs.val + rhs.val, max(lhs.err, rhs.err))

    def __sub__(lhs, rhs):
        return measurement(lhs.val - rhs.val, max(lhs.err, rhs.err))

    def __mul__(lhs, rhs):
        val = lhs.val * rhs.val
        return measurement(val, max(lhs.percent_err(), rhs.percent_err()) * val)

    def __truediv__(lhs, rhs):
        val = lhs.val / rhs.val
        return measurement(val, max(lhs.percent_err(), rhs.percent_err()) * val)

def constant(x):
    return measurement(x, 0)

base_inertia = constant(0.5) * measurement(178, 0.1) / constant(1000) * \
        (measurement(95.39, 0.1) / constant(2000) * 
        measurement(95.39, 0.1) / constant(2000)) 

data = np.loadtxt("data.csv", delimiter=",", skiprows=1, unpack=True)

mass = []
rinner = []
router = []
icm = []
w1 = []
w2 = []
d = []
inertia = []
l1 = []
l2 = []
k1 = []
k2 = []
for i in range(len(data[0])):
    mass += [measurement(data[0][i], data[1][i]) / constant(1000)]
    rinner += [measurement(data[2][i], data[3][i]) / constant(2000)]
    router += [measurement(data[4][i], data[5][i]) / constant(2000)]
    icm += [constant(0.5) * mass[i] * (rinner[i] * rinner[i] + router[i] * router[i])]

    w1 += [measurement(data[6][i], data[7][i])]
    w2 += [measurement(data[8][i], data[9][i])]

    t = router[i] if rinner[i].val == 0 else rinner[i]
    d += [t - measurement(data[10][i], data[11][i]) / constant(1000)]

    inertia += [icm[i] + mass[i] * d[i] * d[i]]

    l1 += [base_inertia * w1[i]]
    l2 += [(base_inertia + inertia[i]) * w2[i]]

    k1 += [constant(0.5) * base_inertia * w1[i] * w1[i]]
    k2 += [constant(0.5) * (base_inertia + inertia[i]) * w2[i] * w2[i]]

print(w1)

lratio = []
kratio = []
for i in range(len(l1)):
    lratio += [l2[i] / l1[i]]
    kratio += [k2[i] / k1[i]]

l1 = np.transpose([m.to_list() for m in l1])
lratio = np.transpose([m.to_list() for m in lratio])
kratio = np.transpose([m.to_list() for m in kratio])

def reject_outliers(data, m=1.5):
    stdev = np.std(data, axis=1)[0]
    mean = np.mean(data, axis=1)[0]
    maskMin = mean - stdev * m
    maskMax = mean + stdev * m
    mask = np.ma.masked_outside(data[0], maskMin, maskMax)
    print('Masking values outside of {} and {}'.format(maskMin, maskMax))
    return [mask, data[1]]

print(lratio)
print(l1)

def linear(x, m, b):
    return m * x + b


def fitfunction(x):
    return popt[0] * x + popt[1]

popt, pcov = np.ma.polyfit(l1[0], lratio[0], 1, cov=True)
start = min(l1[0])
stop = max(l1[0])
xs = np.arange(start, stop, (stop - start) / 1000)
curve = fitfunction(xs)

plt.errorbar(l1[0], lratio[0], xerr=l1[1], yerr=lratio[1], fmt="o")
plt.plot(xs, curve)

plt.xlabel("$L_1$ [kg m^2 s^-1]")
plt.ylabel("$\\frac{L_f}{L_i}$")
plt.title("Best linear fit of $L_i$ to $\\frac{L_f}{L_i}$")

print("Slope: ", popt[0], "+/-", pcov[0, 0] ** 0.5)
print("Intercept: ", popt[1], "+/-", pcov[1, 1] ** 0.5)

plt.savefig("lratio.svg")
plt.show()

popt, pcov = np.ma.polyfit(l1[0], kratio[0], 1, cov=True)
start = min(l1[0])
stop = max(l1[0])
xs = np.arange(start, stop, (stop - start) / 1000)
curve = fitfunction(xs)

plt.errorbar(l1[0], kratio[0], xerr=l1[1], yerr=kratio[1], fmt="o")
plt.plot(xs, curve)

plt.xlabel("$L_1$ [kg m^2 s^-1]")
plt.ylabel("$\\frac{K_f}{K_i}$")
plt.title("Best linear fit of $L_1$ to $\\frac{K_f}{K_i}$")

plt.savefig("kratio.svg")

print("Slope: ", popt[0], "+/-", pcov[0, 0] ** 0.5)
print("Intercept: ", popt[1], "+/-", pcov[1, 1] ** 0.5)

plt.show()