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

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]

{PPrum,PFict,PReal,PKonst,TReal,TPrum}=<<"vykony";

ClearAll[zpracuj];
zpracuj[fceP_,fceT_,popis_]:=Module[{data,ndat=100,Δt,ta,tb,tStart,interval,fit,pl1,pl2,ploty1,errPlot},
interval=0.95*den;
Δt=0.95*(den/(ndat-1));
{ta,tb}=fceP[[1,1]];
tStart=0.5*(ta+tb);
data={fceT[#],fceP[#]}&/@Range[tStart,tStart+interval,Δt];
fit[T_]=Fit[data,{1,T},T];
pl1=Plot[fceP[t],{t,tStart,tStart+interval},AxesLabel->{"t [s]","P [W]"},PlotLabel->popis];
pl2=Plot[fit[fceT[t]],{t,tStart,tStart+interval},AxesLabel->{"t [s]","P [W]"},PlotLabel->popis,PlotStyle->Red];
ploty1=Show[pl1,pl2,PlotRange->All];
errPlot=Plot[((fceP[t]-fit[fceT[t]])/(fceP[t]))*100,{t,tStart,tStart+interval},AxesLabel->{"t [s]","δ [%]"},PlotLabel->popis,PlotStyle->Red];
{ploty1,errPlot}


];
{ga,gb}=zpracuj[#,TReal,"PReal"]&@PReal
{gc,gd}=zpracuj[#,TReal,"PFict"]&@PFict
{ge,gf}=zpracuj[#,TReal,"PPrum"]&@PPrum
{gg,gh}=zpracuj[#,TReal,"PKonst"]&@PKonst
{gi,gk}=zpracuj[#,TPrum,"PKonst"]&@PKonst

exp[{ga,gb,gc,gd,ge,gf,gg,gh,gi,gk}]

ClearAll[exp];
SetAttributes[exp,{HoldAll,Listable}];
exp[co_]:=Module[{a,graf},a=ToString[Unevaluated[co]];graf:=Evaluate[co];Export[Evaluate[a<>".png"],graf,ImageResolution->300]];

ClearAll[zpracuj1,T];
zpracuj1[fceP_,fceT_,popis_]:=Module[{data,ndat=100,Δt,ta,tb,tStart,interval,fit,pl1,pl2,ploty1,errPlot},
interval=0.95*den;
Δt=0.95*(den/(ndat-1));
{ta,tb}=fceP[[1,1]];
tStart=0.5*(ta+tb);
data={fceT[#],fceP[#]}&/@Range[tStart,tStart+interval,Δt];
fit=Fit[data,{1,T},T]
];
fitReal[T_]=zpracuj1[PReal,TReal,"PReal"]
fitFict[T_]=zpracuj1[PFict,TReal,"PFict"]
fitPrum[T_]=zpracuj1[PPrum,TReal,"PPrum"]
fitRealKonst[T_]=zpracuj1[PKonst,TReal,"PKonst"]
fitKonstKonst=zpracuj1[PKonst,TPrum,"PKonst"]

ClearAll[vydej2];
vydej2[{znamyStav:{Ptop_,T1znamy_,ΔT_,Tin_},{tStart_,tEnd_},{Texterni_,Ptopny_}}]:=Module[{dosZnamyStav=identifikaceStav[znamyStav],rce,Tvzduch,dos,rn,dosss,datPel,Tmin=-30.,Tmax=45,intPel,fpom,plP,plt1t2,sole,ple,Tmean,EmeanEl,ErealEl,Etopna,Pftopny,COPTpum,COPReal},
rce[Ta_]:={m*c*(T1-T2)==Ptopny[Ta],Ptopny[Ta]==k*(((T1+T2)/2)-Tin),COP[Ta,T1]*Pel==Ptopny[Ta]}/.dosZnamyStav;
dos[Taa_]:=Union[FindRoot[rce[Taa],{Pel,1000},{T2,40},{T1,25}],{Ta->Taa}]/;Taa≤Tmez;
dos[Taa_]:={Ta->Taa,Pel->0}/;Taa>Tmez;
rn=Range[Tmin,Tmax,0.3];
dosss=dos/@rn;
fpom[{Ta_,Pel_}]:=If[Ta≤Tmez,{Ta,Pel},{Ta,0.}];
datPel=fpom/@({Ta,Pel}/.dosss);
intPel=Interpolation[datPel,InterpolationOrder->2];
plP=Plot[((intPel[Ta])/1000.),{Ta,Tmin,Tmax},AxesLabel-> {"ta [°C]","Pel [kW]"},PlotStyle->{Red,Thick},GridLines->Automatic];
plt1t2=Plot[((intPel[Texterni[t]])/1000.),{t,tStart,tEnd},AxesLabel-> {"ta [°C]","Pel [kW]"},PlotStyle->{Magenta,Thick},GridLines->Automatic,AxesOrigin->{tStart,0}];
Pftopny[t_]:=If[Ptopny[Texterni[t]]≥0,Ptopny[Texterni[t]],0];
sole=NDSolve[{e[tStart]==0,e'[t]==intPel[Texterni[t]],TforMean[tStart]==0,TforMean'[t]==Texterni[t],eTopna[tStart]==0,eTopna'[t]==Pftopny[t]},{e,TforMean,eTopna},{t,tStart,tEnd},MaxSteps->(10^6)][[1]];
ple=Plot[((e[t]/.sole)/(10^6)),{t,tStart,tEnd},AxesLabel-> {"t [s]","E [MJ]"},PlotStyle->{Black,Thick},GridLines->Automatic,AxesOrigin->{tStart,0}];

Tmean=((TforMean[tEnd]-TforMean[tStart])/(tEnd-tStart))/.sole;
EmeanEl=intPel[Tmean]*(tEnd-tStart);

ErealEl=(e[tEnd]-e[tStart])/.sole;
Etopna=(eTopna[tEnd]-eTopna[tStart])/.sole;
{COPTpum,COPReal}=(({((EmeanEl+1)/Etopna),((ErealEl+1)/Etopna)})^(-1));
{Tmean,{ErealEl,Etopna,COPReal}}

];

{ta,tb}=PReal[[1,1]];
ClearAll[vydej3];(*overeni pro Ptopny jako funkce casu*)
vydej3[{znamyStav:{Ptop_,T1znamy_,ΔT_,Tin_},{tStart_,tEnd_},{Texterni_,Ptopny_}}]:=Module[{dosZnamyStav=identifikaceStav[znamyStav],rce,Tvzduch,dos,rn,dosss,datPel,Tmin=-30.,Tmax=45,fpom,plP,plt1t2,sole,ple,Tmean,EmeanEl,ErealEl,Etopna,Pftopny,COPTpum,COPReal,vykon},
rce={m*c*(T1-T2)==Ptopfcet,Ptopfcet==k*(((T1+T2)/2)-Tin),COP[Ta,T1]*Pel==Ptopfcet}/.dosZnamyStav;
vykon[t_]=(Pel/.Quiet@Solve[rce,{T1,T2,Pel}][[1]])/.{Ta->Texterni[t],Ptopfcet->Ptopny[t]};

sole=NDSolve[{e[tStart]==0,e'[t]==vykon[t],TforMean[tStart]==0,TforMean'[t]==Texterni[t],eTopna[tStart]==0,eTopna'[t]==Ptopny[t]},{e,TforMean,eTopna},{t,tStart,tEnd},MaxSteps->(10^6)][[1]];
ple=Plot[((e[t]/.sole)/(10^6)),{t,tStart,tEnd},AxesLabel-> {"t [s]","E [MJ]"},PlotStyle->{Black,Thick},GridLines->Automatic,AxesOrigin->{tStart,0}];
Tmean=((TforMean[tEnd]-TforMean[tStart])/(tEnd-tStart))/.sole;
ErealEl=(e[tEnd]-e[tStart])/.sole;
Etopna=(eTopna[tEnd]-eTopna[tStart])/.sole;

COPReal=((((ErealEl/Etopna)))^(-1));
{Tmean,{ErealEl,Etopna,COPReal}}

];
shift=0*den;
aa=vydej3[{parsDosStav,{ta+shift,ta+den+shift},{TReal,PReal}}]
cc=vydej2[{parsDosStav,{ta+shift,ta+den+shift},{TReal,fitReal}}]
((aa[[2]]-cc[[2]])/(aa[[2]]))*100(*chyby v procentach*)

dd=vydej3[{parsDosStav,{ta+shift,ta+den+shift},{TReal,PReal}}]
bb=vydej2[{parsDosStav,{ta+shift,ta+den+shift},{TReal,fitFict}}]
((dd[[2]]-bb[[2]])/(dd[[2]]))*100(*chyby v procentach*)

((aa[[2]]-bb[[2]])/(aa[[2]]))*100