#! /usr/bin/env python3
# -*- coding: UTF-8 -*-
# TFY4104/TFY4115 Fysikk 2019.
# Cubic spline approximation to roller coaster on the lab
# Length unit used: mm
import matplotlib.pyplot as plt
import numpy as np
# (x,y)-data in file splinebane.txt (1 header row):
# y = measured height of fixed points of roller coaster
# shown in splinebane.jpg
data=np.loadtxt('splinebane.txt',skiprows=1)
print(data)
xfast=data[:,0]
yfast=data[:,1]
h = 200
b = np.zeros(6)
print(b)
for i in range(0,5):
b[i] = (6/h**2)*(yfast[i+2]-2*yfast[i+1]+yfast[i])
#b[0] = (6/h**2)*(yfast[2]-2*yfast[1]+yfast[0])
#b[1] = (6/h**2)*(yfast[3]-2*yfast[2]+yfast[1])
#b[2] = (6/h**2)*(yfast[4]-2*yfast[3]+yfast[2])
#b[3] = (6/h**2)*(yfast[5]-2*yfast[4]+yfast[3])
#b[4] = (6/h**2)*(yfast[6]-2*yfast[5]+yfast[4])
#b[5] = (6/h**2)*(yfast[7]-2*yfast[6]+yfast[5])
print(b)
A = np.zeros((6,6))
A[0,:] = (4,1,0,0,0,0)
A[1,:] = (1,4,1,0,0,0)
A[2,:] = (0,1,4,1,0,0)
A[3,:] = (0,0,1,4,1,0)
A[4,:] = (0,0,0,1,4,1)
A[5,:] = (0,0,0,0,1,4)
z = np.linalg.solve(A,b)
print(z)
print(A)
ddy = np.zeros(8)
ddy[0] = 0
ddy[7] = 0
ddy[1] = z[0]
ddy[2] = z[1]
ddy[3] = z[2]
ddy[4] = z[3]
ddy[5] = z[4]
ddy[6] = z[5]
print(ddy)
x0=np.linspace(0,200,201)
x1=np.linspace(201,400,200)
x2=np.linspace(401,600,200)
x3=np.linspace(601,800,200)
x4=np.linspace(801,1000,200)
x5=np.linspace(1001,1200,200)
x6=np.linspace(1201,1400,200)
x = np.concatenate([x0,x1,x2,x3,x4,x5,x6])
# x = x0+x1+x2+x3+x4+x5+x6
S0 = ddy[0]*(xfast[1]-x0)**3/(6*h) + ddy[1]*(x0-xfast[0])**3/(6*h) \
+ (yfast[1]/h-ddy[1]*h/6)*(x0-xfast[0]) \
+ (yfast[0]/h-ddy[0]*h/6)*(xfast[1]-x0)
S1 = ddy[1]*(xfast[2]-x1)**3/(6*h) + ddy[2]*(x1-xfast[1])**3/(6*h) \
+ (yfast[2]/h-ddy[2]*h/6)*(x1-xfast[1]) \
+ (yfast[1]/h-ddy[1]*h/6)*(xfast[2]-x1)
S2 = ddy[2]*(xfast[3]-x2)**3/(6*h) + ddy[3]*(x2-xfast[2])**3/(6*h) \
+ (yfast[3]/h-ddy[3]*h/6)*(x2-xfast[2]) \
+ (yfast[2]/h-ddy[2]*h/6)*(xfast[3]-x2)
S3 = ddy[3]*(xfast[4]-x3)**3/(6*h) + ddy[4]*(x3-xfast[3])**3/(6*h) \
+ (yfast[4]/h-ddy[4]*h/6)*(x3-xfast[3]) \
+ (yfast[3]/h-ddy[3]*h/6)*(xfast[4]-x3)
S4 = ddy[4]*(xfast[5]-x4)**3/(6*h) + ddy[5]*(x4-xfast[4])**3/(6*h) \
+ (yfast[5]/h-ddy[5]*h/6)*(x4-xfast[4]) \
+ (yfast[4]/h-ddy[4]*h/6)*(xfast[5]-x4)
S5 = ddy[5]*(xfast[6]-x5)**3/(6*h) + ddy[6]*(x5-xfast[5])**3/(6*h) \
+ (yfast[6]/h-ddy[6]*h/6)*(x5-xfast[5]) \
+ (yfast[5]/h-ddy[5]*h/6)*(xfast[6]-x5)
S6 = ddy[6]*(xfast[7]-x6)**3/(6*h) + ddy[7]*(x6-xfast[6])**3/(6*h) \
+ (yfast[7]/h-ddy[7]*h/6)*(x6-xfast[6]) \
+ (yfast[6]/h-ddy[6]*h/6)*(xfast[7]-x6)
S = np.concatenate([S0,S1,S2,S3,S4,S5,S6])
# S = S0+S1+S2+S3+S4+S5+S6
# Her kan analytiske uttrykk for dS0/dx og d2S0/dx2 legges inn,
# og tilsvarende for S1, ..., S6
dS0 = -ddy[0]*(xfast[1]-x0)**2/(2*h) \
+ ddy[1]*(x0-xfast[0])**2/(2*h) \
+ (yfast[1]/h-ddy[1]*h/6) \
- (yfast[0]/h-ddy[0]*h/6)
dS1 = -ddy[1]*(xfast[2]-x1)**2/(2*h) \
+ ddy[2]*(x1-xfast[1])**2/(2*h) \
+ (yfast[2]/h-ddy[2]*h/6) \
- (yfast[1]/h-ddy[1]*h/6)
dS2 = -ddy[2]*(xfast[3]-x2)**2/(2*h) \
+ ddy[3]*(x2-xfast[2])**2/(2*h) \
+ (yfast[3]/h-ddy[3]*h/6) \
- (yfast[2]/h-ddy[2]*h/6)
dS3 = -ddy[3]*(xfast[4]-x3)**2/(2*h) \
+ ddy[4]*(x3-xfast[3])**2/(2*h) \
+ (yfast[4]/h-ddy[4]*h/6) \
- (yfast[3]/h-ddy[3]*h/6)
dS4 = -ddy[4]*(xfast[5]-x4)**2/(2*h) \
+ ddy[5]*(x4-xfast[4])**2/(2*h) \
+ (yfast[5]/h-ddy[5]*h/6) \
- (yfast[4]/h-ddy[4]*h/6)
dS5 = -ddy[5]*(xfast[6]-x5)**2/(2*h) \
+ ddy[6]*(x5-xfast[5])**2/(2*h) \
+ (yfast[6]/h-ddy[6]*h/6) \
- (yfast[5]/h-ddy[5]*h/6)
dS6 = -ddy[6]*(xfast[7]-x6)**2/(2*h) \
+ ddy[7]*(x6-xfast[6])**2/(2*h) \
+ (yfast[7]/h-ddy[7]*h/6) \
- (yfast[6]/h-ddy[6]*h/6)
dSanal = np.concatenate([dS0,dS1,dS2,dS3,dS4,dS5,dS6])
d2S0 = ddy[0]*(xfast[1]-x0)/(h) \
+ ddy[1]*(x0-xfast[0])/(h)
d2S1 = ddy[1]*(xfast[2]-x1)/(h) \
+ ddy[2]*(x1-xfast[1])/(h)
d2S2 = ddy[2]*(xfast[3]-x2)/(h) \
+ ddy[3]*(x2-xfast[2])/(h)
d2S3 = ddy[3]*(xfast[4]-x3)/(h) \
+ ddy[4]*(x3-xfast[3])/(h)
d2S4 = ddy[4]*(xfast[5]-x4)/(h) \
+ ddy[5]*(x4-xfast[4])/(h)
d2S5 = ddy[5]*(xfast[6]-x5)/(h) \
+ ddy[6]*(x5-xfast[5])/(h)
d2S6 = ddy[6]*(xfast[7]-x6)/(h) \
+ ddy[7]*(x6-xfast[6])/(h)
d2Sanal = np.concatenate([d2S0,d2S1,d2S2,d2S3,d2S4,d2S5,d2S6])
#Plotter de ulike funksjonene
plt.figure('S(x)',figsize=(12,3))
# plt.plot(x0,S0,x1,S1,x2,S2,x3,S3,x4,S4,x5,S5,x6,S6)
plt.plot(x,S)
plt.title('S(x)')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$S(x)$',fontsize=20)
plt.ylim(0,350)
plt.show()
# dS = 1. derivert av S
dS = np.gradient(S,x)
plt.figure('dS/dx')
plt.plot(x,dS,x,dSanal)
plt.title('dS/dx')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$dS/dx$',fontsize=20)
# plt.ylim(0,350)
plt.show()
# Calculate slope angle beta (unit: degrees)
#beta = np.arctan(dSanal)
beta = np.arctan(dSanal)*180/np.pi
#print(len(beta))
#print(beta[0])
# plot angle beta
plt.figure('anglebeta')
plt.plot(x,beta)
plt.title('beta')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel(r'$\beta$',fontsize=20)
# plt.ylim(0,350)
plt.show()
# d2S = 2. derivert av S = 1. derivert av dS
d2S = np.gradient(dS,x)
plt.figure('d2S/dx2')
plt.plot(x,d2S,x,d2Sanal)
plt.title('d2S/dx2')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$d^2 S/dx^2$ (1/mm)',fontsize=20)
# plt.ylim(0,350)
plt.show()
# K = krumningen
K = d2S/(1+dS**2)**(1.5)
Kanal = d2Sanal/(1+dSanal**2)**(1.5)
plt.figure('K(x)')
plt.plot(x,K,x,Kanal)
plt.title('K(x)')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$K(x)$ (1/mm)',fontsize=20)
# plt.ylim(0,350)
plt.show()
# R = 1/K = krumningsradius
R = 1/Kanal
plt.figure('R(x)')
plt.plot(x,R)
plt.title('R(x)')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$R(x)$ (mm)',fontsize=20)
plt.ylim(-5000,5000)
plt.show()
# Check difference between numerical and analytical
# dS and d2S
dSreldiff = np.absolute((dS-dSanal)/dS)
d2Sreldiff = np.absolute((d2S-d2Sanal)/(d2S+1E-12))
plt.figure('dScheck')
plt.plot(x,dSreldiff)
plt.title('Checking dS')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$(dS-dSanal)/dS$',fontsize=20)
plt.show()
plt.figure('d2Scheck')
plt.plot(x,d2Sreldiff)
plt.title('Checking d2S')
plt.xlabel('$x$ (mm)',fontsize=20)
plt.ylabel('$(d2S-d2Sanal)/d2S$',fontsize=20)
plt.show()