% LINEARSPECTRA: spectrum of the linearized coupled-mode equations for optical gratings
% 
% The photonic grating is located between z = 0 and z = L
% i( A^+_t + 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^-_t - 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
%
% Linearization:
% u(z,t) = u_st(z) + u_lin(z) exp(lambda t)
% w(z,t) = w_st(z) + w_lin(z) exp(lambda t)
% v(z,t) = v_st(z) + v_lin(z) exp(lambda t)
% y(z,t) = y_st(z) + y_lin(z) exp(lambda t)
% where u_st(z) - stationary solution, u_lin(z) - eigenfunction, lambda - eigenvalue
% Boundary conditions: 
% u_lin(0) = 0, w_lin(0) = 0, v_lin(L) = 0, y_lin(L) = 0
%
% the numerical scheme is based on central differences in space z 

disp('***************************************************************************');
disp('numerical calculator for spectra of linearized coupled-mode equations in 1D');
disp('***************************************************************************');

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

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 < 10
	fprintf('WARNING: results may be inaccurate, not too many points\n');
elseif N > 500 
	fprintf('WARNING: program may require long running time, too many points\n');
end
N2 = 2*N+1;
% a half of the step size (used for Richardson extrapolation)

K = input('enter number of points for the transmitted intensity, K = ');
I_start = input('enter the starting value for transmitted intensity, I_start = ');
I_fin = input('enter the final value for transmitted intensity, I_fin = ');
if I_fin < I_start
   fprintf('ERROR: the starting value of I_out exceeds the final value of I_out \n');
   break;
end
dI = (I_fin-I_start)/K;
% step size in the mesh for transmitted intensity I_out

figure(1); figure(2); 
set(1,'Position',[50 300 400 300]); set(2,'Position',[500 300 400 300]); 
% positioning the figure windows

iD = 0; % index for unstable eigenvalues
Ipair = zeros(2,K); % pair of values for (I_in,I_out)

numberofUnstable = 10;
Leig = zeros(numberofUnstable,K); % we reserve space up to ten unstable eigenvalues for all possible values of I_out

for iK = 1 : K
    I_out = I_start + dI*iK;

    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

    if n_nl==0
   
        if (n_2k*I_out >= abs(n_0k))
            limpar = n_2k*I_out*cos(2*sqrt(n_2k^2*I_out^2-n_0k^2)*L)-n_0k;
        elseif (n_0k > 0)
            limpar = n_0k - n_2k*I_out*cosh(2*sqrt(n_0k^2-n_2k^2*I_out^2)*L);
        else
            limpar = n_2k*I_out*cosh(2*sqrt(n_0k^2-n_2k^2*I_out^2)*L) - n_0k;
        end
    
        if limpar <= 0
	        fprintf('ERROR: the output intensity exceed the limiting intensity!\n');
	        break;   
        end
        [u0,w0,v0,y0,I] = steadyexact(z,I_out,n_0k,n_2k);
    else
        [u0,w0,v0,y0,I,key] = steadysol(z,I_out,n_0k,n_2k,n_nl);
        if key == 0
            fprintf('ERROR: the output intensity exceed the limiting intensity!\n');
	        break; 
        end
    end

    fprintf('incident intensity at the left, I_in = %.4f\n', I);
    fprintf('transmitted intensity at the right, I_out = %.4f\n',I_out);
    visualLin(u0,w0,z,0,1,'incident');
    visualLin(v0,y0,z,0,2,'reflected');
    pause(1);
    % steady-state configuration corresponding to the incident intensity I = I_in

    k = 2*n_0k*dz;
    % parameter kappa in the numerical scheme

    [C] = linearizedmatrix(u0,w0,v0,y0,N,k,n_nl,n_2k,dz);
    % computing the linearized matrix

    [V,D] = eig(C);
    D = D/(2*dz);
    Dv = diag(D);
    % computing proper eigenvalues

    [DR,IR] = sort(real(Dv));
    i1 = 0;
    toler = 0.25*abs(max(imag(Dv)));
    % sorting eigenvalues and designing a tolerance interval to eliminate large imaginary "arti" eigenvalues

    for i2 = length(IR) : -1 : 1
        if DR(i2) <= 0
            break;
        elseif abs(imag(D(IR(i2),IR(i2)))) < toler 
            i1 = i1 + 1;
            L1(i1) = D(IR(i2),IR(i2));
        end
    end

    figure(1); 
    plot(real(Dv),imag(Dv),'*r'); pause(1);
    % plotting the real versus imaginary part for eigenvalues

    if i1 ~= 0
      
        iD = iD + 1;
        Ipair(1,iD) = I;
        Ipair(2,iD) = I_out;
        % unstable eigenvalues exist for the given value of (I,I_out)
   
        if i1 > numberofUnstable
                i1 = numberofUnstable;
        end
    
        for ii1 = 1 : i1
                Leig(ii1,iD) = L1(ii1);
        end
        
    else
        fprintf('No unstable eigenvalues for I_out = %.4f\n\n',I_out);
        % unstable eigenvalues do not exist!
    end
end

figure(1);
plot(Ipair(2,:),real(Leig(1,:)),'*b',Ipair(2,:),real(Leig(2,:)),'*g');
figure(2);
plot(Ipair(2,:),imag(Leig(1,:)),'*b',Ipair(2,:),imag(Leig(2,:)),'*g');

id = 0; i = 0;
while (id ~= -1) 
     i = i + 1;
     data1 = strcat('../DATA/SpecData.',num2str(i));
     id = fopen(data1);
end
id = i; fclose('all'); 
% dynamical search of a number to write the data files

data1 = strcat('../DATA/SpecData.',num2str(id));
data2 = strcat('../DATA/SpecEigen.',num2str(id));

idf1 = fopen(data1,'w');
idf2 = fopen(data2,'w');

fwrite(idf1,L,'float32');
fwrite(idf1,n_0k,'float32');
fwrite(idf1,n_nl,'float32');
fwrite(idf1,n_2k,'float32');
fwrite(idf1,iD,'int32');
fwrite(idf1,Ipair(1,1:iD),'float32');
fwrite(idf1,Ipair(2,1:iD),'float32');
fwrite(idf2,real(Leig(:,1:iD)),'float32');
fwrite(idf2,imag(Leig(:,1:iD)),'float32');

fprintf('\nNORMAL TERMINATION: data are written in data files with id = %d\n\n',id);
fclose(idf1); fclose(idf2); 
% saving data

