% m-file to simulate non-linear effects due to gain in microrings
% Kuldeep Amarnath , 14 Dec 2005.
function Iout = nonlin

% ************ define constants **************

k = 0.4; t = sqrt(1-k^2); % coupling parameter
Gamma = 0.05; % overlap factor
g0 = 550; % unsaturated material gain (/cm)
alpha = 15; % unsaturable loss (scattering)
phi0 = -0.3; % round-trip phase for unsaturated case (rad)
L  = 2*pi*20e-04; % round-trip length (cm)
b = 3; % linewidth-enhancement factor
lambda = 1.55e-04; % wavelength (cm)
% *******************************************

% set input intensity and plot scaling
%Iin = logspace(-4,0,100);
Iin = linspace(0,.1,100);
plotscale='lin';

% solve for ring intensity on the way up: Ir(1,:)

for i = 1:length(Iin)
    % starting guess
    if i == 1
        Ir0 = 0;
    else 
        Ir0 = Ir(1,i-1); % use the previous value
    end
    % solve for Ir(i)
    [Ir(1,i),FVAL,EXITFLAG] = fzero(@residue,Ir0,[],Iin(i),k,t,Gamma,g0,alpha,phi0,L,b);
end

% solve for ring intensity on the way down: Ir(2,:)

for i = length(Iin):-1:1
    % starting guess
    if i == length(Iin)
        Ir0 = Ir(1,i);
    else 
        Ir0 = Ir(2,i+1); % use the previous value
    end
    % solve for Ir(i)
    [Ir(2,i),FVAL,EXITFLAG] = fzero(@residue,Ir0,[],Iin(i),k,t,Gamma,g0,alpha,phi0,L,b);
end


% calculate all relevant parameters for ring intensity

for i = 1:2
    gm(i,:) = g0./(1+Ir(i,:));
    a(i,:) = exp((Gamma*gm(i,:)-alpha)*L/2);
    dn(i,:) = g0*b*lambda*Gamma/(4*pi) * Ir(i,:)./(1+Ir(i,:));
    dphi(i,:) = g0*L*b*Gamma/2 * Ir(i,:)./(Ir(i,:)+1);
    Iout(i,:) = Iin.*(a(i,:).^2 + t^2 - 2*t*a(i,:).*cos(phi0+dphi(i,:)))./(1 + t^2*a(i,:).^2 - 2*t*a(i,:).*cos(phi0+dphi(i,:)));
end

% plot em all

switch plotscale
    case 'lin'
        % linear plot
        subplot(2,2,1), plot(Iin,Ir(1,:),'b-',Iin,Ir(2,:),'r-'), title('I_{R}') ;
        subplot(2,2,2), plot(Iin,a(1,:),'b-',Iin,a(2,:),'r-'), line([min(Iin) max(Iin)], [1/t 1/t]), title('a') ;
        subplot(2,2,3), plot(Iin,phi0+dphi(1,:),'b-',Iin,phi0+dphi(2,:),'r-'), title('\phi');
        subplot(2,2,4), plot(Iin,Iout(1,:),'b-',Iin,Iout(2,:),'r-'), title('I_{T}');
        
    case 'log'
        % log plot
        subplot(2,2,1), loglog(Iin,Ir(1,:),'b-',Iin,Ir(2,:),'r-'), title('I_{R}') ;
        subplot(2,2,2), semilogx(Iin,a(1,:),'b-',Iin,a(2,:),'r-'), line([min(Iin) max(Iin)], [1/t 1/t]), title('a') ;
        subplot(2,2,3), semilogx(Iin,phi0+dphi(1,:),'b-',Iin,phi0+dphi(2,:),'r-'), title('\Delta \phi');
        subplot(2,2,4), semilogx(Iin,Iout(1,:),'b-',Iin,Iout(2,:),'r-'), title('I_{T}');
end

figure(1);

% residue function calculates [Ir - f(Ir,Iin)]
function val = residue(Ir,Iin,k,t,Gamma,g0,alpha,phi0,L,b)

gm = g0/(1+Ir); % material gain
a = exp((Gamma*gm-alpha)*L/2); % round-trip intensity gain factor
phi = phi0 + g0*L*b*Gamma/2 * Ir/(Ir+1); % net round-trip phase
Irf = Iin * k^2/(1+a^2*t^2-2*a*t*cos(phi)); % expected intensity in ring
val = Ir - Irf; % difference between expected and assumed intensity