NL pulse prop.nb

This program calculates the pulse shape after it propagates through an optical fiber, including  self-phase modulation (SPM)  and group velocity dispersion (GVD). It accounts for all dispersion and loss.  

Important parameters set here:
lambda0: 	center wavelength of pulse
tau: 		transform-limited pulsewidth of input
npts:		number of points for FFT, keep to power of 2
dt:		smallest timestep, usually around tau/10.
		If using grating functions, making dt too small will put omega range outside
		the range of the gratings and give nasty imaginary numbers.
		Check output spectrum and pulse duration to make sure there is no energy 
		at the edges of the array. If so, npts and/or dt must be adjusted.
nzpts:		number of z-slices. 
		Program does not have adaptive step sizes.
		Choose nzpts so that answer doesn't change with smaller dz. 

there are various useful functions defined here that are not directly used by the code later.

coding notes:

generalize code:
- derive/correct motion of time grid at group velocity OK
- normalize field to power units as in Agrawal OK
- generalize dispersion selection: materials, fiber/Taylor terms OK
- redefine scaling lengths, compare code to Agrawal
- restore option for input chirp
- make streamlined version 
later:
- restore option for gas dispersion
- intensity: allow avg of Gaussian or a waveguide mode, single intensity value
- loop over intensity values to calculate wavefront distortion due to NL prop
- restore self-steepening, Raman terms
- adapt to spatial NL prop
- work towards spatio-temporal code

other changes:
update ADK calculation
make sure wavelength scaling of ionization spm is right: put plasma dispersion into time domain?

constants and unit definitions

ps = (10^(-12)); fs = (10^(-15));
cm = (10^(-2));  cm3 = (10^6); mm = (10^(-3)); um = (10^(-6)); 
    nm = (10^(-9)); km = (10^3); 
    angs = (10^(-10));
mJ = (10^(-3)); uJ = (10^(-6)); 
    nJ = (10^(-9));        
c = 2.9979 (10^8);     
eps0 = 8.8542 (10^(-12)); 
\(\(re = 2.82\ 10\^\(-15\);\)	\)
\(\(e = 1.6\ 10\^\(-19\);\)	\)
Natm = 760 3.3 (10^16)/cm^3; 
Torr = 1/760;
me = (10^(-30));
W = 
      1;GW = (10^9) W;

refractive index and dispersion defs

Grating phase functions
Lcom: separation of gratings along surface normal
γcom: angle of incidence of gratings
dcom: groove spacing 
Notes:
- large spectral arrays may lead to imaginary phase because of (Sqrt[□]). Physically, wavelengths with imaginary phase have diffraction angles beyond 90deg. Make sure spectral range of grid doesn't produce imaginary phases.
- Negative separations correspond to a stretcher that gives positive (normal) dispersion.

ϕGraCom[λ_, Lcom_, γcom_, 
          dcom_] := ((2 π)/λ) Lcom (Sqrt[1 - ((((λ/dcom) - Sin[γcom Degree]))^2)]);
ϕ1GraCom[λ_, Lcom_, γcom_, 
          dcom_] = -((λ^2)/(2 π c)) (∂_lambda) ϕGraCom[λ, Lcom, γcom, 
              dcom];

define taylor expanded phase

Δω[λ_] := 
        2 π c (1/lambda0 - 1/λ);
ϕTaylor[λ_, kk2_, kk3_, kk4_, kk5_] := 
        kk2 Δω[λ]^2 + 
          kk3 Δω[λ]^3 + 
          kk4 Δω[λ]^4 + 
          kk5 Δω[λ]^5;
ϕ1Taylor[λ_, kk2_, kk3_, kk4_, 
          kk5_] = -((λ^2)/(2 π c)) D[ϕTaylor[λ, kk2, kk3, kk4, kk5], λ];
ϕ2Taylor[lam_, kk2_, kk3_, kk4_, 
          kk5_] = -((λ^2)/(2 π c)) D[ϕ1Taylor[λ, kk2, kk3, kk4, kk5], λ];
ϕ3Taylor[lam_, kk2_, kk3_, kk4_, 
          kk5_] = -((λ^2)/(2 π c)) D[ϕ2Taylor[λ, kk2, kk3, kk4, kk5], λ];

