Off[NIntegrate::"slwcon",InverseFunction::"ifun"];

c=4186;
rce1=Qout==m1*c*(T3-T4);
rce2=Qout==Pel*ϵ;
rce3=m1*c*(T3-T5)==Pk;
rce4=m1*c*(T5-T4)==Qtop;
rce5=Qtop==kTop*(((T5+T4)/2)-Tin);
rce6=Qtop==kDum*(Tin-Tout);

rceStav={rce4,rce5,rce6};
Stav={Qtop->7000,Tout->-12,Tin->22,T5->60,T4->55};
dos=Solve[rceStav/.Stav][[1]]//N

dosrce={rce1,rce2,rce3,rce4,rce5,rce6}/.Union[dos,{Pk->0,ϵ->0.5*((T3+273)/(T3-Tout)),Tin->22}];
Cases[dosrce,_Symbol,∞]//Union

reseni[0]

reseni[ToutMoje_]:=Solve[dosrce/.Tout->ToutMoje][[1]];
epsilon[ToutMoje_]:=0.5((T3+273)/(T3-Tout))/.Union[reseni[ToutMoje],{Tout->ToutMoje}];
vykon[ToutMoje_]:=Qout/.reseni[ToutMoje];
Plot[{epsilon[Tout],((vykon[Tout])/1000)},{Tout,-18,15},AxesLabel->{"Tout[°C]","ϵ[1],Qout[kW]"},AxesOrigin->{-5,0},PlotStyle->Thick]

TemperatureData={{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[month_]=Interpolation[TemperatureData][month];
Temperature[m_]:=pomfce[m]/;1≤m<13;
Temperature[m_]:=Temperature[m-12]/;m≥13;
Temperature[m_]:=Temperature[m+12]/;m<1;
teplsecpom[t_]:=Temperature[(t/(30.5*24*3600))];
rok=12*30.5*24*3600;
Plot[{teplsecpom[t],13},{t,0.2*rok,2.5rok},AxesLabel->{"t[s]","T[°C],13°C"},AxesOrigin->{0,0},PlotStyle->Thick]

strhodn=(1/rok)NIntegrate[teplsecpom[t],{t,0,rok}];
zadstrhodn=6.8;
teplsec[t_]:=teplsecpom[t]-strhodn+zadstrhodn;

touttopim=13;
Plot[{teplsec[t],touttopim},{t,0,2*12*30.5*24*3600},AxesLabel->{"t[s]","T[°C]"},PlotLabel->"  teplota behem roku",PlotStyle->{{Thickness[0.008],GrayLevel[0]},{Thickness[0.008],RGBColor[1,0,0]}}]
castopenipoc=t/.FindRoot[teplsec[t]==touttopim,{t,2*(10^7),3.5*(10^7)}];
Print["zacatek topne sezony =  ",castopenipoc,"  s"]
castopenikonec=t/.FindRoot[teplsec[t]==touttopim,{t,4*(10^7),5*(10^7)}];
Print["konec topne sezony =  ",castopenikonec,"  s"]

epsilon[0]

(*Plot[vykon[teplsec[t]],{t,castopenipoc,castopenikonec}];*)
dataTopnyVykon=Table[{T,vykon[T]},{T,-20,20,0.5}];
vvyk[T_]=Interpolation[dataTopnyVykon][T];
dataElektrickyVykon=Table[{T,((vykon[T])/(epsilon[T]))},{T,-20,20,0.5}];
elvvyk[T_]=Interpolation[dataElektrickyVykon][T];

energienatopeni=NIntegrate[vvyk[teplsec[t]],{t,castopenipoc,castopenikonec}];
Print["potrebne teplo= ",(10^(-9))*%,"  GJ"]
energieProCerpadlo=NIntegrate[elvvyk[teplsec[t]],{t,castopenipoc,castopenikonec}];
Print["koupena elektrika= ",(10^(-9))*%,"  GJ"]
Print["prumerny topny faktor= ",(energienatopeni/energieProCerpadlo),"  -"]