%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% maxwell.m:
% Øving 5 FY1005/TFY4165 Termisk fysikk 2014:
% Maxwells hastighetsfordeling for skiver på luftputebord.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nullstiller variable og lukker figurer
clear all;
close all;
% Leser eksperimentelle data fra 49 filer mass8_1.txt, mass8_2.txt osv 
% Hver fil inneholder 3 kolonner: 
%      tid (i sekunder), horisontal posisjon, vertikal posisjon
% Posisjonsverdiene multipliseres med faktoren 0.22 for
% at de skal angi verdier i cm.
% Det er 2 linjer med overskrift, derfor 'headerlines',2
% Heltallet j holder styr på totalt antall hastighetsmålinger (920)
j=0;
% Tidssteget er 1/15 sekund
dt = 1.0/15.0;
xtot = [];
ytot = [];

for ifil = 1:49
   filnavn = ['mass8_',int2str(ifil),'.txt'];
   [t,x,y] = textread(filnavn,'%f %f %f','headerlines',2);
   x = x*0.22;
   xtot = [xtot; x];
   y = y*0.22;
   ytot = [ytot; y];
   n = length(t);
   for i = 1:n-1, % For i=n kan det ikke beregnes forflytning.
      j = j+1;
      dx(j) = x(i+1) - x(i);
      dy(j) = y(i+1) - y(i);      
%%% HER kan det med fordel lages tabeller xpos og ypos med
%%% samtlige x- og y-verdier, med tanke på plotting av skiveposisjoner
   end
end
% Beregner forflytningene (absoluttverdi)
ds = (dx.^2+dy.^2).^(1/2);
% Beregner fart v og hastighetskomponenter vx og vy
v  = ds./dt;
vx = dx./dt;
vy = dy./dt;
%%% HER kan diverse middelverdier beregnes
vmiddel = sum(v)/length(v);
vrms = sqrt(sum(v.^2)./length(v));
vxmiddel = sum(vx)/length(vx);
vxrms = sqrt(sum(vx.^2)./length(vx)); %er det meningen å lage dette?
vymiddel = sum(vy)/length(vy);
vyrms = sqrt(sum(vy.^2)./length(vy)); %er det meningen å lage dette?


%%% og evt <vx>, <vy>
%%% HER kan midlere kinetisk energi beregnes og ekvipartisjonsprinsippet
%%% Antar at 1/2mv^2 fremdeles virker.
%Kmiddel = 0.5*m *vrms^2;

%%% benyttes til å besvare oppgave b).
%%% HER sorteres v, vx og vy med tanke på å lage histogram
v_sortert  = sort(v);
vx_sortert = sort(vx);
vy_sortert = sort(vy);
%%% HER bestemmes bins = antall intervaller; bruk for eksempel
%%% "The Rice Rule"; kommandoen ceil(x) returnerer første heltall
%%% som er større enn argumentet x
bins = ceil(length(v).^(1/3).*2);
%%% HER beregnes intervallbredder; lik hele intervallet fra minimal [På en
%%% måte "binhøydeforskjell", hvor my høydeforskjell det er innad i hver
%%% bin.
%%% til maksimal verdi dividert med antall intervaller (bins)
v_binbredde  = (max(v)-min(v))/bins;
vx_binbredde = (max(vx)-min(vx))/bins;
vy_binbredde = (max(vy)-min(vy))/bins;
%%% HER nullstilles tabeller for eksperimentelle og teoretiske fordelinger
f = zeros(1,bins);
%eksperimentelle:
gx = zeros(1,bins);
gy = zeros(1,bins);
vbin = zeros(1,bins);
vxbin = zeros(1,bins);
vybin = zeros(1,bins);
%%% osv; du trenger tabeller for f, gx, gy (de eksperimentelle) og tilhørende
%%% teoretiske, samt tabeller for vbin, vxbin og vybin

%%% HER lages eksperimentelle histogrammer og tilhørende teoretiske fordelinger
% Teoretisk fartsfordeling: 
%B = 1/(vrms^2); % :)
%fteo = 2.*B*.v*.exp(-B.*v.^2);
% med B = 1/(vrms^2). Siden integralet av fteo fra 0
% til uendelig er lik 1, må fteo multiplisert
% med intervallbredden, summert over de bins=20
% histogramintervallene, bli lik 1.
B = 1/(vrms^2);
k = 1;
ibin = 0;
while (ibin < bins)
    ibin = ibin + 1;
    while (k <= j && v_sortert(k) <= v_binbredde*ibin)
        f(ibin) = f(ibin)+1;
        k=k+1;
    end
    f(ibin) = f(ibin)/j; %j er totalt antall hastighetsmålinger. :)
    vbin(ibin) = v_binbredde*(ibin-0.5);
    fteo(ibin) = 2*B*vbin(ibin)*exp(-B*vbin(ibin)*vbin(ibin))...
                  *v_binbredde;
