# -*- coding: utf-8 -*-
"""
@author: Vaclav Miklas

Nazev: FERMENTACE
===============================================================================
UVOD

Tento skript obsahuje materialovou a energetickou bilanci fermentace v BPS.
Pokud je cely tento skript naimportovan, obsazene funkce jsou volany zpusobem:
fermentace.MatBilance(argumenty) resp. fermentace.EnerBilance(argumenty)
–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
POKYNY K FUNKCI MatBilance

Argumenty:
df ... pandas dataframe obsahujici informace o substratech (nezbytne hmotnostni 
       toky, obsahy susiny, vyteznost bioplynu) - vstupni xlsx soubor nemusi mit
       ustalenou formu, ale musi byt dodrzeny a v nazvu sloupce uvedeny
       jednotky stezejnich velicin (kg/h, % resp. m3/kg suš.)
recykl ... procenta recyklovaneho digestatu zpet pred fermentor (ockovani
           substratu)
zadano ... vstupni parametr zavisejici na zvolenem modu, jednotky viz nize
mod ... argument na poslednim miste, ve tvaru mod=hodnota; hodnota 1 znamena 
        vypocet dle zadaneho obsahu susiny pred fermentorem (zadano v %), 
        hodnota 2 vychazi ze zadaneho obsahu susiny za fermentorem (zadano v %),
        hodnota 3 pocita na zaklade zadaneho pridavku vody (zadano v kg/h),
        neni-li zadano, je uvazovan mod c. 1.

Priklad volani:
    MatBilance(df, recykl, zadano, mod=X)

Pozn. Vysledky jsou ve forme 14-clenneho pole, prvnich sedm jsou hmotnostni 
      toky [kg/h], dalsich sedm obsahy susiny [–].
–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
POKYNY K FUNKCI EnerBilance 

Argumenty spolecne:
vysledkyMBF ... 14-clenne pole vysledku funkce MatBilance
hod ... pocet hodin, po kterou je denne fermentor michan; na zaklade tohoto
        tohoto udaje se stanovi prumerny prikon michaciho zarizeni (prikon, 
        ktery by pri setrvalem michani daval totoznou spotrebu)
mod ... argument na predposlednim ci poslednim miste, ve tvaru mod=hodnota; 
        hodnota 1 znamena vypocet tepelnych ztrat dle prenosove rovnice,
        hodnota 2 rovnez vyuziva prenosove rovnice, avsak plocha
        vyplyva ze zadane stredni zdrzne doby a pomeru H/D fermentoru,
        hodnota 3 pocita dle procenta vyrobeneho tepla; neni-li zadano, pocita  
        jako 20 % z vyrobeneho tepla a dalsi argumenty nejsou brany v potaz
mich ... argument na predposlednim ci poslednim miste, ve tvaru mich=hodnota;
         hodnota 1 znamena hydraulicke michani, hodnota 2 znamena michani
         mechanickymi (vrtulovymi) michadly, hodnota 3 znamena nemichany
         fermentor; neni-li zadano, je uvazovano hydraulicke michani

Zapis argumentu pro mod = 1:   
    EnerBilance (vysledkyMBF, hod, k, A, T_in, T_out, mod=1, mich=X), kde:
        k ... (prumerny) soucinitel prostupu tepla [W/(m^2·K)]
        A ... plocha fermentoru vystavena vnejsimu prostredi [m^2]
        T_in ... teplota fermentace [°C]
        T_out ... teplota vnejsiho prostredi [°C]
        
Zapis argumentu pro mod = 2:   
    EnerBilance (vysledkyMBF, hod, k, H-D, strDob, T_in, T_out, mod=2, mich=X) kde:
        k ... (prumerny) soucinitel prostupu tepla [W/(m^2·K)]
        H_D ... pomer vyska/prumer fermentoru [–]
        strDob ... stredni zdrzna doba substratu ve fermentoru [den]
        T_in ... teplota fermentace [°C]
        T_out ... teplota vnejsiho prostredi [°C]

Zapis argumentu pro mod = 3:
    EnerBilance (vysledkyMBF, hod, vysledkyKGJ, tep_proc, mod=3, mich=X), kde:
        vysledkyKGJ ... 5-clenne pole vysledku funkce KGJ.EnerBilance
        tep_proc ... procenta vyrobeneho tepla vyuzivane na vytapeni fermentoru [%]
        
Pozn. Vysledky EnerBilance jsou ve forme dvou clenu - prvni vyjadruje 
      elektrickou spotrebu v kW, druhy tepelnou spotrebu, taktez v kW.
      Je zrejme, ze optimalni posloupnost vyuziti funkci je MatBilance, KGJ,
      EnerBilance. Ta umozni pouzit vsechny mody funkce. 
    
===============================================================================
"""

