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

Nazev: MSF
===============================================================================
UVOD

Tento skript  obsahuje vypocetni model vicestupnove mzikove odparky (MSF) doplneny
o algoritmus, ktery MSF integruje do BPS na zaklade daneho odpadniho tepla.
Obsahuje jedinou funkci, ktera je v pripade importu celeho skriptu volana
nasledovne: MSF.MSF(argumenty).
–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
POKYNY K FUNKCI MSF:

Argumenty:
Q ... mnozstvi odpadniho tepla z BPS, ktere je k dispozici [kWt]
vysledkyMBS ... 4-clenne pole sledku funkce separace.MatBilance
R ... pomer chladici vykon/spotreba el. energie; predstavuje ucinnost
      cirkulacniho chladice, viz chlazeni.py
      
Priklad volani:
    MSF.MSF(Q, vysledkyMBS, R)

Pozn. Vysledky EnerBilance jsou ve forme 13-clenneho pole. Prvnich osm clenu
      obsahuje postupne hmotnostni tok m [kg/h] a susinu x [hm.%] ctyr 
      hlavnich proudu v MSF (v poradi fugat, recirkulat, koncentrat, destilat), 
      dale nasleduje tepelna spotreba, velikost teplosmenne plochy ohrevneho
      vymeniku, velikost teplosmenne plochy kondenzacnich sekci, pocet komor 
      chlazenych vodou a nakonec spotreba elektricke energie.
