B=10e9;
alpha=8*B^2;
l=1550;
lo=1310;
So=0.092;
DD=So/4*(l-lo^4/l^3);
D=DD*1e-6;
lambda=l*1e-9;

a=D*lambda^2/(4*pi*3e8)
vg=3e8/1.47;
%wo=7e11;
wo=1e11;

t=linspace(-1/B, 2/B, 1001);
E1=exp(-alpha*t.^2);
E2=exp(j*wo*t);

figure(1)
subplot(2,1,1)
plot(t*1e9,real(E1.*E2))
hold
plot(t*1e9,E1,':')
plot(t*1e9,-E1,':')
hold
title('Single Pulse (\omega set lower to see oscillations)')
ylabel('E')

E3=exp(-alpha*(t-1/B).^2);
E4=exp(j*wo*(t-1/B));

subplot(2,1,2)
plot(t*1e9,E1)
hold
plot(t*1e9,E3,'r')
plot(t*1e9,E1+E3,'k')
hold
xlabel('t (ns)')
ylabel('|E|')
title('Pulse envelope of two adjacent pulses')


figure(2)
z=1000e3;
%z=20e3;
t2=linspace(-2.2e-9+z/vg,2.2e-9+z/vg, 10001)

tmp2=exp(-(t2-z/vg).^2./(1/alpha+16*a^2*z^2*alpha));
tmp2a=exp(-(t2-z/vg-1/B).^2./(1/alpha+16*a^2*z^2*alpha));
tmp2b=exp(-(t2-z/vg+1/B).^2./(1/alpha+16*a^2*z^2*alpha));

phi=wo*t2+(4*a*z./(1/alpha^2+16*a^2*z^2)).*(t2-z/vg).^2;

E2=tmp2;
E2a=tmp2a;
E2b=tmp2b;

subplot(3,1,1)
plot((t2-z/vg)*1e9,real(E2))
hold
plot((t2-z/vg)*1e9,real(E2a),'r')
plot((t2-z/vg)*1e9,real(E2b),'r')
plot((t2-z/vg)*1e9,real(E2a+E2+E2b),'k')
hold
ylabel('|E|')
axis([min(t2-z/vg)*1e9, max(t2-z/vg)*1e9, -.1 3.1])
title('Length increased to show frequency chirp')

to=2*sqrt(log(2)/alpha);
T=to*sqrt(1+(4*log(2)*D*z*1550e-9^2/(4*3e8*to^2))^2)

ww=wo+2*(4*a*z./(1/alpha^2+16*a^2*z^2)).*(t2-z/vg);

subplot(3,1,2)
plot((t2-z/vg)*1e9, ww/wo)
ylabel('f/f_o')
axis([min(t2-z/vg)*1e9, max(t2-z/vg)*1e9, min(ww/wo) max(ww/wo)])

subplot(3,1,3)
plot((t2-z/vg)*1e9, real(tmp2.*exp(j*phi)))
xlabel('t (ns)')
ylabel('E(t)')
axis([min(t2-z/vg)*1e9, max(t2-z/vg)*1e9, -1, 1])