% FORMREGIME: simulating coupled-mode equations for optical gratings under variable boundary condition
% 
% 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
%
% Boundary conditions: 
% u(0,t) = I_in^{1/2}, w(0,t) = 0, v(L,t) = 0, y(L,t) = 0; I_in = I_in(t)
% Initial conditions:
% u(z,0) = 0
% w(z,0) = 0
% v(z,0) = 0
% y(z,0) = 0
%
% the numerical scheme is based on central differences in space z and time t
% the numerical scheme is implicit
% the first time step is performed by using forward difference
% all later time steps are computed by inverting the coefficient matrices

disp('***************************************************');
disp('numerical simulator for coupled-mode equations in 1D');
disp('***************************************************');

L = 20;
%L = 40;
TT = L/2;
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 < 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

I = input('enter intensity of the incident pulse, I_in = ');
%I_in = 'n_0k/cosh(2*n_0k*(t-0.5*L))';
I_in = 'I*exp(-(t-TT)^2)';

u0 = zeros(N+2,1); w0 = zeros(N+2,1); v0 =zeros(N+2,1); y0 = zeros(N+2,1); 
u = zeros(N+2,1); w = zeros(N+2,1); v =zeros(N+2,1); y = zeros(N+2,1); 
figure(1); set(1,'Position',[50 300 400 300]); 
figure(2); set(2,'Position',[500 300 400 300]); 
fprintf('\nDISPLAY: red - steady-state solution, blue - wave amplitude\n\n');
% setting data for figures

t = 0;
iInten = 1;
Iin(1) = eval(I_in);
Iref(1) = 0;
Itrans(1) = 0;
Time(1) = t;
% initialization for reflected and tranmitted intensity data

dt = 0.5*dz;
% the time step is normalized by stepsize unless the nonlinearity is too big 
if (dt * sqrt(I^3) > 0.001)
    dt = 0.001/sqrt(I^3);
end
% this renormalization ensures that the nonlinearity terms remain small in the numerical scheme (order of 0.001)
r = dt/(2*dz);
k = n_0k*dt;
% parameters r and kappa in the numerical scheme
fprintf('time step between two iterations dt = %.3f\n', dt);
fprintf('parameter of the numerical scheme r = %.3f\n\n', r);

UMAX = I; UMIN = I; VMAX = 0; VMIN = 0;
[u_max, u_min] = visual3(u,w,z,t,1,'incident');
[v_max, v_min] = visual3(v,y,z,t,2,'reflected');
UMAX = max(u_max,UMAX); UMIN = min(u_min,UMIN);
VMAX = max(v_max,VMAX); VMIN = min(v_min,VMIN);
% visualization at the initial time instance

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

try 
% the giantic try-catch loop ensures here 
% that the data files in directory "/DATA" contain only complete numerical runs
data1 = strcat('../DATA/incident.',num2str(id));
data2 = strcat('../DATA/reflected.',num2str(id));
data3 = strcat('../DATA/rundata.',num2str(id));
idf1 = fopen(data1,'w');
idf2 = fopen(data2,'w');
idf3 = fopen(data3,'w');
fwrite(idf1,u0,'float32');
fwrite(idf1,w0,'float32');
fwrite(idf2,v0,'float32');
fwrite(idf2,y0,'float32');
fwrite(idf3,L,'float32');
fwrite(idf3,N,'int32');
fwrite(idf3,n_0k,'float32');
fwrite(idf3,n_nl,'float32');
fwrite(idf3,n_2k,'float32');
fwrite(idf3,I,'float32');
data_saved = 0;
% saving initial data and configuration

terminator = 1;
Total = 50;
Inter = 20;
fwrite(idf3,Inter,'int32');
fprintf('ITERATIONS: Press CTRL-C to stop and lose the data or \n            wait until termination in %d iterations\n',Total);
% Inter - number of time instances between two updates of the figure windows
% Total - number of figure window updates before an intermediate termination
% Total * Inter - number of time instance before an intermediate termination

