On Wed, Mar 23, 2016 at 11:35 AM, Emeline Lépine
<[email protected]> wrote:
> Hello,
>
> For my students I have to create a code resolving non-linear Schrodinger
> Equation with Julia.
> This is is one :
>
> #split step fourier method
> #16/03/2016
>
> using Gadfly
>
> ####################################################################
>
> cputime=0
> tic();
Don't do these in global scope.
> ln=1;
> Po=.00064; #input power in watts
> L = 12 # lenght of the fiber
> Α=0;# fiber loss value in dB/km
> α =Α/4.343;#ref page %55 GP Agrawal
> γ=0.003;#fiber non linearity in /W/m
> to=125e-12; #initial pulse width in second
> C=-2; # Input chirp parameter for first calculation
> b2=20e-27 # 2nd order disp. (s2/m)
> Ld=(to^2)/(abs(b2)); # dispersion length in meter
> Ao=sqrt(Po); # Amplitude
> i = complex(1, pi/2);
> start = 0.1;
> step =0.1;
> stop = 1.5;
> nb_points = (stop-start)/step;
> iim = linspace(start, stop, nb_points);
If you just paste the code into the repl you'll see the issue is that
nb_points is not a integer. Round it.
Even if it is not in the repl, the linspace should appear in the backtrace.
> ###################################################################
>
> τ=- 4096e-12:1e-12:4095e-12; #dt=t/to
> dt=1e-12;
> rel_error=1e-5;
> h=1000;#step size
>
> u=Ao*exp(-((1+im.*(-C))/2).*(τ/to).^2);
> plot(x=τ,y=abs(u),Geom.line,Guide.xlabel("Time"), Guide.ylabel("Amplitude"),
> Guide.title("Input pulse"));
>
> for k = 1:nb_points # the various fiber lengths can be varied and this
> vector can be changed
> ii = iim[k];
> l=maximum(size(u));
> fwhm1=find(abs(u).>(maxabs(u)./2));
> fwhm1=length(fwhm1);
> dw = 1/1/dt*2*pi;
> w = (-1*1/2:1:1/2-1)*dw;
> u = fftshift(u);
> w = fftshift(w);
> spectrum = fft(fftshift(u)); #Pulse spectrum
> spectrum = spectrum.*exp(-α*(h/2) + i*b2/2*w.^2*(h/2));
> f = ifft(spectrum);
> f = f.*exp(i*γ*((abs(f).^2)*h));
> spectrum = fft(f);
> spectrum = spectrum.*exp(-α*(h/2)+i*b2*w.^2*(h/2));
> f = ifft(spectrum);
> #op_pulse(ln)= abs(f);
> fwmh = find((abs(f)).>(abs(maxabs(f)./2)));
> fwmh = length(fwmh);
> ratio = fwmh./fwhm1;
> pbratio[k]= ratio;
> dd = atand((abs(imag(f)))/(abs(f)));
> phadisp[k] = dd;
> end
> toc;
This does **NOT** implicitly call the function, this is **NOT** matlab.
Also, don't put `;` at the end of every line.
> cputime = toc;
> plot(x=ln,y=pbratio,Geom.line,Guide.xlabel("Number of steps"),
> Guide.ylabel("Pulse Broadening ratio"));
> plot(x=ln,y=phadisp,Geom.line,Guide.xlabel("Distance travelled"),
> Guide.ylabel("Phase change"));
> draw(PDF("split_step_fourier_method.pdf", 6inch, 6inch), vstack(p1,p2,p3));
>
>
> My problem is a conversion but I don't know where ?
>
> Could you help me ?
>
> Your sincerly,
>
> EL
