%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% SPLIT-STEP FOURIER TRANSFORM METHOD FOU NLSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%輸入光場(chǎng)部分調(diào)制%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             

%%%%%%%%%%參數(shù)選擇lamda=1.054e-4(cm);n0=1.5;n2=2e-7(cm^2/GW);I0=5(GW/cm^2);%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%k0=2*pi*n0/lamda;g=k0*n2/2/n0*I0=0.0298(/cm);%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  B=gz=0.0298*40=1.192%%%%%%%%%%%%%%%  gama=0.0861cm,qx=72.9754%%%%%%%%%%%%

function Transform_NLSE_SPFFT
%%%%%%%%%%%%%%%%%%%%%%
N=512;
%%%%%%%%%%%%%%
DELTA_X=0.001; 
%%%%%%%%%%%%%%
DELTA_Z=0.01;
%%%%%%%%%%%%%%
N_STEP=8000;
%%%%%%%%%%%%%%%%%%%%%%
U=INPUT_FUN;
INTENSITY_OUT_INI(U) ;%%%%%%%output the intensity of the initial%%%%%%%%

LINEAR_FACTOR=LINEAR_FUN;%%%%%%%%%%線性傳輸因子%%%%%%%%%%%%%%%%%%%%%%%%
% CHECK(U);
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
q=0;
%%%%%%%%%%%%%%%%%%%%%%%
while(1)
      
    if mod(q,20)==0
        
        P=POWER_OUT(U);
        fprintf('\n  %d of %d',q,N_STEP);
        fprintf('\n P= %f',P);

    end
     
    U_TEP1=U;

    U_FT=U; 
    U_FT=LINEAR_FT_FUN(LINEAR_FACTOR,U_FT);
    for k=1:2
        U_TEP2=NONLINEAR_FUN(U_FT,U,U_TEP1);
        U_TEP1=LINEAR_FT_FUN(LINEAR_FACTOR,U_TEP2);
    end
    
    U=U_TEP1;
    
    if mod(q,1000)==0
     FILE_ID=q/1000;

     INTENSITY_OUT(U,FILE_ID);
   end
    
    if(q>N_STEP)
        break;
    end
    
    q=q+1;
end

%%%%%%%%%%%輸出最后強(qiáng)度值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid1=fopen('ENDINTENSITY_real.dat','w');
fid2=fopen('ENDINTENSITY_imag.dat','w');
for i=1:N
  fprintf(fid1,'%10.6f',real(U(i,:)));
  fprintf(fid1,'\n');
  fprintf(fid2,'%10.6f',imag(U(i,:)));
  fprintf(fid2,'\n');
