cavity.m

Using spectra to calculate strength of reflections and losses in passive cavities.The wavelength, re

function [cavityvacuum,trans]=cavity(filename, neff, plotoption)
%Running without parameters uses a default of neff=3.26 and plotoption=2
%(display wavelength spectrum and cavity lengths)
%It also prompts the user for the excel filename.
%
%NB: Make sure that xlSheetName, xlWaveRef, and xlPwrdBmRef are set
%correctly below.

%Excel setup
%For HP Spectrum Analyser using TimsVB ims
%xlSheetName='Data';
%xlWaveRef='A16:A2016';
%xlPwrdBmRef='B16:B2016';

%For Nan's tunable laser sweep
xlSheetName='Sheet1';
xlWaveRef='A3:A1503';
xlTrandBmRef='B3:B1503'; %Reflection
xlRefldBmRef='C3:C1503'; %Transmission


%Fixed values
c=3e8;
nplots=4;

%Get filename if none supplied
if (exist('filename')==0)
    [filename filepath]=uigetfile('*.xls','Load excel file');
else
    filepath='';
end;
disp([filepath filename]);

%Set refractive index if none supplied
if (exist('neff')==0)
    neff=3.26;
end;

%Set default plot option
if (exist('plotoption')==0)
    plotoption=2;
end;

%Read spectrum from excel file and pad to 2^n with zeros
xlWave=xlsread([filepath filename],xlSheetName,xlWaveRef);
xlTrandBm=xlsread([filepath filename],xlSheetName,xlTrandBmRef);
xlRefldBm=xlsread([filepath filename],xlSheetName,xlRefldBmRef);

npoints=2^(ceil(log10(length(xlWave))/log10(2)-1e-3));
tranlin=zeros(1,npoints);
refllin=zeros(1,npoints); 
T11=zeros(1,npoints);
S11=zeros(1,npoints);
T21=zeros(1,npoints);
T11win=zeros(1,npoints);
T21win=zeros(1,npoints);


TrandBm=ones(1,npoints)*min(xlTrandBm);
TrandBm(1:length(xlWave))=xlTrandBm;
RefldBm=ones(1,npoints)*min(xlRefldBm);
RefldBm(1:length(xlWave))=xlRefldBm;
lambda=((0:npoints-1)*(xlWave(2)-xlWave(1))+xlWave(1))*1e-9;


%Calculate frequency from wavelength using df=-c.dh/h0^2
lambda0=(lambda(1)+lambda(length(lambda)))/2;
freq0=c/lambda0;
freq=freq0-(c/lambda0^2).*(lambda-lambda0);
%Calculate linear power
Tranlin(1:length(xlWave))=10.^(TrandBm(1:length(xlWave))/10)*1e-3;
Refllin(1:length(xlWave))=10.^(RefldBm(1:length(xlWave))/10)*1e-3;

T11(1:length(xlWave))= 1./(10.^(TrandBm(1:length(xlWave))/10)*1e-3);
S11(1:length(xlWave))= 10.^(RefldBm(1:length(xlWave))/10)*1e-3;
T21 = S11 ./ T11;

T11win(1:length(xlWave)) = hann(length(xlWave)).' .* T11(1:length(xlWave));
T21win(1:length(xlWave)) = hann(length(xlWave)).' .* T21(1:length(xlWave));

figure(1);
    subplot(2,1,1)
    h=plot(lambda*1e6,T11);
    xlabel('Wavelength (um)');
    ylabel('Power (W)');
    title('T11');
    zoomx(h,0.01);

    subplot(2,1,2)
    h=plot(lambda*1e6,T21);
    xlabel('Wavelength (um)');
    ylabel('Power (W)');
    title('T21');
    zoomx(h,0.01);
    
figure(2);
    subplot(2,1,1)
    h=plot(lambda*1e6,T11win);
    xlabel('Wavelength (um)');
    ylabel('Power (W)');
    title('T11');
    zoomx(h,0.01);

    subplot(2,1,2)
    h=plot(lambda*1e6,T21win);
    xlabel('Wavelength (um)');
    ylabel('Power (W)');
    title('T21');
    zoomx(h,0.01);


%Do the FFT and calcualte the x-axis
trans=fft(T11win);
refl=fft(T21win);

transnorm = abs(trans) / abs(trans(1));
reflnorm = abs(refl) / abs(refl(1));

invf=(0:npoints/2)/npoints/(freq(1)-freq(2));
%Convert the 1/df into cavity length 
%(NB: divide by 2 as the cavity resonates at 1/2 periods)
cavityvacuum=c*invf/2;

figure(3)

%figure
subplot(2,1,1)
h=plot(cavityvacuum/neff*1e6,abs(transnorm(1:npoints/2+1)));
xlabel(['Cavity length in neff=' num2str(neff) ' (um)']);
ylabel('Cavity strength (relative units)');
title('Resonant cavity lengths L=c/2f');
xlim([0, 2000]);
ylim([0, max(transnorm(20:1024))])
%ylim([0, 10e-3])


subplot(2,1,2)
h=plot(cavityvacuum/neff*1e6,abs(reflnorm(1:npoints/2+1)));
xlabel(['Cavity length in neff=' num2str(neff) ' (um)']);
ylabel('Cavity strength (relative units)');
title('Resonant cavity lengths L=c/2f');
xlim([0, 2000]);
ylim([0, max(reflnorm(100:1024))]);
%ylim([0, 10e-3])

maxtran = max(transnorm(100:1024))
maxrefl = max(reflnorm(20:1024))

a=1;

function zoomx(hand,fraction)
data=get(hand,'ydata');
over=find(data>max(data(2:length(data)))*fraction);
data=get(hand,'xdata');
ha=get(hand,'parent');
xlim=[data(over(1)) data(over(length(over)))];
if xlim(1)>xlim(2),
    xlim=fliplr(xlim);
end
set(ha,'xlim',xlim);