(*Tabulka parametrů páry s veličinami:
teplota [°C], tlak [MPa], entropie [J/(kg*K)], entalpie [kJ/kg]*)
tab={{540.75,16.817,6420,3405},{535.7,15.551,6452,3405},{347.8,4.372,6523,3078},{540.69,4.115,7195,3538},{538.34,3.912,7214,3534},{282.8,0.593,7316,3026},{39.4,0.00715,7674,2387},{38.9,0.00715,558,162.9},{39,1.2,555,163.3},{70.2,1,955,293.9},{89.5,0.9,1184,374.7},{110.2,0.8,1419,462.2},{129.7,0.7,1630,545.1},{151.1,0.5,1853,636.8},{155.7,20,1881,669.9},{211.8,19.5,2413,912.7},{250.5,19,2764,1089}};

odbery={{450.4,1.969,7295,3359},{220.9,0.3197,7371,2908},{162.3,0.1746,7405,2796},{103.1,0.0853,7455,2684},{78.5,0.0446,7497,2593}};

h=Flatten[Take[tab,All,-1]];
(*vyjádření pouze entalpií*)

ho=Flatten[Take[odbery,All,-1]];
(*vyjádření pouze entalpií odběrů*)

hk={929.5(*216.8°C*),678.3(*160.7°C*),483.4(*115.2°C*),395.9(*94.5*),205.2(*49°C*)};
hOHREV=458(*109°C*);
(*entalpie kondenzátu při daných teplotách*)

(*rovnice pro jednotlivá zařízení:*)

rceVTO2=oVTO2*(h[[3]]-hk[[1]])==h[[17]]-h[[16]];
rceVTO1=oVTO1*(ho[[1]]-hk[[2]])+oVTO2*(hk[[1]]-hk[[2]])==h[[16]]-h[[15]];
rceTO=oTO*(h[[6]]-h[[14]])+(oVTO2+oVTO1)*(hk[[2]]-h[[14]])==(1-oVTO2-oVTO1-oTO)*(h[[14]]-h[[13]]);
rceNTO4=oNTO4*(ho[[2]]-hk[[3]])==(1-oVTO2-oVTO1-oTO)*(h[[13]]-h[[12]]);
rceNTO3=oNTO3*(ho[[3]]-hk[[4]])+oNTO4*(hk[[3]]-hk[[4]])==(1-oVTO2-oVTO1-oTO)*(h[[12]]-h[[11]]);
rceNTO2=oNTO2*(ho[[4]]-h[[11]])+(oNTO4+oNTO3)*(hk[[4]]-h[[11]])==(1-oVTO2-oVTO1-oTO-oNTO4-oNTO3-oNTO2)*(h[[11]]-h[[10]]);
rceNTO1=oNTO1*(ho[[5]]-hk[[5]])+oOHREV*(hOHREV-hk[[5]])==(1-oVTO2-oVTO1-oTO-oNTO4-oNTO3-oNTO2)*(h[[10]]-h[[9]]);
rceTBN=0.75*oTBN*(h[[6]]-h[[7]])==(tab[[15]][[2]]-tab[[14]][[2]])*1.11;
rceOHREV=oOHREV==29200/1567000;
(*bez TBN: rceTBN=oTBN==0;*)
(*bez VTO1: rceVTO2=oVTO2*(h[[3]]-hk[[2]])==h[[17]]-h[[15]];
rceVTO1=oVTO1==0;*)

rce={rceVTO2,rceVTO1,rceTO,rceNTO4,rceNTO3,rceNTO2,rceNTO1,rceTBN,rceOHREV};

promenne={oVTO2,oVTO1,oTO,oNTO4,oNTO3,oNTO2,oNTO1,oTBN,oOHREV};

res=NSolve[rce,promenne][[1]](*řešení soustavy rovnic*)

doNT=oVTO2+oVTO1+oTO;

aVT=h[[2]]-h[[3]];
aST=(1-oVTO2)*(h[[5]]-ho[[1]])+(1-oVTO2-oVTO1)*(ho[[1]]-h[[6]]);
aNT=(1-doNT-oTBN-oOHREV)*(h[[6]]-ho[[2]])+(1-doNT-oTBN-oOHREV-oNTO4)*(ho[[2]]-ho[[3]])+(1-doNT-oTBN-oOHREV-oNTO4-oNTO3)*(ho[[3]]-ho[[4]])+(1-doNT-oTBN-oOHREV-oNTO4-oNTO3-oNTO2)*(ho[[4]]-ho[[5]])+(1-doNT-oTBN-oOHREV-oNTO4-oNTO3-oNTO2-oNTO1)*(ho[[5]]-h[[7]]);
a=aVT+aST+aNT/.res;(*určení vykonané práce*)
(*bez TBN: aELN=(tab[[15]][[2]]-tab[[14]][[2]])*1.11/0.75/0.9;
a=aVT+aST+aNT-aELN/.res;*)

q=h[[1]]-h[[17]]+(1-oVTO2)*(h[[4]]-h[[3]])/.res;(*určení dodaného tepla*)

η=100*(a/q)(*výpočet účinnosti*)

ηmax=100*((tab[[1]][[1]]-tab[[8]][[1]])/(tab[[1]][[1]]+273.15))(*výpočet maximální možné účinnosti*)

(η/(ηmax))*100(*výpočet poměru dosažené a maximální účinnosti*)

P1=1195000.;

P1/q(*výpočet celkového hmotnostního toku*)

(promenne*P1/q)/.res(*výpočet jmenovitých odběrů*)

P1*η/100000(*výpočet vyráběného výkonu*)