–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
"""

from iapws import IAPWS97
from math import log as ln
from scipy.optimize import newton
from scipy.optimize import minimize_scalar

def MSF(Q,vysledkyMBS,R):
#============================ DEKLARACE PROMENNYCH ============================
    Qdane = Q # [kWt] ... mnozstvi disponibilniho tepla
    nCh = 12 # [-] pocet komor je dan
    
    cpW = 4.18 # [kJ/kg·K] ... merna tepelna kapacita vody
    cpDM = 1 # [kJ/kg·K] ... uvazovana merna tepelna kapacita susiny
    
    mFug = vysledkyMBS[1] # [kg/h] ... tok fugatu 
    TFug = 40 # [°C] ... teplota fugatu
    susFug = vysledkyMBS[3]*100 # [%] ... obsah susiny fe fugatu
    susMax = 15 # [%] ... maximalni obsah susiny v koncentratu
    
    TRecOut = 40 # [°C] ... zvolena vystupni teplota recirkulatu z posledni komory, 
                 # odpovida absolutnimu tlaku priblizne 7 kPa          
    ThwIn = 90 # [°C] ... teplota vstupni topne vody do ohrevneho vymniku
    ThwOut = 70 # [°C] ... teplota vystupni topne vody z ohrevneho vymeniku
    Uhw = 0.8 # [kW/m2·K] ... soucinitel prostupu tepla v ohrevnem vymeniku
    Ucond = 1.2 # [kW/m2·K] ... soucinitel prostupu tepla pri kondenzaci (fugat-voda)
    BPE = 0.5 # [°C] ... odhad teplotni deprese
    TCWIn = 25 # [°C] ... vstupni teplota chladici vody
    TCWOut = 33 # [°C] ... vystupni teplota chladici vody
    
    YPump = 300 # [J/kg] ... merna energie cerpadel
    etaPump = 0.5 # [–] ... ucinnost cerpadel
    PucpAtm = 0.2 # [kW] ... dodatecny prikon na ucpavku pri atmosfericke aplikaci
    PucpVak = 0.5 # [kW] ... dodatecny prikon na ucpavku pri vakuove aplikaci
    
    cVzdFug = 0.00005 # [–] ... odhadovana koncentrace vzduchu ve fugatu
    etaVyv = 0.5 # [–] ... ucinnost vyvevy
    n = 1.3 # [–] ... polytropicky koeficient
    mVod = 1.5 # [kg/h] ... prutok vody vyvevou
    
    MAir = 28.96 # [kg/kmol] ... molarni hmotnost vzduchu
    Vn = 22.414 # [Nm3/kmol] ... objem, ktery zaujima idealni plyn za NP
    pAtm = 101.325 # [kPa] ... normalni atmosfericky tlak
    
    
    # Definice vektoru promennych pro jednotlive komory:
    A = [0]*nCh # [m2] ... plocha kondenzacnich sekci
    OHTC = [0]*nCh # [kW/m2·K] ... soucinitel prostupu tepla pri kondenzaci
    LMTD = [0]*nCh # [K] ... stredni logaritmicky teplotni rozdily
    Qcond = [0]*nCh # [kW] ... vykon kondenzacni sekci
    TCW = [0]*nCh # [°C] ... teplota chladici vody/recirkulatu
    mCW = [0]*nCh # [kg/h] ... hmotnostni tok chladici vody/recirkulatu
    susCW = [0]*nCh # [%] ... obsah susiny v chladici vode/recirkulatu
    mdCum = [0]*nCh # [kg/h] ... kumulativni tok destilatu
    mdRevap = [0]*nCh # [kg/h] ... znovuodpareni odtekajiciho destilatu v dalsi komore
    mdTot = [0]*nCh # [kg/h] ... celkove mnozstvi par ke zkondenzovani
    mAir = [0]*nCh # [kg/h] ... tok vzduchu vlivem netesnosti komor (Ing. Matulka)
    TRec = [0]*nCh # [°C] ... teplota recirkulatu
    pCh = [0]*nCh # [kPa] ... absolutni tlak v komore
    TCh = [0]*nCh # [°C] ... saturacni teplota pri danem tlaku
    mdBot = [0]*nCh # [kg/h] ... hmotnostni tok par odparenych ve spodni casti
    mRec = [0]*nCh # [kg/h] ... hmotnostni tok recirkulatu
    cpRec = [0]*nCh # [kJ/kg·K] ... merna tepelna kapacita recirkulatu
    susRec = [0]*nCh # [%] ... obsah susiny v recirkulatu
    P = [0]*6 # [kW] ... prikon zarizeni    
    
#=============================== DEFINICE FUNKCI ==============================
    
    # Funkce na stanoveni merne tepelne kapacity suspenze:
    def cpfug(sus):
        cpFug = cpW*(1-sus/100)+cpDM*sus/100
        return cpFug
    
    # Funkce na stanoveni soucinitele prostupu tepla pri kondenzaci (voda-voda):
    Uwat = lambda x: (1617.5+0.1537*x+0.1825*x**2-0.00008026*x**3)/1000
    
    # Funkce pro stanoveni tlaku odpovidajici dane saturacni teplote:
    def sat_tlak(t):
        water=IAPWS97(T=t+273.15,x=0)
        p = water.P*1000
        return p
    
    # Funkce pro stanoveni vyparneho tepla pri danem tlaku:
    def vyp_tep(p):
        waterV=IAPWS97(P=p/1000,x=1)
        waterL=IAPWS97(P=p/1000,x=0)
        l = waterV.h-waterL.h
        return l
    
    # Funkce pro stanoveni stredniho logaritmickeho rozdilu:
    def fun_LMTD(T1i,T1o,T2i,T2o):
        dT1 = T1i-T2o
        dT2 = T1o-T2i
        dTln = (dT1-dT2)/ln(dT1/dT2) 
        return dTln
    
    def prirazeni(dT):
        global TRecIn
        # Napocitani teplot v komorach a recirkulatu dle daneho dT, se
        # zohlednenim teplotni deprese. Dale jsou stanoveny odpovidajici
        # tlaky v komorach.
        for i in range(nCh):
            TRec[i] = TRecOut +(nCh-i-1)*dT
            TCh[i] = TRec[i]-BPE
            pCh[i] = sat_tlak(float(TCh[i]))
        # Podminka pro zajisteni, ze v dohrevnem vymeniku nedojde k prekrizeni, ani
        # prilisnemu priblizeni teplotniho profilu
        TRecIn = TRecOut + nCh*dT
        if TRecIn > ThwIn-5:
            TRecIn = ThwIn-5
            dT = (TRecIn - TRecOut)/nCh
    
    # Toto je cast vypoctu, pres kterou se bude iterovat k nalezeni ucinnosti
    # vyparovani takove, pro kterou se bude tok destilatu z vypoctu odparovani
    # rovnat toku destilatu z materialove bilance:
    def prirazeni1(etaEv,susConc):
        # Nejprve materialova bilance. Z ucinnosti lze stanovit susinu v recirkulatu:
        TSRec = susConc*(100-etaEv)/100
        # Nasledne take tok recirkulatu pred vstupem do prvni komory:
        mRecIn = mFug*(susConc-susFug)/(susConc-TSRec)
        # A ze samotne definice ucinnosti take tok destilatu:
        mDest = mRecIn*etaEv/100
        for i in range(nCh):
            # Pro prvni komoru se vyuziva okrajovych podminek:
            if i == 0:
                cpRec[i] = cpfug(TSRec)
                mdBot[i] = cpRec[i]*mRecIn*(TRecIn-TRec[i])/vyp_tep(pCh[i])
                mRec[i] = mRecIn - mdBot[i]
                mdCum[i] = mdBot[i]
            # Pro dalsi komory uz vysledku z komory predchozi:
            else:
                # Merna tepelna kapacita se stanovuje dle susiny v predchozi komore:
                cpRec[i] = cpfug(susRec[i-1])
                # Cast entalpie proudu se uvedenim do podtlaku spotrebuje na
                # mzikove odpareni nasledujiciho mnozstvi:
                mdBot[i] = cpRec[i]*mRec[i-1]*(TRec[i-1]-TRec[i])/vyp_tep(pCh[i])
                # Prave o toto mnozstvi potece do dalsi komory mene recirkulatu:
                mRec[i] = mRec[i-1] - mdBot[i]
                # Destilat je odvaden kumulovane, spolecnym odtokem:
                mdCum[i] = mdBot[i] + mdCum[i-1]
            # Tok susiny je v kazde komore totozny, odparem se meni pouze podil na
            # celkovem mnozstvi recirkulatu - dochazi k zahusteni:
            susRec[i] = mRecIn*TSRec/mRec[i]
        # Tato funkce vraci nasledujici veliciny. Veliciny, ktere maji charakter
        # pole, jsou meneny na globalni urovni, zahrnout je do vystupu funkce
        # proto neni zapotrebi:
        return TSRec,mRecIn,mDest
    
    # Tato cast je "nadstavbou" vyse uvedeneho a vede k vypoctu tepelneho vykonu
    # dohrevu, pro jehoz danou hodnotu bude iterativne hledan obsah susiny v
    # koncentratu:
    def prirazeni2():
        # Pro vsechny komory je mozno napocitat nasledujici:
        for i in range(nCh):
            if i != 0:
                # Znovuodpareni odtekajiciho destilatu v dalsi komore (pro
                # prvni komoru je nulove):
                mdRevap[i] = mdCum[i-1]*cpW*(TCh[i-1]-TCh[i])/vyp_tep(pCh[i])
            # Celkove odparene mnozstvi, ktere je nutno zkondenzovat, je pak
            # sumou brydovych par vzniklych dole a znovuodpareneho destilatu:
            mdTot[i] = mdBot[i] + mdRevap[i]
            # Z tohoto celkoveho mnozstvi se pak pocita tepelny vykon
            # jednotlivych kondenzacnich sekci:
            Qcond[i] = mdTot[i]*vyp_tep(pCh[i])/3600
        
        # Hmotnostni prutok koncentratu lze stanovit v vnejsi MB jako rozdil
        # toku fugatu a destilatu (kumulativni odtok z posledni komory):
        mConc = mFug-mdCum[-1]
        # Mnozstvi recurkulatu po odebrani koncentratu z procesu:
        mRecOut = mRec[-1]-mConc
        # Po smiseni s fugatem se jedna o "chladici vodu" do sekci s procesni
        # vymenou tepla:
        mRecCW = mRecOut+mFug
        # Obsah susiny v tomto proudu lze stanovit na zaklade smesovaci rovnice:
        susRecCW = (susRec[-1]*mRecOut+susFug*mFug)/mRecCW
        # Tomuto obsahu susiny odpovida jista merna tepelna kapacita:
        cpRecCW = cpfug(susRecCW)
        # Teplotu lze pak stanovit z entalpicke bilance:
        TRecCW = (cpRec[-1]*mRecOut*TRec[-1]+cpfug(susFug)*mFug*TFug)/(cpRecCW*mRecCW)
        # Pro sekce s procesni vymenou tepla (HRS) recirkulat-brydove pary
        # jsou napocitany nasledujici veliciny:
        for i in range(nCh-nCW-1,-1,-1):
            # Soucinitel prostupu tepla je zde uvazovan nemenny:
            OHTC[i] = Ucond
            # Nemeni se ani obsah susiny a hmotnostni prutok:
            susCW[i] = susRecCW
            mCW[i] = mRecCW
            if i == nCh-nCW-1:
                # Pro nejteplejsi komoru s HRS jsou pouzity okrajove podminky,
                # nejprve je pocitana teplota, na kterou se recirkulat ohreje:
                TCW[i] = TRecCW+Qcond[i]*3600/(cpRecCW*mRecCW)
                # A na zaklade ni lze stanovit stredni logaritmicky teplotni
                # rozdil v dane kondenzacni sekci:
                LMTD[i] = fun_LMTD(TCh[i],TCh[i],TRecCW,TCW[i])
            else:
                # Pro zbylou cast se pak vyuzivaji vysledky z predchozi komory,
                # jinak je postup totozny:
                TCW[i] = TCW[i+1]+Qcond[i]*3600/(cpRecCW*mRecCW)
                LMTD[i] = fun_LMTD(TCh[i],TCh[i],TCW[i+1],TCW[i])
            # Na zaklade prenosove rovnice a tepelnych vykonu lze stanovit
            # velikost teplosmennych ploch:
            A[i] = Qcond[i]/(LMTD[i]*OHTC[i])
        # Po spocitani teploty za predehrevem recirkulatu je pak mozne stanovit
        # vykon ohrevneho vymeniku:
        Qohr = mRecCW*cpRecCW*(TRecIn-TCW[0])/3600
        return mConc,mRecOut,mRecCW,susRecCW,cpRecCW,TRecCW,Qohr
    
    # Tato cast slouzi k minimalizaci velikosti teplosmenne plochy v MSF, slouzi
    # tedy k dopocitani zbylych ploch a vraci prave tyto veliciny:
    def prirazeni3(Qohr):
        # Prutok chladici vody je spocitan na zaklade daneho teplotniho spadu,
        # uhrnneho vykonu sekci chlazenych vodou a merne tepelne kapacity:
        mCWv= sum(Qcond[nCh-nCW:nCh])/(cpW*(TCWOut-TCWIn))*3600
        # Pro tyto sekce (HJS) jsou pak proveden nasledujici vypocty:
        for i in range(nCh-1,nCh-nCW-1,-1):
            # Prutok chladici vody je ve vsech kond. sekcich totozny:
            mCW[i] = mCWv
            # Pro posledni (nejstudenejsi komoru) jsou uvazovany okrajove
            # podminky:
            if i == nCh-1:
                # Teplotu chladici vody lze stanovit z entalpicke bilance:
                TCW[i] = TCWIn+Qcond[i]*3600/(mCW[i]*cpW)
                # A na zaklade ni pak take stredni teplotni logaritmicky rozdil:
                LMTD[i] = fun_LMTD(TCh[i],TCh[i],TCWIn,TCW[i])
            # Pro ostatni komory jsou pocitany stejne veliciny:
            else:
                TCW[i] = TCW[i+1]+Qcond[i]*3600/(mCW[i]*cpW)
                LMTD[i] = fun_LMTD(TCh[i],TCh[i],TCW[i+1],TCW[i])
            # Pro vsechny komory je pak na zaklade aproximacniho vztahu
            # spocitany soucinitele prostupu tepla:
            OHTC[i] = Uwat(TCh[i])
        # Ziskanim vsech strednich teplotnich logaritmickych rozdilu se dospeje
        # do faze, kdy lze dopocitat teplosmenne plochy kondenzacnich sekci:
        for i in range(nCh):
            A[i] = Qcond[i]/(LMTD[i]*OHTC[i])
        LMTDohr = fun_LMTD(ThwIn,ThwOut,TCW[0],TRecIn)
        Aohr = Qohr/(LMTDohr*Uhw)
        Acelkem = Aohr + sum(A)
        return LMTDohr,Aohr,Acelkem

    # Trojstupnova funkce, ktera je volana pro minimalizaci teplosmenne plochy:
    def funkceMinA(dT,*args):
        # Globalni deklarace odhadu susiny v koncentratu umoznuje snadnou cestou
        # vkladat do dalsiho kroku iteracniho vypoctu vysledek predchazejiciho:
        global susKoncOdh
        # Napocitani teplotnich pomeru na zaklade dT, ktere je vstupem funkceMinA:
        prirazeni(dT)
            
        # Iterativni vypocet je dvoji - vnitrni (hleda ucinnost odparovani)
        # a vnejsi (hleda susinu v koncentratu).
        def funkceMid(susConc,*args):
            # I zde se ze stejnych duvodu jako vyse uplatnuje globalni deklarace:
            global etaEvOdh
            # Pro dany odhad susiny v koncentratu je nejprve nalezena ucinnost odparovani.
                
            # Samotna iteracni funkce - v prvni casti prirazeni promennych s uziti vyse
            # definovane funkce, v druhe casti kriterium, ktere se ma rovnat nule (rozdil
            # toku destilatu).
            def funkceIn(etaEv):
                TSRec,mRecIn,mDest = prirazeni1(etaEv,susConc)
                difIn = mdCum[-1]-mDest
                return difIn 
            # Volani iterativni metody secen:
            etaEvRes = newton(funkceIn,etaEvOdh,tol = 1e-4)
            # Vlozeni noveho odhadu do globalni promenne (jinak by po provedeni
            # instance funkce zanikla)
            etaEvOdh = etaEvRes
            # Pro dalsi reseni jsou prirazeny vsechny potrebne veliciny:
            TSRec,mRecIn,mDest = prirazeni1(etaEvRes,susConc)
            mConc,mRecOut,mRecCW,susRecCW,cpRecCW,TRecCW,Qohr = prirazeni2()
            # Na konci je vycislen rozdil skutecneho tepelneho vykonu dohrevu a jeho
            # pozadovane hodnoty, coz je rozhodne kriterium pro vnejsi iteracni vypocet:
            difOut = Qohr - Q
            # Volitelny argument umoznuje, aby vystupem funkce nebylo kriterium, ale
            # ucinnost odpareni:
            if len(args)>0:
                return etaEvRes
            # Pro potreby iterace je pak vracen vyse uvedeny rozdil:
            else:
                return difOut        
        # Volani dvoustupnove iterace:        
        susKoncRes = newton(funkceMid,susKoncOdh,tol = 1e-3)
        # ...a ulozeni noveho odhadu (viz vyse):
        susKoncOdh = susKoncRes
        # Za pouziti libovolneho druheho argumentu lze ziskat ucinnost odpareni:
        etaEvRes = funkceMid(susKoncRes,1)
        # Prirazeni velicin odpovidajicich nalezenemu reseni, vcetne stanoveni
        # velikosti teplosmennych ploch, ktere jsou ucelovou funkci k minimalizaci:
        TSRec,mRecIn,mDest = prirazeni1(etaEvRes,susKoncRes)
        mConc,mRecOut,mRecCW,susRecCW,cpRecCW,TRecCW,Qohr = prirazeni2()    
        LMTDohr,Aohr,Acelkem = prirazeni3(Qohr)
        # Stejne jako vyse volitelny argument umoznuje z funkce "vytahnout"
        # ucinnost odpareni a vysledny obsah susiny v koncenratu:
        if len(args)>0:
            return etaEvRes,susKoncRes
        # Pro potreby minimalizace je vsak vracena teplosmenna plocha:
        else:
            return Acelkem
#=================================== VYPOCET ==================================
    # Na pocatku nastavena takova susina v koncentratu, aby v prvnim kroku
    # nebyla splnena podminka cyklu while nize:
    susConc = 30 # [%] ... obsah susiny v koncentratu
    nCW = int(1/4*nCh) # [%] ... startovaci pocet chladicich komor je uvazovan
                       #         jako ctvrtina komor (3)
    up_bound = (ThwIn-TRecOut)/nCh-5/nCh # [K] ... horni mez dT s ohledem na
                                         #         kvalitu tepla z KGJ
    
    # Tento cyklus zajistuje, ze vysledna susina v koncentratu neprekroci
    # stanovene maximum:
    while susConc > susMax:
        # Cyklus, ktery zabranuje tomu, aby vypocet selhal (z duvodu prekrizeni
        # teplotniho profilu v kondenzacnich sekcich nebo extremnimu obsahu
        # susiny v koncentratu/recirkulatu)
        while nCW <= nCh:
            # Algoritmus zkusi s danymi pocatecnimi odhady provest zjednoduseny
            # vypocet (bez minimalizace):
            try:
                global susKoncOdh, etaEvOdh
                susKoncOdh = 1
                etaEvOdh = 5
                etaEv,susConc = funkceMinA((1+up_bound)/2,1) 
            # Pokud selze, je pridan jeden chladici stupen, coz, jak bylo
            # otestovano, problem spolu s cyklem while nCW <= nCh resi:
            except:
                nCW += 1
                if nCh == nCW:
                    return [0]*13
            # Jakmile je uspesne proveden zjednoduseny vypocet, cyklus konci:
            else:
                break
            
        # Je spocitan tok destilatu a koncentratu:    
        _,_,mDest = prirazeni1(etaEv,susConc)
        mConc,_,_,_,_,_,_ = prirazeni2() 
        
        # Nasledne je z MB stanoveno, jaky prutok koncentratu odpovida dane
        # maximani susine:
        m1 = mFug*susFug/susMax
        # Aktualni merne teplo na kg destilatu:
        sTh = Q/mDest
        # Umoznuje odhadnout potrebu tepla, pokud by se pri danem poctu
        # komor chlazenych vodou mel koncentrat zahustit az na maximalni
        # pripustnou miru:
        Qnove = sTh*(mFug-m1)
        # Pokud toto teplo je mensi nez 80 % dostupneho, pak je pridana sekce
        # chlazena vodou. Tim dojde ke skokovemu narustu potreby tepla:
        if Qnove < 0.8*Qdane:
            # Pokud je pocet komor chlazenych vodou jiz o jedno nizsi nez celkovy
            # pocet komor (na nCh = nCW vypocet neni uzpusoben), pak bude cyklus
            # opusten a kompletni vypocet nasledne proveden pro mensi nez 80%
            # vyuziti tepla:
            if nCh-1 == nCW:
                Q = Qnove
                susKoncOdh = 1
                etaEvOdh = 5
                break
            nCW += 1
        # Pokud je teplo potrebne na maximalni stanovene zahusteni naopak vetsi
        # nez disponibilni teplo, je cyklus opusten (a zahusteni bude nizsi
        # nez maximalni)
        elif Qnove > Qdane:
            break
        # Pokud teplo pro cilove zahusteni lezi teplo mezi 80 a 100 % 
        # disponibilniho tepla, je pouzito pro novy zjednoseny vypocet (viz vyse):
        else:
            Q = Qnove
            
    # Jakmile je pomoci zjednodusenych kalkulaci nalezena mira vyuziti tepla
    # a pocet stupnu chlazenych vodou, je provedena kompletni kalkulace vc.
    # minimalizace teplosmenne plochy:
    res = minimize_scalar(funkceMinA,bounds=(1,up_bound),method='bounded')
    etaEv,susConc = funkceMinA(res.get('x'),1)
    TSRec,mRecIn,mDest = prirazeni1(etaEv,susConc)
    mConc,mRecOut,mRecCW,susRecCW,cpRecCW,TRecCW,Qohr = prirazeni2() 
    LMTDohr,Aohr,Acelkem = prirazeni3(Qohr)

    # Vypocet vyvevy (zavisi na netesnostech a vzduchu vnesenem fugatem):
    for i in range (nCh):
        mAir[i] = (0.011-pCh[i]/10000)*mdTot[i]
    mAirCelk = sum(mAir)+cVzdFug*mFug
    VAir = mAirCelk/MAir*Vn*pAtm/pCh[-1]*(TCh[-1]+273.15)/273.15
    Pkomp = n/(n-1)*pCh[-1]*VAir/3600*((100/pCh[-1])**((n-1)/n)-1)/etaVyv
    Pvod = 90*mVod/1000
    
    # Vypocet prutoku topne vody ohrevnym vymenikem:
    mOhr = Q*3600/cpW/(ThwIn-ThwOut)        
    # Po kterem je k dispozici jiz vse pro spocitani vykonu vsech cerpadel:
    P[0] = Pkomp + Pvod + PucpVak
    P[1] = YPump*mFug/(3600*1000*etaPump)+PucpAtm
    P[2] = YPump*mRec[-1]/(3600*1000*etaPump)+PucpVak
    P[3] = YPump*mdCum[-1]/(3600*1000*etaPump)+PucpVak
    P[4] = YPump*mOhr/(3600*1000*etaPump)+PucpAtm
    P[5] = YPump*mCW[-1]/(3600*1000*etaPump)+PucpAtm
    # Celkova el. spotreba pak jeste zahrnuje nutnost chladit vodu z
    # kondenzacnich sekci:
    E = sum(Qcond[nCh-nCW-1:-1])/R+sum(P)
    # Vystup funkce:
    return mFug, susFug, mRecIn, TSRec, mConc, susConc, mDest, 0, -Q, Aohr, sum(A), nCW, -E