ClearAll["Global`*"];

data={{-5,35,3.7},{-5,50,2.9},{-5,65,2.2},{0,35,4.3},{0,50,3.3},{0,65,2.5},{5,35,5},{5,50,3.7},{5,65,2.8},{10,35,5.8},{10,50,4.2},{10,65,3.1},{15,35,6.8},{15,50,4.8},{15,65,3.5}};
model2[t1_,t4_]=Fit[data,{((t4+273)/(t4-t1))},{t1,t4}]

dos={eps->model2[t1,t4],m1->(1300./3600),m2->(600./3600),cp->4186,Tin->20};
dos1={Qtop->4000,t4->40,Tout->-5};

rovnice={m1*cp*t1+m2*cp*t3+Pel==m1*cp*t2+m2*cp*t4,eps==(Qtop/Pel),Qtop==m2*cp*(t4-t3),Qtop==k1*(Tin-Tout),Qtop==k2*(((t3+t4)/2)-Tin),Tout==t1};

rce=rovnice//.dos/.dos1;
nezname=Union[Cases[rce,_Symbol,{0,∞}]];
res1=Solve[rce,nezname][[1]]//Quiet

datatepl={{1,-1.4},{2,-.2},{3,3.7},{4,8.9},{5,13.8},{6,17.2},{7,19},{8,18.5},{9,14.8},{10,9.1},{11,4.1},{12,.3},{13,-1.4}};
pomfce[mes_]=Interpolation[datatepl][mes];
(*Plot[pomfce[mes],{mes,1,13}]*)
Tepl[m_]:=pomfce[m]/;1≤m<13;
Tepl[m_]:=Tepl[m-12]/;m≥13;
Tepl[m_]:=Tepl[m+12]/;m<1;
Plot[Tepl[mes],{mes,0,24}];
teplsecpom[t_]:=Tepl[(t/(30.5*24*3600))];
rok=365*24*3600;
str=(1/rok)NIntegrate[teplsecpom[t],{t,0,rok}];
strZadana=5.3;
teplsec[t_]:=teplsecpom[t]-str+strZadana;
Plot[teplsec[t],{t,0,2rok}];

touttopim=13;
Plot[{teplsec[t],touttopim},{t,0,2*12*30.5*24*3600},AxesLabel->{"t[s]","T[°C]"},PlotLabel->"Roční průběh teploty"]
castopenipoc=t/.FindRoot[teplsec[t]==touttopim,{t,2*(10^7),3.5*(10^7)}]
castopenikonec=t/.FindRoot[teplsec[t]==touttopim,{t,4*(10^7),5*(10^7)}]

ClearAll[rceCas,resCas];
dos2={k1->(k1/.res1),k2->(k2/.res1)}
rceCas[Tven_]:=rovnice/.Union[dos2,dos]/.Tout->Tven
resCas[Tven_]:=Solve[rceCas[Tven]][[1]]//Quiet;
resCas[5]

datPel=Quiet[{#,Pel/.resCas[#]}&/@Range[-5,15]];
datQtop=Quiet[{#,Qtop/.resCas[#]}&/@Range[-5,15]];

prikon=Interpolation[datPel];
vykon=Interpolation[datQtop];
Plot[{vykon[tepl],prikon[tepl]},{tepl,-5,13},AxesLabel->{"T[°C]","P[kW]"},PlotLabel->"Příkon\!\(\*
StyleBox[\(\*
StyleBox[" ",
FontColor->RGBColor[0, 0, 1]]\\ \\ \\ \)]\)Výkon",PlotStyle->{{Thick,Red},{Thick,Blue}},GridLines->Automatic]

prikonvcase[t_]:=prikon[teplsec[t]];
vykonvcase[t_]:=vykon[teplsec[t]];
Plot[vykonvcase[t],{t,castopenipoc,castopenikonec}];
energieTopna=NIntegrate[vykonvcase[t],{t,castopenipoc,castopenikonec}];
energieElektricka=NIntegrate[prikonvcase[t],{t,castopenipoc,castopenikonec}];
Print["Spotreba elektriny je: ",(energieElektricka/(3600*1000))," kWh, prumerny top. faktor je: ",(energieTopna/energieElektricka)]