# -*- coding: utf-8 -*-
"""
Created on Mon Nov 20 09:49:10 2017

@author: Eben Myrddin
"""


#Q1

import numpy as np
import matplotlib.pyplot as plt
import scipy

plt.close("all")

#myarray = np.fromfile('dat1.dat',dtype=float, count=-1) #this hasn't worked at all.
data=np.loadtxt("dat1.dat",skiprows=1)


#plt.plot(myarray)

X=data[:,0]
Y=data[:,1]
sig=data[:,2]

plt.plot(X,Y,'x')
plt.errorbar(X,Y,yerr=sig, linestyle="None")

#graph formatting
plt.title("plot of data", fontsize=20)
plt.xlabel("X", fontsize=20)
plt.ylabel("Y", fontsize=20)





"""
I thought it would be a good idea to use an off the shelf function to
find a value for the fit with which to compare my own method.

polyfit uses a least squares polynomial fit

"""
ComparisonCoeffs=np.polyfit(X,Y,1)


#define function for chai^2

def chaisquared(x,y,A,B,sigy):
    return sum(((y-A-B*x)**2)/(sigy**2))

#we want to minimise chai sq

N=len(Y) #100 data points

#now define function for finding separate coeffs

def findA(x,y,N):
    numerator=sum(x**2)*sum(y)-sum(x)*sum(x*y)
    denominator=N*sum(x**2)-(sum(x))**2
    return numerator/denominator #checked

def findB(x,y,N):
    num=N*sum(x*y)-sum(x)*sum(y)
    denom=N*sum(x**2)-(sum(x))**2
    return num/denom #checked

#get on with finding them

A=findA(X,Y,N)
B=findB(X,Y,N)

"""
Values compare very favourably with the polyfit least squares method
(8 decimal places of agreement)
"""

#next, must find error in y value


#######################################################################


def findsigY(N,Y,X,A,B):
    return(np.sqrt(sum((Y-A-B*X)**2)/(N-2)))

sigy=findsigY(N,Y,X,A,B)
#######################################################################





print "The A and B values given by my equations are", round(A,5), round(B,5), "respectively"


print "This compares to off-the-shelf least squares fit values of", round(ComparisonCoeffs[1],5),"and",round(ComparisonCoeffs[0],5),"respectively."

print "They clearly use the same method.", "\n", "The error in y (sigma y) was found to be equal to", sigy

#error in A
sigA=sigy*np.sqrt(sum(X**2)/(N*sum(X**2)-(sum(X)**2)))

#error in B
sigB=sigy*np.sqrt(N/(N*sum(X**2)-(sum(X)**2)))

#F*ck it let's plot


x=np.arange(0,1,0.01)
y=B*x+A
plt.plot(x,y, label="line plotted based on least squares fit") #seems to fit quite nicely.

plt.legend(loc=0)


print "\n \n \n Finally, errors for the paramters A and B are", round(sigA,5), round(sigB,5), "and the error in y is", round(sigy,5)


"""
B-B-B-BINGO!
"""