from scipy.optimize import fsolve
# from scipy import interpolate # odkomentovat, pokud budou pouzity radky 215-236
import numpy as np

def MatBilance(df, zadano, recykl, **kwargs):
    # DEKLARACE PROMENNYCH:
    mod = kwargs.get('mod')
    rhoBP = 1.2 # kg/m3
    indexList=[]
    tokyList=[]
    susinyList=[]
    vytezkyList=[]
    if not mod:
        print ('Nebyl zadán mód MatBilance, druhá proměnná je proto brána ' + 
               'jako obsah sušiny na vstupu do fermentoru.')
    
    # ZPRACOVANI DAT:
    # Pro nenulove hmotnostni toky jsou vyctena data z vstupniho dataframe:
    for i in range(df.shape[0]):
        if df.filter(like='kg/h').values[i,0] > 0:
            # Indexy substratu, ktere vstupuji do vypoctu:
            indexList.append(i)
            # Hmotnostni toky jednotlivych substratu [kg/h]:
            tokyList.append(df.filter(like='kg/h').values[i,0])
            # Hmotnostni toky susiny v jednotlivych substratech [kg/h]:
            susinyList.append(df.filter(like='%').values[i,0]/100*tokyList[-1])
            # Vytezky bioplynu z jednotlivych substratu [m3/h]:
            vytezkyList.append(df.filter(like='m3/kg suš.').values[i,0]*susinyList[-1])
    
    # FUNKCE K RESENI MATERIALOVE BILANCE:
    # Pomocna funkce, ktera casti promennych priradi konstanty (ze zadani), 
    # nezname promenne jsou pak prirazeny k poli z. Hlavni vyhoda je, ze
    # vycteni vysledku je pak elegantne ucineno pomoci teto funkce, nikoliv
    # jednotlive:
    def Prirazeni(z):
        # m_i...hmotnostni tok i-teho proud, x_i...odpovidajici obsah susiny
        m1 = sum(tokyList) # vstupni substrat
        m3 = z[0] # vstup do fermentoru
        m4 = sum(vytezkyList)/0.62*rhoBP # produkovany bioplyn
        m5 = z[1] # digestat z fermentoru
        m6 = z[2] # vystupni digestat
        m7 = z[3] # recykl digestatu
        x1 = sum(susinyList)/sum(tokyList)
        x2 = 0
        x4 = 1 # Predpoklad, ze veskery ubytek susiny prechazi do bioplynu.
        # Nezname a zadane parametry se lisi dle zvoleneho modu:
        if mod == 2: # 2, Vstupem susina za fermentorem         
            m2 = z[4]
            x3 = z[5]
            x5 = zadano/100
        elif mod == 3: # 3, Vstupem mnozstvi vody
            m2 = zadano # pridavek vody
            x3 = z[4]
            x5 = z[5]

        else: # 1, (ci nezadano) Vstupem susina pred fermentorem
            m2 = z[4]
            x3 = zadano/100
            x5 = z[5]
        x6 = x5
        x7 = x5
        return m1,m2,m3,m4,m5,m6,m7,x1,x2,x3,x4,x5,x6,x7
    # Funkce, ktera zapisuje soustavu rovnic k reseni. Nejprve se uplatni
    # nomenklatura z predchozi funkce tak, aby rovnice mohly byt zapsany
    # srozumitelnym zpusobem. Tyto rovnice jsou zapsany v homogennim tvaru
    # a jejich vycisleni (formou pole F) je vystupem funkce:
    def SoustavaRovnic(z):
        m1,m2,m3,m4,m5,m6,m7,x1,x2,x3,x4,x5,x6,x7 = Prirazeni(z)
        F = [0,0,0,0,0,0]
        F[0] = m1 + m2 + m7 - m3
        F[1] = m3 - m4 - m5
        F[2] = m5 - m6 - m7
        F[3] = m7 - recykl/100*m5
        F[4] = m1*x1 + m2*x2 + m7*x7 - m3*x3
        F[5] = m3*x3 - m4*x4 - m5*x5
        return F
    
    # POCATECNI ODHADY (lisi se dle modu):
    if mod == 3:
        # Odhad vsech proudu od vstupu do fermentoru dale je proveden na zaklade
        # toku substratu na vstupu a predepsaneho mnozstvi pridavku vody:
        odh = sum(tokyList) + zadano
        # Odhad susiny vychazi z podilu vstupniho toku susiny a vyse uvedeneho
        # odhadu:
        odh_sus = sum(susinyList)/odh
        # Celkovy zapis odhadu. Uplatnuje se zde i vstupni parametr % recyklu. 
        # Odhad mnozstvi susiny za fermentorem je ztotozneno s odhadem pred
        # fermentorem:
        zOdhad = [odh,odh,odh,odh*recykl/100,odh_sus,odh_sus]
        # Odhady koresponduji s m3,m5,m6,m7,x3,x5
    else:
        # Odhad vsech proudu od vstupu do fermentoru dale je proveden na zaklade
        # toku susiny na vstupu a predepsaneho mnozstvi susiny pro cerpani:
        odh = sum(susinyList)/zadano*100
        # Odhad pridavku vody je rozdilem vyse uvedeneho odhadu a vstupniho toku
        # substratu:
        odh_voda = odh - sum(tokyList)
        # Celkovy zapis odhadu. Plati prakticky totez co vyse, je vsak odhadovana
        # pouze jedna ze susin pred a za fermentorem, druha je zadana.
        zOdhad = [odh,odh,odh,odh*recykl/100,odh_voda,zadano/100]
        # Odhady koresponduji s m3,m5,m6,m7,m2,x3/5  

    # RESENI MATERIALOVE BILANCE:    
    # Resitel fsolve z package scipy nalezne takove reseni, ze F je polem nul,
    # tedy plati vsechny rovnice zapsane vyse:
    zReseni = fsolve(SoustavaRovnic,zOdhad)
    # Zapis vysledku do pole celkoveho prehledu materialove bilance:
    vysledkyList = Prirazeni(zReseni)
    # Pokud bude zaporny tok vody (odpovida dehydrataci substratu), bude
    # nezavisle na hodnotu argumentu 'zadano' a zvoleny mod provedena
    # rekalkulace uvazujici nulovy pridavek vody:
    if vysledkyList [1] < 0:
        print ('Pro dané parametry vychází dehydratace substrátu (záporný tok vody). '+
               'Zadané hodnoty budou ignorovány a tok vody vzatý jako nulový.')
        mod = 3
        zadano = 0
        vysledkyList = Prirazeni(fsolve(SoustavaRovnic,zOdhad))
    # ...ktere je zaroven vystupem funkce:
    return vysledkyList

