Wave Form ( WDM+PMD )( pwdm ) 
( T=7.5ps,  40Gb/s,  1000/50km )
 1995.10.6  by K.Shimoura

1998. 8. 21   last  revised
2001. 11. 19  slightly revised

(*-----Initial Conditions-----------------------------*)
(* LD parameters *)
tfw=7.5    ;PdB=+3     ;GHz=40     ;
wl =1552   ;dGz=200    ;ch =+0.0   ;

(* Fiber parameters *)
s  =50.0   ;aa =0.20   ;Sl =0.00   ;Rt =5.0    ;
Dav=+0.03  ;dDp=0.0    ;kd =5      ;PMD=0.0    ;
Dc =+400   ;Nc =1      ;n2 =2.24 10^-20        ;

(* Amplification/Receiver parameters *)
ka =50     ;NF =+5.0   ;fw =100    ;
Rfil=1.4   ;Rthr=0.3   ;Rdec=0.6   ;

(* Calculation parameters *)
h  =1000   ;dt =0.5    ;NM =16     ;ym =1.0    ; 
kmax=1000  ;kw =100    ;kgt=25     ;geta=0     ;

(* Graphics Array *)
Array[   pwdm1,Round[kmax/kgt]];
Array[ snpwdm1,Round[kmax/kgt]];
Array[ pwpwdm1,Round[kmax/kgt]];

Print[ Date[ ]];
Pm = 10 ^ (0.1 PdB);
Lmax= kmax + geta;
dm  = Round[1000/GHz/dt];
M   = 2 Quotient[NM dm,2];
dl  = 3.34 10^-6 ( wl^2 /(M dt) );

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

(* Chirped Gaussian *)  
kt= 1.665 ; t0= tfw /kt ;m= 1.0    ;  
dth=(2 Pi dt dGz)/1000  ;dn =dm/NM ;

Ex0=Table[(Exp[ -.5(1 + ch I)((t+7.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+6.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t+5.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+4.5dm) dt/t0) ^(2m)]  *1
	      +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)]  *1
	      +Exp[ -.5(1 + ch I)((t-2.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t-3.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-4.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-5.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-6.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t-7.5dm) dt/t0) ^(2m)]  *0)
	      
	      *(+Exp[ +4t dth I ]  *1
            +Exp[ +3t dth I ]  *0
	        +Exp[ +2t dth I ]  *1
            +Exp[ +1t dth I ]  *0
            +      1           *1
            +Exp[ -1t dth I ]  *0
            +Exp[ -2t dth I ]  *1
            +Exp[ -3t dth I ]  *0
            +Exp[ -4t dth I ]  *1)	                                              
                                        ,{t,-M/2,M/2-1} ];
Xch= 5;

Ey0=Table[(Exp[ -.5(1 + ch I)((t+7.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+6.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t+5.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t+4.5dm) dt/t0) ^(2m)]  *1
	      +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)]  *1
	      +Exp[ -.5(1 + ch I)((t-2.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t-3.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-4.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-5.5dm) dt/t0) ^(2m)]  *1
	      +Exp[ -.5(1 + ch I)((t-6.5dm) dt/t0) ^(2m)]  *0
	      +Exp[ -.5(1 + ch I)((t-7.5dm) dt/t0) ^(2m)]  *0)
	      
	      *(+Exp[ +4t dth I ]  *0
            +Exp[ +3t dth I ]  *1
	        +Exp[ +2t dth I ]  *0
            +Exp[ +1t dth I ]  *1
            +      1           *0
            +Exp[ -1t dth I ]  *1
            +Exp[ -2t dth I ]  *0
            +Exp[ -3t dth I ]  *1
            +Exp[ -4t dth I ]  *0)
                                        ,{t,-M/2,M/2-1} ];
Ych=4;

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


(* Dispersion length -----------------------*)
Lt = kmax h/1000;
La = ka h/1000;
att= 10^( -aa h /40000);
Dp = Dav - ( 1000 Dc )/( 2 h ka Nc );

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["Dc, Nc      = ",Dc," ( ps/nm ) ",Nc];
Print["Dav,        = ",Dav," ( 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} ];
EEX= Interpolation[ Re[ Ex0 ] ];
Plot[ EEX[t],{t,1,M},PlotRange-> {-ym*Xch,ym*Xch},
                  AxesLabel ->{"time","Ex-field"}];
EEY= Interpolation[ Re[ Ey0 ] ];
Plot[ EEY[t],{t,1,M},PlotRange-> {-ym*Ych,ym*Ych},
                  AxesLabel ->{"time","Ey-field"}];
                  
AE0= (Abs[Ex0])^2 + (Abs[Ey0])^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) *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[ Ex0 ] + InverseFourier[ Ey0 ];
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 ------------------------*)
Dpl= Table[ Dav +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 ] ];

