% Euler simulering av falldata
% modell 1 utan luftmotstånd
% modell 2 linjärt-
% modell 3 kvadratiskt luftmotstånd.
%--------------------------------------------------------
clear
clf
% Experimentella data
% data är tiden (t) i s för olika fallsträckor (y) i m [t,y]
% OBS ta alltid med punkten 0,0
%
data=[
    0.000	0.000   % 0,0 adderas manuellt - en bra punkt!!
	0.154	0.100
	0.187	0.150
	0.215	0.200
	0.241	0.250
	0.263	0.300
	0.304	0.400
	0.340	0.500
	0.374	0.600
	0.403	0.700
	0.432	0.800
	0.459	0.900
	0.486	1.000
	0.547	1.250
	0.605	1.500
	0.660	1.750
	0.710	2.000
	0.813	2.500
	0.906	3.000
	0.992	3.500
	1.078	4.000
	1.244	5.000
	1.406	6.000
	1.573	7.000
	1.727	8.000
	1.889	9.000];
 %
 y=data(:,2);
 t=data(:,1);
 [N,dum]=size(t);
 %
 figure(1)
 % Data punkter ska ALLTID ritas som "+" (helst med felstaplar men det gör vi inte här)
 plot(t,y,'+')
 axis([-.1 2 -.2 10])

 % Axlarna ska ALLTID betecknas med "Storhetsbeteckning / enhet" där Storhetsbeteckning ska vara kursiv (\it)
 % Kommandot \rm återställer till icke kursiv
 xlabel('\itt\rm / s')
 ylabel('\ity\rm / m')
 title('Fallsträckan som funktion av tiden')

 % 3 momentanhastighet 2.
 hold on
 polyf = polyfit(t,y,4);
 xx = linspace(0,20, 1000);
 val = polyval(polyf, xx);
 plot(xx, val);
 maxhast2 = (polyval(polyf,1.727) - polyval(polyf,1.573))/(1.727-1.573);
 polyderiv = polyder(polyf);
 val2 = polyval(polyderiv, xx);
 plot(xx, val2);


 %
 %
 % fortsättning följer.....
 %
 deriv1 = zeros(26,1);
 deriv2 = zeros(25,1);
 n = 1;
 while n < 26
  deriv1(n,1) = (y(n+1)-y(n))/(t(n+1)-t(n))
  n = n + 1;

 end

 k = 1;
 while k < 25
  deriv2(k,1) = deriv1(k+1)-deriv1(k)
  k = k + 1;
 end
 deriv1
 deriv2
 t1 = t(1:end-1)
 %plot (t, deriv1, '-')
 vmax = 6.17024;
 vmax2 = 6.2
 % Accelerationen skall gå ner till 0!
 b = 9.82/6.2
 B = 0.005*9.82

%% euler
t = 0;
dt = 0.01;
v1 = 0;
v0 = 0;
x0 = 0;
x = 0;
x1 = 0;
g = 9.82;
tt = 1;

dist = zeros(200,1);
hast = zeros(200,1);
acc = zeros(200,1);

while t<2
x1 = x0 + v0*dt;
v1 = v0 + g.*dt;
hast(tt) = v1;
dist(tt) = x1;
acc(tt) = g;
x0 = x1;
v0 = v1;
t = t + dt;
tt = tt + 1;
end
X = linspace(0,2,200);

plot(X,dist, 'g-');
hold on


%% euler 2 v2
t = 0;
dt = 0.01;
v1 = 0;
v0 = 0;
x0 = 0;
x = 0;
x1 = 0;
g = 9.82;
a = g;

tt = 1;

dist = zeros(200,1);
hast = zeros(200,1);
acc = zeros(200,1);

while t<2
x1 = x0 + v0*dt;
v1 = v0 + a*dt;
a = g - g*v1/(6.2);
hast(tt) = v1;
dist(tt) = x1;
acc(tt) = a;

x0 = x1;
v0 = v1;
t = t + dt;
tt = tt + 1;
end
X = linspace(0,2,200);

plot(X,dist, 'r-');
hold on

%% euler 3 v1
t = 0;
dt = 0.01;
v1 = 0;
v0 = 0;
x0 = 0;
x = 0;
x1 = 0;
g = 9.82;
a = g;

tt = 1;

dist = zeros(200,1);
hast = zeros(200,1);
acc = zeros(200,1);

while t<2
x1 = x0 + v0*dt;
v1 = v0 + a*dt;
a = g - g*v1^2;
hast(tt) = v1;
dist(tt) = x1;
acc(tt) = a;

x0 = x1;
v0 = v1;
t = t + dt;
tt = tt + 1;
end
X = linspace(0,2,200);

plot(X,dist, 'b-');
hold on