reflect.sci

clear;
epsilon0 = 1/(36*%pi)*10^-9;
mu0 = 4*%pi*10^-7;
epsilon1 = 1;  mu1 = 1; sigma1 = 0;
epsilon2 = 10; mu2 = 1; sigma2 = 0;
epsilon3 = 80; mu3 = 1; sigma3 = 0;
f1 = 300*10^6; f2 = 3*10^9;
omega = 2*%pi*f1;
epsilon1complex = epsilon0*epsilon1*(1-%i*sigma1/(omega*epsilon0*epsilon1));
epsilon2complex = epsilon0*epsilon2*(1-%i*sigma2/(omega*epsilon0*epsilon2));
epsilon3complex = epsilon0*epsilon3*(1-%i*sigma3/(omega*epsilon0*epsilon3));

k1=omega*sqrt(mu0*mu1*epsilon1complex);
k2=omega*sqrt(mu0*mu2*epsilon2complex);
k3=omega*sqrt(mu0*mu3*epsilon3complex);
W1=sqrt(mu0*mu1/epsilon1complex);
W2=sqrt(mu0*mu2/epsilon2complex);
W3=sqrt(mu0*mu3/epsilon3complex);

phipad = 0:%pi/2/1000:%pi/2-%pi/2/1000;
phipr12 = asin(k1/k2*sin(phipad));
phipr13 = asin(k1/k3*sin(phipad));

R12norm = (W2*cos(phipad) - W1*cos(phipr12)) ./ (W2*cos(phipad) + W1*cos(phipr12));
R12paral = (W1*cos(phipad) - W2*cos(phipr12)) ./ (W1*cos(phipad) + W2*cos(phipr12));
R13norm = (W3*cos(phipad) - W1*cos(phipr13)) ./ (W3*cos(phipad) + W1*cos(phipr13));
R13paral = (W1*cos(phipad) - W3*cos(phipr13)) ./ (W1*cos(phipad) + W3*cos(phipr13));
figure(1);
subplot(221);
plot(phipad,abs(R12norm),phipad,abs(R13norm));
subplot(222);
plot(phipad,abs(R12paral),phipad,abs(R13paral))
subplot(223);
plot(phipad,atan(imag(R12norm),real(R12norm)),phipad,atan(imag(R13norm),real(R13norm)));
subplot(224);
plot(phipad,atan(imag(R12paral),real(R12paral)),phipad,atan(imag(R13paral),real(R13paral)));
//figure(1);
//plot(phipad,phipr12,phipad,phipr13,phipad,phipad);