def EnerBilance (vysledkyMBF, hod, *args, **kwargs):  
    # DEKLARACE PROMENNYCH A FUNKCI:
    # Nacteni promennych, ktere se poji s klicovym slovem (mod tepelne bilance
    # a mod michani).
    mod = kwargs.get('mod')
    mich = kwargs.get('mich')
    cpW = 4.18
    cpDM = 1
    
    # Funkce pro tepelnou bilanci v pripade vyuziti prenosove rovnice:
    def Teplo1 (k, A, T_in, T_out):
        # Prenosova rovnice:
        tepelna_ztrata =  k * A * (T_in - T_out)/1000
        # Z DUVODU NADHODNOCENYCH CISEL NENI NASLEDUJICI POSTUP POUZIT:
        """
        # Data pro vypocet tepla potrebneho pro ohrev substratu:
        xdata = np.array([3, 4, 5, 6, 7, 8, 9, 10, 11, 12]) # % susiny
        ydata = np.array([1.48, 1.11, 0.89, 0.74, 0.64, 0.56, 0.49, 0.445, 0.4, 0.37]) # kWh/kg susiny
        # Interpolace dat kubickym splinem:
        f = interpolate.InterpolatedUnivariateSpline(xdata,ydata)
        # Pomocna promenna. V podstate odpovida velicine sus v MB:
        sus = vysledkyMBF[9]*100
        # Osetreni stavu nad a pod rozsahem platnosti interpolace. Algoritmus
        # neextrapoluje, nybrz pri hodnotach mimo vyuziva krajnich hodnot.
        # Zaroven do konzole vypise hlasku, aby byl uzivatel srozumen.
        if sus < 3:
                sus = 3
                print ('Obsah sušiny leží pod rozsahem platnosti interpolace energie na ohřev. ')
                print ('Spotřeba tepla na ohřev bude uvažována pro 3 % sušiny.')
        if sus > 12:
                sus = 12
                print ('Obsah sušiny leží nad rozsahem platnosti interpolace energie na ohřev. ')
                print ('Spotřeba tepla na ohřev bude uvažována pro 12 % sušiny.')
        # Ohrevu substratu [kW] = tok susiny [kg/h] * merna potreba [kWh/kg sus]        
        ohrev_substratu = vysledkyMBF[2] * vysledkyMBF[9] * f(sus)
        """
        # MISTO TOHO JE OHREV POCITAN JAKO ROZDIL ENTALPII 
        # (predpoklad T_subs = T_out):
        cp = vysledkyMBF[9]*cpDM + (1-vysledkyMBF[9])*cpW
        ohrev_substratu = cp * vysledkyMBF[2]/3600 * (T_in-T_out)
        
        # Celkova tepelna spotreba a vystup funkce:
        teplo = tepelna_ztrata + ohrev_substratu
        return teplo
    
    def Teplo2 (KGJ, tep_proc):
        # Vypocet spotreby jako daneho procenta vyrobene tepelne energie:
        teplo = tep_proc * KGJ[1] / 100
        # Vystup funkce:
        return teplo

    # Funkce pro bilanci elektricke energie:
    def Elektrina ():
        # Denni tok susiny do fermentoru v tunach:
        sus_den = vysledkyMBF[2] * vysledkyMBF[9] * 24 / 1000
        # Mod mechanickeho michani:
        if mich == 2:
            # Aproximacni funkce pro mechanicke michani:
            h = lambda x: -0.00563764*x**2 + 1.61550759*x + 0.95956729
            # Osetreni stavu nad a pod rozsahem platnosti aproximace. Algoritmus
            # pri hodnotach mimo vyuziva krajnich hodnot. Zaroven do konzole
            #  vypise hlasku, aby byl uzivatel srozumen.
            if sus_den <= 1:
                sus_den = 1
                print ('Tok sušiny leží pod rozsahem platnosti aproximace energie na míchání. ')
                print ('Spotřeba elektrické energie míchání bude uvažována pro 1 t VL/den.')
            elif sus_den >= 100:
                sus_den = 100
                print ('Tok sušiny leží nad rozsahem platnosti aproximace energie na míchání. ')
                print ('Spotřeba elektrické energie míchání bude uvažována pro 100 t VL/den.')
            # Totez co vyse, avsak s uplatnenim odlisne aproximacni funkce.
            E = h(sus_den)*hod/24
            return E
        # Mod bez michani:
        elif mich == 3:
            E = 0
            return E
        # Pokud nebude zvoleno mechanicke ci zadne michani, vzdy bude pouzito hydraulicke:
        else:
            # Pokud vsak mod nebude zvolen vubec ci spatne, do konzole bude
            # vypsana hlaska:
            if mich != 1:
                print ('Nebyla zadána validní metoda míchání, je proto uvžováno hydraulické.')
            # Aproximacni funkce pro hydraulicke michani:
            g = lambda x: 8.83331056*np.log(x) - 0.05624127*x + 3.3710468
            # Osetreni stavu nad a pod rozsahem platnosti interpolace. Viz Teplo1.
            if sus_den <= 1:
                sus_den = 1
                print ('Tok sušiny leží pod rozsahem platnosti aproximace energie na míchání. ')
                print ('Spotřeba elektrické energie míchání bude uvažována pro 1 t VL/den.')
            elif sus_den >= 100:
                sus_den = 100
                print ('Tok sušiny leží nad rozsahem platnosti aproximace energie na míchání. ')
                print ('Spotřeba elektrické energie míchání bude uvažována pro 100 t VL/den.')
            # Vypocet prumerne spotreby elektricke energie pri uvazeni, ze
            # denne jsou pouze "hod" hodin v chodu a vystup funkce.
            E = g(sus_den)*hod/24
            return E

    # ZASTRESUJICI PROCEDURA:
    # Tato procedura ma za ukol overit, zda byl zadan spravny pocet argumentu
    # pro dany mod a nasledne provest kalkulaci odpovidajici danemu modu. 
    
    # Prvni mod vyuziva prenosovu rovnici, pri cemz A je primo vstupem funkce.
    if mod == 1:
        # Chybova hlaska pro pripad, ze nesedi pocet argumentu:
        if len(args) != 4:
            print ('Pro daný mód je vyžadováno 6 argumentů pro EnerBilance.')
        # Pokud je pocet argumentu v poradku, je zavolana f-ce Teplo1:
        else:
            # Prirazeni vstupu funkce EnerBilance:
            k = args[0]
            A = args[1]
            T_in = args[2]
            T_out = args[3]
            Q = Teplo1(k, A, T_in, T_out)
    # Druhy mod rovnez vyuziva prenosovu rovnici, avsak plochu dopocitava
    # ze stredni zdrzne doby substratu a geometricke charakteristiky fermentoru:
    elif mod == 2:
        # Totez, co vyse.
        if len(args) != 5:
            print ('Pro daný mód je vyžadováno 7 argumentů pro EnerBilance.')
        else:
            # Prirazeni vstupu funkce EnerBilance:
            k = args[0]
            H_D = args[1]
            strDob = args[2]
            T_in = args[3]
            T_out = args[4]
            # Objem je v projekcni praxi pocitavan jako nasobek denniho toku
            # substratu a stredni zdrzne doby:
            V = vysledkyMBF[2]/1000*24*strDob
            # Z rovnic pro povrcha objem valce lze, s prihlednutim ke geom.
            # charakteristice H/D, vyjadrit nasledujici:
            A = V**(2/3)*np.pi**(1/3)*((16*H_D)**(1/3)+(2*H_D)**(-2/3))
            # Pozn. Je uvazovan pouze plast valce a horni podstava. K jistym 
            # tepelnym ztratam vsak samozrejme dochazi i pres dno fermentoru. 
            # Zaroven je potreba mit na pameti, ze soucinitel prostupu tepla u 
            # obvodovych sten a u stropu (plynojemu) jsou rozdilne. Duvodem 
            # jsou jak odlisne materialy, tak rozdilna media, ktera jsou v 
            # kontaktu se stenou (bioplyn ci substrat).
            # Volani funkce Teplo1:
            Q = Teplo1(k, A, T_in, T_out)
    # Treti mod pouziva zjednoduseny vypocet Teplo2:
    elif mod == 3:
        if len(args) != 2:
            print ('Pro daný mód jsou vyžadovány 4 argumenty pro EnerBilance.')
        else:
            # Prirazeni vstupu funkce EnerBilance:
            KGJ = args[0] # doplnit dle charakteru vystupu KGJ!!!
            tep_proc = args[1]
            # Volani funkce:
            Q = Teplo2(KGJ, tep_proc)    
    # Pokud mod neni zadan, nebo je zadan chybne, je pouzit v podstate mod c.3,
    # spotreba vyrobeneho tepla vsak cini 20 % a neni promennou.
    else:
        if len(args) != 1:
            print ('Není-li (správně) zadán mód, jsou vyžadovány 3 argumenty pro EnerBilance.')
        else:
            print ('Nebyl (správně) zvolen mód EnerBilance. Spotřeba tepelné energie je počítána jako 20 % vyrobeného tepla.')
            KGJ = args[0] # doplnit dle charakteru vystupu KGJ!!!
            tep_proc = 20
            Q = Teplo2(KGJ, tep_proc)
    # Volani funkce Elektrina.
    E = Elektrina ()
    # Vystupem funkce jsou spotreby elektriny a tepla. Zaporne vzate jsou prave
    # z duvodu, ze se jedna o spotreby.
    return -E, -Q