(* 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,(Xch+Ych)*ym},
               AxesLabel ->{"time","Intensity"}];
Plot[ FFU[t],{t,1,M},PlotRange-> {-30,20},
               Frame -> True,GridLines -> Automatic];
               

(*---- Pulse propagation by S.S.F.------------------*)

Ei = Ex0;
Ej = Ey0;
Do[

(* If[ Mod[k,kgt] == 1 && k > geta,
   pwdm1[(k-geta-1)/kgt +1]= Ei ]; *)
      
   Dpx= Dp + If[ EvenQ[ Quotient[k,kd]],+dDp,-dDp];
   
   gm = -1.05 10^-8 wl^2 Dpx fr^2
         +1.16 10^-14 wl^4 Sl fr^3; 
   A  =  att Exp[ -(gm h/2)I ];
   B  = -( 3.1416*10^18 (n2/wl) (Pm/s) )I;

   Fx0 = InverseFourier[Ei];
   Ex1 = Fourier[A * Fx0];
   Fy0 = InverseFourier[Ej];
   Ey1 = Fourier[A * Fy0];

(* PMD *) 
   Pmod= Random[ Integer,{-1,1} ];
   Dx1 = Append[ Rest[Ex1],First[Ex1] ];
   Px1 = PMD * Sqrt[ h/1000 ] * (Dx1 - Ex1)/dt;
   Ex1 = Ex1 + Pmod * Px1;
   Dy1 = Append[ Rest[Ey1],First[Ey1] ];
   Py1 = PMD * Sqrt[ h/1000 ] * (Dy1 - Ey1)/dt;
   Ey1 = Ey1 - Pmod * Py1;
      
   Bx1 = Exp[ h B ( (Abs[Ex1])^2 + 0.667(Abs[Ey1])^2 )];
   Ex2 = Bx1 * Ex1;
   Fx2 = InverseFourier[Ex2];
   Fx3 = A * Fx2;
   Ex3 = Fourier[Fx3];

   By1 = Exp[ h B ( (Abs[Ey1])^2 + 0.667(Abs[Ex1])^2 )];
   Ey2 = By1 * Ey1;
   Fy2 = InverseFourier[Ey2];
   Fy3 = A * Fy2;
   Ey3 = Fourier[Fy3];

(*------------Dispersion Compensation--------------*)
If[ Mod[k,ka*Nc] == 0 && OddQ[Quotient[k,ka*Nc]],
   A3 = Exp[ 1.05 10^-5 wl^2 fr^2 Dc I ];
   Ex3 = Fourier[A3 * Fx3];
   Fx3 = InverseFourier[Ex3];
   Ey3 = Fourier[A3 * Fy3];
   Fy3 = InverseFourier[Ey3];     
   ];

(*------------Amplification and Filtering----------*)
If[ Mod[k,ka] == 0,
  
(* Filter *)
   ktf= 0.6931    ;fm = 1      ;
   fc = dGz M dt/1000;
   fwt= fw / dl;
      
   FIL= Table[ 
               +Exp[-ktf Abs[2(f+4fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f+3fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f+2fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f+ fc) /fwt]^(2fm)] *0
               
               +Exp[-ktf Abs[2f/fwt]^(2fm)] *1
               
               +Exp[-ktf Abs[2(f- fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f-2fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f-3fc) /fwt]^(2fm)] *0
               +Exp[-ktf Abs[2(f-4fc) /fwt]^(2fm)] *0
               
                                      ,{f,-M/2,M/2-1}];
   FIL= FIL / Max[FIL];
                                 
(*
(* 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)" ];
*)      
 
   Fx4 = Fx3 * Sqrt[ Part[FIL,tt] ];
   Ex4 = Fourier[ Fx4 ];
   Fy4 = Fy3 * Sqrt[ Part[FIL,tt] ];
   Ey4 = Fourier[ Fy4 ];

(* Amplifier *)
   dAp= 0.0;
   GE4= Interpolation[ (Abs[Ex4])^2 + (Abs[Ey4])^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;                        
   Fx5 = Fx4 * Sqrt[ AMP ];
   Ex5 = Fourier[ Fx5 ];
   Fy5 = Fy4 * Sqrt[ AMP ];
   Ey5 = Fourier[ Fy5 ];

(* Noise *)
   ASE= 0.994 10^-4 * 10^(NF/10) * AMP * (GHz/wl); 
   WNE= Table[ ( ASE/Pm ) * Exp[6.28 Random[] I],{M}];
   Ex3 = Ex5 + Sqrt[ WNE * Xch ];
   Fx3 = InverseFourier[ Ex3 ];
   Ey3 = Ey5 + Sqrt[ WNE * Ych ];
   Fy3 = InverseFourier[ Ey3 ];  	     
   ];   


(*----Pulse Width----------------*)
If[ Mod[k,kgt] == 0 && k > geta,      
   
   AE3 = (Abs[Ex3])^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 ]; *)      
   pwpwdm1[ (k-geta)/kgt ]= Wfwh /tfw;
   ];


