# -*- coding: utf-8 -*-
"""
Created on Tue Jan  5 15:49:07 2021

@author: bjarne
"""

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

#Initialbedingungen
phi0=0.25                   #Anfangswinkel in rad
l=1.45                      #Länge des Stabes in m
g=9.81                      #erdbeschleunigung in m/s^2

#Definition Parameter tau

def tau():
    
    tau = math.sqrt(2*l/(3*g))
    
    return tau

#Solver
    
def solve(phi0,tmax,steps):
    
    sol = list()
    
    sol.append([0,phi0,0,math.sin(phi0) / (tau()**2)])
    
    tstep = tmax / steps
    
    for e in range(steps):
        ddotphi = math.sin(sol[len(sol)-1][1]) / (tau()**2)
        
        dotphi = sol[len(sol)-1][2] + tstep * sol[len(sol)-1][3]
        
        phi = sol[len(sol)-1][1] + tstep * sol[len(sol)-1][2]
        
        sol.append([sol[len(sol)-1][0]+tstep,phi,dotphi,ddotphi])
        
        if phi>=math.pi/2:
            
            break
        
    return np.array(sol)
    


#Berechnung der Fallzeit für verschiedene Anfangswinkel
philist=np.linspace(0.2,1.2,100,endpoint=False)       
def falltime():
    
    timelist = list()
    
    for i in range(len(philist)):
        
        timelist.append([philist[i],solve(philist[i],10,10000)[len(solve(philist[i],10,10000))-1][0]])
        
    return np.array(timelist)

    

    
#Erstellung des Plots        
path='D:/Studium/Praktikum2python'   	#Speicherpfad     
#Achseneinstellungen
xstart=0		#Anfang der x-Achse
xend=1.4		#Ende der x-Achse
xticks=xend/7		#Unterteilung der Achsen
ystart=0		#Anfang der y-Achse
yend=1			#Ende der y-Achse
#Bildeinstellungen
plt.figure(figsize=(20,10), dpi=100)
plt.grid(True)
plt.axis([xstart,xend,ystart,yend])
plt.title('Lösung für den Winkel')                              
plt.xticks(np.arange(xstart,xend,xticks))


plt.plot(falltime()[:,0], falltime()[:,1], markersize=0.5)	#zu plottendes Array
plt.legend(['numerische Lösung'])        			#Legende
plt.xlabel('$\phi$ in rad')                              	#Beschriftung x-Achse
plt.ylabel('Zeit in s')                              		#Beschriftung y-Achse
plt.savefig(path+'/timeplot.jpg')               		#Dateiname
plt.close