%Chapter9 Figure 7(b)
%Fabry-Perot optical resonator
clear;
clf;
err=0.00001;
c=3e8;							%speed of light in vacuum
nr=3.3;							%effective refractive index in cavity
lambda0=1310e-9;				%center emission wavelength
Lc=3e-6;							%Cavity length
r1=0.95;							%Reflectivity of mirror1
r2=0.95;							%Reflectivity of mirror2
r=r1*r2;							%Optical loss per photon round-trip
%r=.8*.4;
F=pi*r^0.5/(1-r)				%Optical Finesse
Fconst=(2*F/pi)^2;
Imax=1/((1-r)^2);				%Peak intensity inside cavity normalized to I0
f0=c*1e-12/lambda0			%center frequency in THz
deltaf=c*1e-12/(2*Lc*nr)	%THz
deltawavelength=lambda0^2/(2*Lc*nr)  %m;
gamma_opt=deltaf/F
Q=f0/gamma_opt

for i=[1:1:10000]
   %x=f0*i*0.001;
   Frequency(i)=(f0+i*0.003);
   Fwavelength(i)=c/Frequency(i);
   x=sin(pi*Frequency(i)/deltaf);
   x2=x*x;
   Intensity(i)=Imax/(1+(Fconst*(x2)));
%y(i)=sin(x*deltax)
end
figure(1);
 plot(Frequency,Intensity);
 %axis([f0,f0+1,0,5]); 
 hold on;
 xlabel('Frequency, \nu (THz)'),ylabel('Intensity');
 title(['r_1 = ',num2str(r1),', r_2 = ',num2str(r2),', n_r = ',...
     num2str(nr),', f_0 = ',num2str(f0),' THz, L_C = ',num2str(Lc),...
     ' m, \lambda_0 = ',num2str(lambda0),' m']);
 hold off;
figure(2);
 plot(Fwavelength*10^-12,Intensity);
xlabel('Wavelength, \lambda (m)'),ylabel('Intensity');
title(['r_1 = ',num2str(r1),', r_2 = ',num2str(r2),', n_r = ',...
    num2str(nr),', f_0 = ',num2str(f0),' THz, L_C = ',num2str(Lc),...
    ' m, \lambda_0 = ',num2str(lambda0),' m']);
