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]];

{data3,denniData}=N[<<"teploty"];
intTBrno=Interpolation[data3];
intPrumBrno=Interpolation[denniData];
{ta,tb}=intPrumBrno[[1,1]];


ClearAll[nahr1]
nahr1[t_]=Fit[denniData,{1,Sin[2*Pi*(t/rok)],Cos[2*Pi*(t/rok)]},t]
Plot[{intPrumBrno[t],nahr1[t]},{t,ta,tb},PlotStyle->{{Blue,Thick,Opacity[0.5]},{Red,Thick,Opacity[0.9]}},AxesLabel->{"t [s]"," \!\(\*
StyleBox[\"Tdenni\",\nFontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\"Tfit\",\nFontColor->RGBColor[1, 0, 0]]\) [°C]"}]

Ttopime=13.;
{t1,t2}=tTopim/.FindRoot[nahr1[tTopim]==Ttopime,{tTopim,#*rok}]&/@{0.7,1.4};
Plot[{intPrumBrno[t],nahr1[t]},{t,t1,t2},PlotStyle->{{Blue,Thick,Opacity[0.5]},{Red,Thick,Opacity[0.9]}},AxesLabel->{"t [s]"," Tdenni Tfit [°C]"}];
nDni=Round[((t2-t1)/den)]+1;
Print["topna sezona trva: ",nDni," dni"]

ClearAll[dejAmplitudu];
dejAmplitudu[indexDne_]:=Module[{dataDen,Δt=15*60,amplituda},
dataDen=intTBrno/@Range[(indexDne-1)*den,indexDne*den-Δt,Δt];
amplituda=0.5*(Max[dataDen]-Min[dataDen])
];
amplitudaDen=Mean[dejAmplitudu/@Range[nDni]];
Ttrial[t_]:=amplitudaDen*Sin[2.*Pi*(t/den)]+nahr1[t];
Plot[{intTBrno[t],Ttrial[t]},{t,t1,t2},PlotStyle->{{Blue,Thick,Opacity[0.5]},{Red,Thick,Opacity[0.9]}},AxesLabel->{"t [s]"," Tdenni Tfit [°C]"}];

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+10,0.1];
Pelint=Interpolation[datPel,InterpolationOrder->1];
(*Plot[Pelint[T],{T,-10,14}]*)
datPost={#,Post[#]}&/@Range[Tmin,Tmax+10,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[denFit,denReal];
(*{Ecelk,Eost,ErealEl,EmeanEl,EostMean,Etopna}*)
denFit[indexDne_]:=vydej2[{parsDosStav,{t1+(indexDne-1)*den,t1+indexDne*den},{Ttrial,Ptopny}}][[2]];
denReal[indexDne_]:=vydej2[{parsDosStav,{t1+(indexDne-1)*den,t1+indexDne*den},{intTBrno,Ptopny}}][[2]];

(*dataFitDen=denFit/@Range[1,nDni,1];
dataFitDen>>"dataFitDen";
dataRealDen=denReal/@Range[1,nDni,1];
dataRealDen>>"dataRealDen";*)

mesicFit[indexMesice_]:=vydej2[{parsDosStav,{t1+(indexMesice-1)*den,t1+indexMesice*den},{Ttrial,Ptopny}}][[2]];
mesicReal[indexMesice_]:=vydej2[{parsDosStav,{t1+(indexMesice-1)*den*30,t1+indexMesice*den*30},{intTBrno,Ptopny}}][[2]];

nMesicu=(nDni/30);
(*dataFitMesic=mesicFit/@Range[1,nMesicu,1];
dataFitMesic>>"dataFitMesic";
dataRealMesic=mesicReal/@Range[1,nMesicu,1];
dataRealMesic>>"dataRealMesic";*)

dataRealDen=<<"dataRealDen";
dataFitDen=<<"dataFitDen";
dataRealMesic=<<"dataRealMesic";
dataFitMesic=<<"dataFitMesic";

{EcelkRD,EostRD,ErealElRD,EmeanElRD,EostMeanRD,EtopnaRD}=celkRealDen=Plus@@dataRealDen;
{EcelkFD,EostFD,ErealElFD,EmeanElFD,EostMeanFD,EtopnaFD}=celkFitDen=Plus@@dataFitDen;

{EcelkRM,EostRM,ErealElRM,EmeanElRM,EostMeanRM,EtopnaRM}=celkRealMesic=Plus@@dataRealMesic;

{EcelkFM,EostFM,ErealElFM,EmeanElFM,EostMeanFM,EtopnaFM}=celkRealMesic=Plus@@dataFitMesic;

Print["EelReal= ",ErealElRD,"  EostReal= ",EostRD," EtopnaRD= ",EtopnaRD," J"]
Print["EelRealFit= ",ErealElFD,"  EostFit= ",EostFD," EtopnaFit= ",EtopnaFD," J"]

Print["EelReal= ",ErealElRM,"  EostReal= ",ErealElRM," EtopnaRD= ",EtopnaRM," J"]
Print["EelRealFit= ",ErealElFM,"  EostFit= ",ErealElRM," EtopnaFit= ",EtopnaFM," J"]

Print["EelMean= ",EmeanElRD,"  EostMean= ",EmeanElRD," EtopnaRD= ",EtopnaRD," J"]
Print["EelMeanFit= ",EmeanElRD,"  EostFit= ",EmeanElRD," EtopnaFit= ",EtopnaRD," J"]

Print["EelMean= ",ErealElRM,"  EostMean= ",ErealElRM]
Print["EelMeanFit= ",ErealElFM,"  EostMeanFit= ",EmeanElFM]

Print["Chyba EelRealFit= ",((ErealElRD-ErealElFD)/ErealElRD) *100," Chyba  EostFit= ",((EostRD-EostFD)/EostRD) *100]

Print["Chyba EelRealFit= ",((ErealElRD-ErealElFD)/ErealElRD) *100," Chyba  EostFit= ",((EostRD-EostFD)/EostRD) *100]