summaryrefslogtreecommitdiff
path: root/examples/unstable_cavity.m
blob: cf7afb619151c455533f6bec712f77bd6b6fb2c3 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
% Illustrates the use of cav.m by modeling a marginally stable optical cavity.
%------------------------------------------------------------------------------------
% PROGRAM NAME: unstable_cavity
% AUTHOR: Andri M. Gretarsson
% DATE: Aug. 24, 2004
%
% SYNTAX:  unstable_cavity
% 
% Uses cav.m to model the fields inside and outside a marginally unstable linear cavity 
% The example is based on the LIGO Livingston recycling cavity in its cold state.
% In the 1D mode, the output is a set of 9 graphs showing the 
% the field inside the cavity and in reflection as a function of length around a resonance. 
% Due to the marginally stable nature of the cavity in this example, we get interesing 
% structure in the fields. The program runs as a script rather than a function and 
% numerous parameters can be set near the beginning of the script.  These are: 
%
% R1, R2 - the radii of curvature of the input and end cavity mirrors
% r1, r2 - the amplitude reflection coefficients of the two mirrors
% l1, l2 - the amplitude loss coefficients of the two mirrors
% L0 - the cavity length
% lambda - the optical wavelength
% l_in, m_in - Hermite-Gaussian mode numbers for the input light (referred to 
%              the basis set by qin (see next line).
% qin - the complex radius of curvature of the input light
% npts - for a 1D calculation, the number of points in the calculation, for a 2D
%        calculation the number of grid points is npts^2
% twoD - Whether the calculation should be done in 1D (cross section) or 2D. The 2D
%        form is significantly slower.  To improve running speeds it may be a good
%        idea to reduce npts.  NOTE: It will help to maximize the figure window and 
%        rerun the program to improve visibility in the 1D mode.
%
% INPUT ARGUMENTS: none
% OUTPUT ARGUMENTS: none
% Last Updated: June 30, 2007 by AMG
%
%------------------------------------------------------------------------------------
% SYNTAX:  unstable_cavity
%------------------------------------------------------------------------------------

% Mirror properties
Ritmx=-14.76e3;             % ITMx radius of curvature [meters]
Ritmy=-14.52e3;             % ITMy radius of curvature [meters]
Rrm=20.78e3;                % RM radius of curvature [meters]
R1=Rrm;             
R2=(Ritmx+Ritmy)/2/1.45;    % 1.45 is due to the index of refraction
r1=sqrt(0.9725);            % Measured reflectance coefficients assuming 150ppm loss in ITMs.
r2=sqrt(0.9731);            % see: http://www.ligo-la.caltech.edu/ilog/pub/ilog.cgi?group=detector&date_to_view=08/19/2003&anchor_to_scroll_to=2003:08:25:10:47:28-rana
%r1=sqrt(1-0.027);          % nominal from COC (Helena Armandula) webpage
%r2=sqrt(1-0.0288);
%r1=0.96;                    % test
l1=0;%sqrt(150e-6);            % ITM loss (nominal)
l2=0;%sqrt(150e-6);            % ITM loss (nominal)
L0=9.2;

% Incident light
lambda=1.064e-6;
l_in=0;
m_in=0;
[qin,pin]=prop(q_(0.03875,7.1e3,lambda),free(L0),[l_in,m_in]); % Input beam is matched to arms so use arm cavity beam at ITMx to get q value.
%[qin,pin]=prop(q_(0.03685,-6.9564e3,lambda),free(L0),[l_in,m_in]);
win=w_(qin,lambda);

% Plot domain
twoD=0;
if twoD
    npts=9; 
    xseed=[-2*win:win/npts:2*win];    
    yseed=xseed';
    [x,y]=meshgrid(xseed,yseed);    
    thecolormap='bone';
else 
    npts=200; 
    xseed=[-2.5*win:win/npts:2.5*win];
    x=xseed;
    y=zeros(size(x));    
end


zin=HermiteGaussianE([l_in,m_in,qin,lambda(1),pin],x,y);


n_bounces=70;
n_lengths=9;
lowerdL=-15e-9; upperdL=15e-9;
deltaL=[lowerdL:(upperdL-lowerdL)/(n_lengths-1):upperdL];
Pcoherent_out=zeros(n_lengths,1);
for s=1:n_lengths
    L=round(L0/lambda)*lambda+deltaL(s);  % Cavity length

    [qrefl,qtrans,qcav,prefl,ptrans,pcav]=...
        cav(qin,l_in,m_in,lambda,L,R1,R2,r1,r2,l1,l2,n_bounces,[1.46,1,1,],[0,0]);
    
    zrefl=HermiteGaussianE([l_in * ones(size(qrefl)),m_in * ones(size(qrefl)),...
            qrefl,lambda(1) * ones(size(qrefl)),prefl],x,y);
    zcav=HermiteGaussianE([l_in * ones(size(qcav)),m_in * ones(size(qcav)),...
            qcav,lambda(1) * ones(size(qcav)),pcav],x,y);
    ztrans=HermiteGaussianE([l_in * ones(size(qtrans)),m_in * ones(size(qtrans)),...
            qtrans,lambda(1) * ones(size(qtrans)),ptrans],x,y);    
    if length(size(zcav))==3
        zsumrefl=sum(zrefl,3);
        zsumcav=sum(zcav,3);
        zsumtrans=sum(ztrans,3);        
    else
        zsumrefl=sum(zrefl,2);
        zsumcav=sum(zcav,2);
        zsumtrans=sum(ztrans,2);
    end
    
    if ~twoD
        interpdom=[min(xseed):(max(xseed)-min(xseed))/999:max(xseed)];
        dx_interp=interpdom(2)-interpdom(1);
        if s==1
            inputintens_interp=interp1(x,abs(zin).^2,interpdom,'spline');
            inputpower=sum(inputintens_interp*dx_interp.*abs(interpdom)*pi);
        end
        cavintens_interp=interp1(x,abs(zsumcav).^2,interpdom,'spline');
        cavpower=sum(cavintens_interp*dx_interp.*abs(interpdom)*pi);
        reflintens_interp=interp1(x,abs(zsumrefl).^2,interpdom,'spline');
        reflpower=sum(reflintens_interp*dx_interp.*abs(interpdom)*pi);
        transintens_interp=interp1(x,abs(zsumtrans).^2,interpdom,'spline');
        transpower=sum(transintens_interp*dx_interp.*abs(interpdom)*pi);
        Pcoherent_out(s)=reflpower+transpower;
    end

    if twoD
        figure(1);
        subplot(3,3,s); h=pcolor(x,y,abs(zsumcav));  % h=pcolor(x,y,abs(zsumcav).*cos(angle(zsumcav)));
        set(h,'EdgeColor','none'); axis square;
        colormap(thecolormap); %colorbar;  
        set(h,'EdgeColor','none'); set(h,'FaceColor','interp');
        set(gca,'Visible','off'); set(gcf,'Color','black');        
        set(gca,'Zlim',[-35 35]);
        drawnow; 
        h=figtext3D(3,1,10,['\Delta{L}: ',num2str(deltaL(s)*1e9,'%0.3g'),' nm'],10);
        set(h,'color','white');
        shg;
    else    
        figure(1); subplot(3,3,s);  
        set(gcf,'Color','default');
        plot(x*100,zin.*conj(zin),'m--','linewidth',2); hold on;
        h=plot(x*100,zsumrefl.*conj(zsumrefl),'-',x*100,zsumcav.*conj(zsumcav),'-');
        set(h,'linewidth',2);
        legend('input','refl','cav');
        drawnow; hold off;     
        xlabel('cm'); ylabel('intensity');
        figtext(0.5,9,['\Delta{L}: ',num2str(deltaL(s)*1e9,'%0.3g'),' nm'],8);
        figtext(0.5,7.5,['P_{cav}/P_{in} = ',num2str(cavpower,'%0.2g')],8);
        shg;
    end

end