function FresnelR
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Fresnel reflectivity characteristics

figure; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----- Intensity reflectivity for different values of b_mu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.35 0.7 0.3 0.2])

q=0.01:0.001:2.5;

b=0.1; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
Rq=rq.*conj(rq); Rn=Rq.*(q.^4)*(2^4);
iq1=find(q<1.5); iq2=find(q>1.5);
[ax,h1,h2]=plotyy(q(iq1),Rq(iq1),q(iq2),Rn(iq2)); 
set(h1,'color','r'); set(h2,'color','r')

axis([0 2.6 0 1.1])
yl=str2mat('    R(q)    ','  Fresnel  ','reflectivity');
text(-1.15,0.5,yl,'FontName','Times','rotation',90,'horizontalalignment','center');

b=0.05; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
Rq=rq.*conj(rq); Rn=Rq.*(q.^4)*(2^4);
line(q(iq1),Rq(iq1),'color','g')
axes(ax(2)); line(q(iq2),Rn(iq2),'color','g')

b=0.01;  qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
Rq=rq.*conj(rq); Rn=Rq.*(q.^4)*(2^4);
axes(ax(1)); line(q(iq1),Rq(iq1),'color','b')
axes(ax(2)); line(q(iq2),Rn(iq2),'color','b')

b=0.001; axes(ax(1)); qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
Rq=rq.*conj(rq); Rn=Rq.*(q.^4)*(2^4);
line(q(iq1),Rq(iq1),'color','m')
axes(ax(2))
text(3.2,1.6,'R(q).(2Q/Q_c)^4','FontName','Times','rotation',90,...
     'horizontalalignment','center')
axis([0 2.6 1 2.1])
set(ax,'Ycolor',[0 0 0],'Ytick',[0.5 1.0 1.5 2.0 2.5],'Xticklabels',[])
set(ax,'FontName','Times','Fontsize',12,'box','on');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----- Penetration length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.35 0.5 0.3 0.2])

q=0.01:0.01:1.4;

b=0.1; qp=sqrt(q.^2-1+2*sqrt(-1)*b); 
line(q,1./imag(qp),'color','r')

b=0.05; qp=sqrt(q.^2-1+2*sqrt(-1)*b); 
line(q,1./imag(qp),'color','g')

b=0.01; qp=sqrt(q.^2-1+2*sqrt(-1)*b);
line(q,1./imag(qp),'color','b')

b=0.001; qp=sqrt(q.^2-1+2*sqrt(-1)*b);
line(q,1./imag(qp),'color','m')

axis([0 2.6 0 1000]);
yl=str2mat('  \Lambda Q_c',' Penetration ','   length   ');
text(-1.1,30,yl,'FontName','Times','rotation',90,'horizontalalignment','center');
set(gca,'Xticklabels',[],'Yscale','log','Ytick',[1 10 100])
set(gca,'FontName','Times','FontSize',12,'box','on')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----- Evanescent intensity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.35 0.3 0.3 0.2])

q=0.01:0.01:2.5;

b=0.1; qp=sqrt(q.^2-1+2*sqrt(-1)*b); ttq=2*q./(q+qp);
h1=line(q,ttq.*conj(ttq),'color','r');

b=0.05; qp=sqrt(q.^2-1+2*sqrt(-1)*b); ttq=2*q./(q+qp);
line(q,ttq.*conj(ttq),'color','g')

b=0.01; qp=sqrt(q.^2-1+2*sqrt(-1)*b); ttq=2*q./(q+qp);
line(q,ttq.*conj(ttq),'color','b')

b=0.001; qp=sqrt(q.^2-1+2*sqrt(-1)*b); ttq=2*q./(q+qp);
line(q,ttq.*conj(ttq),'color','m')

axis([0 2.6 0 4]);
yl=str2mat('   T(q)   ','Evanescent','intensity ');
hl=text(-1.15,2,yl);
set(hl,'FontName','Times','rotation',90,'horizontalalignment','center');
set(gca,'Xticklabels',[],'Ytick',[1 2 3],'Yticklabel',['1';'2';'3'])
set(gca,'FontName','Times','FontSize',12,'box','on')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----- Phase shift of reflected wave
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.35 0.1 0.3 0.2])

q=0.01:0.01:2.5;
b=0.1; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
line(q,angle(rq),'color','r')

b=0.05; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
line(q,angle(rq),'color','g')

b=0.01; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
line(q,angle(rq),'color','b')

b=0.001; qp=sqrt(q.^2-1+2*sqrt(-1)*b); rq=(q-qp)./(q+qp);
line(q,angle(rq),'color','m')

axis([0 2.6 -pi pi/4]);
yl=str2mat('Phase shift of','   reflected  ','    wave      ');
text(-1.,-pi/2,yl,'FontName','Times','rotation',90,'horizontalalignment','center');
text(-0.25,0,'0','horizontalalignment','center','FontName','Times','FontSize',12);
text(-0.25,-pi,'-\pi','horizontalalignment','center','FontSize',14);
xlabel('Q/Q_c or  \alpha/\alpha_c')
set(gca,'Yticklabels',[],'Ytick',[-pi -pi/2 0 pi/2],'Xtick',[0.5 1.0 1.5 2.0 2.5])
set(gca,'FontName','Times','FontSize',12,'box','on')