Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

(*Upozorneni: prilia mnoho parametru jde volit, tedy nejde v ramci diplomove prace prostor moznosti "prozkoumat" zamerne je pro rozumnou nemodelovatelnost opominuto osluneni budovy (projevilo by se v okrajove podmince a v solernich ziscich), take tepelne zisky uvnitr atd. Tedy v podstate simulujeme prazdnou budovu za zamracenych dni... *)
Ampl=7;
hodina=3600.
den=hodina*24.;
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}];

cm=0.01;
dCihly=30cm;
dIzol=16cm;
dmax=dCihly+dIzol;
λCihly=0.8;
λIzol=0.04;
rocCihly=3000*200;
rocIzol=50*1000;
αin=5;
αext=15;
U=(1/((1/(αin))+(1/(αext))+(dCihly/(λCihly))+(dIzol/(λIzol))));
ClearAll[prechod];
prechod[vlevo_,vpravo_]:=Module[{res,a,b},res=Solve[{a-b==vlevo,a+b==vpravo}];a+b*Tanh[((x-dCihly)/(0.1cm))]/.res[[1]]];
roc=prechod[rocCihly,rocIzol];
λ=prechod[λCihly,λIzol];
pl[co_,popis_]:=Plot[co/.x->xx,{xx,0,dmax},PlotRange->All,AxesLabel->{"x [m]",popis},AxesOrigin->{0,0},PlotStyle->{Hue[0.63],Thickness[0.004]}];
plRoc=pl[roc,"ρc [(J/(K*(m^3)))]"];
plLambda=pl[λ,"λ [(W/(K*m))]"];

Tin=21.5;
ClearAll[vydej];
vydej[{tStart_,nDni_}]:=Module[{nabeh,rce,poc, okrVlevo,okrVpravo,tmax,res,tokTprum,tokZed,pltok,plDelta,pltokAll,toky,naDisk},
nabeh=Tanh[t-tStart];
rce={roc*D[T[x,t],t]==D[λ*D[T[x,t],x],x]};
poc={T[x,tStart]==Tin};
okrVlevo={-λ*D[T[x,t],x]==αin*(Tin-T[x,t])*nabeh}/.x->0;
okrVpravo={-λ*D[T[x,t],x]==αext*(T[x,t]-teplDenni[t])*nabeh}/.x->dmax;
tmax=tStart+nDni*den;
res=NDSolve[Union[rce,poc,okrVlevo,okrVlevo,okrVpravo],{T},{x,0,dmax},{t,tStart,tmax},MaxStepSize->{.05cm,0.1*den}][[1]];
plot3D=Plot3D[T[x,t]/.res,{x,0,dmax},{t,tStart,tmax},AxesLabel->{"x [m]","t [s]","T [°C]"},PlotPoints->100];
tokTprum[t_]:=U*(Tin-teplsec[t]);
tokZed[t_]=(αin*(Tin-T[x,t]/.res))/.x->0;
pltok=Plot[{tokTprum[t],tokZed[t]},{t,tmax-(nDni-1)*den,tmax},AxesOrigin->{tStart,0},PlotStyle->{{Blue,Thick},{Red,Thick}},AxesLabel->{"t [s]","pstr p[t]  [(W/(m^2))]"},GridLines->Automatic];
plDelta=pltok=Plot[((tokZed[t]-tokTprum[t])/(tokTprum[t]))*100,{t,tmax-(nDni-1)*den,tmax},AxesOrigin->{tStart,0},PlotStyle->{Blue,Thickness[0.006]},AxesLabel->{"t [s]","δp [%]"},GridLines->Automatic];
pltokAll=Plot[{tokTprum[t],tokZed[t],U*(Tin-teplDenni[t])},{t,tmax-(nDni-1)*den,tmax},AxesOrigin->{tStart,0},PlotStyle->{{Blue,Thickness[0.006]},{Red,Thickness[0.006]},{Black,Thickness[0.006]}},PlotLabel->{"vykonZedPrum[t] vykonZedReal[t] vykonZedFic[t] [(W/(m^2))]"},GridLines->Automatic];
toky[t_]:={tokTprum[t],tokZed[t],U*(Tin-teplDenni[t])};
naDisk=Table[{t,toky[t]},{t,tStart,tmax,0.1*hodina}];
naDisk>>"datToky";
{pltokAll,plDelta,plot3D}

];
pom=vydej[{10000,15}]//Quiet

(*abychom odhadli vliv sertvacnosti zdi, vytvorime si fiktivni dum s rozumnym vetranim, pomerem ploch zdi, oken a dveri atd... Samosebou pouzivame obalkovou metodu.
dum je krychle o strane a*)
a=7.5;
nOken=16;
S1okno=1.5*1.5;
S1dvere=2.1*1;
nDveri=2;
Sdveri=nDveri*nDveri;
Soken=nOken*S1okno;
Szed=5*(a^2)-Sdveri-Soken;
Szem=(a^2);
Objem=(((a-2*dmax))^3);