end
fclose(fid1);
fclose(fid2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%the end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% CHECK(U);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% SUB FUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% SUB FUNCTION 1-----INPUT DATA %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=INPUT_FUN

%%%%%%%%%%
N=512;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fid=fopen('RE_U.DAT','r');
% RE_U=fscanf(fid,'%f',[400,inf]);
% fclose(fid);
% 
% fid=fopen('IM_U.DAT','r');
% IM_U=fscanf(fid,'%f',[400,inf]);
% fclose(fid);
% 
% U=complex(RE_U,IM_U);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=0.1;
DELTA_X=0.001; 
X=-N/2*DELTA_X:DELTA_X:N/2*DELTA_X;
Y=X;
qx=72.9754;
for i=1:N
    for j=1:N
        if abs(X(i))<=0.0861/2&&abs(Y(j))<=0.0861/2
        U(i,j)=1+a*abs(sin(qx*X(i))*sin(qx*Y(j)));
       else 
        U(i,j)=1;
       end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%INTENSITY_OUT_INITIAL%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function INTENSITY_OUT_INI(U)
fid=fopen('INIINTENSITY.DAT','W')
N=512
for i=1:N
fprintf(fid,'%9.6f ',abs(U(i,:).^2));
fprintf(fid,'\n');
end
fclose(fid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% LINEAR MULTIPLIER FUNCTION %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function LINEAR_FACTOR=LINEAR_FUN

%%%%%%%%%%%%%%
MAX=0.1;
%%%%%%%%%%%%%%
MAX_PS=0;
lamda=1.054e-4;
%%%%%%%%%%%%%%
PI=3.141592654;
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%
N=512;
%%%%%%%%%%%%%
DELTA_X=0.001;
%%%%%%%%%%%%%%
DELTA_Z=0.01;
%%%%%%%%%%%%%%
ii=sqrt(-1);
%%%%%%%%%%%%%
for i=1:N
    if (i<=N/2+1)
       FX(i)=((i-1)/(N*DELTA_X))^2;
   else
       FX(i)=((i-1-N)/(N*DELTA_X))^2;
   end
end
FACTOR=-lamda/2*DELTA_Z;
% FACTOR=-2*PI^2*DELTA_Z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MAX_PS=abs(FACTOR)*(FX(N/2)+FX(N/2));

if MAX_PS>MAX
    disp('The chosen DELTA_Z may be too large');
end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:N
    for j=1:N
        LINEAR_FACTOR(i,j)=exp(ii*FACTOR*(FX(i)+FX(j)));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% OUTPUT POWER FUNCTION %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function P=POWER_OUT(U)

%%%%%%%%%%%%%%%%
DELTA_X=0.001;
%%%%%%%%%%%%%%%%
N=512;
%%%%%%%%%%%%%%%%

I=abs(U).^2;

P=sum(sum(DELTA_X^2*I));

% P1=0;
% P2=0;
% P3=0;
% for i=2:N
%     for j=2:N
%         P1=P1+DELTA_X^2*(I1(i,j)+I1(i-1,j)+I1(i,j-1)+I1(i-1,j-1))/4;
%         P2=P2+DELTA_X^2*(I2(i,j)+I2(i-1,j)+I2(i,j-1)+I2(i-1,j-1))/4;
%         P3=P3+DELTA_X^2*(I3(i,j)+I3(i-1,j)+I3(i,j-1)+I3(i-1,j-1))/4;
%     end
% end
% P=P1+P2+P3;

fid1=fopen('POWER.DAT','w');
fprintf(fid1,'P=%10.6f\n',P);
fclose(fid1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% OUTPUT INTENSITY FUNCTION %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function INTENSITY_OUT(U,FILE_ID)

%%%%%%%%%%%%%%
N=512;
%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if FILE_ID<9 
    FILE_NAME1=['00' char(FILE_ID+48) char(49) 'I' '.' 'DAT']; 
else
    FILE_NAME1=[char(floor(FILE_ID/10)+48)  char(mod(FILE_ID,10)+48) char(49) 'I' '.' 'DAT'];   
end

fid2=fopen(FILE_NAME1,'w');

I=abs(U).^2;

for i=1:N
    fprintf(fid2,'%10.6f',I(i,:));
    fprintf(fid2,'\n');    
end

 fclose(fid2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  I=abs(U).^2;
%  
%  fid3=fopen('I.DAT','a+');
%  
%  for i=1:N
%      fprintf(fid3,'%15.9f',I(i,:));
%      fprintf(fid3,'\n');
%  end
%  
%  fclose(fid3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% THE LINEAR PART’S FOURIER TRANSFORM %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U_FT=LINEAR_FT_FUN(LINEAR_FACTOR,U_FT)

%%%%%%%%%
%%%%%%%%%

U_FT=fft2(U_FT);

U_FT=U_FT.*LINEAR_FACTOR;

U_FT=ifft2(U_FT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% NONLINEAR MULTIPLIER FUNCTION %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U_TEP1=NONLINEAR_FUN(U_FT,U,U_TEP1)

%%%%%%%%%%%%%
N=512;
%%%%%%%%%%%%%
DELTA_Z=0.01;
%%%%%%%%%%%%%
g=0.0298;
I=abs(U).^2;
I1=abs(U_TEP1);

ii=sqrt(-1); 
%  U_TEP=U_TEP+F_FUN(U_TEP2);
% 
% U_TEP=0.5*U_TEP;

U_TEP1=U_FT.*exp(ii*DELTA_Z*g*0.5.*(I+I1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% NONLINEAR TERM IN NLSE %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function f_fun=F_FUN(U)
% 
% % %%%%%%%%%
% % N=400;
% % %%%%%%%%%
% % s=0.5;
% % %%%%%%%%%
% % 
%  I=abs(U).^2;
% % 
%  g=0.0298;
% % 
% f_fun=g.*I;
% % for i=1:N
% %     for j=1:N
% %        f1fun(i,j)=(abs(U(i,j))^2+abs(U2(i,j))^2+abs(U3(i,j))^2)/(1+s*(abs(U(i,j))^2+abs(U2(i,j))^2+abs(U3(i,j))^2));
% %    end
% % end