# compare simulation results with analytical results

# Material parameters

K = 1.38; # conduction coefficient
Cp = 709;  # specific heat
rho = 2203;  # density

# Boundary conditions

T_left = 400;
T_right = 300;
L = 100e-6;   # length of slab

# Temp on right edge is assumed to be zero in analytic expression.  So T0 = T_left - T_right

T0 = T_left - T_right;  

# Initializing vectors and matrices

x = linspace(0,100e-6,51);
t = linspace(0,0.005,51);

T = matrix(51,51);

# Analytic calculation

C = sqrt(K/Cp/rho);

for (i=1:51) {
    
    T(1:51,i) = T0*(1-x/L);    
}

for (j=1:51) {
    
    for (n=1:1000) {
        
        T(1:51,j) = T(1:51,j) - 2*T0/pi * (1/n) * exp(-n^2* (C*pi/L)^2 * t(j)) * sin(n*pi/L*x);        
    }    
}

T = T+300;  # set the base temperature to 300 K

# Plot result

# from simulation

data = getresult('HEAT::monitor','temperature');
T_sim = pinch(data.T);
x_sim = pinch(data.x);
t_sim = pinch(data.t);
N_t = length(t_sim);
N_x = length(x_sim);

kern = find(t_sim>=1e-4);
t_loc_1 = kern(1);
kern = find(t_sim>=1e-3);
t_loc_2 = kern(1);

# Temp profile at t = 0.1, 1, and 5 ms

plotxy(x_sim,T_sim(1:N_x,t_loc_1),x,T(1:51,2),'x (m)','T (K)');
holdon;
plotxy(x_sim,T_sim(1:N_x,t_loc_2),x,T(1:51,11));
plotxy(x_sim,T_sim(1:N_x,N_t),x,T(1:51,51));
holdoff;
legend('t = 0.1 ms (simulation)','t = 0.1 ms (analytic)',
't = 1 ms (simulation)','t = 1 ms (analytic)','t = 5 ms (simulation)','t = 5 ms (analytic)');

# Temp varitaion at x=60 micron w.r.t. time

kern = find(x_sim>6e-5);
x_sim_loc = kern(1)-1;

kern = find(x>6e-5);
x_loc = kern(1)-1;

plotxy(t_sim,T_sim(x_sim_loc,1:N_t),t(1:51),T(x_loc,1:51),'t (s)','T (K)','Temperature at x = 60 micron');
legend('simulation','Analytic');