# Computational physics

1. Apr 3, 2014

### eahaidar

I made this code for a non linear process called Four Wave Mixing in which. Two pulses will enter a fiber to get a new wave called idler and try to amplify the signal but I think it's not working well because of the units so I called them in the figure arbitrary units until the problem is fixed any thoughts how? Here is my code thank you

clear all;

T = 500e-12; % Time window (period)
nt = 2^14; % Number of points
dt = T/nt; % Time resolution
% Make w and t windows
t = ((1:nt)-(nt+1)/2)*dt;
t = fftshift(t);

% Create frequency window, then shift it

w = 2*pi*([-nt/2:nt/2-1])/T;
w = fftshift(w);
%Pump pulse probe
P0 = 100;
FWHM = 1e-12;
c = FWHM/(2*sqrt(2*log(2)));
u0 = sqrt(P0)*exp(-0.5*(t/c).^2);
%Signal pulse probe
Ps = 1;
Ts = 5e-12;
detun = 6.05e12;

Pf = Ps*3.042e4;
detun = detun*2*pi;
fbw = 0.44/Ts;
FWHM_omega = 2*pi*fbw;
c = FWHM_omega/(2*sqrt(2*log(2)));
uw = sqrt(2*Pf)*exp(-0.25*((w-detun)/c).^2);
u1 = (ifft(uw));

u = fftshift(u0+u1);

figure
subplot(1,2,1)
plot(w/(1e12)/(2*pi),abs(fft(u)).^2)
xlabel('Frequency (THz)')
ylabel('Arb. Units')
xlim([-10 10])
ut = ifft(uw);
subplot(1,2,2)
plot(t/(1e-12),fftshift(abs(u).^2))
xlabel('Time(ps)')
ylabel('Power (W)')
xlim([-6 6])

%split step method
beta2=0;
beta3=0; alpha=0;
h=0.001;%h = L/step_num; % step size in z
D=(1i.*beta2.*0.5.*w.^2)-((1./6).*1i.*beta3.*w.^3); % dispersion operator
step_num = 1000; % No. of z steps

gamma=1;
temp2=fft(u);
temp3=temp2.*exp(u.*h.*0.5);
temp4=ifft(temp3);
%%%%%%%%%

for m=1:step_num
Nonlinear=exp(1i.*h*gamma.*u.^2).*temp4;
temp5=fft(Nonlinear);
temp6=temp5.*exp(D.*h./2);
temp4=ifft(temp6);
end
temp8=fft(temp4);
temp9=temp8.*exp(-1.*(h./2).*D);
temp10=ifft(temp9);
% final pulse
tempo=fftshift(fft(temp10)).*(nt*dt)/sqrt(2*pi); %Final spectrum

%----Plot output pulse shape and spectrum
figure;
subplot(2,1,1);
plot(t/(1e-12), abs(temp10).^2,'-k'); hold on;
axis([-20 20 0 inf]);
xlabel(' Time');
ylabel(' Power');
title(' Output Pulse Shape and Spectrum');
subplot(2,1,2);
plot(w/(1e12)/(2*pi), abs(tempo).^2, '-k');
hold on;
axis([-10 10 0 inf]);
xlabel('Normalized Frequency');
ylabel('Spectral Power');