Quiet@Remove["Global`*"];
$HistoryLength=2;
Off[FindRoot::lstol]
SetDirectory[NotebookDirectory[]];
c=4186.;
hodina=3600.;
den=hodina*24.;
rok=365.25*den;
dosStav={Ptop->6000.,T1->45,ΔT->5.,Tin->20};
parsDosStav=#[[2]]&/@dosStav;
Ptopny[T_]=<<"Ptopny";
Tmez=Tt/.Solve[Ptopny[Tt]==0][[1]];

Adenni=7.5;
Arocni=8.1;
meanRocni=8.2;
Trok[t_]:=meanRocni+Arocni*Sin[2.*Pi*(t/rok)+1.4*Pi];
Ttrial[t_]:=Trok[t]+Adenni*Sin[2.*Pi*(t/den)];
Ttopime=13.;
{t1,t2}=tTopim/.FindRoot[Trok[tTopim]==Ttopime,{tTopim,#*rok}]&/@{0.7,1.4};
Plot[Trok[t],{t,t1,t2}];
Plot[Ttrial[t],{t,t1,t2}];
{{t1,t2},Ttrial[t],Trok[t]}>>"oTeplote";

kPower=0.3;
{PelHP[Tvzduch_,Tvoda_],PtopHP[Tvzduch_,Tvoda_],COP[Tvzduch_,Tvoda_]}={kPower,kPower,1}*<<"struct";

{{Tamin,Tamax},{T1min,T1max}}=PelHP[a,b][[2,0,1]];

ClearAll[identifikaceStav];
identifikaceStav[{Ptop_,T1_,ΔT_,Tin_}]:=Module[{rceZnamyStav,T2=T1-ΔT},
rceZnamyStav={m*c*(T1-T2)==Ptop,Ptop==k*(((T1+T2)/2)-Tin)};
Solve[rceZnamyStav][[1]]

];
identifikaceStav[parsDosStav];

ClearAll[vydej2];
vydej2[{znamyStav:{Ptop_,T1znamy_,ΔT_,Tin_},{tStart_,tEnd_},{Texterni_,Ptopny_}}]:=Module[{krokCasu=0.5*hodina,krokTeploty=1.5,dosZnamyStav=identifikaceStav[znamyStav],rce,Tvzduch,dos,rn,dosss,datPel,Tmin=-30.,Tmax,intPel,fpom,plP,plt1t2,sole,ple,Tmean,EmeanEl,ErealEl,Etopna,Pftopny,COPTpum,COPReal,rceHranice,poc,resHranice,TmezBivalence,rceCerpadlo100,res100,rceBivalenti,resBival,res,Pelektricky,Post,Pelint,datPost,Postint,Pcelkint,PfelTC,Pfost,Pfcelk,Eost,Ecelk,EostMean,dilciVys1,teploty},
Tmax=T/.Solve[Ptopny[T]==0,T][[1]];
rceHranice[Tvzduch_]={m*c*(T1-T2)==Ptopny[Tvzduch],Ptopny[Tvzduch]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==Ptopny[Tvzduch],Ptopny[Tvzduch]==PtopHP[Tvzduch,T1]}/.dosZnamyStav;
poc={{Pel,6000},{T1,45},{T2,40},{Tvzduch,-5}};
resHranice=FindRoot[rceHranice[Tvzduch],poc];
TmezBivalence=Tvzduch/.resHranice;
(*Print["hranice bivalence je ",TmezBivalence," °C"]*)
rceCerpadlo100[Tvzduch_]={m*c*(T1-T2)==Ptopny[Tvzduch],Ptopny[Tvzduch]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==Ptopny[Tvzduch]}/.dosZnamyStav;

res100[Tvzduch_]:=Union[FindRoot[rceCerpadlo100[Tvzduch],{{Pel,6000},{T1,45},{T2,40}}],{PtopOstatni->0}];
(*res100[TmezBivalence]*)
rceBivalenti[Tvzduch_]={m*c*(T1-T2)==PtopHP[Tvzduch,T1],PtopHP[Tvzduch,T1]==k*(((T1+T2)/2)-Tin),COP[Tvzduch,T1]*Pel==PtopHP[Tvzduch,T1],PtopOstatni+PtopHP[Tvzduch,T1]==Ptopny[Tvzduch]}/.dosZnamyStav;
resBival[Tvzduch_]:=FindRoot[rceBivalenti[Tvzduch],{{Pel,6000},{T1,45},{T2,40},{PtopOstatni,1000}}];
(*resBival[TmezBivalence]*)
res[Tvzduch_]:=Quiet@res100[Tvzduch]/;Tvzduch≥TmezBivalence;
res[Tvzduch_]:=Quiet@resBival[Tvzduch]/;Tvzduch<TmezBivalence;
Pelektricky[Tvzduch_]:=Pel/.res[Tvzduch];
Post[Tvzduch_]:=PtopOstatni/.res[Tvzduch];
datPel={#,Pelektricky[#]}&/@Range[Tmin,Tmax,0.1];
Pelint=Interpolation[datPel,InterpolationOrder->1];
(*Plot[Pelint[T],{T,-10,14}]*)
datPost={#,Post[#]}&/@Range[Tmin,Tmax,0.1];
Postint=Interpolation[datPost,InterpolationOrder->1];
(*Plot[Postint[T],{T,-10,14}]*)
Pcelkint:=Pelint[#]+Postint[#]&;
(*Plot[Pcelkint[T],{T,-10,14}];*)
Pftopny[t_]:=If[Ptopny[Texterni[t]]≥0,Ptopny[Texterni[t]],0.];
PfelTC[t_]:=If[Ptopny[Texterni[t]]≥0,Pelint[Texterni[t]],0.];
Pfost[t_]:=If[Ptopny[Texterni[t]]≥0,Postint[Texterni[t]],0.];
Pfcelk[t_]:=If[Ptopny[Texterni[t]]≥0,Pcelkint[Texterni[t]],0.];

sole=NDSolve[{e[tStart]==0,e'[t]==PfelTC[t],TforMean[tStart]==0,TforMean'[t]==Texterni[t],eTopna[tStart]==0,eTopna'[t]==Pftopny[t],eOst[tStart]==0,eOst'[t]==Pfost[t],eCelk[tStart]==0,eCelk'[t]==Pfcelk[t]},{e,TforMean,eTopna,eOst,eCelk},{t,tStart,tEnd},MaxSteps->(10^6),MaxStepSize->0.1*hodina][[1]];

ple=Plot[((e[t]/.sole)/(10^6)),{t,tStart,tEnd},AxesLabel-> {"t [s]","E [MJ]"},PlotStyle->{Black,Thick},GridLines->Automatic,AxesOrigin->{tStart,0}];
(*Plot[#[t],{t,tStart,tEnd}]&/@{Pftopny,PfelTC,Pfost,Pfcelk}*)
Tmean=((TforMean[tEnd]-TforMean[tStart])/(tEnd-tStart))/.sole;
EmeanEl=If[Ptopny[Tmean]≥0,Pelint[Tmean]*(tEnd-tStart),0.];
ErealEl=(e[tEnd]-e[tStart])/.sole;
Etopna=(eTopna[tEnd]-eTopna[tStart])/.sole;
Eost=(eOst[tEnd]-eOst[tStart])/.sole;
Ecelk=(eCelk[tEnd]-eCelk[tStart])/.sole;
EostMean=If[Ptopny[Tmean]≥0,Postint[Tmean]*(tEnd-tStart),0.];
dilciVys1={Ecelk,Eost,ErealEl,EmeanEl,EostMean,Etopna};
{{tStart,tEnd},dilciVys1,{Pelint,Postint}}

];
{{tStart,tEnd},dilciVys1,{Pelint,Postint}}=vydej2[{parsDosStav,{0.5*(t1+t2),0.5*(t1+t2)+1*den},{Ttrial,Ptopny}}];

ClearAll[teplarenskyVypocet];
teplarenskyVypocet[{krokT_,krokt_},{Texterni_,Ptopny_},{tStart_,tEnd_}]:=Module[{Tmax,znuluj,Pel,Post,teploty,obory,oborySCetnostmi},
Tmax=T/.Solve[Ptopny[T]==0,T][[1]];
Pel[T_]:=If[T<Tmax,Pelint[T],0.];
Post[T_]:=If[T<Tmax,Postint[T],0.];
(*Plot[{Pel[T],Post[T]},{T,-10,20}]*)
teploty=Texterni/@Range[tStart,tEnd,krokt];
obory=Partition[Range[-30,20,krokT],2,1];
oborySCetnostmi=obory/.{T1_,T2_}:>{{T1,T2},Length[Select[teploty,T1≤#<T2&]]};
Plus@@Quiet[oborySCetnostmi/.{{T1_,T2_},n_}:>n*krokt*((Length[teploty])/(Length[teploty]+1))*{Pelint[0.5*(T1+T2)],Postint[0.5*(T1+T2)]}]
];
teplarenskyVypocet[{0.5,hodina/60.},{Ttrial,Ptopny},{0.5*(t1+t2),0.5*(t1+t2)+1*den}];
pom=dilciVys1[[2;;3]]//Reverse;
((%-%%)/(%))*100;
ClearAll[vydej4];
vydej4[{krokT_,krokt_}]:=Module[{zasrTepl,errs},
zasrTepl=teplarenskyVypocet[{krokT,krokt},{Ttrial,Ptopny},{0.5*(t1+t2),0.5*(t1+t2)+1*den}];
errs=((pom-zasrTepl)/zasrTepl)*100;
errs

];
vydej4[{1.2,0.6*hodina}];
ClearAll[lajnaCasu];
maxdT=5;
mindT=0.5;
lajnaCasu[krokT_]:=Module[{errs,dats,lplPel,lplPost},
errs={#,vydej4[{krokT,#}]}&/@Range[0.1*hodina,2*hodina,0.1*hodina];
lplPel=ListPlot[errs/.{dt_,{Pel_,Post_}}:>{dt,Pel},PlotRange->All,AxesLabel->{"Δt","δPel  [%]"},PlotStyle->{Thickness[0.004],Hue[(krokT/maxdT)]},Joined->True];
lplPost=ListPlot[errs/.{dt_,{Pel_,Post_}}:>{dt,Post},PlotRange->All,AxesLabel->{"Δt"," δPost [%]"},PlotStyle->{Thickness[0.004],Hue[(krokT/maxdT)]},Joined->True];
{lplPel}

];
datsLajna=lajnaCasu/@Range[0.5,5,0.2];
Show[Flatten[datsLajna],PlotRange->All]

ClearAll[lajnaTeploty];
maxdt=2*hodina;
mindt=0.1*hodina;
lajnaTeploty[krokCasu_]:=Module[{errs,dats,lplPel,lplPost},
errs={#,vydej4[{#,krokCasu}]}&/@Range[0.1,maxdT,0.1];
lplPel=ListPlot[errs/.{dt_,{Pel_,Post_}}:>{dt,Pel},PlotRange->All,AxesLabel->{"ΔT","δEel  [%]"},PlotStyle->{Thickness[0.004],Hue[(krokCasu/maxdt)]},Joined->True];
lplPost=ListPlot[errs/.{dt_,{Pel_,Post_}}:>{dt,Post},PlotRange->All,AxesLabel->{"ΔT"," δPost [%]"},PlotStyle->{Thickness[0.004],Hue[(krokCasu/maxdt)]},Joined->True];
{lplPel}

];
datsLajnaT=lajnaTeploty/@Range[mindt,maxdt,mindt];
Show[Flatten[datsLajnaT],PlotRange->All]

Plot[{Pelint[T],Postint[T]},{T,-10,13.5},AxesLabel->{"T [C]","Pel  Post [W] "},PlotStyle->{{Thick,Blue},{Thick,Red}},GridLines->Automatic];