(*o vetrani*)
ucinnostRekuperace=0.6;
casVymenyVzduchu=4*hodina;
cp=1000.;
hustotaVzduchu=(349./(Tin+273));
vykonVetrani[t_]:=(1-ucinnostRekuperace)*(Tin-teplDenni[t])*hustotaVzduchu*cp*(Objem/casVymenyVzduchu);
vykonVetraniPrum[t_]:=(1-ucinnostRekuperace)*(Tin-teplsec[t])*hustotaVzduchu*cp*(Objem/casVymenyVzduchu);

(*o dverich a oknech a podlaze*)
Uokno=1.5;
Udvere=0.5;
vykonOknaDvere[t_]:=(Tin-teplDenni[t])*(Udvere*Sdveri+Uokno*Soken);
vykonOknaDverePrum[t_]:=(Tin-teplsec[t])*(Udvere*Sdveri+Uokno*Soken);
dIzolPodlaha=10cm;
Upodlaha=(1/((1/(αin))+(dIzolPodlaha/(λIzol))));
Tzem=0;
vykonZem=(Tin-Tzem)*Szem*Upodlaha;

(*vyuzijeme reseni DFourier-Kirchhoffovy rovnice*)
datToky=<<"datToky";
tok1mkvTprum=Interpolation[datToky/.{t_,{Prum_,Preal_,Pfict_}}:>{t,Prum}];
tok1mkvTreal=Interpolation[datToky/.{t_,{Prum_,Preal_,Pfict_}}:>{t,Preal}];
tok1mkvFict=Interpolation[datToky/.{t_,{Prum_,Preal_,Pfict_}}:>{t,Pfict}];
{tapom,tb}=tok1mkvTreal[[1,1]];
Plot[tok1mkvTreal[t],{t,tapom,tb}];(*slouzi pro zvoleni ta tak, abychom byli za prechodnym dejem*)
ta=6*(10^5);
Plot[tok1mkvTreal[t],{t,ta,tb}];
vykonZedReal[t_]:=Szed*tok1mkvTreal[t];
vykonZedFict[t_]:=Szed*tok1mkvFict[t];
vykonZedPrum[t_]:=Szed*tok1mkvTprum[t];

(*ted vyvorime tri celkove toky odpovidajici situaci s uvazovanim prumerne teploty pro prostup zdi (situace "Prum"), situaci s uvazovanim skutecneho toku zdi podle reseni FKR (situace "Real") a situaci, kdy zed povazujeme za tepelny odpor bez tepelne kapacity (situace "Fict"), v techto pripadech je uvazovano skutecne vetrani a prostup okny, dalsi situaci znacime "Konst", kde bereme prumerne denni teploty i pro vetrani a prostup dvermi a okny*)
vykonOstatni[t_]:=vykonOknaDvere[t]+vykonVetrani[t]+vykonZem;
vykonOstatniPrum[t_]:=vykonOknaDverePrum[t]+vykonVetraniPrum[t]+vykonZem;
PReal[t_]:=vykonZedReal[t]+vykonOstatni[t];
PPrum[t_]:=vykonZedPrum[t]+vykonOstatni[t];
PFict[t_]:=vykonZedFict[t]+vykonOstatni[t];
PKonst[t_]:=vykonZedPrum[t]+vykonOstatniPrum[t];
pltokAll=Plot[{PReal[t],PPrum[t],PFict[t],PKonst[t]},{t,ta,tb},AxesOrigin->{ta,0},PlotStyle->{{Blue,Thickness[0.006]},{Red,Thickness[0.006]},{Black,Thickness[0.006]},{Magenta,Thickness[0.006]}},AxesLabel->{"t [s]","pReal pPrum pFict pKonst [W]"},GridLines->Automatic]

(*pro prehlednost si radeji vytvorime data a ulozime na disk*)
Δt=0.1*hodina;
dPrum=Table[{t,PPrum[t]},{t,ta,tb,Δt}];
dFict=Table[{t,PFict[t]},{t,ta,tb,Δt}];
dReal=Table[{t,PReal[t]},{t,ta,tb,Δt}];
dKonst=Table[{t,PKonst[t]},{t,ta,tb,Δt}];
dTReal=Table[{t,teplDenni[t]},{t,ta,tb,Δt}];
dTPrum=Table[{t,teplsec[t]},{t,ta,tb,Δt}];
{dPrum,dFict,dReal,dKonst,dTReal,dTPrum}>>"dataZed"

(*{dPrum,dFict,dReal,dKonst,dTReal,dTPrum}=<<"dataZed";*)

Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
hodina=3600.;
den=hodina*24.;
dataZed={dPrum,dFict,dReal,dKonst,dTReal,dTPrum}=<<"dataZed";
{PPrum,PFict,PReal,PKonst,TReal,TPrum}=Interpolation/@dataZed;
{ta,tb}=PReal[[1,1]]
Plot[Evaluate@Through[{PPrum,PFict,PReal,PKonst,TReal,TPrum}[t]],{t,ta,tb}];  
{PPrum,PFict,PReal,PKonst,TReal,TPrum}>>"vykony"