Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

λConcrete=1;
λizol=0.04;
dConcr=0.2;
dIzol=0.15;
αin=5;
αout=15;
Δ=dConcr+dIzol;
a=10;
Stotal=5*(((a+2Δ))^2);
nWindows=24;
S1window=1.5*2;
Swindows=nWindows*S1window;
Swalls=Stotal-Swindows;
Tin=20;
Uwindows=1;
RthWalls=(1/Swalls)*((1/(αin))+(dConcr/(λConcrete))+(dIzol/(λizol))+(1/(αout)));
powerWalls=((Tin-Tout)/RthWalls);
powerWindows=Uwindows*Swindows*(Tin-Tout);
cAir=1000;
massFlow=densityAir*VolumeFlow;
timeOfAirExchange=3*3600;
Volume=(a^3);
VolumeFlow=(Volume/timeOfAirExchange);

densityAir=(349./(273+Tin));
powerVentilation=massFlow*cAir*(Tin-Tout);
powerTotal[Tout_]=powerVentilation+powerWalls+powerWindows;
powerTotal[T]

(*known state of heating system*)
powerKnown=powerTotal[-18];
T2Known=60;
T1Known=T2Known-5;
cH2O=4186;
eqnsOfKnownState={powerKnown==mH2O1*cH2O*(T2Known-T1Known),
powerKnown==kS1*((T2Known-Tin+T1Known-Tin)/2)
};
{kS,mH2O}={kS1,mH2O1}/.Solve[eqnsOfKnownState][[1]]

COP[T2_,T1_]:=0.5*((T2+273)/(T2-T1));
ClearAll[eqns]
eqns[Tout_]:={powerTotal[Tout]==mH2O*cH2O*(T2-T1),
powerTotal[Tout]==kS*((T2-Tin+T1-Tin)/2),
powerTotal[Tout]==Pel*COP[T2,T1]

};
Solve[eqns[-18]][[1]]
Pelectrical[Tout_]:=Pel/.Solve[eqns[Tout]][[1]];
Pelectrical[-13]

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=Interpolation[TemperatureData];
Temperature[m_]:=pomfce[m]/;1≤m<13;
Temperature[m_]:=Temperature[m-12]/;m≥13;
Temperature[m_]:=Temperature[m+12]/;m<1;
secondsInMonth=365.*24*3600/12;
teplsecpom[t_]:=Temperature[(t/secondsInMonth)];
year=12*secondsInMonth;
meanValueOfTemperature=(1/year)NIntegrate[teplsecpom[t],{t,0,year}];
desiredMeanValue=4.5;
teplsec[t_]:=teplsecpom[t]-meanValueOfTemperature+desiredMeanValue;
touttopim=13;
Plot[{teplsec[t],touttopim},{t,0,2*12*30.5*24*3600},AxesLabel->{"t[s]","T[°C]"},PlotLabel->"  Průměrná denní teplota jako funkce času",PlotStyle->{{Thickness[0.005],GrayLevel[0]},{Thickness[0.007],RGBColor[1,0,0]}},GridLines->Automatic];
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)}];
Ampl=7;
den=3600*24;
teplDenni[t_]:=Ampl*Sin[2*Pi*(t/den)]+teplsec[t];
Plot[teplDenni[t],{t,castopenipoc,castopenikonec}];

dataForInterpolation={#,Pelectrical[#]}&/@Range[-25,14,0.5];
intPel=Interpolation[dataForInterpolation];
energyElectrical=Quiet@NIntegrate[intPel[teplDenni[t]],{t,castopenipoc,castopenikonec}];
heat=Quiet@NIntegrate[powerTotal[teplDenni[t]],{t,castopenipoc,castopenikonec}];
averageCOP=(heat/energyElectrical)