lab4_code.m

%s = 500;
s = [1.5, 2, 3, 4, 6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 80, 100, 120, 140, 160, 180, 200, 250, 300, 350, 400, 500, 600, 800, 1000];
a= [-0.980685 -0.771995 -0.563305 -0.354615 -0.145925 0.062765 0.271455  0.480145  0.688835  0.897525  1.106215  1.314905  1.523595  1.732285  1.940975  2.149665  2.358355  2.567045  2.775735  ];
f=[0.00097112 -0.00102152     0.00906965      0.01404316      0.09012         0.30171582      0.99627084      1.3690832       -2.99681171 1.65463068      -0.59399277     0.22329813      -0.10119309     0.05186135 -0.02748647 0.01384932 -0.00599074     0.00190463      -0.0003216];
for jj=1:length(s)
lambda=[];
for j=1:19
lambda(j)=10^( a(j)-log10(s(jj)) );
end


rho=[1000,200,1000];
h=[5,10];
%T=[];

k=length(rho);
T(k,:)=rho(k);
j=1;
while((k>1) & (j<=19))
    T(k-1,j)=( T(k,j)+rho(k-1)*tanh(lambda(j)*h(k-1)) ) / ( 1+(T(k,j)*tanh(lambda(j)*h(k-1)))/rho(k-1) );
    k=k-1;
    j=j+1;
end
k=length(rho);
endrho=[1000,200,1000];
h=[5,10];

rho1(jj)=0;
for j=1:19
rho1(jj)=rho1(jj)+f(j)*T(1,j);
end

end
figure
loglog(s,rho1)