% COUPLEMODE: simulating 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
%
% Boundary conditions: 
% u(0,t) = I_in^{1/2}, w(0,t) = 0, v(L,t) = 0, y(L,t) = 0
% Initial conditions:
% u(z,0) = u_st(z) + u_pert(z)
% w(z,0) = w_st(z) + w_pert(z)
% v(z,0) = v_st(z) + v_pert(z)
% y(z,0) = y_st(z) + y_pert(z)
%
% 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;
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_out = input('enter output intensity of the incident wave, I_out = ');
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
% steady-state configuration corresponding to incident intensity I = I_in

fprintf('incident intensity at the left, I_in = %.3f\n\n', I);

pert = input('enter a perturbation factor to the reflected wave, pert = ');
if (pert == 0)
   fprintf('WARNING: the initial condition is unperturbed steady-state regime!\n');
end

u = u0; w = w0; v = v0 + pert*v0.*exp(-z); y = y0 + pert*y0.*exp(-z); 
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;
Iref(1) = v(1)^2+y(1)^2;
Itrans(1) = u(N+2)^2+w(N+2)^2;
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);

t = 0;
UMAX = I; UMIN = I; VMAX = 0; VMIN = 0;
[u_max, u_min] = visual2(u,w,u0,w0,z,t,1,'incident');
[v_max, v_min] = visual2(v,y,v0,y0,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 = 10;
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

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,I,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;
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,
          
	        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,I,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 
		    
            t = t + dt;
            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;
            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] = visual2(uu,ww,u0,w0,z,t,1,'incident');
        [v_max, v_min] = visual2(vv,yy,v0,y0,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,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,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


