%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file shows the propagation of a trapezoidal pulse
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%(1) Calculate the frequency spectrum of the pulse

w=linspace(-10e11,10e11,401);   %calculated frequency range

%trapezoidal specs
a=40e-12;
b=60e-12;

% plot the trapezoidal pulse
figure(1)
t=linspace(-1.1*b,1.1*b,1001);
f=zeros(size(t));
I=find(abs(t)>a & abs(t)<b);
f(I)=(b-abs(t(I)))/(b-a);
J=find(abs(t)<=a);
f(J)=ones(size(J));
subplot(2,1,1)
plot(t*1e12,f)
xlabel('t (ps)')
ylabel('f(t)')

%calculate the frequency spectrum numerically
clear F1
for lp_w=1:length(w)
    ww=w(lp_w);
    tmp=f.*exp(j*ww*t);
    F1(lp_w)=trapz(t,tmp);
end

%the analytic calculation of the frequency spectrum
F2=(cos(a.*w)-cos(b*w))./(pi*w.^2*(b-a));

%plot the frequency spectrum
subplot(2,1,2)
plot(w/(2*pi)*1e-9,F1./max(F1))  %numeric
hold
plot(w/(2*pi)*1e-9,F2./max(F2),'r')  %analytic
hold
ylabel('F(f)')
xlabel('frequency (GHz)')
axis([-40 40 -.2 1])

%(2) Now calculate the pulse shape after propagating
% we will plot four different ranges
figure(2)
clf
D=17e-6;    %dispersion in s/m^2
lambda=1550e-9; %wavelength in m
c=3e8;      %speed of light in m/s
vg=c/1.46;  %group velocity
aa=-lambda^2*D/(4*pi*c);    

%remove the w=0 so that there is not a divide by zero
I=find(w==0);
w=[w(1:I-1),w(I+1:end)];
F2=[F2(1:I-1),F2(I+1:end)];

%the three lengths
LL=[1e3, 10e3, 100e3, 200e3];
for lp=1:length(LL)
    L=LL(lp);    %the propagation distance

    tt=linspace(L/vg-4*b,L/vg+4*b,1001); %the time range
    
    %the numerical inverse Fourier transform
    for lp_t=1:length(tt)
        t=tt(lp_t);
        ff=F2.*exp(j*w.*((t-L/vg)-aa*w*L));
        val1=trapz(w,real(ff));
        val2=trapz(w,imag(ff));
        val(lp_t)=val1+j*val2;    
    end

    %get an equivalent pulse if there was no dispersion
    fo=zeros(size(tt));
    I=find(abs(tt-L/vg)>a & abs(tt-L/vg)<b);
    fo(I)=(b-abs(tt(I)-L/vg))/(b-a);
    J=find(abs(tt-L/vg)<=a);
    fo(J)=ones(size(J));

    %plot the dispersed pulse
    subplot(length(LL),1,lp)
    hold 'on'
    plot((tt-L/vg)*1e9,abs(val)./max(abs(val)))
    plot((tt-L/vg)*1e9,fo./max(fo),'r')
    title(['L= ', num2str(L*1e-3),' km'])
end
hold 'off'
xlabel('time (ps)')