I1 = eval(I_in);
NonU = 2*dt*nonlinearity(u,w,v,y,n_nl,n_2k);
NonY = 2*dt*nonlinearity(y,v,w,u,n_nl,n_2k);
[uu,yy] = step1(u,y,r,k,I1,NonU,NonY);
NonW = 2*dt*nonlinearity(w,u,y,v,n_nl,n_2k);
NonV = 2*dt*nonlinearity(v,y,u,w,n_nl,n_2k);
[ww,vv] = step1(w,v,r,-k,0,-NonW,-NonV);
t = dt;
% obtaining the first time step 

iInten = 2;
Time(2) = t;
Iin(2) = eval(I_in);
Iref(2) = vv(1)^2+yy(1)^2;
Itrans(2) = uu(N+2)^2+ww(N+2)^2;
% data after first step

A1 = coefmatrix(N,r,k);
AA1 = inv(A1);
B1 = coefmatrix(N,-r,-k);
% coefficient matrix for the system of u and y
A2 = coefmatrix(N,r,-k);
AA2 = inv(A2);
B2 = coefmatrix(N,-r,k);
% coefficient matrix for the system of w and v

while terminator > 0
   
   for icircle = 1 : Total, 
      
      for icircleinter = 1 : Inter,
          
            t = t - dt; I1 = eval(I_in);
            t = t + 2*dt; I2 = eval(I_in);
	        NonU = 2*dt*nonlinearity(uu,ww,vv,yy,n_nl,n_2k);
   		    NonY = 2*dt*nonlinearity(yy,vv,ww,uu,n_nl,n_2k);
	  		[uuu,yyy] = timestep(u,y,r,k,[I1 I2],NonU,NonY,AA1,B1);
			NonW = 2*dt*nonlinearity(ww,uu,yy,vv,n_nl,n_2k);
			NonV = 2*dt*nonlinearity(vv,yy,uu,ww,n_nl,n_2k);
			[www,vvv] = timestep(w,v,r,-k,0,-NonW,-NonV,AA2,B2);
			% obtaining a time step 

            u = uu; w = ww; v = vv; y = yy;
   		    uu = uuu; ww = www; vv = vvv; yy = yyy;
      	    % re-assignment of the values
            
            iInten = iInten + 1;
            Time(iInten) = t;
            Iin(iInten) = I2;
            Iref(iInten) = vv(1)^2 + yy(1)^2;
            Itrans(iInten) = uu(N+2)^2 + ww(N+2)^2;
            
        end
        
		fwrite(idf1,uuu,'float32');
   		fwrite(idf1,www,'float32');
   		fwrite(idf2,vvv,'float32');
   		fwrite(idf2,yyy,'float32');
   		data_saved = data_saved + 1;
   		% saving data for the next time step
      
  		[u_max, u_min] = visual3(uu,ww,z,t,1,'incident');
        [v_max, v_min] = visual3(vv,yy,z,t,2,'reflected');
        UMAX = max(u_max,UMAX); UMIN = min(u_min,UMIN);
		VMAX = max(v_max,VMAX); VMIN = min(v_min,VMIN);
   	    % visualization at the current time instance	
   
	end
    
    figure(3); 
    plot(Time,Iin,'y',Time,Iref,'b',Time,Itrans,'r');
    
	terminator = input('1 - to continue or 0 - to terminate, answer = ');
end

fprintf('\nNORMAL TERMINATION: data are written in data files with id = %d\n\n',id);
fwrite(idf3,UMAX,'float32');
fwrite(idf3,UMIN,'float32');
fwrite(idf3,VMAX,'float32');
fwrite(idf3,VMIN,'float32');
fwrite(idf3,data_saved,'int32');

fwrite(idf3,iInten,'int32');
fwrite(idf3,Time,'float32');
fwrite(idf3,Iin,'float32');
fwrite(idf3,Iref,'float32');
fwrite(idf3,Itrans,'float32');
fclose(idf1); fclose(idf2); fclose(idf3);

catch 
    
% all incomplete data files will be removed
fclose('all');
fprintf('\nABNORMAL TERMINATION: all data files with id = %d will be removed',id);
eval(strcat('delete ../DATA/incident.',num2str(id)));
eval(strcat('delete ../DATA/reflected.',num2str(id)));
eval(strcat('delete ../DATA/rundata.',num2str(id)));
close(1); close(2); 

end


