% Koeffizienten
h = 12.2;

% copper conductivity
kf = 400; % Monier-Vinard
% dielectric material conductivity
km = 0.3; % Monier-Vinard
% copper fill factor
fillfactor = 0:0.01:1;

% Transistorfläche: 8*8 mm² = 64 mm²
% Induktorfläche: 23.5 mm * 24 mm = 564 mm²
% Ac = 23.5 * 24;
Ac = 8*8;
r1 = sqrt(Ac/pi);

% Kupferfläche (heat sink) von 64 mm² bis 9754 mm² (9750 ist maximale
% Fläche für heat spreading)
As = 64:96.9:9754;
r2 = sqrt(As/pi);
% Dicke von 35 mikrometer bis 1035 mikrometer (also 0.035 mm bis 1.035 mm)
t = 0.035:0.010:1.035;

[tt, r2r2] = meshgrid(t,r2);
[ff] = meshgrid(fillfactor);

tau = tt./r2r2;

% effective axial thermal conductvity of each layer 
gamma = (3*ff- 1)*kf + (3*(1-ff)-1)*km; % array
kk = 0.25*(gamma + sqrt(gamma.^2 + 8*kf*km)); %array

epsilon = r1./r2r2;
Bi = h*r2r2./kk;
lambda = pi + 1./(epsilon*sqrt(pi));
Phi = (tanh(lambda.*tau) + (lambda./Bi)) ./ (1 + (lambda .* (tanh(lambda.*tau))./Bi));

Rsp = (1-epsilon)*Phi ./ (pi * kk * r1);
Rcond = tt ./ (pi * kk .* r2r2); 
Rconv = 1 ./ (pi * h * r2r2);

Rtot = Rcond + Rsp + Rconv;

surf(tau(:), ff(:), Rtot)

zlabel('thermal spreading resistance in K/W')
xlabel('tau = t/r2 (dimensionless)')
ylabel('copper fill factor in decimals')