# -*- coding: utf-8 -*-
"""
Created on Mon Oct 20 09:31:43 2014

@author: mrazny
"""

import matplotlib as mpl
import numpy as np
import math
from scipy.interpolate import interp1d

#norma pierwsza, to całka fun od a do b
#norma druga, to całka fun^2 od a do b

def caleczkaSimps(a, b, pktWylicz, pktWyliczDokl, dokl):
    sumaT = 0;
    suma = 0;
    h = (b-a)/dokl;
    i = a;
    j = 0;
    while i<=b-h:
        sumaT += (np.abs((pktWylicz[j]-pktWyliczDokl[j])) + np.abs((pktWylicz[j+1]-pktWyliczDokl[j+1])))/2.;
        j+=1;
        i+=h;
        if i<h:
            suma += np.abs(pktWylicz[j] - pktWyliczDokl[j]);
    return h/6.*(pktWylicz[0]+pktWylicz[j]+2.*suma+4*sumaT);

def caleczkaSimps2(a, b, pktWylicz, pktWyliczDokl, dokl):
    sumaT = 0;
    suma = 0;
    h = (b-a)/dokl;
    i = a;
    j = 0;
    while i<=b-h:
        sumaT += (((pktWylicz[j]-pktWyliczDokl[j]) + (pktWylicz[j+1]-pktWyliczDokl[j+1]))/2.)*(((pktWylicz[j]-pktWyliczDokl[j]) + (pktWylicz[j+1]-pktWyliczDokl[j+1]))/2.);
        j+=1;
        i+=h;
        if i<h:
            suma += (pktWylicz[j] - pktWyliczDokl[j])*(pktWylicz[j] - pktWyliczDokl[j]);
    return math.sqrt(h/6.*(pktWylicz[0]+pktWylicz[j]+2.*suma+4*sumaT));

def funkcja(x):
    return np.exp(x/(x-1))

def pochodna(x):
    return -(np.exp(x/(x-1))/((x-1)*(x-1)))

################################################
############## TUTAJ MASZ SPLAJNA ##############
################################################
def szukaj_wart(x, xs, ys, ks): #szukana wart, tabx, taby, tabpochodna
    i = 1
    while xs[i]<x:
        i+=1

    t = (x - xs[i-1]) / (xs[i] - xs[i-1])

    a =  ks[i-1]*(xs[i]-xs[i-1]) - (ys[i]-ys[i-1])
    b = -ks[i  ]*(xs[i]-xs[i-1]) + (ys[i]-ys[i-1])

    q = (1-t)*ys[i-1] + t*ys[i] + t*(1-t)*(a*(1-t)+b*t)
    return q
################################################
############## TUTAJ MASZ SPLAJNA ##############
################################################

n1 = 5
n2 = 9
n3 = 17
n4 = 33
n5 = 65
dokl = 257

a = 1.075
b = 1.08

# tworzenie siatek z wat x i y danych pkt
siatka1 = np.linspace(a, b, n1)
siatka2 = np.linspace(a, b, n2)
siatka3 = np.linspace(a, b, n3)
siatka4 = np.linspace(a, b, n4)
siatka5 = np.linspace(a, b, n5)
siatkaDokl = np.linspace(a, b, dokl)

pkt1 = funkcja(siatka1)
pkt2 = funkcja(siatka2)
pkt3 = funkcja(siatka3)
pkt4 = funkcja(siatka4)
pkt5 = funkcja(siatka5)
pktDokl = funkcja(siatkaDokl)

# i liczenie pochodnych w tych pkt
pochodna1 = pochodna(siatka1)
pochodna2 = pochodna(siatka2)
pochodna3 = pochodna(siatka3)
pochodna4 = pochodna(siatka4)
pochodna5 = pochodna(siatka5)
pochodnaDokl = pochodna(siatkaDokl)

i = 0
pktWylicz1 = []
while i<dokl:
    pktWylicz1.append(szukaj_wart(siatkaDokl[i], siatka1, pkt1, pochodna1))
    i+=1
i = 0
pktWylicz2 = []
while i<dokl:
    pktWylicz2.append(szukaj_wart(siatkaDokl[i], siatka2, pkt2, pochodna2))
    i+=1
i = 0
pktWylicz3 = []
while i<dokl:
    pktWylicz3.append(szukaj_wart(siatkaDokl[i], siatka3, pkt3, pochodna3))
    i+=1
i = 0
pktWylicz4 = []
while i<dokl:
    pktWylicz4.append(szukaj_wart(siatkaDokl[i], siatka4, pkt4, pochodna4))
    i+=1
i = 0
pktWylicz5 = []
while i<dokl:
    pktWylicz5.append(szukaj_wart(siatkaDokl[i], siatka5, pkt5, pochodna5))
    i+=1
i = 0
pktWyliczDokl = []
while i<dokl:
    pktWyliczDokl.append(szukaj_wart(siatkaDokl[i], siatkaDokl, pktDokl, pochodnaDokl))
    i+=1
i = 0


mpl.pyplot.axis([1.075, 1.08, 0, 2500000])
mpl.pylab.plot(siatka1, pkt1, 'ro')
mpl.pylab.plot(siatkaDokl, pktWylicz1)
#mpl.pylab.plot(siatka2, pkt2, 'ro')
#mpl.pylab.plot(siatkaDokl, pktWylicz2)
#mpl.pylab.plot(siatka3, pkt3, 'ro')
#mpl.pylab.plot(siatkaDokl, pktWylicz3)
#mpl.pylab.plot(siatka4, pkt4, 'ro')
#mpl.pylab.plot(siatkaDokl, pktWylicz4)
#mpl.pylab.plot(siatka5, pkt5, 'ro')
#mpl.pylab.plot(siatkaDokl, pktWylicz5)
#mpl.pylab.plot(siatkaDokl, pktWyliczDokl)

#cube = interp1d(siatka3, pkt3, kind='cubic')
#mpl.pylab.plot(siatkaDokl, cube(siatkaDokl))

mpl.pylab.show()

print "norma1               norma2"
print caleczkaSimps(a, b, pktWylicz1, pktWyliczDokl, n1), caleczkaSimps2(a, b, pktWylicz1, pktWyliczDokl, n1)
print caleczkaSimps(a, b, pktWylicz2, pktWyliczDokl, n2), caleczkaSimps2(a, b, pktWylicz2, pktWyliczDokl, n2)
print caleczkaSimps(a, b, pktWylicz3, pktWyliczDokl, n3), caleczkaSimps2(a, b, pktWylicz3, pktWyliczDokl, n3)
print caleczkaSimps(a, b, pktWylicz4, pktWyliczDokl, n4), caleczkaSimps2(a, b, pktWylicz4, pktWyliczDokl, n4)
print caleczkaSimps(a, b, pktWylicz5, pktWyliczDokl, n5), caleczkaSimps2(a, b, pktWylicz5, pktWyliczDokl, n5)