Wave Form Calculation Code  ( wave ) 
( T=10ps,  20Gb/s,  10000/50km )
 1995.10.6  by K.Shimoura  

1998. 8. 21   last  revised


(*-----Initial Conditions-----------------------------*)
(* LD parameters *)
wl =1552   ;tfw=10     ;Pm =12.59  ;GHz=20     ;
ch =+0.4   ;dnm=2.0    ; 

(* Fiber parameters *)
s  =50.0   ;aa =0.20   ;Sl =0.07   ;Rt =5.0    ;
Dp =-0.473 ;dDp=0.0    ;kd =5      ;
Dc =+80    ;Nc =3      ;n2 =2.24 10^-20        ;

(* Amplification/Receiver parameters *)
ka =5      ;NF =6      ;Neff=0.1   ;fw =3      ;
dAp=0.0    ;Rfil=1.4   ;Rthr=0.3   ;Rdec=0.6   ;

(* Calculation parameters *)
h  =10000  ;dt =2      ;NM =10     ;ym =1.0    ; 
kmax=1000  ;kw =100    ;kgt=20     ;geta=0     ;

(* Graphics Array *)
Array[   wave1,Round[kmax/kgt]];
Array[ snwave1,Round[kmax/kgt]];
Array[ pwwave1,Round[kmax/kgt]];
Print[ Date[ ]];

Lmax= kmax + geta;
dm  = Round[1000/GHz/dt];
M   = 2 Quotient[NM dm,2];

(*----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;
*)

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

E0= Table[ Sech[(t+4.5dm) dt/t0]  *0
		  +Sech[(t+3.5dm) dt/t0]  *1
          +Sech[(t+2.5dm) dt/t0]  *0
		  +Sech[(t+1.5dm) dt/t0]  *1
          +Sech[(t+0.5dm) dt/t0]  *1
          +Sech[(t-0.5dm) dt/t0]  *1
          +Sech[(t-1.5dm) dt/t0]  *0
          +Sech[(t-2.5dm) dt/t0]  *1
          +Sech[(t-3.5dm) dt/t0]  *1
          +Sech[(t-4.5dm) dt/t0]  *0
                                       ,{t,-M/2,M/2-1}];
*)

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

EE0=Table[ ( Sech[(t-1.5dm) dt/t0]*Exp[ phd Random[] I ] *1
            +Sech[(t-0.5dm) dt/t0]*Exp[ phd Random[] I ] *0
            +Sech[ t dt/t0]  *0
            +Sech[(t+0.5dm) dt/t0]*Exp[ phd Random[] I ] *1
            +Sech[(t+1.5dm) dt/t0]*Exp[ phd Random[] I ] *1)
                      
          *( Exp[-2t dth I ]  *1
            +Exp[-1t dth I ]  *1
            +      1          *0
            +Exp[+1t dth I ]  *1
            +Exp[+2t dth I ]  *1 )
                                     ,{t,-M/2,M/2-1} ];

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

(* Chirped Gaussian *)  
kt= 1.665 ; t0= tfw /kt ; m= 1.0;  

E0= Table[ Exp[ -.5(1 + ch I)((t+4.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+3.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t+2.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+1.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t+0.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-0.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-1.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t-2.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-3.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-4.5dm) dt/t0) ^(2m)]  *0			
                                        ,{t,-M/2,M/2-1} ];


(* Chirped Gaussian 2  
kt= 1.665 ; t0= tfw /kt ; m= 1.0;  
E0= Table[ Exp[ -.5(1 + ch I)( t dt/t0) ^(2m) ]		
                                        ,{t,-M/2,M/2-1} ];
*)

(* 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];
      
(* 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);
trms= Sqrt[ (ITV/IEU) - (ITU/IEU)^2 ] * dt;

