%% board measurements (height Hb, length Lb, width Wb)
Hb = 1.6;
Lb = 150;
Wb = 100;

%% layer parameters %% 

% number of layers 
nl = 3; 
i = 1:nl;

% (array of) layer thicknesses
% bottom layer
ti(1) = 0.035;
ti(2) = 1.53; %prepreg
% top layer
ti(nl) = 0.035;
z = ti;

% copper fill factor of layers 
Phi(1) = 1; % bottom layer
Phi(2) = 0; % no copper fill factor for prepreg layers! :)
Phi(nl) = 0.5; % top layer

%% components / heat sources

% T1H: high-sider power semiconductor of first half-bridge
% T1L: low-sider power semiconductor of first half-bridge
% T2H: high-sider power semiconductor of second half-bridge
% T2L: low-sider power semiconductor of second half-bridge
% L1: inductor of first half-bridge
% L2: inductor of first half-bridge

% number of heat sources
ns = 6;
% j = 1:ns;

% power losses (different operatng points tbd)
qsT1H = 1;
qsT2H = 1;
qsT1L = 1;
qsT2L = 1;
qsL1 = 2;
qsL2 = 2;
qsj = [qsT1H qsT2H qsT1L qsT2L qsL1 qsL2];

% lenght of components / heat sources
LsT1H = 5; LsT2H = 5; LsT1L = 5; LsT2L = 5; LsL1 = 18.2; LsL2 = 18.2;
Lsj = [LsT1H LsT2H LsT1L LsT2L LsL1 LsL2];

% width of components / heat sources
WsT1H = 5; WsT2H = 5; WsT1L = 5; WsT2L = 5; WsL1 = 18.3; WsL2 = 18.3;
Wsj = [WsT1H; WsT2H; WsT1L; WsT2L; WsL1; WsL2];

% x-, y- and z- coordinates for heat source / components center positon
xcT1H = 70; ycT1H  = 80; zcT1H = Hb;
xcT1L = 70; ycT1L = 60; zcT1L = Hb;
xcT2H = 70; ycT2H = 40; zcT2H = Hb;
xcT2L = 70; ycT2L = 20; zcT2L = Hb;
xcL1 = 80; ycL1 = 70; zcL1 = Hb;
xcL2 = 80; ycL2 = 30; zcL2 = Hb;
xcj = [xcT1H; xcT2H; xcT1L; xcT2L; xcL1; xcL2];
ycj = [ycT1H; ycT2H; ycT1L; ycT2L; ycL1; ycL2];
zcj = [zcT1H; zcT2H; zcT1L; zcT2L; zcL1; zcL2];

%% heat transfer parameters
% uniform heat transfer coefficient for convecton and radiaton
% top layer
ht = 12.2;
% bottom layer
hr = 12.2;

%ambient temperature
Tamb = 70;

% copper conductvity
kf = 400; % Monier-Vinard
% dielectric material conductvity
km = 0.3; % Monier-Vinard
% effectve axial thermal conductvity of each layer 
gamma = (3*Phi - 1)*kf + (3*(1-Phi)-1)*km; % array
ke = 0.25*(gamma + sqrt(gamma.^2 + 8*kf*km)); %array

% conduction heat transfer coefficients in x-, y- and z-direction
kxi = ke;
kyi = ke;
kzi = ke;

%% calculation parameters / solution of partial differential equation
% upper limit values for truncated Fourier series
M = 20; %225
syms m;
N = 30; %300 
syms n;

% solution of differential equation for temperature 

for j = 1:ns   
    for m = sym(0):sym(M)
        for n = sym(0):sym(N)

            % Fourier coefficients dedicated to heat sources j
            Amj = 4*(sin(m*pi*Lsj(j)/Lb*2)*cos(m*pi*xcj(j)/Lb))/(m*pi + kroneckerDelta(m)) + Lsj(j)*kroneckerDelta(m)/Lb;
            Bnj = 4*(sin(n*pi*Wsj(j)/Wb*2)*cos(n*pi*ycj/Wb))/(n*pi + kroneckerDelta(n)) + Wsj*kroneckerDelta(n)/Wb;

            % parameters
            D_mni = alpha_mni(nl) + ht*chi_mni(nl)*beta_mni(nl) + (chi_mni(nl) + ht)*gamma_mni(nl);

                for p = 1:(nl-1)
                    N_mni = N_mni * ((exp(r(m,n,nl)*(z(nl-1)-Hb)))/(exp(r(m,n,i)*(z(i-1)-Hb)))) * ((2*exp(-ti(p)*r_mni(p)))/(gamma_mni(p)+chi_mni(p)*beta_mni(p)))
                end

            % summation
            theta_sum = theta_sum + theta_mni(x,y,z);
            Ti = Ti(x,y,z)

        end
    end
end

% functions (must appear at the end of a file)

function theta = theta_mni(x,y,z)
    theta = (qsj(j)/(Lsj(j)*Wsj(j))) * Amj * Bmj * cos(m*pi*x/Lb)*cos(n*pi*y/Wb) * omega_mni(z);
end

function temperature = Ti(x,y,z)
    temperature = Tamb + theta_sum;
end

function r = r_mni(i)
    r = sqrt((m*pi/Lb)^2*kxi(i)/kzi(i)+(n*pi/Wb)^2*kyi(i)/kzi(i));
end

function omega = omega_mni(z) 
    omega = N_mni/D_mni*(omegac_mni(z-z(i-1)) + chi(m,n,i)/kz(i)*(omegas_mni(z-z(i-1))))*exp((z-Hb)*r(m,n,i));
end

function omegac = omegac_mni(z) 
    omegac = 1 + exp(-2*z*r_mni(i));
end

function omegas = omegas_mni(z) 
    omegas = (1-exp(-2*z*r_mni(i)))/r_mni(i); 
end

function alpha = alpha_mni(i)
    alpha = kzi(i) * r_mni(i)^2 * omegas_mni(ti(i));
end

function beta = beta_mni(i)
    beta = omegas_mni(ti(i)) / kzi(i);
end

function gamma = gamma_mni(i)
    gamma = omegac_mni(ti(i));
end

function chi = chi_mni(i) 
    if(i==1)
        chi = hr;
    else
        chi = (alpha_mni(i) + chi_mni*gamma_mni(i))/(gamma_mni(i) + chi_mni(i)*beta_mni(i));
    end
end