%Eigenschaften Glas
lambda_glas=1.38;
rho_glas=2201;
r_glas=1052;
l_glas=0.0002;

%Eigenschaften Luft
lambda_luft=0.0261;
rho_luft=1.234;
c_luft=1005;
l_luft=0.0004;

%Eigenschaften Polyimid Mantel
lambda_PI=0.34;
rho_PI=1400;
c_PI=1088.568;
l_PI=0.00005;

%Eigenschaften CFK-Wickel-Dings
lambda_CFK=2;4;
rho_CFK=1400;
c_CFK=1088.568;
l_CFK=0.0002;

%Temperaturen Defimieren
T_0=23;
T_rand=100;

%Zeitschritte und Gesamtzeit 
M=101;
t=2;

%Zeitschritt
DELTA_t=t/(M-1);
for j=1:M 
    time(j)=(j-1)*DELTA_t;
end

%Diskrete Punkte pro Schicht
N1=10;
N2=10;
N3=10;
N4=11;

L=l_glas+l_luft+l_PI+l_CFK;
N=N1+N2+N3+N4;

%Längenschritte Länge
DELTA_r1=l_CFK/(N1); 
DELTA_r2=l_PI/(N2); 
DELTA_r3=l_glas/(N3);
DELTA_r4=l_luft/(N4-1); 

%Position jedes Punktes
r1=1:DELTA_r1:l_CFK; 
r2=l_CFK+DELTA_r2:DELTA_r2:l_CFK+l_PI; 
r3=l_CFK+l_PI+DELTA_r3:DELTA_r3:l_CFK+l_PI+l_glas;
r4=l_CFK+l_PI+l_glas+DELTA_r4:DELTA_r4: l_CFK+l_PI+l_glas+l_luft;

r=[r1 r2 r3 r4];

%Anfangstemperatur definieren
for i=1:N   
   T(i,1)=T_0;       
end 



%Schleife
for j=1:(M-1)

    for i=2:(N-1) 
       
%Dirichlet-Randbedingen (Sprunghafte Temperaturänderung am Rand)
 T(1,j) =0.5*( T_rand-T(2,j));

 %Eigenschaften sind ortsabhängig
 if i<=N1 
         rho=rho_CFK; 
         c=c_CFK; 
         k=lambda_CFK; 
         DELTA_r=DELTA_r1;
  r=r1; 
        elseif i<=N2 and i>N1
         rho=rho_PI; 
         c=c_PI; 
         k= lambda_PI; 
         DELTA_r=DELTA_r2; 
  r=r2;
        elseif i<=N3 and i>N2
         rho=rho_glas; 
         c=c_glas; 
         k= lambda_glas; 
         DELTA_r=DELTA_r3; 
  r=r3;
        elseif i<=N4 and i>N3
         rho=rho_luft; 
         c=c_luft; 
         k= lambda_luft; 
         DELTA_r=DELTA_r4; 
  r=r4;

        end 
        %Fourrier-Gleichung-in-Zylinderkoordinaten
         T(i,j+1)=T(i,j)+(k*DELTA_t./(rho*c*DELTA_r^2))*(T(i-1,j)+T(i+1,j)-2*T(i,j))+(T(i-1,j)-T(i+1,j))*(1./r)* (DELTA_t./(rho*c.*DELTA_r^2)); 
    end 
        
%Neumann-Randbedingung im Zentrum =0, für Symmetrie
    T(N-1,j)=T(N-2,j); 
    
end
plot(time,T) 
figure 
plot(r,T) 
figure 
g=T(:,M);   
plot(r,g)