Print["W(FWHM),chirp = ",tfw,"( ps )  ",ch ];
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[ Reverse[ 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[ Reverse[ Ddf ] ];

(* Dispersion compensation parameters-------*)
Dav= ( Dp h ka Nc + 1000 Dc ) / ( h ka Nc );
Print["Dc, Nc      = ",Dc," ( ps/nm ) ",Nc];
Print["Dav         = ",Dav," ( ps/nm/km )"];

(* Initial shape drawing--------------------*)
Plot[ DDT[t],{t,1,M},PlotRange-> {-2,+2},
               AxesLabel ->{"W.L.","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.------------------*)

Ei = E0;
Do[

If[ Mod[k,kgt] == 1 && k > geta,
   wave1[(k-geta-1)/kgt +1]= Ei ];
      
   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];


(*------Dispersion Compensation before Amplicication------*)
If[ Mod[k,ka*Nc] == 0,
   A3 = Exp[ 25.2 fr^2 Dc I ];
   E3 = Fourier[A3 * F3];
   F3 = InverseFourier[E3]     
   ];
   
   
(*------------Amplification and Filtering-----------------*)
If[ Mod[k,ka] == 0,
  
(* Filter *)
   ktf= 0.6931;
   fm = 1     ;
   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); 
   gau= 0.441 ( M dt /tfw );
   
   WNE= Table[ ( ASE/Pm ) * ( Neff * gau ) 
              * 2 Random[] * Exp[6.28 Random[] I],{M}];
   E3 = E3 + Sqrt[ WNE ];
   F3 = InverseFourier[ E3 ];  
         
   F4 = F3 * Sqrt[ Part[FIL,tt] ];
   E4 = Fourier[ F4 ];
   
   GE4= Interpolation[ (Abs[E4])^2 ];
   IE4= NIntegrate[ GE4[t],{t,1,M},AccuracyGoal->3];
   
   dmp= If[ EvenQ[ Quotient[k,ka] ]
                         ,10^(-dAp/10),10^(dAp/10)];
   AMP= IEU/IE4 * dmp;                        
(* AMP= Ga; *)
(* Print["AMP  = ",10 Log[10,AMP]," ( dB )"]; *)

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

(*----Pulse Width----------------*)
If[ Mod[k,kgt] == 0 && k > geta,      
   
   AE3 = (Abs[E3])^2;
   Pmax= Max[ AE3 ];
   GE3 = Interpolation[ AE3 ];
      
(* 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 ];

 (*Print[ "Wfwh   =  ", Wfwh ]; *)      
   pwwave1[ (k-geta)/kgt ]= Wfwh /tfw;
   ];


(*---- Q-factor Calculation ----------------------------*)
If[ Mod[k,kgt] == 0 && k > geta,      
   
(* E-filter *)
   ktf= 0.6931;
   fwR= Rfil * GHz ( wl^2 10^-8 /2.998 );
   fwQ= fwR / dl; 
   FIE= Table[ Exp[ -ktf Abs[ 2f/fwQ ]^(2fm) ]
                                      ,{f,-M/2,M/2-1}];
   F4 = F3 * Sqrt[ Part[ FIE,tt ] ];
   E4 = Fourier[ F4 ];   
   
   AE4 = (Abs[ E4 ])^2;
   Ymax= Max[ AE4 ];
   GE4 = Interpolation[ AE4 ];
   
(* Auto Search ------------*)   
   ww  = AE4;   
   AEF = 0;
   Do[
      AEF = AEF + Take[ ww, dm ];
      ww  = Drop[ ww, dm ]
   ,{t,1,NM}];
   
   GEF = Interpolation[ AEF ];
   IX  = NIntegrate[ GEF[t],{t,1,dm},AccuracyGoal->3];
   AMAX= Max[ AEF ];
 (*Plot[ GEF[t]/Pmax,{t,1,dm},PlotRange->{0,NM*ym},
                    AxesLabel ->{"time","Intensity"}];*)
   Do[
      CEN = t;
      If[ Part[ AEF,t ] > 0.99 AMAX, Break[ ]]
   ,{t,1,dm} ];   
      
   ww  = AE4;   
   If[ CEN > dm/2, ww = Drop[ ww, Round[CEN -dm/2] ]
                 , ww = Drop[ ww, Round[CEN +dm/2] ]];
 
 (* eye-pattern -----------*)
   nn  = Table[ t,{t,1,M} ];
   cn1 = Table[ Mod[(t +dm -CEN), 2dm],{t,1,M} ];
   yey1= { Part[cn1,nn], Part[AE4,nn] };
   eye1= Transpose[ yey1 ];
   
   cn2 = Table[ Mod[(t -CEN), 2dm],{t,1,M} ];
   yey2= { Part[cn2,nn], Part[AE4,nn] };
   eye2= Transpose[ yey2 ];
   
   eye = Join[ eye1,eye2 ];
 (*ListPlot[ eye, Prolog->AbsolutePointSize[3]
                , PlotRange->{0,Ymax} 
                (*PlotJoined -> True*) ]; *)
 
 (* Q-factor --------------*)     
   iww= Interpolation[ ww ];
   
 (* energy  
   WD = Rdec * dm;
   PWD= Table[ NIntegrate[ iww[t],{t, (n-.5)dm -WD/2,
                       (n-.5)dm +WD/2}] ,{n,1,NM-1}];*)
                       
 (* level *)
   PWD= Table[ iww[(n-.5)dm ] ,{n,1,NM-1}];  
   PWE = Sort[ PWD ];
 (*Print[ PWE ];*)
   thr = Rthr * Max[ PWE ];
  
   For[ i=1, PWE[[i]] <thr, i++ ];
   zer = Take[ PWE, i-1];
   one = Drop[ PWE, i-1];
   n0  = Part[ Dimensions[ zer ],1] +0.001;
   n1  = Part[ Dimensions[ one ],1] +0.001;

   zrm= Apply[Plus,zer] /n0;
   zrd= (zer-zrm)^2;
   zrs= Sqrt[ Apply[Plus,zrd] /n0 ];
   
   oem= Apply[Plus,one] /n1;
   oed= (one-oem)^2;
   oes= Sqrt[ Apply[Plus,oed] /n1 ];
   
   soz = zrs + oes + 0.0001;
   Qfac= Abs[ (oem-zrm)/soz ];
   
 (*Print[ n0,"  ",zrm,"  ",zrs ];
   Print[ n1,"  ",oem,"  ",oes ];
   Print[ "Qfac=  ", Qfac ];
   Print[ "Rthr,thr =  ", Rthr,"  ",thr ];*)
 
   snwave1[ (k-geta)/kgt ]= Qfac;
   ];

(*
(*---Dispersion Compensation after Amplicication---------*)
If[ Mod[k,ka*Nc] == 0,
   E3 = Fourier[ A3 * F3 ]
   ];
*)

(*---- Pulse Shape Drawing ------------------------------*)
If[ Mod[k,kw] == 0 && k > geta,
   
   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);

   AF3 = Part[ Abs[F3],tt ]^2;
   GF3 = Interpolation[ AF3 ];
   LF3 = Interpolation[ Reverse[ 10 Log[10,AF3] ] ];
   
   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 )"];*)
      
   Plot[ GE3[t]/Pmax,{t,1,M},PlotRange->{0,ym},
                    AxesLabel ->{"time","Intensity"}];
   Plot[ LF3[t],{t,1,M},PlotRange->{-30,20},
                 Frame -> True,GridLines -> Automatic];
   ListPlot[ eye, Prolog->AbsolutePointSize[3]
                , PlotRange->{0,Ymax} 
                (*PlotJoined -> True*) ];              
  ];  

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

Array[ wave1,kmax/kgt]   >> A:\soliton\wave1.m;
Array[ snwave1,kmax/kgt] >> A:\soliton\snwave1.m;
Array[ pwwave1,kmax/kgt] >> A:\soliton\pwwave1.m;
Print[ Date[ ]]
(*----------------------------------------------------*)



(* T=10ps, 10000/50km, fw=3nm, NF=6dB,
	  Dp=-0.473, Pm=11dBm, Dc=-80, Nc=3, ch=+0.4 *)

wave1 =<<A:\soliton\wave1.m;
Power3D =Abs[ wave1 ^2];
ListPlot3D[ Power3D,  PlotRange->{0,1.1},
	           BoxRatios->{1,2,0.4},ViewPoint->{1,-2,.4} ];
ListPlot3D[ Log[10,Power3D],  PlotRange->{-5,1},
	           BoxRatios->{1,2,0.4},ViewPoint->{1,-2,.4} ];



pwwave1 =<<A:\soliton\pwwave1.m;

ListPlot[ pwwave1, PlotRange->{0.8,1.5},
	         PlotJoined ->True, AxesLabel ->{"Lt","t(FWHM)"} ];



snwave1 =<<A:\soliton\snwave1.m;

Qwave1 = Sqrt[ snwave1 ];
ListPlot[ Qwave1, PlotRange->{0,20},
	                  PlotJoined ->True, AxesLabel ->{"Lt","Q-fac"} ];
(*ListPlot[ Log[10,.5Erfc[Qwave1/1.414] ], PlotRange->{0,-30},
	                  PlotJoined ->True, AxesLabel ->{"Lt","BER"} ];
ListPlot[ 20 Log[10,Qwave1] , PlotRange->{10,30},
	        PlotJoined ->True, AxesLabel ->{"Lt","20 LogQ"} ];*)



wave1 =<<A:\soliton\wave1.m;

Power3D =Abs[ wave1 ^2];
ListContourPlot[ -Power3D, PlotRange->{0,-1.4},
	Contours->{0,-0.2,-0.4,-0.6,-0.8,-1,-1.2,-1.4} ] ;
ListContourPlot[ -Log[10,Power3D], PlotRange->{-1,+5},
	Contours->{5,4,3,2,1,0,-1} ];       



(* Filter Shapes *) 
FI = Interpolation[ FIL ];
Plot[ FI[t],{t,1,M},PlotRange->{0,1},
             Frame-> True, GridLines-> Automatic ];
FIP= Interpolation[ 10 Log[10,FIL] ];
Plot[ FIP[t],{t,1,M},PlotRange->{-50,0},
             Frame-> True, GridLines-> Automatic ];
dlr= dl * 50;  Print[ dlr,"  (nm/div)" ]