% STATREGIME: Matlab script for computing the transmittance curve I_out as function of I_in 
% for the given optical grating
% according to the numerical solution of steady-state equations
%
% z - column vector of (N+2) mesh points in z between z=0 and z=L
% I_in - row vector of intensity of the incident wave at the left
% I_out - row vector of output intensity of the incident wave at the right
% n_0k, n_2k, n_nl - parameters of the optical grating

% Steady-state equations:
% i A^+_z + n_0k A^- + n_nl ( |A^+|^2 + 2 |A^-|^2 ) A^+ 
%			+ n_2k [ ( 2 |A^+|^2 + |A^-|^2 ) A^- + A^+^2 bar{A^-} ] = 0
% - i A^-_z + n_0k A^+ + n_nl ( 2 |A^+|^2 + |A^-|^2 ) A^- 
%			+ n_2k [ ( |A^+|^2 + 2 |A^-|^2 ) A^+ + A^-^2 bar{A^+} ] = 0
% A^+ = u + i w
% A^- = v + i y

disp('********************************************');
disp('computing transmittance for optical gratings');
disp('********************************************');

L = 20;
fprintf('\nthe length of the optical grating is: L = %d', L);

% n_2k = 1;
% fprintf('\nvariance of the nonlinearity index n_2k = %d\n', n_2k);

n_2k = input('\nenter variance of the Kerr nonlinearity, n_2k = ');
n_0k = input('enter variance of the linear grating index, n_0k = ');
n_nl = input('enter average nonlinearity index, n_nl = ');

if (n_2k > 0.75*abs(n_nl))
    if ( n_0k*abs(n_nl) > - pi * n_2k/(3*L) ) 
        disp('the optical grating is in LIMITING regime');
    else
        disp('the optical grating is in locally BISTABLE LIMITING regime');
    end
else
    disp('the optical grating is in MULTISTABLE regime');
end


N = input('\nenter number of mesh points for the optical gratings, N = ');
if N < 100
	fprintf('WARNING: results may be inaccurate, too few points\n');
elseif N > 10000 
	fprintf('WARNING: program may require long running time, too many points\n');
end

dz = L/(N + 1);
% the space step between two mesh points
z = 0 : dz : L; z = z';
% the vector for the mesh points between z = 0 and z = L

answer = 1; close all;

while (answer == 1)

    clear u; clear w; clear v, clear y; clear uu; clear vv; clear ww; clear yy; clear I_in; clear I_out;
    K = input('enter number of points for output intensity, K = ');
    I_start = input('enter the starting value for transmitted intensity, I_start = ');
    I_fin = input('enter a final value for output intensity, I_fin = ');
    dI = (I_fin-I_start)/K;
    I_out = (I_start + dI) : dI : I_fin;

    u(N+2,:) = sqrt(I_out);
    w(N+2,:) = zeros(1,K);
    v(N+2,:) = zeros(1,K);
    y(N+2,:) = zeros(1,K);
    % boundary conditions at the right boundary (output intensity)

    for n = (N+1):-1:1
        % matrix computation for the whole range of output intensities
        uu = u(n+1,:); ww = w(n+1,:); vv = v(n+1,:); yy = y(n+1,:);
    
        % steady-state equations:
        u(n,:) = uu + dz*( n_0k*yy + n_nl*(uu.^2+ww.^2+2*vv.^2+2*yy.^2).*ww + n_2k*((uu.^2+3*ww.^2+vv.^2+yy.^2).*yy+2*uu.*ww.*vv));
        w(n,:) = ww - dz*( n_0k*vv + n_nl*(uu.^2+ww.^2+2*vv.^2+2*yy.^2).*uu + n_2k*((3*uu.^2+ww.^2+vv.^2+yy.^2).*vv+2*uu.*ww.*yy));
        v(n,:) = vv - dz*( n_0k*ww + n_nl*(2*uu.^2+2*ww.^2+vv.^2+yy.^2).*yy + n_2k*((uu.^2+ww.^2+vv.^2+3*yy.^2).*ww+2*uu.*vv.*yy));
        y(n,:) = yy + dz*( n_0k*uu + n_nl*(2*uu.^2+2*ww.^2+vv.^2+yy.^2).*vv + n_2k*((uu.^2+ww.^2+3*vv.^2+yy.^2).*uu+2*ww.*vv.*yy));
    
    end

    I_in = u(1,:).^2 + w(1,:).^2;
    
    plot(I_in,I_out,'m-');
    figure(1);

    answer = input('Enter 1 to try again or 0 to terminate, answer = ');

end













      
      
      
      
      
      