% 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