Question

The Language is MATLAB

Modify Software Write an interactive computer program that prompts for the following inputs:

relative permittivity

relative permeability

conductivity in S/m

frequency in MHz

and generates as outputs the following variables

Propagation constant (gamma)

attenuation constant (alpha)

phase constant (beta)

intrinsic impedance (eta) to include the option to analyze reflection and transmission of plane waves between two lossy materials.

The program should prompt the user for the parameters of medium 1, the parameters of medium 2 and the frequency.

The program should output:

A plot of the magnitude of the reflection coefficient versus incidence angle (from 0 to 90 degrees) (include both TE and TM cases on the same plot)

A plot of the magnitude of the transmission coefficient versus incidence angle (from 0 to 90 degrees) (include both TE and TM cases on the same plot).

Submit with your code the plots for the following two cases:

Case 1: medium 1 free space, medium 2: epsilonr=2.25, sigma=0, mur=1

Case 2: medium 1: epsilonr=2.25, sigma=0, mur=1 medium 2: free space

, the transmission coefficient can be greater than 1, the reflection coefficient cannot be greater than 1.The transmitted power however relative to the incident power cannot be greater than 1. same for the relative reflected power.

For the assignment you need to plot the magnitude (absolute value) of the coefficients. For extra credit you can also if you want plot the relative reflected and transmitted powers.

Answer 1

clear all close all

thetai=0:0.1:90; %angle of incidence (degrees) lambda=670; %vacuum wavelength (nm) h=[NaN,50,NaN,NaN]; %film thicknesses in nm, equal in length to n, start and end with NaN pol=1; %polarization, 1 for p and 0 for s n=[1.52,sqrt(Au(lambda)),1.4,1.33]; %refractive index data, NaN for frequency dependence ds=0:10:100; %film thicknesses

for b=1:length(ds) h(3)=ds(b); for a=1:length(thetai) [FR(a,b),FT(a,b),FA(a,b)]=Fresnel(lambda,thetai(a),h,n,pol); end disp([num2str(b/length(ds)*100) '% done...']) end

%plot results: figure hold on plot(thetai,FT,'b') plot(thetai,FR,'r') plot(thetai,FA,'g') xlim([thetai(1),thetai(length(thetai))]) ylim([0 1]) xlabel('incident angle (degrees)') ylabel('Fresnel coefficient') title('Fresnel coefficients for transmission (blue), reflection

and absorption (green)')

end

function [FR,FT,FA]=Fresnel(lambda,thetai,h,n,pol)

%Snell's law: theta(1)=thetai*pi/180; for a=1:length(n)-1 theta(a+1)=real(asin(n(a)/n(a+1)*sin(theta(a))))-1i*abs(imag(asin(n(a)/n(a+1)*sin(theta(a))))); end

%Fresnel coefficients: if pol==0 %formulas for s polarization for a=1:length(n)-1 Fr(a)=(n(a)*cos(theta(a))-n(a+1)*cos(theta(a+1)))/(n(a)*cos(theta(a))+n(a+1)*cos(theta(a+1))); Ft(a)=2*n(a)*cos(theta(a))/(n(a)*cos(theta(a))+n(a+1)*cos(theta(a+1))); end elseif pol==1 %formulas for p polarization for a=1:length(n)-1 Fr(a)=(n(a)*cos(theta(a+1))-n(a+1)*cos(theta(a)))/(n(a)*cos(theta(a+1))+n(a+1)*cos(theta(a))); Ft(a)=2*n(a)*cos(theta(a))/(n(a)*cos(theta(a+1))+n(a+1)*cos(theta(a))); end end

%phase shift factors: for a=1:length(n)-2 delta(a)=2*pi*h(a+1)/lambda*n(a+1)*cos(theta(a+1)); end

%build up transfer matrix: M=[1,0;0,1]; %start with unity matrix for a=1:length(n)-2 M=M*1/Ft(a)*[1,Fr(a);Fr(a),1]*[exp(-1i*delta(a)),0;0,exp(1i*delta(a))]; end M=M*1/Ft(length(n)-1)*[1,Fr(length(n)-1);Fr(length(n)-1),1];

%total Fresnel coefficients: Frtot=M(2,1)/M(1,1); Fttot=1/M(1,1);

%special case of single interface: if length(n)==2 Frtot=Fr(1); Fttot=Ft(1); end

%total Fresnel coefficients in intensity: FR=(abs(Frtot))^2; FT=(abs(Fttot))^2*real(n(length(n))*cos(theta(length(n))))/real(n(1)*cos(theta(1))); FA=1-FR-FT;

end

function epsilon=Au(lambda)

%analytical formula for gold based on wavelength in nm, fits J&C data: epsiloninf=1.54; lambdap=143; gammap=14500; A1=1.27; lambda1=470; phi1=-pi/4; gamma1=1900; A2=1.1; lambda2=325; phi2=-pi/4; gamma2=1060;

%other parameters, worse fit to J&C but seems more accurate often: %epsiloninf=1.53; %lambdap=155; %gammap=17000; %A1=0.94; %lambda1=468; %phi1=-pi/4; %gamma1=2300; %A2=1.36; %lambda2=331; %phi2=-pi/4; %gamma2=940;

for a=1:length(lambda) epsilon(a)=epsiloninf-1/(lambdap^2*(1/lambda(a)^2+1i/(gammap*lambda(a))))... +A1/lambda1*(exp(phi1*1i)/(1/lambda1-1/lambda(a)-1i/gamma1)+exp(-phi1*1i)/(1/lambda1+1/lambda(a)+1i/gamma1))... +A2/lambda2*(exp(phi2*1i)/(1/lambda2-1/lambda(a)-1i/gamma2)+exp(-phi2*1i)/(1/lambda2+1/lambda(a)+1i/gamma2)); end

end