%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% losning11.m:
% Oving 11 TFY4115 Fysikk 2013:
% Maxwells hastighetsfordeling for skiver paa 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 maa multipliseres med faktoren 0.22 for
% at de skal angi enheten cm.
% Det er 2 linjer med overskrift, derfor 'headerlines',2
% Heltallet j holder styr paa totalt antall hastighetsmaalinger
j=0;
% Tidssteget er 1/15 sekund
dt = 1.0/15.0;
for ifil = 1:49
   filnavn = ['mass8_',int2str(ifil),'.txt'];
   [t,x,y] = textread(filnavn,'%f %f %f','headerlines',2);
   x = x*0.22;
   y = y*0.22;
   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);
      % Lager tabeller xpos og ypos med tanke paa plotting
      % av skiveposisjoner
      xpos(j) = x(i);
      ypos(j) = y(i);
   end
end
% Beregner forflytningene (absoluttverdi)
ds = sqrt(dx.^2+dy.^2);
disp('Totalt antall beregnede hastigheter');
disp(j);
% Beregner fart v og hastighetskomponenter vx og vy
v  = ds/dt;
vx = dx/dt;
vy = dy/dt;
middelfart = sum(v)/j;
vrms = sqrt(sum(v.^2)/j);
kinetiskenergi = 0.5*0.032*(vrms/100)^2;
temperatur = kinetiskenergi/1.38e-23;
%Sorterer v, vx og vy med tanke paa aa lage histogram
v_sortert = sort(v);
vx_sortert = sort(vx);
vy_sortert = sort(vy);
% binantall = antall intervaller for v
binantall = 21;
% binkomp = antall intervaller for vx og vy
% velger her binkomp = binantall
binkomp = binantall;
% beregner intervallbredder
v_binbredde = max(v)/binantall;
vx_binbredde = 2*max(abs(max(vx)),abs(min(vx)))/binkomp;
vy_binbredde = 2*max(abs(max(vy)),abs(min(vy)))/binkomp;
% nullstiller eksperimentelle og teoretiske fordelinger
f = zeros(1,binantall);
gx = zeros(1,binkomp);
gy = zeros(1,binkomp);
f_teo = zeros(1,binantall);
gx_teo = zeros(1,binkomp);
gy_teo = zeros(1,binkomp);
vbin = zeros(1,binantall);
vxbin = zeros(1,binkomp);
vybin = zeros(1,binkomp);
% Teoretisk fartsfordeling: f_teo = Ateo*v*exp(-Bteo*v^2)
% Velger Bteo = 1/(vrms^2) og tilpasser deretter Ateo
% slik at f_teo summert blir lik f summert, dvs 920
Ateo = 5.498525;
Bteo = 1/(vrms^2);
k = 1;
ibin = 1;
while (ibin <= binantall)
    while (k <= j && v_sortert(k) <= v_binbredde*ibin)
        f(ibin) = f(ibin)+1;
        k=k+1;
    end
    vbin(ibin) = v_binbredde*(ibin-0.5);
    f_teo(ibin) = Ateo*vbin(ibin)*exp(-Bteo*vbin(ibin)*vbin(ibin));
    ibin=ibin+1;
end
% Teoretisk komponentfordeling, gx_teo = Axteo*exp(-Bxteo*vx^2), 
% og tilsvarende for gy. Bruker samme verdi for B i alle fordelinger.
% Tilpasser Axteo og Ayteo slik at gx_teo og gy_teo summert blir lik 920
Axteo = 124.770706;
Bxteo = Bteo;
k = 1;
ibin = 1;
while (ibin <= binkomp)
    while (k <= j && vx_sortert(k) <= min(vx)+vx_binbredde*ibin)
        gx(ibin) = gx(ibin)+1;
        k=k+1;
    end
    vxbin(ibin) = min(vx)+vx_binbredde*(ibin-0.5);
    gx_teo(ibin) = Axteo*exp(-Bxteo*vxbin(ibin)*vxbin(ibin));
    ibin=ibin+1;
end
Ayteo = 102.608471;
Byteo = Bteo;
k = 1;
ibin = 1;
while (ibin <= binkomp)
    while (k <= j && vy_sortert(k) <= min(vy)+vy_binbredde*ibin)
        gy(ibin) = gy(ibin)+1;
        k=k+1;
    end
    vybin(ibin) = min(vy)+vy_binbredde*(ibin-0.5);
    gy_teo(ibin) = Ayteo*exp(-Byteo*vybin(ibin)*vybin(ibin));
    ibin=ibin+1;
end
% Sjekker avvik mellom maalte og teoretiske fordelinger
f_avvik = sum(sqrt((f-f_teo).^2));
gx_avvik = sum(sqrt((gx - gx_teo).^2));
gy_avvik = sum(sqrt((gy - gy_teo).^2));
fsum = sum(f);
fteosum = sum(f_teo);
gxsum = sum(gx);
gxteosum = sum(gx_teo);
gysum = sum(gy);
gyteosum = sum(gy_teo);
% Plotting. Lager ett panel med 2 x 3 subplott.
figure;
subplot(2,3,1);
hold on;
bar(vbin,f);
plot(vbin,f_teo,'ro');
xlabel('v (cm/s)');
ylabel('Antall fartsmålinger');
title('Fordeling v');
axis([0 120 0 120]);
subplot(2,3,2);
hold on;
bar(vxbin,gx);
plot(vxbin,gx_teo,'ro');
xlabel('vx (cm/s)');
ylabel('Antall fartsmålinger');
title('Fordeling vx');
axis([-120 120 0 150]);
subplot(2,3,3);
hold on;
bar(vybin,gy);
plot(vybin,gy_teo,'ro');
xlabel('vy (cm/s)');
ylabel('Antall fartsmålinger');
title('Fordeling vy');
axis([-120 120 0 150]);
subplot(2,3,4);
plot(v);
xlabel('Nr j');
ylabel('v (cm/s)');
title('Skivenes fart');
subplot(2,3,5);
plot(vx,vy,'bo');
xlabel('vx (cm/s)');
ylabel('vy (cm/s)');
title('Skivenes hastigheter');
axis([-120 120 -120 120]);
subplot(2,3,6);
plot(xpos,ypos,'bx');
xlabel('x (cm)');
ylabel('y (cm)');
title('Skivenes posisjoner');
