clear all
tau=5e-11;
N=3;
NN=N+1;
for nn=1:NN;
    n=nn-1;
    E(nn)=1;
    k(nn)=sqrt(E(nn)/tau/factorial(n)/sqrt(2*pi));
end;
it=0;
for t=-5e-10:2.5e-12:5e-10;
    it=it+1;
    for nn=1:NN;
        n=nn-1;
        hh=0;
        for i=0:fix(n/2);
            hh=hh+(-.5)^i*(t/tau)^(n-2*i)/factorial(n-2*i)/factorial(i);
        end;
        h(nn,it)=k(nn)*exp(-t^2/4/tau^2)*factorial(n)*hh;
    end,
end;
figure (1);t=-5e-10:2.5e-12:5e-10;
for nn=1:NN;
    if nn==1;co='r';end;if nn==2;co='g';end;if nn==3;co='b';end;if nn==4;co='c';end;if nn==5;co='k';end;if nn==6;co='m';end;if nn==7;co='y';end;
    plot(t,h(nn,:),co);hold on;
end;
legend('n=0','n=1','n=2','n=3');
xlabel('Time, sec');ylabel('Amplitude');grid;
% e_1=sum(h(1,:).^2)*2.5e-12,e_2=sum(h(2,:).^2)*2.5e-12,e_3=sum(h(3,:).^2)*2.5e-12,e_4=sum(h(4,:).^2)*2.5e-12,%test for energy
 
figure (2);f=1e8:1e8:4e11;
for nn=1:NN;
    if nn==1;co='r';end;if nn==2;co='g';end;if nn==3;co='b';end;if nn==4;co='c';end;if nn==5;co='k';end;if nn==6;co='m';end;if nn==7;co='y';end;
    H=fft(h(nn,:),10*(it-1))/4e11;
    plot(f,abs(H),co);hold on;
end;
legend('n=0','n=1','n=2','n=3');
xlabel('Frequency, Hz');ylabel('Magnitude');grid;
axis([0 1.2e10 0 1.6e-5]);
 
iif=0;
for f=0:2e6:1.2e10;
    iif=iif+1;
    H_0(iif)=k(1)*tau*2*sqrt(pi)*exp(-(2*pi*tau*f)^2);
    H_1(iif)=k(2)*tau*(-j*4*pi*f*tau)*2*sqrt(pi)*exp(-(2*pi*tau*f)^2);
    H_2(iif)=k(3)*tau*(1-16*pi^2*f^2*tau^2)*2*sqrt(pi)*exp(-(2*pi*tau*f)^2);
    H_3(iif)=k(4)*tau*(-j*12*pi*f*tau+j*64*pi^3*f^3*tau^3)*2*sqrt(pi)*exp(-(2*pi*tau*f)^2);
end,
figure (3);f=0:2e6:1.2e10;plot(f,abs(H_0),'r');hold on;plot(f,abs(H_1),'g');hold on;plot(f,abs(H_2),'b');hold on;plot(f,abs(H_3),'c');
legend('n=0','n=1','n=2','n=3');
xlabel('Frequency, Hz');ylabel('Magnitude');grid;
% ef_0=sum(abs(H_0).^2)*2e6,ef_1=sum(abs(H_1).^2)*2e6,ef_2=sum(abs(H_2).^2)*2e6,ef_3=sum(abs(H_3).^2)*2e6    %test for energy (.5)
% of_01=sum(h(1,:).*h(2,:))*2.5e-12, of_02=sum(h(1,:).*h(3,:))*2.5e-12, of_12=sum(h(2,:).*h(3,:))*2.5e-12, %test for orthogonality (0)