simple soliton simulation ( N=3 soliton ) 
( T=10ps,  20Gb/s,  400/50km )
 1995.10.6      by K.Shimoura

1997.8.15     last  revised

(*-----Initial Conditions-----------------------------*)
(* LD parameters *)
wl =1550   ;tfw=10     ;Pm =39.42  ;GHz=20     ;
ch =-0.0   ;dnm=2.0    ; 

(* Fiber parameters *)
s  =50.0   ;aa =0.00   ;Sl =0.00   ;
Dp =0.2    ;dDp=0.0    ;kd =5      ;
n2 =2.24 10^-20        ;Rt =0.0    ;

(* Amplification parameters *)
ka =25     ;NF =-100   ;fw =100    ;fm =1      ;
 
(* Calculation parameters *)
h  =1000   ;dt =0.5    ;NM =1      ; 
kmax=400   ;kw =100    ;kgt=5      ;ym =1.2    ;
dm =1000/GHz/dt        ;M =2 Quotient[NM dm,2] ;

(* Graphics Array *)
Array[ nse1,kmax/kgt];
Print[Date[ ]];

(*----Pulse shape definition--------------------------*)

(* Hyperbolic-secant *)  
kt =1.763 ; t0= tfw /kt ;  

EE0=Table[ Sech[t dt/t0] ,{t,-M/2,M/2-1} ];
Chp=Table[ Exp[-((ch/2)(t dt/t0)^2 )I],{t,-M/2,M/2-1}];
E0 =EE0 * Chp;


(* WDM of Hyperbolic-secant  
kt =1.763 ; t0= tfw /kt ; 
dth=(2 Pi dt dnm)/8.04  ;
dn =dm/NM ;

EE0=Table[   Sech[t dt/t0] 
           + Sech[(t-1dn) dt/t0] * Exp[-1t dth I ] 
           + Sech[(t-2dn) dt/t0] * Exp[-2t dth I ] 
           + Sech[(t+1dn) dt/t0] * Exp[+1t dth I ] 
           + Sech[(t+2dn) dt/t0] * Exp[+2t dth I ] 
                                 ,{t,-M/2,M/2-1} ];

SLOT=Table[ If[Abs[t dt]<=dm/2,1,1 ],{t,-M/2,M/2-1} ];
E0 =EE0;
*)

(* Super Gausian 
kt =1.665 ; t0= tfw /kt ;  m= 1.436 ;
EE0=Table[ Exp[ -0.5(t dt/t0)^(2m)],{t,-M/2,M/2-1} ];
Chp=Table[Exp[-((ch/2)(t dt/t0)^(2m))I],{t,-M/2,M/2-1}];
E0 =EE0 * Chp; 
*)

(* Dark soliton 
kt =1.763 ; t0= tfw /kt ; 
EE0=Table[ Tanh[t dt/t0] * 
                Exp[-(t dt/(4t0) )^4] ,{t,-M/2,M/2-1}];
Chp=Table[ Exp[-((ch/2)(t dt/t0)^2 )I],{t,-M/2,M/2-1}];
E0 =EE0 * Chp; 
*)

(* Dispersion length -----------------------*)
Lt = kmax h/1000;
La = ka h/1000;
att= 10^( -aa h /40000);
Ld = t0^2 /Abs[ 5.31 10^-7 wl^2 Dp ];
Lq = t0^3 /Abs[ 2.82 10^-13 wl^4 Sl ];
Ls = 1.571 * Ld;
Print["Dp, dDp      = ",Dp," +- ",dDp,"( ps/nm/km )"];
Print["Ld, Lq       = ",Ld,"( km )  ",Lq ,"( km ) "];
Print["Lt, Lt/Ls    = ",Lt,"( km )  ",Lt/Ls];
Print["La, La/Ld    = ",La,"( km )  ",La/Ld];
Print[" "];
      
(* Initial width & power -------------------*)
uu = Table[ t ,{t,-M/2,M/2-1} ];
EET= Interpolation[ Re[ E0 ] ];
    Plot[ EET[t],{t,1,M},PlotRange-> {-ym,ym},
                  AxesLabel ->{"time","E-field"}];
AE0= (Abs[E0])^2;
EEU= Interpolation[ AE0 ];
TEU= Interpolation[ uu * AE0 ];
TEV= Interpolation[ uu^2 * AE0 ];

IEU= NIntegrate[ EEU[t],{t,1,M},AccuracyGoal->3];
ITU= NIntegrate[ TEU[t],{t,1,M},AccuracyGoal->3];
ITV= NIntegrate[ TEV[t],{t,1,M},AccuracyGoal->3];
Pow= IEU Pm dt GHz /(1000 NM);
ktE= tfw / (Sqrt[ (ITV/IEU) - (ITU/IEU)^2 ] * dt);

Print["W(FWHM,rms)   = ",tfw,"( ps )  ",tfw/ktE ];
Print["Peak Power    = ",Pm,"( mW ) ",
                     N[ 10 Log[10,Pm]],"( dBm )"];
Print["Average Power = ",Pow,"( mW ) ",
                        10 Log[10,Pow],"( dBm )"];

