theta = linspace(10,80)*pi/180;
da=1e-9:3e-7;
na = [1.30 1.31 1.32 ] ;
P = zeros(length(na),length(theta)) ;
for i = 1:length(na)
P(i,:) = myfunc(theta,na(i),da) ;
end
plot(theta,P)
legend({'na = 1.30', 'na = 1.31', 'na = 1.32'})
% to plot(da,p)
function y=myfunc(theta,na,da)
l=6330*(10)^(-10);
k0=2*pi/l;
%na=1.340;
nc=1.33;
ec=nc^2;
ea=na^2;
nf=1.59;
ef=nf^2;
nm=0.064+1i*4;
em=nm^2;
ns=1.50;
es=ns^2;
df=0.5e-6;
dm=0.03e-6;
%da=0.06e-6;
np=1.340;
x=np*sin(theta);
kc=sqrt(x.^2-ec.*k0.^2);
ka=sqrt(ea.*k0.^2-x.^2);
kf=sqrt(ef.*k0.^2-x.^2);
km=sqrt(em.*k0.^2-x.^2);
ks=sqrt(x.^2-es.*k0^.2);
rfm=(kf-km)./(kf+km);
rms=(km-ks)./(ks+km);
rfa=(kf-ka)./(kf+ka);
rac=(ka-kc)./(ka+kc);
rsm=(ks-km)./(ks+km);
rmf=(km-kf)./(kf+km);
rfac=((rfa+rac.*exp(2*1i.*ka.*da))./(1+rfa.*rac.*exp(2*1i.*ka.*da)));
rmfac=((rmf+rfac.*exp(2*1i.*kf.*df))./(1+rmf.*rfac.*exp(2*1i*kf.*df)));
r=((rsm+rmfac.*exp(2*1i.*km.*dm))./(1+rsm.*rmfac.*exp(2*1i.*km.*dm)));
r1=conj(r);
y=r.*r1;
end
