import numpy as np
import matplotlib.pyplot as plt

#define the constants
l = 1.45 #length of the broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#define the second derivation of phi
def dd_phi(phi):
    a = np.sin(phi) / (tau ** 2)
    return a

#define the algorythmen
def iteration(phi0, t0, d_phi0, dd_phi, steps_count, steps_wide):
    values = [[phi0, t0, d_phi0]]
  
    i = 0
    t = t0
    while i <= steps_count:
        
        next_d_phi = values[i][2] + dd_phi(values[i][0]) * steps_wide
        
        next_phi = values[i][0] + values[i][2] * steps_wide

        #break, when stick reaches the ground
        if next_phi >= (np.pi / 2.0):
            break
        
        t = values[i][1] + steps_wide
        
        new_values = [next_phi, t, next_d_phi]
        values.append(new_values)
        
        i += 1 
    return(values) 

#define function to reproduce graphic 3
def time_ankle():
    c = 0.2 #starting ankel
    ankel = [] #empty list for ankels
    time = [] #empty list for times
    while c <= 1.4:
        iterate_values = iteration(c, 0, 0, dd_phi, 10 ** 7, 10 ** -2)
        ankel.append(c)
        
        new_time = [i[1] for i in iterate_values]
        last_new_time = new_time[-1] #last time in list is time needed to fall
        time.append(last_new_time)
        
        c = c + 0.05
    return(ankel, time)

x, y = time_ankle()

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_title('Fallzeit')#title of the plot
ax.set_xlabel('Startwinkel in rad')#label of x axes
ax.set_ylabel('Zeit in s')#label of y axes
ax.set_xlim(0, 1.4) #interval of x value
ax.set_ylim(0, 1) #inerval of y value
ax.grid()
ax.plot(x, y)
plt.show()


#values for plotting the numerical solution and experimental values
#gives the needed time to fall for starting ankle
def diagram_points(phi0):
    values = iteration(phi0, 0, 0, dd_phi, 10 ** 7, 10 ** -2)
    
    new_time = [i[1] for i in values]
    last_new_time = new_time[-1] 
    diagramm = [phi0, last_new_time]

    return (diagramm)

#define the new constants for first broomstick
l = 1.285 #length of the first broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
d = diagram_points(0.1308996938995747)
e = diagram_points(0.2617993877991494)
f = diagram_points(0.4363323129985824)
g = diagram_points(0.6981317007977318)
h = diagram_points(0.8726646259971648)
i = diagram_points(1.0471975511965976)
j = diagram_points(1.2217304763960306)

print("Messwerte erster Besen:")
print(d, e, f, g, h, i, j)

#uncertainties for first broomstick
l = 1.284 #minimal length of the first broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
dmin = diagram_points(0.1308996938995747)
emin = diagram_points(0.2617993877991494)
fmin = diagram_points(0.4363323129985824)
gmin = diagram_points(0.6981317007977318)
hmin = diagram_points(0.8726646259971648)
imin = diagram_points(1.0471975511965976)
jmin = diagram_points(1.2217304763960306)

print("Minimum erster Besen:")
print(dmin, emin, fmin, gmin, hmin, imin, jmin)


l = 1.286 #maximal length of the first broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
dmax = diagram_points(0.1308996938995747)
emax = diagram_points(0.2617993877991494)
fmax = diagram_points(0.4363323129985824)
gmax = diagram_points(0.6981317007977318)
hmax = diagram_points(0.8726646259971648)
imax = diagram_points(1.0471975511965976)
jmax = diagram_points(1.2217304763960306)

print("Maximum erster Besen:")
print(dmax, emax, fmax, gmax, hmax, imax, jmax)


#define the new constants for second broomstick
l = 1.397 #length of the second broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
k = diagram_points(0.1308996938995747)
l = diagram_points(0.2617993877991494)
m = diagram_points(0.4363323129985824)
n = diagram_points(0.6981317007977318)
o = diagram_points(0.8726646259971648)
p = diagram_points(1.0471975511965976)
q = diagram_points(1.2217304763960306)

print("Messwerte zweiter Besen:")
print(k, l, m, n, o, p, q)

#uncertainties for second broomstick
l = 1.396 #minimal length of the second broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
kmin = diagram_points(0.1308996938995747)
lmin = diagram_points(0.2617993877991494)
mmin = diagram_points(0.4363323129985824)
nmin = diagram_points(0.6981317007977318)
omin = diagram_points(0.8726646259971648)
pmin = diagram_points(1.0471975511965976)
qmin = diagram_points(1.2217304763960306)

print("Minimum zweiter Besen:")
print(kmin, lmin, mmin, nmin, omin, pmin, qmin)

l = 1.398 #maximal length of the second broomstick
g = 9.81 
tau = np.sqrt((2 * l) / (3 * g))

#same ankels as in the experiment
kmax = diagram_points(0.1308996938995747)
lmax = diagram_points(0.2617993877991494)
mmax = diagram_points(0.4363323129985824)
nmax = diagram_points(0.6981317007977318)
omax = diagram_points(0.8726646259971648)
pmax = diagram_points(1.0471975511965976)
qmax = diagram_points(1.2217304763960306)

print("Maximum zweiter Besen:")
print(kmax, lmax, mmax, nmax, omax, pmax, qmax)