(* Dynamic soliton parameters --------------*)
P0 = 2.63 10^-28 wl^3 s Dp/(n2 tfw^2);
Ga = 10^(aa ka h/10000);
gam= Ga Log[Ga]/(Ga-1);
Pam= P0 * gam;
Ne = Sqrt[ Pm /P0 ];
Nt = Sqrt[ gam ];
Print["P0, Pa, gam   = ",P0,"( mW )  ",Pam,"( mW )  "
      ,gam];
Print["N-exp,N-th    = ",Ne,"   ",Nt ];
Print[" "];

(* Spectrum parameters ---------------------*)
F0 = InverseFourier[ E0 ] ;
tt = Table[ If[ t <= M/2,(t+M/2),(t-M/2)],{t,1,M}];
AF0= Part[ Abs[F0],tt ]^2 ; 
FFU= Interpolation[ 10 Log[10,AF0] ];
ff = Table[ If[ t <= M/2,(t-1),(t-M-1) ],{t,1,M} ];
fr = ff /(M dt);

(* Dispersion slope ------------------------*)
dl = 3.34 10^-6 ( wl^2 /(M dt) );
Dpl= Table[ Dp +Sl(t dl),{t,-M/2,M/2-1} ];
ltf= Table[ -t+M+1,{t,1,M} ];    
Ddf= Part[ Dpl,ltf ];
DDT= Interpolation[Ddf];

(* Initial shape drawing--------------------*)
Plot[ DDT[t],{t,1,M},PlotRange-> {-2,+2},
               AxesLabel ->{"freq","D( ps/nm/km )"} ];
Plot[ EEU[t],{t,1,M},PlotRange-> {0,ym},
               AxesLabel ->{"time","Intensity"}];
Plot[ FFU[t],{t,1,M},PlotRange-> {-30,20},
               Frame -> True,GridLines -> Automatic];