fused silica refractive index nSiO2[λ]
n2 value from Boyd table 4.1.2

as1 =  .6961663;
as2 =  .4079426;
as3 =  .8974794;
bs1 =  .0684043^2;
bs2 =  .1162414^2;
bs3 = 9.896161^2;
nSiO2[lambda_] := 
        Sqrt[1 + (as1 (lambda/um)^2)/((lambda/um)^2 - 
                  bs1) + (as2 (lambda/um)^2)/((lambda/
                        um)^2 - 
                  bs2) + (as3 (lambda/um)^2)/((lambda/
                        um)^2 - bs3)];
n2SiO2 = 3.2 (10^(-16)) (cm^2)/W;

set up for calculation

define FFT functions

FFT[fLst_, npts_, dt_] := (Sqrt[((npts (dt^2))/(2 π))]) 
          RotateLeft[Fourier[RotateRight[fLst, npts/2]], 
            npts/2];
InverseFFT[FLst_, npts_, 
          dt_] := (Sqrt[((2 π)/(npts (dt^2)))]) 
          RotateLeft[InverseFourier[RotateRight[FLst, npts/2]], 
            npts/2];

define input pulse functions
I = \(TraditionalForm\`\(1\/2\) n\ ϵ0\ c\  | E\( | \^2\)\  = \ 
        Ipk\ \(f(t)\)\ \(g(r)\)\)
note correction for refr index: this may affect value of n2, check calc of E-field from intensity
choose f(t),  g(r), integrate t, r to get total energy, calculate Ipk
calculate Epk from Ipk

lambda[w_] := 2 π c/w;
fwhmFactor = 4 Log[2];
gaussT[t_, tauTmp_] :=  
        Exp[-fwhmFactor ((((t/tauTmp)))^2)];
pkFactor = 1;
pkInten[eplsTmp_, spotradTmp_, tauTmp_, pkFactorTmp_, n_] := 
        n (( pkFactorTmp eplsTmp)/(π (spotradTmp^2) tauTmp));
ef0[
          intenD_] := (Sqrt[((2 intenD)/(eps0 c))]);
inten[ef_]  :=  
        Re[nEff[lambda0] eps0 c Conjugate[ef]*ef/2]  // 
          N;
A0[intenD_, area_] := (Sqrt[intenD area]);

gaussR[r_, wTmp_] := 
        Exp[-((((r/wTmp)))^2)];
areaEffGauss[
          wTmp_] = ((((2 π Integrate[r gaussR[r, wTmp]^2, {r, 0, ∞}, GenerateConditions -> False]))^2)/(2 π Integrate[r gaussR[r, wTmp]^4, {r, 0, ∞}, GenerateConditions -> False]));
x0 = 2.4048;
besselJ0R[r_, wTmp_] := 
        BesselJ[0, x0 r/wTmp];
areaEffBesselJ0[
          wTmp_] = ((((2 π Integrate[wTmp rp BesselJ[0, x0 rp]^2, {rp, 0, 1}, GenerateConditions -> False]))^2)/(2 π NIntegrate[rp BesselJ[0, x0 rp]^4, {rp, 0, 1}]));

areaEffGauss[w]
areaEffBesselJ0[w]

Define characteristic lengths for dispersion and nonlinearities
Note: the soliton number works out to be = 1 when an e-field in the form of a Sech[τ/T] is used. But this T gives a pulse width that is larger by about 1/0.56 than the T in a gaussian Exp[-(t^2)/(T^2)]. So if you are choosing parameters to generate solitons, you must account for this extra factor.

dispLen[ω_, 
          tauD_] := ((tauD^2)/fwhmFactor) (1/Abs)[β2[ω]] // 
          N;
nonlinLen[γD_, P0_] := 
        (1/(γ P0));
solitonN[ω_, tauD_, γ_, P0_] := 
        Sqrt[dispLen[ω, tauD 0.565]/nonlinLen[γ, P0]] // 
          N;

initialization for calculation

initialization and input

notes

Coding: 
change internal field units to Sqrt[power] as in Agrawal, use γ instead of n2
incorporate different types of dispersion: Taylor (β1,β2), glasses, gases 

Important adjustable parameters:
		0 = off, 1 = on:
		lossesOn	for fiber losses
		sstOn		for self-steepening
		RamanOn	for raman nonlinearities (air only, so far)
		IonizeOn	for ionization
		
	capillary parameters (use units! internal calculations are done mostly in SI units)
		a		fiber radius
		zMax		fiber length
		pr		gas pressure
		pEl		equivalent pressure of static free electrons
		epulse		input pulse energy
		n2factor	this multiplies the n2 for the gas chosen to scale it. Usually need 
				it to be around 1/3 to match measured spectra (defined in spm init)
	input chirp parameters: code first generates transf. limited pulse, FWHM tau, then 
	applies input chirp with either a grating pair or Taylor exp. terms
		gamacom	grating incident angle, in degrees
		dcom		groove density
		lggcom	grating separation relative to the compressed position, 
				negative value puts on positive chirp
		k2,k3,k4,k5	Taylor coeffic., see notes in that section
		lsFrog		amount of material (fused silica) before the fiber

Calculation notes: 
try to subtract off energy that goes into ionization

input parameters
a = beam/fiber radius

a = 5 um; 
zMax = 20 cm;
epulse = 4 3/(1.69^2) nJ;
lambda0 = 800 nm;
k0 = 2 π/lambda0;
ω0 = k0 c;
tau = 50 fs;

x compressor and input chirp settings

gamacom=20;			(*compressor incident angle*)
dcom=1/600 mm;		(*groove spacing*)
lggcom=-15 mm 0;
(*set orders of spectral phase on input pulse*)
k2=0 300 fs^2;  		(*positive is positive chirp*)
k3=0 fs^3; 		(*negative puts ringing on back (i.e.pos.third-order)*)
k4=0 fs^4;
k5=0 fs^5;
(*subtract material of frog*)
lsFrog=0 mm;

select dispersion for medium
for now, using fused silica parameters
later, add gas, free electrons, modes

no dispersion

nEff[ω_] = 
        nSiO2[lambda[ω0]];
β[ω_] := (ω/c) 
          nEff[ω];
n2 = n2SiO2;
areaEff := areaEffGauss[a]
γ := ((k0 n2)/(areaEff ));

bulk fused silica

nEff[ω_] = 
        nSiO2[lambda[ω]];
β[ω_] := (ω/c) 
          nEff[ω];
n2 = n2SiO2;
areaEff := areaEffGauss[a]
γ := ((k0 n2)/(areaEff ));

fiber with taylor series dispersion

nFib = nSiO2[lambda0];
β1SiO2 = 
        D[(ω/c) 
            nSiO2[lambda[ω]], ω];
β1Fib = β1SiO2 /. ω -> ω0;
β2Fib = -20 (ps^2)/km;
β[ω_] = ω(0/c) 
            nFib + (ω - ω0) β1Fib  + (1/2) (((ω - ω0))^2) β2Fib;
nEff[ω_] = β[ω] 
          c/ω;
γ = 2/W/km;
areaEff := areaEffGauss[a]
γ := ((k0 n2)/(areaEff ));

notes on dispersion and moving frame:

NLS 
\(TraditionalForm\`∂\_z\ 
          A\  = \ \(-\ β\_1\) ∂\_t\ 
                A - \ i\ \(1\/2\) \(β\_2\) ∂\_\(t, t\)A\  - \(α\/2\) A\  + i\ γ\  | 
          A\( | \^2\)A\)
γ = (((Subscript[k, 0]) (Subscript[n, 2]))/(Subscript[A, eff])), where (Subscript[A, eff]) is the spatial mode effective area
This expression uses a Taylor expansion for dispersion:
\(TraditionalForm\`β(ω) = \(ω\/c\) \(\(n\_eff\)(ω)\) ≈ β\_0 + \((ω - ω\_0)\) β\_1 + \(1\/2\) \(\((ω - ω\_0)\)\^2\) β\_2 + \(1\/6\) \(\((ω - ω\_0)\)\^3\) β\_3 + ⋯\)
(Subscript[n, eff]) is the effective index, which includes any refractive index dispersion as well as any modal dispersion. 
(Subscript[β, n]) = ((∂_ω)^n) β  
change to a frame of reference moving at the pulse group velocity: 
T = t - z/(Subscript[v, g]) = t - z (Subscript[β, 1])
(∂_t) 
          A  =  ∂(T/∂)t (∂_T) 
              A  = (∂_T) A 
(∂_z) A(z,  
            T(t, z))  =  (∂_z) 
              A + ∂(T/∂)z (∂_t) 
                A = (∂_z) 
              A - (Subscript[β, 1]) (∂_t) A
If we had no nonlinearity, we would work in the spectral domain by doing a FT of the equation: 
\(TraditionalForm\`∂\_z\ 
          A\_ω\  = \ \(-\ β\_1\) ℱ {∂\_t\ 
                A} - \ i\ \(1\/2\) \(β\_2\) ℱ {∂\_\(t, t\)A}\  - \(\(α\/2\) \(A\_ω\)\(\ \)\)\)
\(TraditionalForm\`ℱ {∂\_t\ 
              A}\  = \ \(∫\((∂\_t\ 
                  A)\)\ \(e\^\(i\ ω\ t\)\) dt\  = \ \(A\ e\^\(i\ ω\ t\)\( | \_\(-∞\)\)\^∞\ \(-i\)\ ω \(∫\ 
                  A\ \(e\^\(i\ ω\ t\)\) dt\) = \(\(-i\)\(\ \)\(ω\)\(\ \)\(A\_ω\)\(\ \)\)\)\)\)
\(TraditionalForm\`ℱ {∂\_\(t, t\)\ 
              A}\  -> \ \(\(-ω\^2\)\(\ \)\(A\_ω\)\(\ \)\)\)
\(TraditionalForm\`∂\_z\ 
          A\_ω\  = \ \(i\ ω\ \(β\_1\) 
              A + i\ \(1\/2\) \(ω\^2\) \(β\_2\) 
              A\_ω\  - \(α\/2\) 
              A\_ω\  = \ i\ \(β(ω)\)\ A\)\)

dispersion parameters
βR  = 
        k0 ((Subscript[n, eff])(ω)  -  (Subscript[n, eff])((Subscript[ω, 0]))) real part of propagation constant
β1 = dk/dω = 1/vgroup
nGroup = c β1

β[ω] dz (ω/c) nEff dz
-> (β[ω] - β[ω0] - (ω-ω0)β1[ω0]) dz = \(TraditionalForm\`\((\(ω\/c\) nEff[ω] - \(ω0\/c\) 
                nEff[ω0] - \((ω - ω0)\) β1[ω0])\)\ dz\  = \ \(1\/c\) \((ω\ nEff[ω] - ωc\ β1[ω0] + ω0\ \((c\ β1[ω0] - 
                    nEff[ω0]\ )\)\ )\)\ dz\)

β0 = β[ω0];
β1[ω_] = 
        D[β[ω], ω];
β2[ω_] = 
        D[β1[ω], ω];
βR[ω_] := β[ω] - β0 - (ω - ω0) β1[ω0];   (* 
      grid moves at group velocity *) 

 (* βR[ω_] := β[ω];  
      stationary grid *) 
βI[ω_] = 
        0 ;  (* no loss *) 
\(\(Print[β2[ω0]/\((ps\^2/
              km)\), \*"" \!\(ps\^2\)/km""];\)\)

Ipk = 
        pkInten[epulse, a, tau, pkFactor, 1];
\(\(Print["input peak intensity: ", 
        Ipk\ \ cm\^2, \*"" W/\!\(cm\^2\)""];\)\)
Print["nominal B-integral: ", 
        k0 n2 Ipk zMax]; 
Print["dispersion length: ",  
        dispLen[ω0, tau]/cm, " cm"];
Print["nonlinear length: ", 
        nonlinLen[γ, Ipk areaEff]/
          cm, " cm"];
Print["soliton number: ", 
        solitonN[ω0, tau, γ, 
          Ipk areaEff]];

calculation/grid parameters

npts = (2^14);  
dt = tau/35;
nzpts = 100;
dw := 
        2 π/(npts dt);

Print["npts: ",  npts];
Print["time range of array: +/-",  
        npts dt/2 /fs  // 
          N,  " fs"];
Print["wavelength range of array: ",  
      lambda[ω0 + dw npts/2]/nm, " to ", 
      lambda[ω0 - dw npts/2]/nm, " nm"]

Set up time, frequency, and spatial grids:
    Note that t=0 and w = w0 at npts/2 in array
    need to implement checks on frequency space to make sure grid is
    large enough.

timeLst = 
        Table[t, {t, -(npts/2 - 1) dt, (npts/2) dt, 
            dt}];
omegaLst = 
        Table[w + ω0, {w, -(npts/2 - 1) dw, (npts/2) dw, 
            dw}];
wvlengthLst = 2 π c/omegaLst // N;
dz := zMax/nzpts // N;

Print[dz nEff[ω0]/c /fs];
Print[dz β1[ω0] /fs];
Print[
        nzpts dz (β1[ω0] - nEff[ω0]/c) /
            fs];

set up field arrays in both domains for all z points
At: (Sqrt[power]) in time domain (proportional to E-field)
Aω: (Sqrt[power]) in freq domain

AtLst = 0 Range[nzpts + 3];
AωLst = 0 Range[nzpts + 3];

calculation options
control whether fiber losses and self-steepening terms are used, set to 1 if using these 
   time delay not yet implemented in this model

lossesOn = 0;
inputPhaseOn = 0;

display/plotting functions 
define these so that they will adjust if run parameters (e.g. tau) are changed

<< Graphics`Colors`
$TextStyle  =  {FontFamily -> Times,  
          FontSize -> 12};

define plotting ranges

timeRange := 3 tau/fs;
lamRange = 400;

calculate fwhm from field, using gaussian fit

tsamp = half width of time window around guess position to use in fit (too large a window weights baseline in fit too much)

<< Statistics`NonlinearFit`;
gaussFcn[t_, tauFWHM_, tCtr_, ampl_] := 
        ampl Exp[-fwhmFactor ((((t - tCtr)/tauFWHM))^2)];
\(\(fwhmFit[efieldTmp_, tLstTmp_, ctrGuess_, tsamp_] := 
        Block[{}, 
normFactor = 
            Max[inten[efieldTmp]]; 
tstLst = 
            Take[Transpose[{tLstTmp/fs - ctrGuess, 
                  inten[efieldTmp]/normFactor}], {npts/2 + 
                  Floor[\((ctrGuess - tsamp)\)/\((dt/fs)\)], 
                npts/2 + 
                  Floor[\((ctrGuess + tsamp)\)/\((dt/
                          fs)\)]}]; 
ans = 
            NonlinearRegress[tstLst, 
              gaussFcn[t, tauTst, t0, amplitude], {t}, {tauTst, t0, 
                amplitude}, 
              RegressionReport -> {BestFitParameters}]; 
Print[ans]; 
parms\  = \ 
            BestFitParameters\  /. 
              ans[\([1]\)]; 
Print["FWHM: ", 
            tauTst /. 
              parms\_\(\([[\)\(1\)\(]]\)\), " fs, center: ", \((t0 + ctrGuess)\) /. 
              parms\_\(\([[\)\(2\)\(]]\)\), " fs"]; 
Return[
            tauTst /. parms[\([1]\)]]
];\)\)

Calculation section

select spatial mode and input pulse shape
spatialMode
1: Gaussian (free space, a = 1/(e^2) intensity radius)
2: Bessel (lowest mode, a = fiber radius)
3: fixed intensity
pulseShape
1: Gaussian (tau = fwhm in intensity)
2: parabolic
Also select here any input dispersion.

spatialMode = 1;
pulseShape  =  1;

    
Ipk =  pkInten[epulse, a, tau, pkFactor, 1];
If[
        spatialMode == 1, 
APk  =  
          A0[Ipk, areaEff], 
Null];
If[
        pulseShape == 
          1, 
(Subscript[AtLst, [[1]]]) = 
          APk gaussT[timeLst, 
              tau], 
Null];

    
Aω(Subscript[Lst, [[1]]]) = FFT[(Subscript[AtLst, [[1]]]), 
          npts, dt];

plot input

iPlt = 1;
Print["time domain: "];
ListPlot[
        Transpose[{timeLst/fs, 
            (Subscript[AtLst, [[iPlt]]])}], PlotJoined -> True, 
        PlotRange -> {{-timeRange, timeRange}, 
            All}];
Print["spectral domain: "];
ListPlot[
        Transpose[{wvlengthLst/nm, 
            Abs[Aω(Subscript[Lst, [[iPlt]]])]}], PlotJoined -> True, 
        PlotRange -> {{lambda0/nm - lamRange, lambda0/nm + lamRange}, 
            All}];

x apply input spectral phase: convention is i(k z- w t), spectral phase is like spatial prop phase and is positive
check to see if code properly handles first and last half steps

Aω(Subscript[Lst, [[2]]]) = Aω(Subscript[Lst, [[1]]]);
(Subscript[AtLst, [[2]]]) = 
        (Subscript[AtLst, [[1]]]);

Main calculation loop

split-step gridding arrays: 
1: transform-limited input
2: pre-chirped input
3: iz = 1, freq * half step spectral->time * full step temporal -> freq
...
nzpts+2: iz = nzteps for temporal phase, nzsteps-1/2 spectral
nzpts+3: freq * half step spectral->time = final output

sphase = 
        Exp[(I βR[omegaLst] - 
                lossesOn*βI[
                    omegaLst]) dz];
sphaseHalf = 
        Exp[(I βR[omegaLst] - 
                lossesOn*βI[omegaLst]) dz/
              2];
iz = 0;
\(\(For[i = 3, 
        i ≤ 
          nzpts + 3, \(i++\), 
\(iz++\); 
 (*\(Print["zstep: ", iz];\)*) 
If[
          iz == 
            1, 
ASpect = \(AωLst\_\(\([[\)\(i - 1\)\(]]\)\)\) 
              sphaseHalf; 
AtLst\_\(\([[\)\(i\)\(]]\)\) = 
            InverseFFT[ASpect, npts, dt]; 
tphase = 
            Exp[I\ γ\ Abs[AtLst\_\(\([[\)\(i\)\(]]\)\)]\^2\ dz]; 
AtLst\_\(\([[\)\(i\)\(]]\)\) = \(AtLst\_\(\([[\)\(i\)\(]]\)\)\) 
              tphase; 
AωLst\_\(\([[\)\(i\)\(]]\)\) = 
            FFT[AtLst\_\(\([[\)\(i\)\(]]\)\), npts, dt], 

\(If[
              iz ≤ 
                nzpts, 
AωLst\_\(\([[\)\(i\)\(]]\)\) = \(AωLst\_\(\([[\)\(i - 1\)\(]]\)\)\) 
                  sphase; 
AtLst\_\(\([[\)\(i\)\(]]\)\) = 
                InverseFFT[
                  AωLst\_\(\([[\)\(i\)\(]]\)\), npts, dt]; 
tphase = 
                Exp[I\ γ\ Abs[AtLst\_\(\([[\)\(i\)\(]]\)\)]\^2\ dz]; 
AtLst\_\(\([[\)\(i\)\(]]\)\) = \(AtLst\_\(\([[\)\(i\)\(]]\)\)\) 
                  tphase; 
AωLst\_\(\([[\)\(i\)\(]]\)\) = 
                FFT[AtLst\_\(\([[\)\(i\)\(]]\)\), npts, 
                  dt], 

AωLst\_\(\([[\)\(i\)\(]]\)\) = \(AωLst\_\(\([[\)\(i - 1\)\(]]\)\)\) 
                  sphaseHalf; 
AtLst\_\(\([[\)\(i\)\(]]\)\) = 
                InverseFFT[
                  AωLst\_\(\([[\)\(i\)\(]]\)\), npts, 
                  dt]
];\)
];
];\)\)

iPlt = nzpts + 3;
Print["time domain: "];
ListPlot[
        Transpose[{timeLst/fs, 
            Abs[(Subscript[AtLst, [[iPlt]]])]}], PlotJoined -> True, 
        PlotRange -> {{-timeRange, timeRange}, 
            All}];
fwhmFit[
      (Subscript[AtLst, [[iPlt]]]), 
      timeLst, 0, 200]
Print["spectral domain: "];
ListPlot[
        Transpose[{wvlengthLst/nm, 
            Abs[Aω(Subscript[Lst, [[iPlt]]])]}], PlotJoined -> True, 
        PlotRange -> {{lambda0/nm - lamRange, lambda0/nm + lamRange}, 
            All}];