(*---- Q-factor Calculation ----------------------------*)
If[ Mod[k,kgt] == 0 && k > geta,      

(* E-filter *)
   ktf= 0.6931  ;fm = 1  ;
   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}];
   Fx4 = Fx3 * Sqrt[ Part[ FIE,tt ] ];
   Ex4 = Fourier[ Fx4 ];   
   Fy4 = Fy3 * Sqrt[ Part[ FIE,tt ] ];
   Ey4 = Fourier[ Fy4 ];   
   
   AE4 = (Abs[ Ex4 ])^2;
(* AE4 = (Abs[ Ex4 ])^2 + (Abs[ Ey4 ])^2; *)
   Ymax= Max[ AE4 ];
   GE4 = Interpolation[ AE4 ];
   
   pwdm1[(k-geta)/kgt]= Ex4;
   
(* 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 ];*)
 
   snpwdm1[ (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[ Ex3 ])^2 + (Abs[ Ey3 ])^2;
   Pmax= Max[ AE3 ];
   GE3 = Interpolation[ AE3 ];
   IE3 = NIntegrate[ GE3[t],{t,1,M},AccuracyGoal->3];
   Po3 = IE3 Pm dt GHz /(1000 NM);

   (* single-channel *)
   ktf= 0.6931  ;fm = 1  ;
   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}];
   Fx4 = Fx3 * Sqrt[ Part[ FIE,tt ] ];
   Ex4 = Fourier[ Fx4 ];   
   Fy4 = Fy3 * Sqrt[ Part[ FIE,tt ] ];
   Ey4 = Fourier[ Fy4 ];   
   
   AE4 = (Abs[ Ex4 ])^2;
   Ymax= Max[ AE4 ];
   GE4 = Interpolation[ AE4 ];
   
   AF3 = Part[ Abs[Fx3],tt ]^2 + Part[ Abs[Fy3],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[ GE4[t]/Ymax,{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= Ex3;
 Ej= Ey3
 
,{k,1,Lmax} ];

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

(* T = 7.5  ps, 1000/50  km, fw = 100  nm, NF = +5 dB,
Dav = +0.03, Pm = 1  dBm, Dc = +400, Nc = 1, ch = +0.0 *)

\(pwdm1\  = << C:\soliton\pwdm1.m;\)
Power3D  = Abs[ pwdm1 ^2] *10;
ListPlot3D[ Power3D,   PlotRange -> {0, 1.1}, 
          BoxRatios -> {1, 2, 0.4}, 
Mesh  ->  False, 
          ViewPoint -> {1, -2,  .4} ];
\(ListPlot3D[\ 
        Log[10, Power3D], \ \ PlotRange -> {\(-5\), 
            1}, 
	\ \ \ \ \ \ \ \ \ \ \ BoxRatios -> {1, 2, 
            0.8}, 
Mesh\  -> \ False, 
        ViewPoint -> {3, 0, 0.4}\ ];\)

\(
\(pwpwdm1\  = << C:\soliton\pwpwdm1.m;\)
\)
\(\(ListPlot[\ pwpwdm1, \ 
        PlotRange -> {0.8, 1.5}, 
	\ \ \ \ \ \ \ \ \ PlotJoined\  -> 
          True, \ AxesLabel\  -> {"Lt", "t(FWHM)"}\ ];\)\)

\(
\(snpwdm1\  = << C:\soliton\snpwdm1.m;\)
\)
Qfactor  =  Sqrt[ snpwdm1 ];
\(\(ListPlot[\ Qfactor, \ 
        PlotRange -> {0, 
            20}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotJoined\  -> 
          True, \ AxesLabel\  -> {"Lt", "Q-fac"}\ ];\)\)
\(\(ListPlot[\ 20\ Log[10, Qfactor]\ , \ 
        PlotRange -> {0, 30}, 
	\ \ \ \ \ \ \ \ PlotJoined\  -> True, \ 
        AxesLabel\  -> {"Lt", "20 LogQ"}\ ];\)
\)

\(pwdm1\  = << C:\soliton\pwdm1.m;\)
Power3D  = 
        Abs[ pwdm1 ^2] *10 ;
\(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)"\ ]
\)