Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
c=4186.;
den=3600*24.;
rok=365.25*den;
kPower=1;
{PelHP[Tvzduch_,Tvoda_],PtopHP[Tvzduch_,Tvoda_],COP[Tvzduch_,Tvoda_]}={kPower,kPower,1}*<<"struct";

TprumDen=-10;
Ampl=7;
teplota[t_]:=TprumDen+Ampl*Sin[2*Pi*(t/den)];

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)*Quiet@NIntegrate[teplsecpom[t],{t,0,year}];
desiredMeanValue=7.5;
teplsec[t_]:=teplsecpom[t]-meanValueOfTemperature+desiredMeanValue;
teplDenni[t_]:=Ampl*Sin[2*Pi*(t/den)]+teplsec[t];
Plot[teplsec[t],{t,0,12*secondsInMonth}];

Quiet@Plot3D[#[Tvzduch,Tvoda],{Tvzduch,-15.`,12.`},{Tvoda,25.`,55.`}]&/@{PelHP,PtopHP,COP}

(*funkce "zpracuj" je vytvori z puvodnich interpolacnich funkci PelHP,PtopHP a COP nove interpolacni funkce tak, aby byla interpolace prvniho radu. Pri vypoctech se totiz muzeme dostat mimo rozsah puvodnich interpolacnich interpolacnich funkci a proklady vyssiho radu, kdyz se dostaneme mimo rozsah, "utikaji" daleko*)
ClearAll[pelhp,ptophp,zpracuj,cop,Tvzduch,Tvoda];
zpracuj[fce_]:=Module[{Tamin,Tamax,Tvmin,Tvmax,data,ΔT=0.2},
{{Tamin,Tamax},{Tvmin,Tvmax}}=fce[T1,T2][[2,0,1]];
data=Flatten[Table[{Ta,Tv,fce[Ta,Tv]},{Ta,Tamin,Tamax,ΔT},{Tv,Tvmin,Tvmax,ΔT}],1];
Interpolation[data,InterpolationOrder->1]
];
{pelhp,ptophp}=zpracuj/@{PelHP,PtopHP};

cop[Tvzduch_,Tvoda_]:=((ptophp[Tvzduch,Tvoda])/(pelhp[Tvzduch,Tvoda]));
Quiet@Plot3D[#[Tvzduch,Tvoda],{Tvzduch,-20.`,20.`},{Tvoda,15.`,70.`}]&/@{pelhp,ptophp,cop};

ClearAll[vydej];
vydej[{{Pel_,m1_,MA_,TStart_,Tzadana_},{tStart_,Tvzduch_}}]:=Module[{rce,tEnd,tmax,sol,c=4186.,ploty,Tmean,tEndpom,soltEnd,T1mean,EmeanEl,ErealEl,COPTprum,COPReal,Etopna,TmeanDenni,EmeanElDenni,COPPrumDenni,Tnadrz},

rce={Ptopny[t]==Pel*cop[Tvzduch[t],T1[t]],Ptopny[t]==m1*c(T1[t]-TA[t]),
MA*c*TA'[t]==Ptopny[t]
,TA[tStart]==TStart,e'[t]==Ptopny[t],e[tStart]==0,eEl'[t]==Pel,eEl[tStart]==0,TforMean'[t]==Tvzduch[t],TforMean[tStart]==0};
tmax=tStart+1*den;
sol=Quiet@NDSolve[rce,{TA,Ptopny,e,eEl,T1,TforMean},{t,tStart,tmax}][[1]];
Tnadrz:=TA/.sol;
soltEnd=FindRoot[Tnadrz[tEndpom]==Tzadana,{tEndpom,tStart+0.05*den}];
tEnd=tEndpom/.soltEnd;
Tmean=((TforMean[tEnd]-TforMean[tStart])/(tEnd-tStart))/.sol;
TmeanDenni=((TforMean[tmax]-TforMean[tStart])/(tEnd-tStart))/.sol;
T1mean=0.5*(TStart+Tzadana);
ErealEl=(eEl[tEnd]-eEl[tStart])/.sol;
Etopna=c*MA*(Tzadana-TStart);
EmeanEl=(Etopna/(COP[Tmean,T1mean]));
EmeanElDenni=(Etopna/(COP[TmeanDenni,T1mean]));
{COPTprum,COPReal,COPPrumDenni}=(({((EmeanEl+1)/Etopna),((ErealEl+1)/Etopna),((EmeanElDenni+1)/Etopna)})^(-1));
ploty=Quiet@{Plot[{COP[Tvzduch[t],T1[t]]/.sol,COPTprum,COPPrumDenni},{t,tStart,tmax},AxesLabel->{"t [s]","COP [-]"},PlotStyle->{{Red,Thick},Black,Blue},PlotRange->{{tStart,tEnd},Automatic},GridLines->Automatic],
Plot[{TA[t]/.sol,Tzadana},{t,tStart,tmax},AxesLabel->{"t [s]","TA [°C]"},PlotStyle->{{Red,Thick},Black},PlotRange->{{tStart,tEnd},{0,Tzadana+10}},GridLines->Automatic]
};
{ploty,{Tmean,{EmeanEl,ErealEl,Etopna,{COPTprum,COPReal},100.*((ErealEl-EmeanEl)/ErealEl),"%  ",100.*((ErealEl-EmeanElDenni)/ErealEl),"%  "}}}

];
pooom=vydej[{{2000,0.5,1500,20,50},{180000,teplDenni}}]

gra=pooom[[1,1]];
Export["gra.png",gra,ImageResolution->300]

grb=pooom[[1,2]];
Export["grb.png",grb,ImageResolution->300]