clear all
clc

%KE2 SIM1 WS 14/15
%Aufgabe 3a Flachdach
%Julian Kammerer 172276
%Emanuel Forko 172061

Dichte=2500;                                            %Rohdichte kg/m³
Lambda=2.3;                                             %Wärmeleitfähigkeit [W/(m*k)]
Dicke= 0.20;                                            %Dicke in [m]
Rsi=0.1;                                                %Übergangswiderstand innen
Rse=0.04;                                                %Übergangswiderstand außen
Hr=4.5;                                                  %radiativer Wärmeübergangskoeffizent W/m²*k
Hcv=20;                                                 %konvektiver Wärmeübergangskoeffizent
Hsi=24.5;                                                %Gesamtwiderstand
c=880;                                                  %spez. Wärmekapazität in K/kg*K
Zaehler=1; 


%Delta x und delta t festelegen

x=10;                                                    %Anzahl der Streckenabschnitte

dx=0.02;         %dx=Dicke/x;
t=1/60;                           

dt1=172;         %dt1=((dx*Dichte*c)/(2*(Lambda/dx)+Hsi));                  %außen
dt2=191;         %dt2=((dx*dx)*Dichte*c)/(2*Lambda);                        %im Bauteeil


dt=150;                                                    %Wert ist durch 3600teilbar und der kleinste von dt1,2
%%%%%%
% Wetterdaten interpoliert;


d = load('wetterdaten.prn');
Monat=d(:,1);
Tag=d(:,2);
Stunde=d(:,3);
direkt=d(:,4);
diffus=d(:,5);
temp=d(:,6);
Wind=d(:,7);
feuchte=d(:,8);

GS=(direkt+diffus);             %Geammtstrahlung ausrechnen
T=(temp+273);                   %Temperaturen in K umrechnen
I=(GS/3.6);                     %Ummrechnung von [Kj/m²h] auf [W/m²]

vqw=[ ]

%Interpol. von SZtemp. und WZtemp. (Sommer- und Winterzeittemperaturen)

%winterwoche                                 
xq=(1:0.043:168);
xw=(1:1:168);                         %letzte Dezember Woche
yw=T(6569:8760);
vqw=interp1(xw,yw,xq,'spline');

Iw=I;
vqiw=interp1(xw,Iw,xq,'spline');

%Sommerwoche
xs=(1:1:168);
ys=T(4345:4512);
vqs=interp1(xs,ys,xq,'spline');

Is=I;
vwis=interp1(xs,Is,xq,'spline');



Talt=[293 293 293 293 293 293 293 293 293 293 293 293];
Tneu=[0 0 0 0 0 0 0 0 0 0 0 0];

hr=4.5;
hcv=20;                  %Stadardwert
alpha=0.55;              %Abstrahlungswert Beton

for Zeitk=1:4016
    
    Tneu(1,1)=(((Talt(1,2)-Talt(1,1))*Lambda/dx+hcv*(vqiw(1,Zeitk)-Talt(1,1))+hr+(0.0552*(vqw(1,Zeitk))^1.5-Talt(1,1))+vqiw(1,Zeitk)*alpha*2*1/60/(0.02*Dichte*c)+Talt(1,1)));
    
    
    Tneu(1,10)=(((hr+hcv)+(293-Talt(1,10))+Lambda/0.02+(Talt(1,9)-Talt(1,10)))*2*1/60)/(0.02*Dichte*c)+Talt(1,10);
    
    for Schichtk=2:9;
        
    Tneu(1,Schichtk)=Lambda/(Dichte*c)*((1/60)/0.02^2)*Talt(1,Schichtk)+Talt(1,Schichtk-1)+Talt(1,Schichtk);
    
    end
    
    Talt=Tneu;
    TGraph1=Talt-273;
    figure=(1);
    subplot(1,2,2)
    plot=TGraph1;
    ylim([0 40])
    xlim([1 10])
    title('Temperaturverlauf im Flachdach "Winterwoche"')
    xlabel('Bauteildicke 0.20 [m]')
    ylabel('Temperatur in [°C]')
    
end  
    
for Zeitk=1:4016
    
    Tneu(1,1)=(((Talt(1,2)-Talt(1,1))*Lambda/0.02+hcv*(vqis(1,Zeitk)-Talt(1,1))+hr+(0.0552*(vqs(1,Zeitk))^1.5-Talt(1,1))+vqis(1,Zeitk)*alpha*2*1/60/(0.02*Dichte*c)+Talt(1,1)));
    
    
    Tneu(1,10)=(((hr+hcv)+(293-Talt(1,10))+Lambda/0.02+(Talt(1,9)-Talt(1,10)))*2*1/60)/(0.02*Dichte*c)+Talt(1,10);
    
    for Schichtk=2:9;
        
    Tneu(1,Schichtk)=Lambda/(Dichte*c)*((1/60)/0.02^2)*Talt(1,Schichtk)+Talt(1,Schichtk-1)+Talt(1,Schichtk);
    
    end
    
    Talt=Tneu;
    TGraph2=Talt-273;
    figure=(1);
    subplot(1,2,1)
    plot=TGraph2;
    ylim([0 40])
    xlim([1 10])
    title('Temperaturverlauf im Flachdach "Sommerwoche"')
    xlabel('Bauteildicke 0.20 [m]')
    ylabel('Temperatur in [°C]')

end