So, what do I need, is to have an ability to calculate body position and velocity at given time. At first I have body position and velocity What do I can calculate:
Angular momentum $$ \textbf{h} = \textbf{r}\times\textbf{v} $$ Eccentricity: $$ \textbf{e} = \frac{\textbf{v}\times\textbf{h}}\mu - \frac{\textbf{r}}{r} $$ Also I have found that I can calculate position: $$ x=a\left(e-\cosh\tau\right) $$ $$ y=a\sqrt{e^2-1}\sinh\tau $$ And also there some relation between eccentricity anomaly and time: $$ t=\sqrt{a^3\over\mu}\left(e\sinh\tau-\tau\right) $$ But I don't understand, how to calculate that eccentricity anomaly? How to calculate velocity? Is there something like mean anomaly and "mean velocity" (like in ecliptic trajectories)? And is there some formula to calculate it?
Solution
Okay, guys, thank you all, finally I got it. So what conclusion I came to:
If we have some body with velocity and position relative to the body around which first one is going to rotate and time, then we should:
Define orbit
Angular momentum: $$ \bar{h} = \bar{r}\times \bar{v} $$ Eccentricity vector: $$ \bar{e} = \frac{\bar{v}\times \bar{h}}\mu - \frac{\bar{r}}{r} $$ Perigee inclination angle: $$ i=\textrm{atan2}(-\bar{e}_y,-\bar{e}_x) $$ Semi major axis: $$ a=(\frac{2}{r}-\frac{v^2}{\mu})^{-1} $$ If we will rotate position vector by -i, then we will have raw position vector(which not affected by trajectory inclination), this is just inverted 2D vector rotation formula: $$ \bar{p}=\binom{\bar{r}_x\textrm{cos}i+\bar{r}_y\textrm{sin}i}{\bar{r}_y\textrm{cos}i-\bar{r}_x\textrm{sin}i} $$ Then by having raw position and by calculating sinh and cosh we can find creation hyperbolic anomaly: $$ H_0=\textrm{atanh}(\frac{\frac{\bar{p}_y}{a\sqrt{e^2-1}}}{e-\frac{\bar{p}_x}{a}}) $$ And then simply find creation mean anomaly: $$ M_0=\textrm{sinh}(H_0)e-H_0 $$ Also there is a need of mean velocity: $$ n=\sqrt{\frac{\mu}{|a^3|}}\textrm{sign}(h) $$ And for future we need some creation time, which should be gotten from simulation time, when orbit is created: $$ t_0 = T\textrm{, where T is currect simulation time.} $$
Any time get orbit status by time
Calculate current mean anomaly: $$ M=M_0 + (T-t_0)n $$ Find current hyperbolic anomaly via newton's method: $$ H_{j+1}=H_j+\frac{M-e\textrm{sinh}(H_j)+H_j}{e\textrm{cosh}(H_j)-1} $$ And get current raw position vector: $$ \bar{p}=\binom{a(e-\textrm{cosh}H)}{a\sqrt{e^2-1}\textrm{sinh}H} $$ And finally to get body relative position at time T we should rotate raw position vector by inclination angle: $$ \bar{r}=\binom{\bar{p}_x\textrm{cos}i-\bar{p}_y\textrm{sin}i}{\bar{p}_x\textrm{sin}i+\bar{p}_y\textrm{cos}i} $$ Have fun!