end

%%% DERETTER lages på tilsvarende vis eksperimentelle histogrammer og
%%% tilhørende teoretiske fordelinger for komponentene vx og vy.
B = 1/(vxrms^2);
k = 1;
ibin = 0;
while (ibin < bins)
    ibin = ibin + 1;
    while (k <= j && vx_sortert(k) <= vx_binbredde*ibin)
        gx(ibin) = gx(ibin)+1;
        k=k+1;
    end
    gx(ibin) = gx(ibin)/j; %j er totalt antall hastighetsmålinger. :)
    vxbin(ibin) = vx_binbredde*(ibin-0.5);
    gxteo(ibin) = 2*B*vxbin(ibin)*exp(-B*vxbin(ibin)*vxbin(ibin))...
                  *vx_binbredde;
end
B = 1/(vyrms^2);
k = 1;
ibin = 0;
while (ibin < bins)
    ibin = ibin + 1;
    while (k <= j && vy_sortert(k) <= vy_binbredde*ibin)
        gy(ibin) = gy(ibin)+1;
        k=k+1;
    end
    gy(ibin) = gy(ibin)/j; %j er totalt antall hastighetsmålinger. :)
    vybin(ibin) = vy_binbredde*(ibin-0.5);
    gyteo(ibin) = 2*B*vybin(ibin)*exp(-B*vybin(ibin)*vybin(ibin))...
                  *vy_binbredde;
end


%%%
%%% HER kan du vurdere å kontrollere at du har fått riktig normering av
%%% de ulike fordelingene: Summen av alle elementene i tabellene f og fteo
%%% skal nå være (omtrent) lik 1. Tilsvarende for gx osv.
%%%
%%% PLOTTING: Her angis essensielt kun bruken av kommandoen bar(x,y). 
%%% Du må selv utvide programmet slik at det totalt sett lager et panel
%%% med 2 x 3 figurer, som vist i øvingsteksten.

subplot(2,3,1),plot(v_sortert);
xlabel('verdinummer');
title('Hastigheter v');
ylabel('v (cm/s)');

subplot(2,3,2,'align'), bar(vbin,f);
hold on;
plot(vbin,fteo,'ro');
xlabel('v (cm/s)');
ylabel('Andel fartsmålinger');
title('Fordeling for fart v');
axis([0 120 0 max(f)*1.1]);

subplot(2,3,3,'align'), bar(vxbin,gx);
hold on;
plot(vxbin,gxteo,'ro');
xlabel('vx (cm/s)');
ylabel('Andel fartsmålinger');
title('Fordeling for fart vx');
axis([0 120 0 max(gx)*1.1]);

subplot(2,3,4,'align'), bar(vbin,gy);
hold on;
plot(vybin,gyteo,'ro');
xlabel('vy (cm/s)');
ylabel('Andel fartsmålinger');
title('Fordeling for fart v');
axis([0 120 0 max(gy)*1.1]);


subplot(2,3,5), scatter(vx, vy);
xlabel('vx (cm/s)');
ylabel('vy (cm/s)');
title('Fartsretninger vx mot vy');


subplot(2,3,6), scatter(xtot, ytot);
xlabel('x (cm)');
title('Posisjoner x mot y');
ylabel('y (cm)');


%b) Boltzmanns plastskivekonstant k_p:
k_p = 0.5 * 0.032*sum((v.*0.01).^(0.5))/length(v)/300  %1/2 * (m i kg) *gjennomsnittet(v^2 i (m/s)^2) / (T i Kelvin)
% = 3.0640*10^-5 = 3.1*10^-5 (kg m^2)/(T s^2)
%c) vmiddel, vrms, vmiddel/vrms = sqrt(pi)/2 ?

vmiddel = sum(v)/length(v); % =36.3132 cm/s
vrms = sqrt(sum(v.^2)./length(v)); % =41.0324 cm/s %rms-hastigheten er større som forventet. :)
a = vmiddel / vrms; % = 0.8850   omtrent det samme :)
b = sqrt(pi)/2; % = 0.8862       omtrent det samme :)

%d)  tegn til drift i eksemlpene under?
%vy går hurtigere mot høyre