GEAkonstanty;
aa=RZ+100000
exc=0.2
tau=0
inc=60
OMEGA=60
omega=30
elementy=[aa,exc,inc,OMEGA,omega,tau,0];
t0=0;
rv0=rvzelementu(t0,elementy,'Ea','m','s',8)
te=86400
pt0=rv0(:,1)
vt0=rv0(:,2)
rtel0=sqrt(dot(pt0,pt0))
vtel0=sqrt(dot(vt0,vt0))
rvt=rvzelementu(te,elementy,'Ea','m','s',8);
ptt=rvt(:,1)
vtt=rvt(:,2)
rtel=sqrt(dot(ptt,ptt))
vtel=sqrt(dot(vtt,vtt))
[t,u]=ode23(@(t,u)satellite3j20(t,u,GMZ,RZ),[t0, te], [pt0(1) vt0(1) pt0(2) vt0(2) pt0(3) vt0(3)],opt); 
m=length(t)
rt45=sqrt(u(m,1)^2+u(m,3)^2+u(m,5)^2)
vt45=sqrt(u(m,2)^2+u(m,4)^2+u(m,6)^2)
deltar=rtel-rt45
deltav=vtel-vt45
for i=1:m
rt=[u(i,1) u(i,3)  u(i,5) ];
vt=[u(i,2) u(i,4)  u(i,6)];
ele=elementyzrv(te,rt,vt,'Ea','m','s');
OO(i)=ele(4);
ww(i)=ele(5);
aaa(i)=ele(1);
eee(i)=ele(2);
iii(i)=ele(3);
endfor
outf=fopen('vystup.dat','wt');
for i=1:m
  fprintf(outf,'%12.3f %15.2f %10.5f %10.4f %10.4f %10.4f\n', t(i),aaa(i),eee(i),iii(i),OO(i),ww(i))
endfor  
%figure;
%u(:,1)
%plot3 (u(:,1),u(:,3),u(:,5));
n=sqrt(GMZ/aa^3);
dotOMEGA=3*n/2*J20/(1-exc^2)^2*(RZ/aa)^2*cos(inc/180*pi);
disp("dot OMEGA =")
dotOMEGA*86400*180/pi