Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
{{Pel,Ptop},info}=<<"struct";
(*Ta vstupni vzduch do TC, Tw vystupni voda z TC*)
Plot3D[{Pel[Ta,Tw],Ptop[Ta,Tw]},{Ta,-40,40},{Tw,20,60}];
pel[Ta_,Tw_]:=1000*Pel[Ta,Tw];
ptop[Ta_,Tw_]:=1000*Ptop[Ta,Tw];
Plot3D[{pel[Ta,Tw],ptop[Ta,Tw]},{Ta,-40,40},{Tw,20,60}];

(*simulace tepelnych ztrat domu*)
ThraniceSezony=13;
Pztraty[Ta_]:=1000*(ThraniceSezony-Ta)/;Ta≤ThraniceSezony;
Pztraty[Ta_]:=0/;Ta>ThraniceSezony;
Plot[Pztraty[Ta],{Ta,-20,20}];

(*identifikace vlastnosti domu a otopne soustavy*)
Tvypoctova=-12;
Pvypoctovy=Pztraty[Tvypoctova];
Print["vypoctova ztrata= ",(Pvypoctovy/(	1000))," kW"]
Tw1Nom=40;
Tw2Nom=35;
c=4186.;
Tin=20.;
rovniceZnamehostavu={Pvypoctovy==c*msteckou1*(Tw1Nom-Tw2Nom),
Pvypoctovy==kS1*(((Tw1Nom+Tw2Nom)/2)-Tin)
};
{kS,msteckou}={kS1,msteckou1}/.Solve[rovniceZnamehostavu][[1]]

(*rovnice hranice bivalence*)
rovniceBivalence={msteckou*c*(Tw1-Tw2)==Pztraty[Ta],
Pztraty[Ta]==kS*(((Tw1+Tw2)/2)-Tin),
Pztraty[Ta]==ptop[Ta,Tw1]
};
solHranice=Quiet@FindRoot[rovniceBivalence,{Tw1,Tw1Nom},{Tw2,Tw2Nom},{Ta,0}];
THraniceBivalence=Ta/.solHranice;
Print["hranice bivalence geft ",THraniceBivalence," °C"]

rovniceNadBivalenci[Ta_]:={msteckou*c*(Tw1-Tw2)==Pztraty[Ta],
Pztraty[Ta]==kS*(((Tw1+Tw2)/2)-Tin),
Pztraty[Ta]==kVyuziti*ptop[Ta,Tw1]
};
solNadBivalenci[Ta_]:=Quiet@Union[FindRoot[rovniceNadBivalenci[Ta],{Tw1,Tw1Nom},{Tw2,Tw2Nom},{kVyuziti,0.5}],{Pk->0}];
solNadBivalenci[5]

(*teplota pod hranicí bivalence, vytápění TČ s pomocí elektrokotle*)rovnicePodBivalenci[Ta_]:={msteckou*c*(Tw1-Tw2)==ptop[Ta,Tw1],
ptop[Ta,Tw1]==kS*(((Tw3+Tw2)/2)-Tin),
Pztraty[Ta]==ptop[Ta,Tw1]+Pk,
Pk==msteckou*c*(Tw3-Tw1)

};
ΔTpoc=1.;
solPodBivalenci[Ta_]:=Quiet@Union[FindRoot[rovnicePodBivalenci[Ta],{Tw1,Tw1Nom-ΔTpoc},{Tw3,Tw1Nom},{Tw2,Tw2Nom-0.8ΔTpoc},{Pk,1050}],{kVyuziti->1}];
solPodBivalenci[-15]

ClearAll[sol]
sol[Ta_]:=solNadBivalenci[Ta]/;Ta≥THraniceBivalence;
sol[Ta_]:=solPodBivalenci[Ta]/;Ta<THraniceBivalence;
sol[5]
sol[-15]

prikonTC[Ta_]:=kVyuziti*pel[Ta,Tw1]/.sol[Ta];
prikonKotle[Ta_]:=Pk/.sol[Ta];
vykonTC[Ta_]:=kVyuziti*ptop[Ta,Tw1]/.sol[Ta];

intPelTC=Interpolation[Table[{Ta,prikonTC[Ta]},{Ta,-20,35,0.5}]];
intVykonTC=Interpolation[Table[{Ta,vykonTC[Ta]},{Ta,-20,35,0.5}]];
intPk=Interpolation[Table[{Ta,prikonKotle[Ta]},{Ta,-20,35,0.5}]];
CHarTČElk=Plot[{intPelTC[Ta],intVykonTC[Ta],intPk[Ta],Pztraty[Ta]},{Ta,-15,20},PlotStyle->{Blue,{Red,Dashed},Magenta,{Black,Opacity[0.5]}},AxesLabel->{"Ta [°C]","P[W]"},GridLines->Full, PlotLegends->Placed["PelTČ VýkonTČ VýkonElektrokotle ZtrátyObjektu", Below]]
Export["CHarTČELK.jpg",CHarTČElk]

teploty=(<<"teplotyBrno");
posun=1*365.25*24*3600;
teplDenniPom=Interpolation[teploty];
tmaxPom=teplDenniPom[[1,1,2]];
(*resamplovani teplot*)
teplotypom={#,teplDenniPom[#]}&/@Range[0.,tmaxPom,1800];
teplDenni1=Interpolation[teplotypom];

teplDenni=teplDenni1[#+posun]&;

teplDenniPom[[1,1]]

posun

tmax=posun;
rcepom=Union[Thread[{intPelTC[Ta],intVykonTC[Ta],intPk[Ta],Pztraty[Ta]}=={eel'[t],etop'[t],eKotel'[t],eZtraty'[t]}],{eel'[t],etop'[t],eKotel'[t],eZtraty'[t]}/.v_'[t]:>v[0]==0];
(*teplota=9+10*Sin[2*((Pi*t)/(tmax=365.25*24*3600))]+8*Sin[2*((Pi*t)/(24*3600))];*)

rce=rcepom/.Ta->teplDenni[t];
nezname={eel'[t],etop'[t],eKotel'[t],eZtraty'[t]}/.v_'[t]:>v;
vysledek=NDSolve[rce,nezname,{t,0,tmax}][[1]];

Print["elektrika TC = ",((eel[tmax])/(10^9))/.vysledek," GJ"]
Print["elektrika elektrokotel = ",((eKotel[tmax])/(10^9))/.vysledek," GJ"]
Print["potreba objektu = ",((eZtraty[tmax])/(10^9))/.vysledek," GJ"]
Print["cerpadlo dodalo = ",((etop[tmax])/(10^9))/.vysledek," GJ"]

Print["COP cerpadlo = ",((etop[tmax])/(eel[tmax]))/.vysledek," -"]
Print["COP objekt = ",((eZtraty[tmax])/(eel[tmax]+eKotel[tmax]))/.vysledek," -"]