(*---- Pulse propagation by S.S.F.M.------------------*)
Ei = E0;
Do[   
If[ Mod[k,kgt] == 1,
   nse1[(k-1)/kgt+1]= Abs[ Ei^2 ] ];
   
   Dpx= Dp + If[ EvenQ[ Quotient[k,kd]],+dDp,-dDp];
   gm = -0.0252 Dpx fr^2 +0.0673 Sl fr^3; 
   A  = att Exp[ -(gm h/2)I ];
   B  = -( 6.283*10^18 (n2/wl) (Pm/s) )I;

   F0 = InverseFourier[Ei];
   E1 = Fourier[A * F0];

(* Raman *)  
   AE1= (Abs[E1])^2;
   DE1= Append[ Rest[AE1],First[AE1] ];
   RAM= ( DE1 - AE1 ) * 0.001 Rt/dt ;
   B1 = Exp[ h B ( AE1 - RAM ) ];
   E2 = B1 * E1;
   F2 = InverseFourier[E2];
(*   
   B1 = Exp[ h B (Abs[E1])^2 ];
   E2 = B1 * E1;
   F2 = InverseFourier[E2];     
*)
   F3 = A * F2;
   E3 = Fourier[F3];

(*------------Amplification and Filtering---------- *)
If[ Mod[k,ka] == 0,
  
(* Filter *)
   ktf= 0.6931;
   fwt= fw / dl; 
   FIL= Table[Exp[ -ktf Abs[2f/fwt]^(2fm) ]
                                      ,{f,-M/2,M/2-1}];
(* Noise *)    
   ASE= 0.1988 * 10^(NF/10) /(wl M dt); 
   WhN= Table[ ( ASE / Pm )
         (* 2 Random[] * Exp[ 6.28 Random[] I ]*),{M}];
   F4 = ( F3 + Sqrt[WhN] ) * Sqrt[ Part[FIL,tt] ];
   E4 = Fourier[ F4 ];
   GE4= Interpolation[ (Abs[E4])^2 ];
   IE4= NIntegrate[ GE4[t],{t,1,M},AccuracyGoal->3];
   
   AMP= IEU/IE4;
(* AMP= Ga; *)
(* Print["AMP  = ",10 Log[10,AMP]," ( dB )"]; *)

   F5 = F4 * Sqrt[ AMP ];
   E3 = Fourier[ F5 ]    
   ];   

(*----Pulse Width----------------*)
If[ Mod[k,kw] == 0,      
   L  = k h /1000;     
      
   AE3 = (Abs[E3])^2;
   Pmax= Max[ AE3 ];
   GE3 = Interpolation[ AE3 ];
   IE3 = NIntegrate[ GE3[t],{t,1,M},AccuracyGoal->3];
   Po3 = IE3 Pm dt GHz /(1000 NM);

(* FWHM *)   
Do[
   x1 = x;
   S1 = GE3[x1]  -Pmax/2;
   S2 = GE3[x1+1]-Pmax/2;
   If[ S1*S2<0 && S1<0, Break[] ]
,{x,1,M-1}];
   x1 = x1 -S1/(-S1+S2);
Do[
   x2 = M-x;
   S1 = GE3[x2]  -Pmax/2;
   S2 = GE3[x2+1]-Pmax/2;
   If[ S1*S2<0 && S1>0, Break[] ]
,{x,1,M-1}];
   x2 = x2 +S1/( S1-S2);
      
   Wid1= If[x2>x1, (x2-x1)*dt, 999];
   Wid2= (M+1-x1+x2) *dt;
   Wfwh= Min[ Wid1,Wid2 ];
   Wfwh= If[ Wfwh>dt, Wfwh, dt ];

(* rms *)
   TE1 = Interpolation[ uu * AE3 ];
   TE2 = Interpolation[ uu^2 * AE3 ];
   IT1 = NIntegrate[ TE1[t],{t,1,M},AccuracyGoal->3];
   IT2 = NIntegrate[ TE2[t],{t,1,M},AccuracyGoal->3];
   Wrms= Sqrt[ (IT2/IE3) - (IT1/IE3)^2 ] *ktE *dt;
   Wrms= If[ Wrms>dt, Wrms, dt ];
      
(*----Spectrum Width----------------*)  
   AF3 = Part[ Abs[F3],tt ]^2;
   GF3 = Interpolation[ AF3 ];
   LF3 = Interpolation[ 10 Log[10,AF3] ];
   Fma= Max[ AF3 ];
(*
(* FWHM *)
Do[
   y1 = y;
   T1 = GF3[y1]  -Fma/2;
   T2 = GF3[y1+1]-Fma/2;
   If[ T1*T2<0 && T1<0, Break[] ]
,{y,1,M-1}];
   y1 = y1 -T1/(-T1+T2);
Do[
   y2 = M-y;
   T1 = GF3[y2]  -Fma/2;
   T2 = GF3[y2+1]-Fma/2;
   If[ T1*T2<0 && T1>0, Break[] ]
,{y,1,M-1}];
   y2 = y2 +T1/(+T1-T2);

   Wf1 = If[y2>y1, (y2-y1)*dl, 999];
   Wf2 = (M+1-y1+y2) *dl;
   Lfwh= Min[ Wf1,Wf2 ];
   Lfwh= If[ Lfwh>dl, Lfwh, dl ];

(* rms *)
   TF1 = Interpolation[ uu * AF3 ];
   TF2 = Interpolation[ uu^2 * AF3 ];
      
   IF3 = NIntegrate[ GF3[t],{t,1,M},AccuracyGoal->3];
   IU1 = NIntegrate[ TF1[t],{t,1,M},AccuracyGoal->3];
   IU2 = NIntegrate[ TF2[t],{t,1,M},AccuracyGoal->3];
   Lrms= Sqrt[ (IU2/IF3) - (IU1/IF3)^2 ] *ktE *dl;
   Lrms= If[ Lrms>dl, Lrms, dl ];

TI86[kp,kd] = tfy ;
WF86[kp,kd] = Wfwh ;
WR86[kp,kd] = Wrms;
LF86[kp,kd] = Lfwh  ;
LR86[kp,kd] = Lrms ;
*)

(*---- Pulse shape -------------------------------- *)
   Print[" "];
   Print["L, L/Ls       = ",L,"( km )  ",L/Ls];
   Print["Peak Power    = ",Pmax * Pm,"( mW )  ",
                     N[ 10 Log[10,Pmax*Pm]],"( dBm )"];
   Print["Average Power = ",Po3,"( mW )  ",
                              10 Log[10,Po3]"( dBm )"];
   Print["Width(rms,fw) = ",Wrms,"( ps )  ",Wfwh];

   Plot[ GE3[t]/Pmax,{t,1,M},PlotRange->{0,ym},
                     AxesLabel ->{"time","Intensity"}];
   Plot[ LF3[t],{t,1,M},PlotRange->{-40,10},
                 Frame -> True,GridLines -> Automatic];
   Print[" "]
   ];  

  Ei= E3
 
,{k,1,kmax} ];

Array[ nse1,kmax/kgt] >> C:\soliton\nse1.m   ;
Print[Date[ ]]
(*----------------------------------------------------*)

\( (*\ T = 10  ps, \ 400/50  km, \ fw = 3  nm, \ 
      NF = \(-100\) dB, 
Dp = 0.2, \ Pm = 39.42  mW, \ 
      sl = 0.00\ *) 

\(nse1\  = << 
          C:\soliton\nse1.m\ \ \ ;\)\)
\(ListPlot3D[\ nse1, PlotRange -> {0, 6}, BoxRatios -> {1, 2, 0.5}, 
	
      ViewPoint -> { .5, \(-2\),  .4}\ ]
	\ \ \ \ \ \)

\(\(nse1\  = << C:\soliton\nse1.m\ \ \ ;\)\)
\(ListContourPlot[\ \(-nse1\), \ PlotRange -> {0, \(-4\)}, 
	
      Contours -> {0, \(-0.25\), \(-0.5\), \(-0.75\), \(-1.0\), \(-1.25\), \(-1.5\), \(-2\), \(-3\), \(-4\)}, 
	
      ColorFunction\  -> \ Hue]\ \ \ \ \ \)