#! /usr/bin/env python3
# -*- coding: UTF-8 -*- 
##############################
# pythontips2.py
# Numerisk beregning av hastighet vx og akselerasjon ax fra t og x.
##############################
# Importerer nødvendige funksjoner fra numpy og matplotlib
from numpy import loadtxt,zeros,diff
from matplotlib.pyplot import figure,title,xlabel,ylabel,plot,hold,show,rc,grid

# Leser data fra bane.txt, se pythontips1.py.
data=loadtxt('bane.txt',skiprows=1)
t=data[:,0]
x=data[:,1]
y=data[:,2]

# Det enkleste er å bruke diff(x). Element nr n
# i tabellen diff(x) er x(n+1)-x(n).
# NB: Når x har N elementer, dvs len(x) har verdien N,
# vil diff(x) ha N-1 elementer, dvs len(diff(x)) har verdien N-1.
# Her er tidssteget konstant, så vi regner ut dt en gang for alle.
dt = t[1] - t[0]
vx = diff(x)/dt
ax = diff(vx)/dt
#
# Eller med for-løkker:
N = len(t)
# Initialisering:
vx2 = zeros(N)
ax2 = zeros(N)
# v(i) = Delta x/Delta t = (x[i+1]-x[i])/dt
for i in range(0,N-1):
    vx2[i] = (x[i+1] - x[i])/dt

# setter v(N-1) = (x(N-1)-x(N-2))/dt, dvs vx2[N-1]=vx2[N-2]
vx2[N-1] = vx2[N-2]
# a(i) = Delta v/Delta t = (vx2[i+1]-vx2[i])/dt
for i in range(0,N-1):
    ax2[i] = (vx2[i+1] - vx2[i])/dt

# setter a(N-1) = (v(N-1)-v(N-2))/dt, dvs ax2[N-1]=ax2[N-2]
# (siden vi satte vx2[N-1]=vx2[N-2], blir nå ax2[N-2]=ax2[N-1]=0)
ax2[N-1] = ax2[N-2]

# Plotter vx(t) og vx2(t)
figure('Figur 1')
title('Figur 1. Klossens hastighet på skråplanet')
hold(True)
xlabel(r'Tid (s)')
ylabel(r'Hastighet (cm/s)')
plot(t[0:len(vx)],vx)
plot(t,vx2)
grid()
show()

# Plotter ax(t) og ax2(t)
figure('Figur 2')
title('Figur 2. Klossens akselerasjon på skråplanet')
hold(True)
xlabel(r'Tid (s)')
ylabel(r'Akselerasjon (cm/s$^2$)')
plot(t[0:len(ax)],ax)
plot(t,ax2)
grid()
show()