% tfyhjennetään ensiksi muisti kaikesta vanhasta roskasta

clear all;
close all;

% ALUSTUS:

% tarkastellaan aaltoja suljetussa L-kokoisessa laatikossa, määritellään 
% koordinaatisto siten, että laatikon vasen reuna on x=-L/2 ja oikea reuna 
% x=L/2. Jaetaan x-akseli pieniin 0.01-kokoisiin väleihin.
% Määrittelemme siis vektorin (listan) xvals, jonka komponentit ovat vain kunkin
% paikkaavaruuden pisteen x-koordinaatti, eli 
%   xvals = (-L/2,-L/2+0.1,...,L/2).

L=100.0;
xvals=-L/2.0:0.01:L/2.0;

% yllä muuttuja x on nyt itseasiassa vektori, jonka elementit saavat arvot
% -L/2, -L/2+0.01, -L/2+0.02, ..., L/2,
% eli siinä on lueteltu kaikki tarkastelemamme x-akselin pisteet

% määritellään seuraavaksi gaussinen aaltopakettimme

y_x = 1.0/sqrt(2.0*pi)*exp(-0.5*xvals.^2);

% myös y_x on nyt vektori, jonka n:s elementti kertoo funktion arvon 
% vastaavassa x-akselin n:ssä pisteessä.

% piirretään vielä kuvaaja

figure % luodaan uusi kuvaaja-ikkuna
plot(xvals,y_x) % piirretään ikkunaan y_x:n arvot xvals:in funktiona
title('Alkuperäinen aaltopaketti paikka-avaruudessa') % annetaan kuvaajalle nimi
xlabel('Paikka x') % annetaan x-akselille nimi
ylabel('Aallon korkeus y(x)') % annetaan y-akselille nimi
text(10,0.2,"std = "+std(xvals,y_x));
% kuvaajasta nähdään siisti gaussinen aaltopaketti (kellokäyrä), jonka
% keskipiste on nollassa (x=0)
pause;


% FOURIER-MUUNNOS

% seuraavaksi tehdään aaltopaketille Fourier-muunnos, tai tarkemmin sanoen
% cos-muunnos. Koska aaltopaketti on parillinen (sama muoto ja merkki
% x:n nollakohdan molemmilla puolilla), katoavat kaikki sin-funktiot
% Fourier muunnoksesta.

% määritellään ensiksi Fourier muunnoksessa laskemamme n-indeksin arvot. 
% summaushan menee periaatteessa äärettömyyteen asti, mutta yhä korkeampia 
% n-indeksejä vastaavat kertoimet \beta_n pienenevät nopeasti. Käytetään
% sataa Fourier-termiä (plus vakiotermi, eli n=0, sillä cos(2*pi/L*0*x) = 1). 
% Määrittelemme siis vektorin nvals, jossa on vain lueteltuna nuo
% n-indeksit, eli nvals = (0,1,2,...,100).

nvals=0:1:100;

% sitten tehdään varsinainen Fourier muunnos, eli jokaista n-indeksin arvoa
% nval kohden lasketaan integraali 
%      \beta_nval = (2/L)*\int_{-L/2}^{L/2} f(x)*cos(2*pi/L*nval*x) dx
%
% integrointi toteutetaan trapz(x,y(x))-funktiolla, joka laskee y(x):n 
% integraalin yli muuttujan x trapezoidi-menetelmällä. Kyseessä on määrätty
% integraali, jossa integroimisalueena on x-muuttujan määrittelyväli.

% tehdään silmukka joka juoksee läpi kaikkien n-indeksien arvojen yli

beta = zeros(1,length(nvals)); %alustetaan sopivan kokoinen vektori, sisällöllä ei väliä

for nindex=1:length(nvals)
    nval = nvals(nindex);
    beta(nindex) = (2.0/L)*trapz(xvals,y_x.*cos(2.0*pi/L*nval*xvals));
end


% piirretään saatu Fourier-muunnos

plot(nvals,beta,'o')
xlabel('n-indeksi')
ylabel('Fourier kerroin beta_n')
title('Alkuperäisen aaltopaketin Fourier-kertoimet')

% sekin on gaussinen! Huomaa, että x-akselina on nyt n-indeksit, ja
% kuvattuna on itseasiassa Fourier-kertoimet \beta_n

pause;

% edellä Fourier kertoimet on esitetty indeksin n funktiona. Mutta indeksillä on fysikaalinen merkityksensäkin, 
% sillä kutakin arvoa n vastaa harmonisen aallon aaltovektori k_n = 2*pi/L*n. Esitetään nyt aaltopaketti
% paikka-avaruudessa sekä sen Fourier muunnos aaltovektorin k_n funktiona.

kvals = 2.0*pi/L*nvals;

tiledlayout(2,1)
nexttile
plot(xvals,y_x) % piirretään ikkunaan y_x:n arvot xvals:in funktiona
title('Alkuperäinen aaltopaketti paikka-avaruudessa') % annetaan kuvaajalle nimi
xlabel('Paikka x') % annetaan x-akselille nimi
ylabel('Aallon korkeus y(x)') % annetaan y-akselille nimi
text(10,0.2,"std = "+std(xvals,y_x));

nexttile
plot(kvals,beta,'o')
xlabel('n-indeksi')
ylabel('Fourier kerroin beta_n')
title('Alkuperäisen aaltopaketin Fourier-kertoimet')
text(5,max(beta)/2.0,"std = "+std(kvals,abs(beta)));

pause;

for leveys=0.5:0.5:5.0
    
    y_x = 1.0/sqrt(2.0*pi)*exp(-0.5*(xvals/leveys).^2);

    for nindex=1:length(nvals)
        nval = nvals(nindex);
        beta(nindex) = (2.0/L)*trapz(xvals,y_x.*cos(2.0*pi/L*nval*xvals));
    end

    tiledlayout(2,1)
    nexttile
    plot(xvals,y_x) % piirretään ikkunaan y_x:n arvot xvals:in funktiona
    title('Aaltopaketti paikka-avaruudessa') % annetaan kuvaajalle nimi
    xlabel('Paikka x') % annetaan x-akselille nimi
    ylabel('Aallon korkeus y(x)') % annetaan y-akselille nimi
    text(5,0.2,"std = "+std(xvals,y_x));
    
    nexttile
    plot(kvals,beta,'o')
    xlabel('n-indeksi')
    ylabel('Fourier kerroin beta_n')
    title('Aaltopaketin Fourier-kertoimet')
    text(5,max(beta)/2.0,"std = "+std(kvals,abs(beta)));
    
    pause
end