function golfovymicek()

%modelovani pohybu golfoveho micku s uvazenim rotace a odporu vzduchu

clear; clc;

global m omega g w rho S Koefodp dir;

%parametry micku
m=0.1; %hmotnost micku
g=9.81; %tihove zrychleni
d=0.05; %prumer micku
S=pi*d^2/4; %prurezova plocha micku
rho=1.29;   %hustota vzduchu 
Koefodp=S*rho/(2*m); %odporovy koeficient

%otaceni micku
dir=1; %topspin / dir=1
frekv=50; %frekvence otaceni micku
omega=2*pi*frekv;
w=d/2*omega;


options = odeset('RelTol',1e-4,'AbsTol',1e-5);

%pocatecni podminky
t0=0; tmax=10;

v0=60; %velikost pocatecni rychlosti micku
fi=pi/4; %uhel vypaleni micku

x0=0; y0=0; vx0=v0*cos(fi); vy0=v0*sin(fi);

%s odporem vzduchu a rotaci micku
[T,Y] = ode45(@vrhfunkce,[t0 tmax],[x0 vx0 y0 vy0],options);


%bez rotace micku
d=0; omega=0; S=pi*d^2/4;
w=d/2*omega;
[T2,Y2] = ode45(@vrhfunkce,[t0 tmax],[x0 vx0 y0 vy0],options);


%bez odporu vzduchu a rotace micku
d=0; omega=0; S=pi*d^2/4;
w=d/2*omega;
Koefodp=0;
[T3,Y3] = ode45(@vrhfunkce,[t0 tmax],[x0 vx0 y0 vy0],options);

figure(1)
plot(Y(:,1),Y(:,3),'r')
hold on
plot(Y2(:,1),Y2(:,3),'b')
hold on
plot(Y3(:,1),Y3(:,3),'g')
hold off
ylim([0 1.05*max(Y3(:,3))])
legend('s odporem a rotaci','s odporem, bez rotace','bez odporu a rotace');
xlabel('X [m]');
ylabel('Y [m]');




function dy = vrhfunkce(t,y)

global m omega g w rho S Koefodp dir;

dy = zeros(4,1);  

v=sqrt(y(2).^2+y(4).^2);

Cd=(0.508+w./(22.503*(w).^2.5+4.196*(v).^2.5).^0.4).*Koefodp;
Cm=dir*w./(2.022*w+0.981*v).*Koefodp;

dy(1) =y(2);
dy(2) = -Cd*v*y(2)+Cm*v*y(4);
dy(3) = y(4);
dy(4) = -g-Cd*v*y(4)-Cm*v*y(2);






