# Author: Alex Gezerlis

# Numerical Methods in Physics with Python (CUP, 2020)

import numpy as np

def helpers(dataxs,datays,datasigs):

S = np.sum(1/datasigs**2)

Sx = np.sum(dataxs/datasigs**2)

Sy = np.sum(datays/datasigs**2)

Sxx = np.sum(dataxs**2/datasigs**2)

Sxy = np.sum(dataxs*datays/datasigs**2)

Del = S*Sxx - Sx**2

return S, Sx, Sy, Sxx, Sxy, Del

def computecs(dataxs,datays,datasigs):

S,Sx,Sy,Sxx,Sxy,Del = helpers(dataxs,datays,datasigs)

cs = np.zeros(2); dcs = np.zeros(2)

cs[0] = (Sxx*Sy - Sx*Sxy)/Del

cs[1] = (S*Sxy - Sx*Sy)/Del

dcs[0] = np.sqrt(Sxx/Del)

dcs[1] = np.sqrt(S/Del)

return cs, dcs

def computechisq(dataxs,datays,datasigs,cs):

chisq = np.sum((datays-cs[0]-cs[1]*dataxs)**2/datasigs**2)

return chisq

dataxs = np.linspace(0,1,6)

datays = np.array([3.085, 3.123, 3.224, 3.360, 3.438, 3.569])

datasigs = np.array([0.048, 0.053, 0.02, 0.005, 0.023, 0.07])

cs, dcs = computecs(dataxs, datays, datasigs)

print(cs); print(dcs)

chisq = computechisq(dataxs, datays, datasigs, cs)

print(chisq/(dataxs.size - cs.size))