import numpy as np   
import matplotlib.pyplot as plt

def dv_dt(v): # definisanje diferencijalne jednačine
    return g-(k/m)*v**2

def v_ex(t): # definisanje egzaktnog resenja dif. jednacine
    return np.sqrt((m*g)/k)*np.tanh(g*t*(np.sqrt((k)/(m*g))))
      
m = 80 
k = 0.24
g = 9.81
t0 = 0 # početni uslovi  
v0 = 0 # početni uslovi  

tf = 30 # krajnja vrednost t, vrednost za koju zelimo da izracunamo vrednost jednačine  
n = 31 # broj tačaka na koje se deli interval od "t0" do "tf"  
l = n - 1 # broj podintervala (koraka) između "t0" i "tf"  

t = np.linspace(t0,tf,n) # interval od "t0" do "tf" deli na "n" tačaka  
v = np.zeros(n) # pravi skup od "n" clanova, pri cemu je svaki clan nula  
  
h = (tf-t0)/l # ovako se definiše "h" na osnovu "np.linspace" komande  
              # a može i da se definiše "h" kao "h=x[1]-x[0]"   
              # aproksimativno rešenje, Ojlerov metod  
for i in range (0,n-1):  
    t[0] = t0  
    v[0] = v0  
    v[i+1]= v[i] + h * dv_dt(v[i])  

v_exact = v_ex(t)
podaci=list(zip(t,v,v_exact))

#------------------------
#Crtanje grafika
#------------------------
plt.plot(t,v,label="Aproksimativno resenje")
plt.plot(t,v_exact,label="Egzaktno resenje")
plt.title("Slobodan pad sa silom otpora sredine")
plt.xlabel("Vreme [s]")
plt.ylabel("Brzina [m/s]")
plt.legend()
plt.show()