import numpy as np
import matplotlib.pyplot as plt


def phi(t,deltat,l,g,startwinkel):
    '''Als erstes definieren wir die Konstante Tau und erstellen ein 
    array, in welchem wir unsere numerischen Werte speichern. Hierbei
    beschreibt len(t) die Anzahl der Elemente im jeweiligen Array'''
    tau = np.sqrt(2*l/(3*g))
    phiarray = np.zeros(len(t))
    
    '''Nun der erste Schritt für das Zeitschrittverfahren. Die Winkelbeschl. ist
    durch die jeweilige Funktion im Skript gegeben. Hierbei gehen wir von einer 
    Anfangswinkelgeschw.=0.'''
    winkelbeschleunigung = np.sin(startwinkel)/tau**2
    winkelgeschwindigkeit = 0
    phi = startwinkel
    phiarray[0] = np.pi/2 - phi

    '''Nun definieren wir eine Schleife, welche die jeweiligen Werte berechnet und 
    eine Termination der Schleife bewirkt, wenn der Winkel über pi/2 (90Grad) fällt.
    Die jeweilige Prozedur übernehmen wir schlichtweg aus dem Skript.'''
    for i in range(len(t)-1):
        winkelgeschwindigkeit = winkelgeschwindigkeit + deltat*winkelbeschleunigung
        phi = phi + deltat*winkelgeschwindigkeit
        winkelbeschleunigung = np.sin(phi)/tau**2
        phiarray[i+1] = np.pi/2 - phi
        if(phi > np.pi/2):
            phiarray[i+1] = 0
        
    return phiarray


    ''' Danach definieren wir einen Endzeitpunkt, welcher eintritt, sobald im Array 
        das erste Element mit einer 0 gefunden wird.
        Falls die der Endzeitpunkt nicht im Intervall bis tmax liegt, erhält man 
        die Fehlermeldung -1000.'''
        
def t_Aufprall(t,phi):
    for i in range(len(t)):
        if(phi[i] == 0):
            return(t[i])
    return(-1000)
    

def main():
    '''Definieren der Konstanten und des Zeitschrittverfahrens:
        Der startwinkel wird in rad angegeben, somit ist die Umrechnung
        Pi/360Grad notwendig. Delta(t) gibt die Zeitschritte an, in welchen
        die Prozedur vollzogen werden soll. Je kleiner die Zeitschritt, desto
        genauer auch das Ergebnis'''
    l = 1.45
    startwinkel = 0.25
    tmin = 0
    tmax = 1.2
    deltat = 0.00001
    g = 9.81
    t = np.linspace(tmin,tmax,int((tmax - tmin)/deltat))
    
    plt.plot(t,phi(t,deltat,l,g,startwinkel))
    plt.ylabel('Fallwinkel [rad]')
    plt.xlabel('t [s]')
    plt.grid(True)
    plt.show()
  
    
    print('t_Aufprall =',t_Aufprall(t,phi(t,deltat,l,g,startwinkel)),'s')
      
  
if __name__=="__main__":
    main()