# -*- coding: utf-8 -*-

"""
Created on Mon Nov 27 13:38:08 2017

@author: Ptaaq
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

"""
### dane wejsciowe ###
"""
beta = 1
gamma = 10
N = 1000

"""
### warunki początkowe ###
"""
x0 = 900
y0 = 90
z0 = 10
init = [x0, y0, z0] #wektor z warunkami początkowymi, co jest przekazywany do funkcji ODEINT

### czas symulacji ###
n = 1001 #liczba punktów
t = np.linspace(0,70,n) #wektor czasu - punkty w których będzie obliczana symulacja

wart_zadana = 10
par = [beta,gamma,N,wart_zadana] #parametry



def moj_PID(wartosc_rzeczywista, wartosc_zadana, Kp, Ki, Kd, wartosc_I_pocz = 0, wartosc_epsilon_pocz = 0):
    """
    argumenty funkcji:
        - wartosc_rzeczywista - aktualna wartosc mierzonej/obliczanej zmiennej
        - wartosc_zadana - wartosc jaka chcemu uzyskac
        - Kp - wspolczynnik wzmocnienia czlonu proporcjonalnego
        - Ki - wspolczynnik wzmocnienia czlonu calkujacego
        - Kd - wspolczynnik wzmocnienia czlonu rozniczkujacego
    """
    k = 0
    #inicjalizacja
    wartosc_I_pop = wartosc_I_pocz
    epsilon_pop = wartosc_epsilon_pocz

    while 1:
        epsilon = (wartosc_zadana - wartosc_rzeczywista) #uchyb regulacji
        wartosc_P = Kp*epsilon                          #odpowiedz czlonu proporcjonlanego
        wartosc_I = wartosc_I_pop + Ki*epsilon          #odpowiedz czlonu calkujacego
        wartosc_D = Kd*(epsilon-epsilon_pop)            #odpowiedz czlonu rozniczkujacego
    #przypisanie poprzednich wartosci
    epsilon_pop = epsilon
    wartosc_I_pop = wartosc_I
    #oczliczanie wyjscia systemu
    PID = wartosc_P  + wartosc_I + wartosc_D
    k += 1
    return PID

def model(wektor_stanu,t,par):       #funkcja model
#   print(t)
   x = wektor_stanu[0]
   y = wektor_stanu[1]
   z = wektor_stanu[2]
   wart_obliczana = (1/gamma)*y
   //*
   regulator = moj_PID(wart_obliczana,wart_zadana,1,1,1)
   uchyb = wart_zadana - wart_obliczana
   if (uchyb > 0):
       gamma = gamma + regulator
   else:
       gamma = gamma - regulator
   *//
   x_prim = (-beta/N)*x*y
   y_prim = (beta/N)*x*y - wart_obliczana
   z_prim = (1/gamma)*y
   return [x_prim, y_prim, z_prim]

solution = odeint(model,init,t,args=(par,))

plt.plot(t, solution[:,0],'b', label='x - narazeni')
plt.plot(t, solution[:,1],'g', label='y - nosiciele')
plt.plot(t, solution[:,2],'r', label='z - uodpornieni')