From 9b43dd5cfda00e61c7969e1cc44ee9f679fceedf Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Tue, 15 Nov 2022 15:57:29 +0100 Subject: [PATCH 01/19] replace old inline math to suit with new Github markdown style --- .../Component/AOCS/Spec_RWJitter.md | 32 ++++----- Specifications/Component/AOCS/Spec_STT.md | 12 ++-- .../Component/Abstract/Spec_ComponentBase.md | 4 +- .../Component/Mission/Spec_Telescope_en.md | 18 ++--- .../Component/Mission/Spec_Telescope_ja.md | 14 ++-- .../Propulsion/Spec_SimpleThruster.md | 8 +-- Specifications/Disturbance/Spec_GGTorque.md | 4 +- .../Disturbance/Spec_GeoPotential.md | 36 +++++----- .../Disturbance/Spec_MagDisturbance.md | 4 +- .../Disturbance/Spec_SurfaceForce.md | 12 ++-- .../Disturbance/Spec_SurfaceForce_AirDrag.md | 46 ++++++------ .../Spec_SurfaceForce_SolarRadiation.md | 20 +++--- .../Disturbance/Spec_ThirdBodyGravity.md | 6 +- .../Dynamics/Spec_AttitudeDynamics.md | 10 +-- .../Dynamics/Spec_ControlledAttitude.md | 24 +++---- Specifications/Dynamics/Spec_EnckeMethod.md | 22 +++--- Specifications/Dynamics/Spec_KeplerOrbit.md | 72 +++++++++---------- Specifications/Dynamics/Spec_RelativeOrbit.md | 34 ++++----- Specifications/Dynamics/Spec_Rk4Orbit.md | 10 +-- Specifications/Environment/Spec_Atmosphere.md | 2 +- .../Environment/Spec_CelestialRotation.md | 38 +++++----- .../Environment/Spec_MagEnvironment.md | 4 +- .../Environment/Spec_SRPEnvironment.md | 48 ++++++------- 23 files changed, 240 insertions(+), 240 deletions(-) diff --git a/Specifications/Component/AOCS/Spec_RWJitter.md b/Specifications/Component/AOCS/Spec_RWJitter.md index 229deeb8..d4461d6f 100644 --- a/Specifications/Component/AOCS/Spec_RWJitter.md +++ b/Specifications/Component/AOCS/Spec_RWJitter.md @@ -22,15 +22,15 @@ 3. how to use - Set the harmonics coefficients in `radial_force_harmonics_coef.csv` and `radial_torque_harmonics_coef.csv` - - The first column is an array of the $`h_i`$($`i`$-th harmonic number). The second column is an array of the $`C_i`$(amplitude of the $`i`$-th harmonic). + - The first column is an array of the $h_i$($i$-th harmonic number). The second column is an array of the $C_i$(amplitude of the $i$-th harmonic). - Set parameters in `RW.ini` - - When only the static imbalance and dynamic imbalance(correspond to $`C_i`$ at $`h_i`$≒1) is known according to the spec sheet, edit the files as follows. + - When only the static imbalance and dynamic imbalance(correspond to $C_i$ at $h_i$≒1) is known according to the spec sheet, edit the files as follows. + `radial_force_harmonics_coef.csv` - * Set $`h_1`$(the line 1 of the first column) as $`1.0`$. - * Set $`C_1`$(the line 1 of the second column) as the static imbalance on the spec sheet. + * Set $h_1$(the line 1 of the first column) as $1.0$. + * Set $C_1$(the line 1 of the second column) as the static imbalance on the spec sheet. + `radial_torque_harmonics_coef.csv` - * Set $`h_1`$(the line 1 of the first column) as $`1.0`$. - * Set $`C_1`$(the line 1 of the second column) as the dynamic imbalance on the spec sheet. + * Set $h_1$(the line 1 of the first column) as $1.0$. + * Set $C_1$(the line 1 of the second column) as the dynamic imbalance on the spec sheet. + `RW.ini` * Set `harmonics_degree = 1`. - Set the jitter update period to an appropriate value. @@ -57,8 +57,8 @@ u(t)=\sum_{i=1}^n C_i\Omega^2\sin(2\pi h_i\Omega t+\alpha_i) ``` - - where $`u(t)`$ is the disturbance force and torque in Newton (N) or Newton-meters (Nm), $`n`$ is the number of harmonics included in the model, $`C_i`$ is the amplitude of the $`i`$th harmonic in $`\mathrm{N^2/Hz}`$ (or $`\mathrm{(Nm)^2/Hz}`$), $`\Omega`$ is the wheel speed in Hz, $`h_i`$ is the $`i`$th harmonic number and $`\alpha_i`$ is a random phase (assumed to be uniform over $`[0, 2\pi]`$) [1]. - - $`\alpha_i`$ is generated as a uniform random number in the constructor. + - where $u(t)$ is the disturbance force and torque in Newton (N) or Newton-meters (Nm), $n$ is the number of harmonics included in the model, $C_i$ is the amplitude of the $i$th harmonic in $\mathrm{N^2/Hz}$ (or $\mathrm{(Nm)^2/Hz}$), $\Omega$ is the wheel speed in Hz, $h_i$ is the $i$th harmonic number and $\alpha_i$ is a random phase (assumed to be uniform over $[0, 2\pi]$) [1]. + - $\alpha_i$ is generated as a uniform random number in the constructor. - When users want to use a more precise model, set `considers_structural_resonance` to ENABLE in `RW.ini` and use a model that takes structural resonance inside the RW into account. + If structural resonances are not taken into account, the RW disturbance will be underestimated, but it is not a significant change in general. + See the description of `AddStructuralResonance()` for the algorithm to calculate the structural resonance. @@ -72,15 +72,15 @@ - output: + jitter force and torque with structural resonance in component frame 3. algorithm - - The transfer function from disturbance by harmonics of RW without resonance ($`u(t)`$) to disturbance with resonance ($`y(t)`$) is modeled as following equation: + - The transfer function from disturbance by harmonics of RW without resonance ($u(t)$) to disturbance with resonance ($y(t)$) is modeled as following equation: ```math G(s)=\frac{s^2+2\zeta\omega_ns+\omega_n^2}{s^2+2d\zeta\omega_ns+\omega_n^2} ``` ```math Y(s)=G(s)U(s) ``` - - where $`\omega_n`$ is the angular frequency on the structural resonance. Other parameters such as $`\zeta`$, $`d`$ are determined by the result of experiments. - - To perform the simulation in discrete time, A bi-linear transformation $`G(s)\rightarrow H(z)`$ is applied. $`T`$ is the jitter update period. + - where $\omega_n$ is the angular frequency on the structural resonance. Other parameters such as $\zeta$, $d$ are determined by the result of experiments. + - To perform the simulation in discrete time, A bi-linear transformation $G(s)\rightarrow H(z)$ is applied. $T$ is the jitter update period. ```math \begin{aligned} @@ -91,13 +91,13 @@ \end{aligned} ``` - - The $`\omega_n`$ should be the fixed value by pre-warping because there is frequency distortion due to bilinear transformation. The formula for calculating $`\omega_n`$ for the true resonant frequency $`\omega_d`$ is as follows: + - The $\omega_n$ should be the fixed value by pre-warping because there is frequency distortion due to bilinear transformation. The formula for calculating $\omega_n$ for the true resonant frequency $\omega_d$ is as follows: ```math \omega_n=\frac{2}{T}\tan(\frac{T\omega_d}{2}) ``` - - The bi-linear transformation transforms the relationship between input $`u`$ and output $`y`$ as follows: + - The bi-linear transformation transforms the relationship between input $u$ and output $y$ as follows: ```math \begin{aligned} @@ -106,12 +106,12 @@ \end{aligned} ``` - - By applying the inverse z-transform, the continuous relationship between $`y(t)`$ and $`u(t)`$ can be expressed as a discrete relationship of a difference equation between $`y[n]`$ and $`u[n]`$, where $`[n]`$ is the current simulation time step. The difference equation is as follows: + - By applying the inverse z-transform, the continuous relationship between $y(t)$ and $u(t)$ can be expressed as a discrete relationship of a difference equation between $y[n]$ and $u[n]$, where $[n]$ is the current simulation time step. The difference equation is as follows: ```math c_0y[n]+c_1y[n-1]+c_2y[n-2]=c_3u[n]+c_4u[n-1]+c_5u[n-2] ``` - - Therefore, $`y[n]`$ is calculated as follows. + - Therefore, $y[n]$ is calculated as follows. ```math y[n]=\frac{(-c_1y[n-1]-c_2y[n-2]+c_3u[n]+c_4u[n-1]+c_5u[n-2])}{c_0} ``` @@ -166,7 +166,7 @@ Simulated RW jitter in time domain - - Both the first-order mode and the structural resonance ($`\omega_n=585\mathrm{Hz}`$) are approximately simulated. + - Both the first-order mode and the structural resonance ($\omega_n=585\mathrm{Hz}$) are approximately simulated. ## 4. References 1. Masterson, R. A. (1999). Development and validation of empirical and analytical reaction wheel disturbance models (Doctoral dissertation, Massachusetts Institute of Technology). diff --git a/Specifications/Component/AOCS/Spec_STT.md b/Specifications/Component/AOCS/Spec_STT.md index 7b18a712..06c1d647 100644 --- a/Specifications/Component/AOCS/Spec_STT.md +++ b/Specifications/Component/AOCS/Spec_STT.md @@ -21,9 +21,9 @@ - TBW 2. `Judgement` 1. `EarthJudgement` - - Calculate the angle $`\theta_{ce}`$ between the earth's center direction $`\boldsymbol{r_{sc}}`$ and the earth's edge direction $`\boldsymbol{r_{se}}`$. $`R_e`$ is the earth's radius. - - Calculate the angle $`\theta_{es}`$ between the sight direction $`\boldsymbol{r_{sight}}`$ and the earth's edge direction $`\boldsymbol{r_{se}}`$. - - Judge the STT error flag by comparing $`\theta_{es}`$ with the earth forbidden angle $`\theta_{efa}`$. If $`\theta_{es} > \theta_{efa}`$, the earth is completely outside the earth forbidden angle. + - Calculate the angle $\theta_{ce}$ between the earth's center direction $\boldsymbol{r_{sc}}$ and the earth's edge direction $\boldsymbol{r_{se}}$. $R_e$ is the earth's radius. + - Calculate the angle $\theta_{es}$ between the sight direction $\boldsymbol{r_{sight}}$ and the earth's edge direction $\boldsymbol{r_{se}}$. + - Judge the STT error flag by comparing $\theta_{es}$ with the earth forbidden angle $\theta_{efa}$. If $\theta_{es} > \theta_{efa}$, the earth is completely outside the earth forbidden angle. ```math @@ -58,15 +58,15 @@ 3. results - The angle between pointing direction and earth center = 15deg - - STT flag is always 1, since the angle $`\theta_{es}`$ between the sight direction and the earth's edge direction is 15 - 8.6 = 6.4deg < 10deg. + - STT flag is always 1, since the angle $\theta_{es}$ between the sight direction and the earth's edge direction is 15 - 8.6 = 6.4deg < 10deg. ![](./figs/stt_flag_15.png) - The angle between pointing direction and earth center = 20deg - - STT flag is always 0, since the angle $`\theta_{es}`$ between the sight direction and the earth's edge direction is 20 - 8.6 = 11.4deg > 10deg. + - STT flag is always 0, since the angle $\theta_{es}$ between the sight direction and the earth's edge direction is 20 - 8.6 = 11.4deg > 10deg. ![](./figs/stt_flag_20.png) - The angle between pointing direction and earth center = 30deg - - STT flag is always 0, since the angle $`\theta_{es}`$ between the sight direction and the earth's edge direction is 30 - 8.6 = 21.4deg > 10deg. + - STT flag is always 0, since the angle $\theta_{es}$ between the sight direction and the earth's edge direction is 30 - 8.6 = 21.4deg > 10deg. ![](./figs/stt_flag_30.png) diff --git a/Specifications/Component/Abstract/Spec_ComponentBase.md b/Specifications/Component/Abstract/Spec_ComponentBase.md index 4d19269a..69d27f16 100644 --- a/Specifications/Component/Abstract/Spec_ComponentBase.md +++ b/Specifications/Component/Abstract/Spec_ComponentBase.md @@ -21,7 +21,7 @@ 2. inputs - `prescaler` + `prescaler` determines the execution cycle of the `MainRoutine` function. - + The period of `MainRoutine` equals to `SimTime::compo_update_interval_sec` $`\times`$ `prescaler`. + + The period of `MainRoutine` equals to `SimTime::compo_update_interval_sec` $\times$ `prescaler`. - `clock_gen` + `clock_gen` is an instance that simulates the clock of a component. + Users do not need to care about this. @@ -29,7 +29,7 @@ + `power_port` is an instance that simulates the power supply - `fast_prescaler` + `fast_prescaler` determines the execution cycle of the `FastUpdate` function. - + The period of `FastUpdate` equals to `SimTime::compo_update_interval_sec` $`\times`$ `fast_prescaler`. + + The period of `FastUpdate` equals to `SimTime::compo_update_interval_sec` $\times$ `fast_prescaler`. + If you don't need to use `FastUpdate`, you don't need to specify this (it is set to 1 by default). 3. algorithm - N/A diff --git a/Specifications/Component/Mission/Spec_Telescope_en.md b/Specifications/Component/Mission/Spec_Telescope_en.md index 05b8bbe1..d29de342 100644 --- a/Specifications/Component/Mission/Spec_Telescope_en.md +++ b/Specifications/Component/Mission/Spec_Telescope_en.md @@ -55,12 +55,12 @@ + false: The celestial body is not in forbidden angle 3. process to judge - - The judging process is calculated in the telescope's component coordinate. $`q_{b2c}`$ is the quaternion to convert from the body coordinate(B) to the component coordinate(C). Specify $`q_{b2c}`$ in `Telescope.ini`. The X-axis of the component coordinate is defined as the line of sight of the telescope. + - The judging process is calculated in the telescope's component coordinate. $q_{b2c}$ is the quaternion to convert from the body coordinate(B) to the component coordinate(C). Specify $q_{b2c}$ in `Telescope.ini`. The X-axis of the component coordinate is defined as the line of sight of the telescope. 2. `Observe` 1. overview - This function judges whether the celestial bodies(provided by `CelesInfo`) are in the field of view and outputs the position of them on the image sensor if they are in the field of view - - If they are not in the field of view, this function outputs $`(-1,-1)`$. + - If they are not in the field of view, this function outputs $(-1,-1)$. 2. input and output - input @@ -76,7 +76,7 @@ - A 2D coordinate on the image sensor is defined to handle the position on the image sensor. The definition is as follows: - The X-axis of the image sensor coordinate corresponds with the Z-axis of the component coordinate. - The Y-axis of the image sensor coordinate corresponds with the Y-axis of the component coordinate. - - Then, the inclination angle from the X-axis of the celestial body's direction in the XZ plane of the component coordinate is calculated using $`(x_c, y_c, z_c)`$ as follows: + - Then, the inclination angle from the X-axis of the celestial body's direction in the XZ plane of the component coordinate is calculated using $(x_c, y_c, z_c)$ as follows: ```math tan^{-1}⁡\frac{z_c}{x_c} ``` @@ -86,7 +86,7 @@ tan^{-1}⁡\frac{y_c}{x_c} ``` - - They are written as `arg_x` and `arg_y` in the code. In this manual, $`\theta_x`$ and $`\theta_y`$ are used for them. If $`\theta_x`$ is within FOV_x and $`\theta_y`$ is within FOV_y, the celestial body is judged to be in the field of view. + - They are written as `arg_x` and `arg_y` in the code. In this manual, $\theta_x$ and $\theta_y$ are used for them. If $\theta_x$ is within FOV_x and $\theta_y$ is within FOV_y, the celestial body is judged to be in the field of view.
@@ -96,7 +96,7 @@
2. Calculate process for calculating the position of the image - - The origin of the sensor coordinate is the corner of the image sensor, so $`x_{imgsensor}`$ and $`y_{imgsensor}`$ have positive values. The unit of them is pixel(pix). In this manual, $`N_x`$ and $`N_y`$ are used for the total number of pixels along x, y axes of the sensor coordinate (They are `x_num_of_pix` and `y_num_of_pix` in the code). In the same way, `X` and `Y` are used for the position of the celestial body on the image sensor (They are `pos_imgsensor[0]` and `pos_imgsensor[1]`). Then, `X` and `Y` are calculated as follows: + - The origin of the sensor coordinate is the corner of the image sensor, so $x_{imgsensor}$ and $y_{imgsensor}$ have positive values. The unit of them is pixel(pix). In this manual, $N_x$ and $N_y$ are used for the total number of pixels along x, y axes of the sensor coordinate (They are `x_num_of_pix` and `y_num_of_pix` in the code). In the same way, `X` and `Y` are used for the position of the celestial body on the image sensor (They are `pos_imgsensor[0]` and `pos_imgsensor[1]`). Then, `X` and `Y` are calculated as follows: ```math X=\frac{N_x}{2}\times\frac{\tan(\theta_x)}{\tan(FOV_x)}+\frac{N_x}{2} ``` @@ -104,7 +104,7 @@ Y=\frac{N_y}{2}\times\frac{\tan(\theta_y)}{\tan(FOV_y)}+\frac{N_y}{2} ``` - If the celestial body is not in the field of view, the output is $`X=Y=-1`$. + If the celestial body is not in the field of view, the output is $X=Y=-1$. 3. `ObserveStars` 1. overview @@ -132,7 +132,7 @@ In this section, the output of the functions when some angular velocity is input 1. Input of angular velocity around X-axis of the body coordinate 1. overview - - input $`ω_b=[0.1~0~0]^T`$ ． + - input $ω_b=[0.1~0~0]^T$ ． 2. conditions for the verification 1. input files - `SampleSimBase.ini` @@ -192,7 +192,7 @@ In this section, the output of the functions when some angular velocity is input
The track of the image of the Moon on the image sensor
- This figure shows the track makes a circle. This result seems to be reasonable, because the angular velocity around the X-axis of the body coordinate correspond with that of the component coordinate, for $`q_{b2c}=[0~0~0~1]^T`$ . The 3D plot of MOON_POS_B for further verification is as follows: + This figure shows the track makes a circle. This result seems to be reasonable, because the angular velocity around the X-axis of the body coordinate correspond with that of the component coordinate, for $q_{b2c}=[0~0~0~1]^T$ . The 3D plot of MOON_POS_B for further verification is as follows:
3D plot of MOON_POS_B @@ -229,7 +229,7 @@ In this section, the output of the functions when some angular velocity is input 2. input of angular velocity around y axis of the body coordinate 1. overview - - The angular velocity input is $`ω_b=[0.1~0~0]^T`$ ．The other condition is the same as the case of 1. Note that the verification of the case around z axis is omitted because y and z are equivalent under this condition. + - The angular velocity input is $ω_b=[0.1~0~0]^T$ ．The other condition is the same as the case of 1. Note that the verification of the case around z axis is omitted because y and z are equivalent under this condition. 1. result 1. judge for forbidden angle - The angle from the line of sight about the direction of the Sun, the Earth, the Moon is as follows: diff --git a/Specifications/Component/Mission/Spec_Telescope_ja.md b/Specifications/Component/Mission/Spec_Telescope_ja.md index 47b4738f..66577fe5 100644 --- a/Specifications/Component/Mission/Spec_Telescope_ja.md +++ b/Specifications/Component/Mission/Spec_Telescope_ja.md @@ -51,7 +51,7 @@ + 禁止角に入っているならばtrue，そうでないならばfalse 3. 判定処理 - - 処理はコンポ座標系で行っている．$`q_{b2c}`$はTelescope.iniで指定する機体座標系(B)→コンポ座標系(C)変換Quaternionである．コンポ座標系は，望遠鏡視線方向をx方向としている． + - 処理はコンポ座標系で行っている．$q_{b2c}$はTelescope.iniで指定する機体座標系(B)→コンポ座標系(C)変換Quaternionである．コンポ座標系は，望遠鏡視線方向をx方向としている． 2. `Observe` @@ -79,7 +79,7 @@
- このとき，コンポ座標系からみた天体の位置 $`(x_c,y_c,z_c)`$ より，コンポ座標系xz平面でみた天体の方向のx軸からの偏角は + このとき，コンポ座標系からみた天体の位置 $(x_c,y_c,z_c)$ より，コンポ座標系xz平面でみた天体の方向のx軸からの偏角は ```math tan^{-1}⁡\frac{z_c}{x_c} @@ -91,10 +91,10 @@ tan^{-1}⁡\frac{y_c}{x_c} ``` - で計算される．これらの値はそれぞれソースコード上でarg_x，およびarg_yである（説明上，それぞれ $`\theta_x`$，$`\theta_y`$と表記する）． $`\theta_x`$ がセンサ座標系x方向の視野角（FOV_x）以下であり，かつ $`\theta_y`$ がセンサ座標系y方向の視野角（FOV_y）以下であれば対象の天体が画角内に入っていると判定する． + で計算される．これらの値はそれぞれソースコード上でarg_x，およびarg_yである（説明上，それぞれ $\theta_x$，$\theta_y$と表記する）． $\theta_x$ がセンサ座標系x方向の視野角（FOV_x）以下であり，かつ $\theta_y$ がセンサ座標系y方向の視野角（FOV_y）以下であれば対象の天体が画角内に入っていると判定する． 4. イメージセンサ上における位置計算処理 - - センサ座標系の原点はイメージセンサの角に位置しており， $`x_{imgsensor}`$ および $`y_{imgsensor}`$ は原点を基準とした正の値で表現される．単位はピクセル（pix）である．ここでは説明のためイメージセンサのx方向全ピクセル数（ソースコード上ではx_num_of_pix），y方向全ピクセル数（ソースコード上ではy_num_of_pix）をそれぞれ $`N_x`$ および $`N_y`$ と表記する．同様に，天体のイメージセンサ上の位置（ソースコード上ではpos_imgsensor）をX，Yと表記することにすれば，これらは + - センサ座標系の原点はイメージセンサの角に位置しており， $x_{imgsensor}$ および $y_{imgsensor}$ は原点を基準とした正の値で表現される．単位はピクセル（pix）である．ここでは説明のためイメージセンサのx方向全ピクセル数（ソースコード上ではx_num_of_pix），y方向全ピクセル数（ソースコード上ではy_num_of_pix）をそれぞれ $N_x$ および $N_y$ と表記する．同様に，天体のイメージセンサ上の位置（ソースコード上ではpos_imgsensor）をX，Yと表記することにすれば，これらは ```math X=\frac{N_x}{2}\times\frac{\tan(\theta_x)}{\tan(FOV_x)}+\frac{N_x}{2} @@ -103,7 +103,7 @@ Y=\frac{N_y}{2}\times\frac{\tan(\theta_y)}{\tan(FOV_y)}+\frac{N_y}{2} ``` - として計算される．なお，天体が画角内に入っていないときは， $`X=Y=-1`$ なる出力を与える． + として計算される．なお，天体が画角内に入っていないときは， $X=Y=-1$ なる出力を与える． 3. `ObserveStars` 1. 概要 @@ -129,7 +129,7 @@ 禁止角判定， `Observe` 関数の動作， `ObserveStars` 関数の動作が正しく行えているかどうかを判定する．ここでは，body座標系において各軸まわりの角速度入力があったときに，どのような出力を得るかを確認し，その妥当性を検証することにする． 1. x軸周りの角速度を与えた場合 1. 概要 - - 角速度入力　$`ω_b=[0.1~0~0]^T`$を与える． + - 角速度入力　$ω_b=[0.1~0~0]^T$を与える． 2. シミュレーション条件 1. 入力ファイル @@ -225,7 +225,7 @@ 2. y軸周りの角速度を与えた場合 1. 概要 - - 角速度入力 $`ω_b=[0 ~0.1~0]^T`$ を与える．その他のシミュレーション条件はx軸周りの場合と同様である．なお，本シミュレーションでは望遠鏡に関してはy軸とz軸が等価であるから，z軸周りの回転を与えた場合の検証は割愛する． + - 角速度入力 $ω_b=[0 ~0.1~0]^T$ を与える．その他のシミュレーション条件はx軸周りの場合と同様である．なお，本シミュレーションでは望遠鏡に関してはy軸とz軸が等価であるから，z軸周りの回転を与えた場合の検証は割愛する． 1. 結果 1. 禁止角判定 - 太陽・地球・月の視線方向からの角度は，以下のようである． diff --git a/Specifications/Component/Propulsion/Spec_SimpleThruster.md b/Specifications/Component/Propulsion/Spec_SimpleThruster.md index 75bed441..4fcb26d2 100644 --- a/Specifications/Component/Propulsion/Spec_SimpleThruster.md +++ b/Specifications/Component/Propulsion/Spec_SimpleThruster.md @@ -48,7 +48,7 @@ \tag{1} ``` - where $`F_{thrust} `$ is thrust magnitude, $`\epsilon`$ is the duty of thruster, $`F_{max}`$ is the maximum thrust magnitude, and $`n_{f}`$ is the error of thrust magnitude. + where $F_{thrust} $ is thrust magnitude, $\epsilon$ is the duty of thruster, $F_{max}$ is the maximum thrust magnitude, and $n_{f}$ is the error of thrust magnitude. Thrust direction can be calculated as follows: @@ -61,7 +61,7 @@ \tag{2} ``` - where $`\boldsymbol{d}_{true}`$ is the thrust vector without errors, $`n`$ is the random angles to rotate the direction of error $`\boldsymbol{d}_{err}`$, $`\boldsymbol{d}_{x}`$ is the vector which is not equal to $`\boldsymbol{d}_{true}`$, $`n_d`$ is the directional error, and $`\boldsymbol{d}_{thrust}`$ is the thrust vector with errors. $`\boldsymbol{q}(\boldsymbol{d},n)`$ is the quaternion which has the rotation axis $`\boldsymbol{d}`$ and the rotation angle $`n`$. + where $\boldsymbol{d}_{true}$ is the thrust vector without errors, $n$ is the random angles to rotate the direction of error $\boldsymbol{d}_{err}$, $\boldsymbol{d}_{x}$ is the vector which is not equal to $\boldsymbol{d}_{true}$, $n_d$ is the directional error, and $\boldsymbol{d}_{thrust}$ is the thrust vector with errors. $\boldsymbol{q}(\boldsymbol{d},n)$ is the quaternion which has the rotation axis $\boldsymbol{d}$ and the rotation angle $n$. Thrust can be calculated as follows: @@ -70,7 +70,7 @@ \tag{3} ``` - where $`\boldsymbol{F}_{thrust}`$ is thrust. + where $\boldsymbol{F}_{thrust}$ is thrust. 2. Torque 1. overview @@ -90,7 +90,7 @@ \tag{4} ``` - where $`\boldsymbol{T}_{thrust} `$ is torque by the thruster, $`\boldsymbol{v}_{thruster}`$ is thruster position and $`\boldsymbol{v}_{SC}`$ is the mass center of spacecraft. + where $\boldsymbol{T}_{thrust} $ is torque by the thruster, $\boldsymbol{v}_{thruster}$ is thruster position and $\boldsymbol{v}_{SC}$ is the mass center of spacecraft. ## 3. Results of verifications 1. Case1 diff --git a/Specifications/Disturbance/Spec_GGTorque.md b/Specifications/Disturbance/Spec_GGTorque.md index 96ea3861..579ed59b 100644 --- a/Specifications/Disturbance/Spec_GGTorque.md +++ b/Specifications/Disturbance/Spec_GGTorque.md @@ -33,7 +33,7 @@ \tag{2} ``` - where $`\mu`$ is the gravitational constant of the Earth, $`R_0`$ is the distance between the Earth center and the satellite, $`\boldsymbol{R_0}`$ is the vector from Earth center to the satellite, $`\boldsymbol{I}`$ is the inertia tensor of the satellite. + where $\mu$ is the gravitational constant of the Earth, $R_0$ is the distance between the Earth center and the satellite, $\boldsymbol{R_0}$ is the vector from Earth center to the satellite, $\boldsymbol{I}$ is the inertia tensor of the satellite. ## 3. Results of verifications @@ -52,7 +52,7 @@ - Disturbance torque: All Disable except GGTorque 3. results - - The order of the gg torque is $`10^7`$, which seems a proper value. + - The order of the gg torque is $10^7$, which seems a proper value. ![](./figs/test_ggtorque.png) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 1178916d..2bcfd143 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -28,11 +28,11 @@ 1. overview - This function reads the geopotential coefficients from the EGM96 file `egm96_to360.ascii`. - The file doesn't have coefficients for `n=0,m=0`, `n=1,m=0`, and `n=1,m=1`. - - All coefficients are completely normalized by following normalization function $`N_{n,m}`$ + - All coefficients are completely normalized by following normalization function $N_{n,m}$ ```math N_{n,m}=\sqrt{\frac{(n+m)!}{(2-\delta_{0m})(2n+1)(n-m)!}}\\ ``` - - where $`\delta_{0m}`$ is the Kronecker delta. + - where $\delta_{0m}$ is the Kronecker delta. 2. inputs and outputs - Input @@ -40,7 +40,7 @@ - maximum degree for geopotential calculation - Output - **Return**: reading is succeeded or not. - - Normalized coefficient $`C_{n,m}`$ and $`S_{n,m}`$ + - Normalized coefficient $C_{n,m}$ and $S_{n,m}$ 3. algorithm - The file format of `egm96_to360.ascii` is `n,m,Cnm,Snm,sigmaCnm,sigmaSnm` in line with space delimiter. In this calculation, the `sigmaCnm` and `sigmaSnm` are not used. @@ -49,7 +49,7 @@ ```math N_{line}=\frac{1}{2}(n+1)(n+2)-3 ``` - where $`n`$ is maximum degree, and -3 is for the coefficients of `n=0,m=0`, `n=1,m=0`, and `n=1,m=1`, which are not in the file. + where $n$ is maximum degree, and -3 is for the coefficients of `n=0,m=0`, `n=1,m=0`, and `n=1,m=1`, which are not in the file. 4. note @@ -60,7 +60,7 @@ 1. overview - We chose to use the recursion algorithm written in chapter 3.2.4 of [Satellite Orbits](https://www.springer.com/jp/book/9783540672807) since the calculation of the Legendre polynomials for spherical harmonics is time-consuming. - However, the original equation in the book is unnormalized form, and it is not suitable with the normalized coefficients. - - For a small degree, users can directly use the normalized function $`N_{n,m}`$ to unnormalize the coefficients or to normalize the functions $`V_{n,m}`$ and $`V_{n,m}`$ . But for a large degree, the factorial calculation in the $`N_{n,m}`$ reaches a huge value, which standard `double` variables cannot handle. + - For a small degree, users can directly use the normalized function $N_{n,m}$ to unnormalize the coefficients or to normalize the functions $V_{n,m}$ and $V_{n,m}$ . But for a large degree, the factorial calculation in the $N_{n,m}$ reaches a huge value, which standard `double` variables cannot handle. - To avoid that, the normalized recursion algorithm was derived as described in Section 3. - There are the following two functions: - `v_w_nn_update` @@ -69,13 +69,13 @@ 2. inputs and outputs - Inputs - Both functions - - Satellite position in ECEF frame $`x, y, z`$ - - degree and order as $`n`$ and $`m`$ - - `v_w_nn_update`: $`V_{n-1,n-1}`$ and $`W_{n-1,n-1}`$ - - `v_w_nm_update`: $`V_{n-1,m}, W_{n-1,m}, V_{n-2,m}`$, and $`W_{n-2,m}`$ + - Satellite position in ECEF frame $x, y, z$ + - degree and order as $n$ and $m$ + - `v_w_nn_update`: $V_{n-1,n-1}$ and $W_{n-1,n-1}$ + - `v_w_nm_update`: $V_{n-1,m}, W_{n-1,m}, V_{n-2,m}$, and $W_{n-2,m}$ - Outputs - - `v_w_nn_update`: $`V_{n,n}`$ and $`W_{n,n}`$ - - `v_w_nm_update`: $`V_{n,m}`$ and $`W_{n,m}`$ + - `v_w_nn_update`: $V_{n,n}$ and $W_{n,n}$ + - `v_w_nm_update`: $V_{n,m}$ and $W_{n,m}$ 3. algorithm @@ -148,34 +148,34 @@ The recurrence relations for normalized V and W are derived as follows 2. inputs and outputs - Input - - normalized Coefficients: $`\bar{C}_{n,m}`$ and $`\bar{S}_{n,m}`$ - - normalized function: $`\bar{V}_{n,m}`$ and $`\bar{W}_{n,m}`$ + - normalized Coefficients: $\bar{C}_{n,m}$ and $\bar{S}_{n,m}$ + - normalized function: $\bar{V}_{n,m}$ and $\bar{W}_{n,m}$ - Output - Gravity acceleration in ECEF frame 3. algorithm For unnormalized algorithms, See chapter 3.2.5 of [Satellite Orbits](https://www.springer.com/jp/book/9783540672807). -When we use the normalized coefficients $`\bar{C}_{n,m}`$ and $`\bar{S}_{n,m}`$ and $`\bar{V}_{n,m}`$ and $`\bar{W}_{n,m}`$ functions, the acceleration calculation is described like follows +When we use the normalized coefficients $\bar{C}_{n,m}$ and $\bar{S}_{n,m}$ and $\bar{V}_{n,m}$ and $\bar{W}_{n,m}$ functions, the acceleration calculation is described like follows ```math \ddot{x}_{n,m}=-\frac{GM}{Re^{2}}\bar{C}_{n,0}\bar{V}_{n+1,1} =-\frac{GM}{Re^{2}} C_{n,0}V_{n+1,1} \frac{N_{n,0}}{N_{n+1,1}}\quad(m=0) ``` -The division of normalized function $`\frac{N_{n,0}}{N_{n+1,1}}`$ should be removed, so we have to multiply following correction factors into the equation. +The division of normalized function $\frac{N_{n,0}}{N_{n+1,1}}$ should be removed, so we have to multiply following correction factors into the equation. When $m=0$, following correction factors are useful for x and y acceleration ```math \frac{N_{n+1,1}}{N_{n,0}}=\sqrt{\frac{1}{2}}\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+2)(n+1)} ``` -When $`m=1`$, following correction factors are useful for x and y acceleration +When $m=1$, following correction factors are useful for x and y acceleration ```math \frac{N_{n+1,0}}{N_{n,1}}=\sqrt{2}\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{n(n+1)}}\quad(m=1)\\ ``` -When $`m>1`$, following correction factors are useful for x and y acceleration +When $m>1$, following correction factors are useful for x and y acceleration ```math \frac{N_{n+1,m+1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+m+1)(n+m+2)}\\ \frac{N_{n+1,m-1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{(n-m+1)(n-m+2)}}\\ ``` -When $`m>=0`$, following correction factors are useful for z acceleration +When $m>=0$, following correction factors are useful for z acceleration ```math \frac{N_{n+1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{n+m+1}{n-m+1}} ``` diff --git a/Specifications/Disturbance/Spec_MagDisturbance.md b/Specifications/Disturbance/Spec_MagDisturbance.md index 6039cc66..2797f96e 100644 --- a/Specifications/Disturbance/Spec_MagDisturbance.md +++ b/Specifications/Disturbance/Spec_MagDisturbance.md @@ -32,7 +32,7 @@ \tag{1} ``` - where $`\boldsymbol{M}`$ is the residual magnetic moment in the body-fixed frame, $`\boldsymbol{B}`$ is the magnetic field in the body-fixed frame. + where $\boldsymbol{M}$ is the residual magnetic moment in the body-fixed frame, $\boldsymbol{B}$ is the magnetic field in the body-fixed frame. 2. `CalcRMM` function 1. overview @@ -53,7 +53,7 @@ \boldsymbol{r}(t_{k+1}) = \boldsymbol{r}(t_{k}) +\boldsymbol{w_2}(t_{k+1}) \tag{2} ``` - where $`\boldsymbol{M}_0`$ is the average residual magnetic moment in the body-fixed frame, $`\boldsymbol{r}`$ is the random walk of RMM, and $`\boldsymbol{w}_i \sim N([0,0,0],\Sigma_i)`$ is the white noise. + where $\boldsymbol{M}_0$ is the average residual magnetic moment in the body-fixed frame, $\boldsymbol{r}$ is the random walk of RMM, and $\boldsymbol{w}_i \sim N([0,0,0],\Sigma_i)$ is the white noise. ## 3. Results of verifications diff --git a/Specifications/Disturbance/Spec_SurfaceForce.md b/Specifications/Disturbance/Spec_SurfaceForce.md index 1ffa5f05..0185372b 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce.md +++ b/Specifications/Disturbance/Spec_SurfaceForce.md @@ -47,27 +47,27 @@ 3. algorithm - Let us consider the following coordinate on a surface for surface force calculation - - $`\boldsymbol{n}`$ is the normal vector of the surface - - $`\boldsymbol{v}`$ is the direction vector of the disturbance source (e.g., sun direction vector or velocity vector) + - $\boldsymbol{n}$ is the normal vector of the surface + - $\boldsymbol{v}$ is the direction vector of the disturbance source (e.g., sun direction vector or velocity vector) SummaryCalculationTime - - $`\boldsymbol{t}`$ is the direction of in-plane force. + - $\boldsymbol{t}$ is the direction of in-plane force. ```math \boldsymbol{t}=\frac{\boldsymbol{v}\times\boldsymbol{n}}{|\boldsymbol{v}\times\boldsymbol{n}|}\times\boldsymbol{n} ``` - Surface force and torque acting on the surface is expressed as following equation - - $`\boldsymbol{r}_{s}`$ is the position vector of the surface - - $`\boldsymbol{r}_{cg}`$ is the position vector of the center of mass + - $\boldsymbol{r}_{s}$ is the position vector of the surface + - $\boldsymbol{r}_{cg}$ is the position vector of the center of mass ```math \boldsymbol{F}=-C_{n}\boldsymbol{n}+C_{t}\boldsymbol{t}\\ \boldsymbol{T}=(\boldsymbol{r}_{s}-\boldsymbol{r}_{cg})\times\boldsymbol{F} ``` - - Detail of the $`C_{n}`$ and $`C_{t}`$ are defined by subclasses by using `CalcCoef` function + - Detail of the $C_{n}$ and $C_{t}$ are defined by subclasses by using `CalcCoef` function 4. note - NA diff --git a/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md b/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md index 326984b7..8a894d5d 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md +++ b/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md @@ -30,15 +30,15 @@ C_{t}=\frac{1}{2}\rho A v^2 C_{t}^{\prime} ``` - - This function mainly calculates the common part of the coefficient calculation. $`C_{n}^{\prime}`$ and $`C_{t}^{\prime}`$ are calculated in `CalCnCt` function, and they will be used in this function. + - This function mainly calculates the common part of the coefficient calculation. $C_{n}^{\prime}$ and $C_{t}^{\prime}$ are calculated in `CalCnCt` function, and they will be used in this function. 2. inputs and outputs - input - - $`\rho`$: air density [kg/m3] - - $`v`$: Relative velocity between the spacecraft and the atmosphere [m/s] - - $`A`$: Area of the surface [m2] (given by ) + - $\rho$: air density [kg/m3] + - $v$: Relative velocity between the spacecraft and the atmosphere [m/s] + - $A$: Area of the surface [m2] (given by ) - output - - coefficients $`C_{n}`$ and $`C_{t}`$ + - coefficients $C_{n}$ and $C_{t}$ 3. algorithm - See above equations. @@ -47,35 +47,35 @@ 2. `CalCnCt` 1. overview - - `CalCnCt` calculates $`C_{n}^{\prime}`$ and $`C_{t}^{\prime}`$. + - `CalCnCt` calculates $C_{n}^{\prime}$ and $C_{t}^{\prime}$. 2. inputs and outputs - input variables - - $`\boldsymbol{v}`$: Relative velocity vector between the spacecraft and the atmosphere [m/s] + - $\boldsymbol{v}$: Relative velocity vector between the spacecraft and the atmosphere [m/s] - Currently, we assume that this value equals spacecraft velocity in the body-fixed frame. - input parameters - - $`\sigma_{d}`$: Diffuse coefficients for air drag - - Ini file provide specularity for air drag $`\sigma_{s}`$, and the diffuse coefficient is derived as $`\sigma_{d}=1-\sigma_{s}`$. + - $\sigma_{d}$: Diffuse coefficients for air drag + - Ini file provide specularity for air drag $\sigma_{s}$, and the diffuse coefficient is derived as $\sigma_{d}=1-\sigma_{s}$. - **Note**: There is no absorption term for air drag. Thus total reflectivity is set as 1. - - $`T_{w}`$: Temperature of the surface [K] - - $`T_{m}`$: Temperature of the atmosphere [K] - - $`M`$: Molecular weight of the thermosphere [g/mol] - - In the default ini file, we use $`M=18`$, and it is a little bit smaller than the molecular weight of atmosphere $`M=29`$. [Structure of the Thermosphere](https://www.sciencedirect.com/science/article/pii/0032063361900368?via%3Dihub) provides information on the molecular weight of the thermosphere. + - $T_{w}$: Temperature of the surface [K] + - $T_{m}$: Temperature of the atmosphere [K] + - $M$: Molecular weight of the thermosphere [g/mol] + - In the default ini file, we use $M=18$, and it is a little bit smaller than the molecular weight of atmosphere $M=29$. [Structure of the Thermosphere](https://www.sciencedirect.com/science/article/pii/0032063361900368?via%3Dihub) provides information on the molecular weight of the thermosphere. - outputs - - $`C_{n}^{\prime}`$ and $`C_{t}^{\prime}`$ + - $C_{n}^{\prime}$ and $C_{t}^{\prime}$ 3. algorithm - - $`C_{n}^{\prime}`$ and $`C_{t}^{\prime}`$ are calculated as following equations + - $C_{n}^{\prime}$ and $C_{t}^{\prime}$ are calculated as following equations ```math C_{n}^{\prime} = \frac{2-\sigma_{d}}{\sqrt{\pi}}\frac{\Pi(S_{n})}{S^{2}}+\frac{\sigma_{d}}{2}\frac{\chi(S_{n})}{S^{2}}\sqrt{\frac{T_{w}}{T_{m}}}\\ C_{t}^{\prime} =\frac{\sigma_{d}}{\sqrt{\pi}}\frac{\chi(S_{n})}{S^{2}}S_{t} ``` - - $`S, S_{n}, S_{t}`$ are defined as follows - - $`k=1.38064852E-23`$ is the Boltzmann constant - - $`\theta`$ is the angle between the normal vector and the velocity vector - - $`\cos{\theta}`$ and $`\sin{\theta}`$ are calculated in `SurfaceForce` base class. + - $S, S_{n}, S_{t}$ are defined as follows + - $k=1.38064852E-23$ is the Boltzmann constant + - $\theta$ is the angle between the normal vector and the velocity vector + - $\cos{\theta}$ and $\sin{\theta}$ are calculated in `SurfaceForce` base class. ```math S=\sqrt{\frac{Mv^{2}}{2kT_{w}}}\\ @@ -83,7 +83,7 @@ S_{t}=S\sin{\theta}\\ ``` - - $`\Pi(x)`$ and $`\chi(x)`$ are defined as follows + - $\Pi(x)$ and $\chi(x)$ are defined as follows - where `erf` is the [Gauss error function](https://en.wikipedia.org/wiki/Error_function). ```math @@ -107,9 +107,9 @@ | parameters/results | Case 1 | Case 2 | Case 3 | | --------------------------- | ------- | ------- | ------- | - | $`\sigma_{d}`$ | 0.8 | 0.6 | 0.4 | - | $`\theta`$ rad | 0.202 | 0.202 | 0.202 | - | $`v`$ m/s | 7420 | 7420 | 7420 | + | $\sigma_{d}$ | 0.8 | 0.6 | 0.4 | + | $\theta$ rad | 0.202 | 0.202 | 0.202 | + | $v$ m/s | 7420 | 7420 | 7420 | | Out-plane force (S2E) | 2.30297 | 2.68680 | 3.07062 | | Out-plane force (reference) | 2.30297 | 2.68680 | 3.07062 | | Out-plane force (S2E) | 0.31514 | 0.23636 | 0.15757 | diff --git a/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md b/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md index 5c59dc67..c81ae525 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md +++ b/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md @@ -26,19 +26,19 @@ 2. inputs and outputs - inputs - - $`v_{s}`$:Direction vector of the sun (s/c -> sun) at the body frame - - $`P`$ Solar pressure at the position of the spacecraft [N/m^2] + - $v_{s}$:Direction vector of the sun (s/c -> sun) at the body frame + - $P$ Solar pressure at the position of the spacecraft [N/m^2] - setting parameters - - $`\nu`$ : Total reflectance - - $`\nu = 1-\alpha`$, where $`\alpha`$ is the absorption of the sun spectrum. - - $`\mu`$ : Specularity. Ratio of specular reflection inside the total reflected light. - - $`A`$ : Area of the surface + - $\nu$ : Total reflectance + - $\nu = 1-\alpha$, where $\alpha$ is the absorption of the sun spectrum. + - $\mu$ : Specularity. Ratio of specular reflection inside the total reflected light. + - $A$ : Area of the surface - outputs - - $`C_{n}`$ and $`C_{t}`$ + - $C_{n}$ and $C_{t}$ 3. algorithm - - $`C_{n}`$ and $`C_{t}`$ are calculated as follows: - - $`\theta`$ is the angle between the normal vector and the sun vector. + - $C_{n}$ and $C_{t}$ are calculated as follows: + - $\theta$ is the angle between the normal vector and the sun vector. ```math C_{n} = AP \left((1+\nu\mu)\cos^{2}{\theta}+\frac{2}{3}\nu(1-\mu)\cos{\theta} \right)\\ @@ -55,7 +55,7 @@ - In the perfect reflection case, the direction of the SRP force will be opposite from the direction of the sun. 2. conditions for the verification - - We assumed that the structure of the spacecraft is a 50-cm cube whose optical property is the perfect specular reflection($`\nu=\mu=1`$). + - We assumed that the structure of the spacecraft is a 50-cm cube whose optical property is the perfect specular reflection($\nu=\mu=1$). 3. results - We confirmed that the direction of the SRP force is opposite from the direction of the sun at the body frame. diff --git a/Specifications/Disturbance/Spec_ThirdBodyGravity.md b/Specifications/Disturbance/Spec_ThirdBodyGravity.md index 82e3e4c9..de1b403e 100644 --- a/Specifications/Disturbance/Spec_ThirdBodyGravity.md +++ b/Specifications/Disturbance/Spec_ThirdBodyGravity.md @@ -38,9 +38,9 @@ 3. algorithm - Definition of the variables - - \boldsymbol{r}`$: the position of the spacecraft (written as `sat_pos_i` in the code) - - $`\boldsymbol{s}`$: the position of the third body (written as `third_body_pos_i` in the code) - - $`\boldsymbol{s}-\boldsymbol{r}`$: the position of the third body seen from the spacecraft (written as `third_body_pos_from_sc_i` in the code) + - \boldsymbol{r}$: the position of the spacecraft (written as `sat_pos_i` in the code) + - $\boldsymbol{s}$: the position of the third body (written as `third_body_pos_i` in the code) + - $\boldsymbol{s}-\boldsymbol{r}$: the position of the third body seen from the spacecraft (written as `third_body_pos_from_sc_i` in the code) ThirdBodyGravity_general diff --git a/Specifications/Dynamics/Spec_AttitudeDynamics.md b/Specifications/Dynamics/Spec_AttitudeDynamics.md index e38e25cd..a2ef38b4 100644 --- a/Specifications/Dynamics/Spec_AttitudeDynamics.md +++ b/Specifications/Dynamics/Spec_AttitudeDynamics.md @@ -58,7 +58,7 @@ - (double) dt: The duration for the attitude propagation 3. algorithm - If the differential equation (1) is given, the state quantity in $`n+1`$ step can be calculated as (2). + If the differential equation (1) is given, the state quantity in $n+1$ step can be calculated as (2). ```math \hat{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x},t) @@ -69,14 +69,14 @@ +2\boldsymbol{k_3}+\boldsymbol{k_4}) \tag{2} ``` - where $`\Delta t`$ is a time step, which meets the equation (3). + where $\Delta t$ is a time step, which meets the equation (3). ```math t_{n+1} = t_{n} + \Delta t \tag{3} ``` - $`k_i \, (i=1,2,3,4)`$, which has the same number of elements, can be calculated as the equations (4). + $k_i \, (i=1,2,3,4)$, which has the same number of elements, can be calculated as the equations (4). ```math \boldsymbol{k_{1}} = \boldsymbol{f}(\boldsymbol{x_n},t_n) \\ @@ -86,7 +86,7 @@ \tag{4} ``` - In this attitude propagation, the quantity of state $`\boldsymbol{x}`$ consists of 7 elements, including `Quaternion_i2b` and angular velocity $`\boldsymbol{\omega}_b`$. + In this attitude propagation, the quantity of state $\boldsymbol{x}$ consists of 7 elements, including `Quaternion_i2b` and angular velocity $\boldsymbol{\omega}_b$. ```math \boldsymbol{\omega}_b = [{\omega}_{bx} \, {\omega}_{by} \, {\omega}_{bz}]^T \\ @@ -118,7 +118,7 @@ \tag{6} ``` - where $`\boldsymbol{\omega}_b`$[rad/s] is angular velocity in the body-fixed coordinate, $`\boldsymbol{I}_b`$[kgm$`^2`$] is inertia tensor of the satellite, $`\boldsymbol{T}_b`$[Nm] is torque in the body-fixed coordinate, $`\boldsymbol{h}_b`$[Nms] is angular momentum of the satellite in the body-fixed coordinate. + where $\boldsymbol{\omega}_b$[rad/s] is angular velocity in the body-fixed coordinate, $\boldsymbol{I}_b$[kgm$^2$] is inertia tensor of the satellite, $\boldsymbol{T}_b$[Nm] is torque in the body-fixed coordinate, $\boldsymbol{h}_b$[Nms] is angular momentum of the satellite in the body-fixed coordinate. Quaternion_i2b is calculated from the kinematics equation (7). This equation is executed in `Omega4Kinematics` function. ```math diff --git a/Specifications/Dynamics/Spec_ControlledAttitude.md b/Specifications/Dynamics/Spec_ControlledAttitude.md index 5938648a..584fbda4 100644 --- a/Specifications/Dynamics/Spec_ControlledAttitude.md +++ b/Specifications/Dynamics/Spec_ControlledAttitude.md @@ -39,7 +39,7 @@ - EARTH_CENTER_POINTING = 2 - VELOCITY_DIRECTION_POINTING = 3 - ORBIT_NORMAL_POINTING = 4 - - orbit normal $`n`$ is defined as $`n=r\times v`$, where $`r`$ is radial direction and $`v`$ is velocity direction. + - orbit normal $n$ is defined as $n=r\times v$, where $r$ is radial direction and $v$ is velocity direction. ## 2. Explanation of Algorithm @@ -102,21 +102,21 @@ 2. inputs and outputs - inputs - - the main direction of the target object in ECI frame $`t_m^i`$ - - the sub direction of the target object in ECI frame $`t_s^i`$ - - the main controlled direction in the body frame $`d_m^b`$ - - the sub controlled direction in the body frame $`d_s^b`$ + - the main direction of the target object in ECI frame $t_m^i$ + - the sub direction of the target object in ECI frame $t_s^i$ + - the main controlled direction in the body frame $d_m^b$ + - the sub controlled direction in the body frame $d_s^b$ - outputs - quaternion_i2b 3. algorithm - - Firstly, the $`DCM_{t2i}`$, which is the frame transformation from the target frame to the inertial frame, is calculated using $`t_m^i`$ and $`t_s^i`$ with `CalcDCM`. - - Next, the $`DCM_{t2b}`$, which is the frame transformation from the target frame to the body-fixed frame, is calculated using $`d_m^b`$ and $`d_s^b`$ with `CalcDCM`. + - Firstly, the $DCM_{t2i}$, which is the frame transformation from the target frame to the inertial frame, is calculated using $t_m^i$ and $t_s^i$ with `CalcDCM`. + - Next, the $DCM_{t2b}$, which is the frame transformation from the target frame to the body-fixed frame, is calculated using $d_m^b$ and $d_s^b$ with `CalcDCM`. - Finally, both DCMs are combined as ```math DCM_{i2b} = DCM_{t2b} \cdot DCM_{t2i}' ``` - and `quaternion_i2b` is calculated from the $`DCM_{i2b}`$. + and `quaternion_i2b` is calculated from the $DCM_{i2b}$. 4. note - NA @@ -129,17 +129,17 @@ 2. inputs and outputs - inputs - - the main direction in the frame $`a`$ : $`d_m^a`$ - - the sub direction in the frame $`a`$ : $`d_s^a`$ + - the main direction in the frame $a$ : $d_m^a$ + - the sub direction in the frame $a$ : $d_s^a$ - outputs - - Coordinate transform matrix from the new frame to the original frame $`a`$ + - Coordinate transform matrix from the new frame to the original frame $a$ 3. algorithm - The first basis vector of the new frame is defined as the main direction. ```math e_1=d_m^a ``` - - The second basis vector needs to be direct to the sub direction, but it should be vertical with $`e_1`$. + - The second basis vector needs to be direct to the sub direction, but it should be vertical with $e_1$. ```math e_2 = \frac{(e_1\times d_s^a) \times e_1}{|(e_1\times d_s^a) \times e_1|} ``` diff --git a/Specifications/Dynamics/Spec_EnckeMethod.md b/Specifications/Dynamics/Spec_EnckeMethod.md index 8931b6ed..a6d988a2 100644 --- a/Specifications/Dynamics/Spec_EnckeMethod.md +++ b/Specifications/Dynamics/Spec_EnckeMethod.md @@ -25,10 +25,10 @@ 2. Inputs and outputs - Input - - $`\mu`$ : The standard gravitational parameter of the central body - - $`t`$ : Time in Julian day - - $`\boldsymbol{r}_{i}`$ : Initial position in the inertial frame - - $`\boldsymbol{v}_{i}`$ : Initial velocity in the inertial frame + - $\mu$ : The standard gravitational parameter of the central body + - $t$ : Time in Julian day + - $\boldsymbol{r}_{i}$ : Initial position in the inertial frame + - $\boldsymbol{v}_{i}$ : Initial velocity in the inertial frame - Output - The reference orbit - The difference is set as zero @@ -40,17 +40,17 @@ 1. Overview - This function is the main algorithm of Encke's method and calculates the orbit of the spacecraft. - The method separates the orbit to the reference and the difference. The reference is calculated with the Kepler orbit method as a two-body problem, and the difference is calculated, including the disturbances. - - $`\boldsymbol{r}_{ref}`$ : Reference orbit - - $`\boldsymbol{\delta}`$ : Difference + - $\boldsymbol{r}_{ref}$ : Reference orbit + - $\boldsymbol{\delta}$ : Difference - Please refer to the references to learn the original idea of Encke's method. 2. Inputs and outputs - Input - - $`\boldsymbol{a}_d`$ : Acceleration - - $`t`$ : Current time + - $\boldsymbol{a}_d$ : Acceleration + - $t$ : Current time - Output - - $`\boldsymbol{r}_{i}`$ : Initial position in the inertial frame - - $`\dot{\boldsymbol{r}}_{i}`$ : Initial velocity in the inertial frame + - $\boldsymbol{r}_{i}$ : Initial position in the inertial frame + - $\dot{\boldsymbol{r}}_{i}$ : Initial velocity in the inertial frame 3. Algorithm 1. Rectification @@ -112,7 +112,7 @@ - The following figure shows the difference between orbit derived with RK4 mode and Encke mode with all disturbances. - - The error is larger than the non-disturbance case, but the $`10^4 [m]`$ error between the Encke method and the Cowell method is compatible with the ref[2] when using the RK4 in LEO. We confirmed that the Encke propagation mode is correct when all disturbances are included. + - The error is larger than the non-disturbance case, but the $10^4 [m]$ error between the Encke method and the Cowell method is compatible with the ref[2] when using the RK4 in LEO. We confirmed that the Encke propagation mode is correct when all disturbances are included.
diff --git a/Specifications/Dynamics/Spec_KeplerOrbit.md b/Specifications/Dynamics/Spec_KeplerOrbit.md index 51a63b01..003689d5 100644 --- a/Specifications/Dynamics/Spec_KeplerOrbit.md +++ b/Specifications/Dynamics/Spec_KeplerOrbit.md @@ -24,12 +24,12 @@ - Defined by the orbital elements - Select `init_mode_kepler = INIT_OE` - Set the value of the following orbital elements - - $`a`$ : Semi major axis [m] - - $`e`$ : Eccentricity - - $`i`$ : Inclination [rad] - - $`\Omega`$ : Right Ascension of the Ascending Node (RAAN) [rad] - - $`\omega`$ : Argument of Perigee [rad] - - $`t_{epoch}`$ : Epoch [julian day] + - $a$ : Semi major axis [m] + - $e$ : Eccentricity + - $i$ : Inclination [rad] + - $\Omega$ : Right Ascension of the Ascending Node (RAAN) [rad] + - $\omega$ : Argument of Perigee [rad] + - $t_{epoch}$ : Epoch [julian day] ## 2. Explanation of Algorithm @@ -42,10 +42,10 @@ 2. Inputs and outputs - Input - Orbital Elements - - $`\mu`$ : The standard gravitational parameter of the central body + - $\mu$ : The standard gravitational parameter of the central body - Output - - $`n`$ : Mean motion - - $`R_{p2i}`$ : Frame conversion matrix from in-plane position to the inertial frame position + - $n$ : Mean motion + - $R_{p2i}$ : Frame conversion matrix from in-plane position to the inertial frame position 3. Algorithm - Mean motion @@ -77,21 +77,21 @@ 2. Inputs and outputs - Input - - $`t`$ : Time in Julian day + - $t$ : Time in Julian day - Orbital Elements - Constants - Output - - $`\boldsymbol{r}_{i}`$ : Position in the inertial frame - - $`\boldsymbol{v}_{i}`$ : Velocity in the inertial frame + - $\boldsymbol{r}_{i}$ : Position in the inertial frame + - $\boldsymbol{v}_{i}$ : Velocity in the inertial frame 3. Algorithm - - Calculate mean anomaly $`l`$[rad] + - Calculate mean anomaly $l$[rad] ```math l = n * (t-t_{epoch}) ``` - - Calculate eccentric anomaly $`u`$[rad] by solving the Kepler Equation + - Calculate eccentric anomaly $u$[rad] by solving the Kepler Equation - Details are described in `KeplerOrbit::SolveKeplerFirstOrder` - - Calculate two dimensional position $`x^{\*}`$, $`y^{\*}`$ and velocity $`\dot{x}^{\*}`$, $`\dot{y}^{\*}`$ in the orbital plane + - Calculate two dimensional position $x^{\*}$, $y^{\*}$ and velocity $\dot{x}^{\*}$, $\dot{y}^{\*}$ in the orbital plane ```math x^* = a(\cos{u}-e)\\ y^* = a\sqrt{1-e^2}\sin{u}\\ @@ -117,22 +117,22 @@ 2. Inputs and outputs - Input - - $`e`$ : eccentricity - - $`l`$ : mean anomaly [rad] - - $`\epsilon`$ : threshold for convergence [rad] + - $e$ : eccentricity + - $l$ : mean anomaly [rad] + - $\epsilon$ : threshold for convergence [rad] - Limit of iteration - Output - - $`u`$ : eccentric anomaly [rad] + - $u$ : eccentric anomaly [rad] 3. Algorithm - Set the initial value of eccentric anomaly as follows - - $`u_0=l`$ - - Calculate $`u_{n+1}`$ with the following equation + - $u_0=l$ + - Calculate $u_{n+1}$ with the following equation ```math u_{n+1} = l + e\sin{u_n} ``` - Iterate the calculation until the following conditions are satisfied - - $`|u_{n+1} - u_{n}| < \epsilon`$ + - $|u_{n+1} - u_{n}| < \epsilon$ - The iteration number over the limit of iteration 4. `OrbitalElements::CalcOeFromPosVel` function @@ -141,52 +141,52 @@ 2. Inputs and outputs - Input - - $`\mu`$ : The standard gravitational parameter of the central body - - $`t`$ : Time in Julian day - - $`\boldsymbol{r}_{i}`$ : Initial position in the inertial frame - - $`\boldsymbol{v}_{i}`$ : Initial velocity in the inertial frame + - $\mu$ : The standard gravitational parameter of the central body + - $t$ : Time in Julian day + - $\boldsymbol{r}_{i}$ : Initial position in the inertial frame + - $\boldsymbol{v}_{i}$ : Initial velocity in the inertial frame - Output - orbital element 3. Algorithm - - $`\boldsymbol{h}_{i}`$ : Angular momentum vector of the orbit + - $\boldsymbol{h}_{i}$ : Angular momentum vector of the orbit ```math \boldsymbol{h}_{i} = \boldsymbol{r}_{i} \times \boldsymbol{v}_{i} ``` - - $`a`$ : Semi-major axis + - $a$ : Semi-major axis ```math a = \frac{\mu}{2\frac{\mu}{r} - v^2} ``` - - $`i`$ : Inclination + - $i$ : Inclination ```math i = \cos^{-1}{h_z} ``` - - $`\Omega`$ : Right Ascension of the Ascending Node (RAAN) - - Note: This equation is not support $`i = 0`$ case. + - $\Omega$ : Right Ascension of the Ascending Node (RAAN) + - Note: This equation is not support $i = 0$ case. ```math \Omega = \sin^{-1}\left(\frac{h_x}{\sqrt{h_x^2 + h_y^2}}\right) ``` - - $`x_{p}, y_{p}`$ : Position in the orbital plane + - $x_{p}, y_{p}$ : Position in the orbital plane ```math x_{p} = r_{ix} \cos{\Omega} + r_{iy} \sin{\Omega};\\ y_{p} = (-r_{ix} \sin{\Omega} + r_{iy} \cos{\Omega})\cos{i} + r_{iz}\sin{i}; ``` - - $`\dot{x}_{p}, \dot{y}_{p}`$ : Velocity in the orbital plane + - $\dot{x}_{p}, \dot{y}_{p}$ : Velocity in the orbital plane ```math \dot{x}_{p} = v_{ix} \cos{\Omega} + v_{iy} \sin{\Omega};\\ \dot{y}_{p} = (-v_{ix} \sin{\Omega} + v_{iy} \cos{\Omega})\cos{i} + v_{iz}\sin{i}; ``` - - $`e`$ : Eccentricity + - $e$ : Eccentricity ```math c_1 = \frac{h}{\mu}\dot{y}_p - \frac{x_p}{r}\\ c_2 = -\frac{h}{\mu}\dot{x}_p - \frac{y_p}{r}\\ e = \sqrt{c_1^2 + c_2^2} ``` - - $`\omega`$ : Argument of Perigee + - $\omega$ : Argument of Perigee ```math \omega = \tan^{-1}\left(\frac{c_2}{c_1}\right) ``` - - $`t_{epoch}`$ : Epoch [Julian day] + - $t_{epoch}$ : Epoch [Julian day] ```math f = \tan^{-1}\left(\frac{y_p}{x_p}\right) - \omega\\ u = \tan^{-1}\frac{\frac{r \sin{f}}{\sqrt{1-e^2}}}{r\cos{f} + ae}\\ diff --git a/Specifications/Dynamics/Spec_RelativeOrbit.md b/Specifications/Dynamics/Spec_RelativeOrbit.md index 6aaaa936..38eb5e84 100644 --- a/Specifications/Dynamics/Spec_RelativeOrbit.md +++ b/Specifications/Dynamics/Spec_RelativeOrbit.md @@ -30,14 +30,14 @@ - Set up the configuration of the `[ORBIT]` section in the `Sat.ini`. + Set `propagate_mode =2` to use the relative orbit propagation + Choose `update_method`. - * `update_method = 0` means the orbit is updated with the propagation of the relative dynamics equation($`\dot{\boldsymbol{x}}=\boldsymbol{Ax}+\boldsymbol{Bu}`$, i.e., Hill equation). - * `update_method = 1` means the orbit is updated with the STM($`\boldsymbol{x}(t)=\boldsymbol{\Phi}(t,t_0)\boldsymbol{x}(t_0)`$, i.e., Clohessy-Wiltshire solution). + * `update_method = 0` means the orbit is updated with the propagation of the relative dynamics equation($\dot{\boldsymbol{x}}=\boldsymbol{Ax}+\boldsymbol{Bu}$, i.e., Hill equation). + * `update_method = 1` means the orbit is updated with the STM($\boldsymbol{x}(t)=\boldsymbol{\Phi}(t,t_0)\boldsymbol{x}(t_0)$, i.e., Clohessy-Wiltshire solution). + When you choose `update_method = 0`, set `relative_dynamics_model_type`. + When you choose `update_method = 1`, set `stm_model_type`. + Set the initial relative position of the satellite in the LVLH frame. LVLH frame definition is: - * $`\boldsymbol{x}`$ is a direction vector from the reference satellite ("chief" in the figure) radially outward. - * The direction of $`\boldsymbol{z}`$ corresponds to the angular momentum vector of the reference satellite orbit. - * The direction of $`\boldsymbol{y}`$ is determined by $`\boldsymbol{z}\times\boldsymbol{x}`$. + * $\boldsymbol{x}$ is a direction vector from the reference satellite ("chief" in the figure) radially outward. + * The direction of $\boldsymbol{z}$ corresponds to the angular momentum vector of the reference satellite orbit. + * The direction of $\boldsymbol{y}$ is determined by $\boldsymbol{z}\times\boldsymbol{x}$.
@@ -79,7 +79,7 @@ * `current_jd` - The initial Julian day * `mu` - - The gravity constant of the reference celestial body $`\mu`$ + - The gravity constant of the reference celestial body $\mu$ * `timestep` - RK4 propagation timestep * `wgs` @@ -105,7 +105,7 @@ * `reference_sat_orbit` - The orbit of the reference satellite * `mu` - - The gravity constant $`\mu`$ + - The gravity constant $\mu$ + output * none @@ -125,7 +125,7 @@ * `reference_sat_orbit` - The orbit of the reference satellite * `mu` - - The gravity constant $`\mu`$ + - The gravity constant $\mu$ * `elapsed_sec` - Elapsed simulation time + output @@ -144,21 +144,21 @@ 2. inputs and outputs + input * `orbit_radius` - - Radius of the reference satellite orbit $`R`$ + - Radius of the reference satellite orbit $R$ * `mu` - - The gravity constant $`\mu`$ + - The gravity constant $\mu$ + output * `system_matrix` - system matrix 3. algorithm - + The mean motion of the reference orbit ($`n`$) is calculated as follows: + + The mean motion of the reference orbit ($n$) is calculated as follows: ```math n=\sqrt{\frac{\mu}{R^3}} ``` - + Then, the system matrix ($`\boldsymbol{A}`$) is calculated as follows: + + Then, the system matrix ($\boldsymbol{A}$) is calculated as follows: ```math \boldsymbol{A}= @@ -183,9 +183,9 @@ 2. inputs and outputs + input * `orbit_radius` - - Radius of the reference satellite orbit $`R`$ + - Radius of the reference satellite orbit $R$ * `mu` - - The gravity constant $`\mu`$ + - The gravity constant $\mu$ * `elapsed_sec` - Elapsed simulation time + output @@ -193,13 +193,13 @@ - system matrix 3. algorithm - + The mean motion of the reference orbit ($`n`$) is calculated as follows: + + The mean motion of the reference orbit ($n$) is calculated as follows: ```math n=\sqrt{\frac{\mu}{R^3}} ``` - + Then, the system matrix ($`\boldsymbol{A}`$) is calculated as follows: + + Then, the system matrix ($\boldsymbol{A}$) is calculated as follows: ```math \boldsymbol{\Phi}(t,t0)= @@ -228,7 +228,7 @@ - `Sat0.ini` - `Sat1.ini` 2. initial values - - The orbit of the reference satellite (Sat0) is GEO. The initial position of the satellite (Sat1) is $`(0\mathrm{m}, 100\mathrm{m}, 0\mathrm{m})^\mathrm{T}`$ in LVLH frame. The orbit was propagated for 86400 sec (the period of GEO). + - The orbit of the reference satellite (Sat0) is GEO. The initial position of the satellite (Sat1) is $(0\mathrm{m}, 100\mathrm{m}, 0\mathrm{m})^\mathrm{T}$ in LVLH frame. The orbit was propagated for 86400 sec (the period of GEO). - `Sat0.ini` ``` diff --git a/Specifications/Dynamics/Spec_Rk4Orbit.md b/Specifications/Dynamics/Spec_Rk4Orbit.md index 5ab99680..4bfdcc73 100644 --- a/Specifications/Dynamics/Spec_Rk4Orbit.md +++ b/Specifications/Dynamics/Spec_Rk4Orbit.md @@ -28,8 +28,8 @@ ## 2. Explanation of Algorithm 1. `Propagate` function - - The position and velocity of the satellite are updated by using RK4. As the input of RK4, the six-state variables are set. These state variables are the three-dimensional position [$`x`$, $`y`$ ,$`z`$] and three-dimensional velocity [$`v_x`$, $`v_y`$, $`v_z`$] at the inertial coordinate. Here, the inertial coordinate is decided by the `PlanetSelect.ini` - - As the force which works to the satellite motion is the external acceleration [$`a_x`$,$`a_y`$,$`a_z`$] calculated from the disturbance class or thruster class and the gravity force from the center planet, which is defined in `PlanetSelect.ini`. As a summary, the orbit is calculated as the following equation. + - The position and velocity of the satellite are updated by using RK4. As the input of RK4, the six-state variables are set. These state variables are the three-dimensional position [$x$, $y$ ,$z$] and three-dimensional velocity [$v_x$, $v_y$, $v_z$] at the inertial coordinate. Here, the inertial coordinate is decided by the `PlanetSelect.ini` + - As the force which works to the satellite motion is the external acceleration [$a_x$,$a_y$,$a_z$] calculated from the disturbance class or thruster class and the gravity force from the center planet, which is defined in `PlanetSelect.ini`. As a summary, the orbit is calculated as the following equation. ```math \dot{x} = v_x\\ \dot{y} = v_y\\ @@ -75,8 +75,8 @@
- - In the cases of `OrbitPropagateStepTimeSec=0.1` and `OrbitPropagateStepTimeSec=1`, the error is kept within $`10^{-6}`$ order. However, once the error grows, it will get bigger and bigger. - - In the case of `OrbitPropagateStepTimeSec=10`, the error quickly grows up to $`10^{-4}`$ order. + - In the cases of `OrbitPropagateStepTimeSec=0.1` and `OrbitPropagateStepTimeSec=1`, the error is kept within $10^{-6}$ order. However, once the error grows, it will get bigger and bigger. + - In the case of `OrbitPropagateStepTimeSec=10`, the error quickly grows up to $10^{-4}$ order.
orbit_steptimesec_1_longterm @@ -99,7 +99,7 @@
- As this figure shows, the initial values in the result are slightly different from the input. - + In the .sa files, the initial values of $`x, y, z, v_x, v_y, v_z`$ are converted into elements of orbit and stored. The error might occur in the process of this conversion. + + In the .sa files, the initial values of $x, y, z, v_x, v_y, v_z$ are converted into elements of orbit and stored. The error might occur in the process of this conversion. ## 4. References diff --git a/Specifications/Environment/Spec_Atmosphere.md b/Specifications/Environment/Spec_Atmosphere.md index 48bb8cb9..057fc128 100644 --- a/Specifications/Environment/Spec_Atmosphere.md +++ b/Specifications/Environment/Spec_Atmosphere.md @@ -53,7 +53,7 @@ ```math \rho \left(h \right) = \rho_0 \exp \left(-\frac{h - h_0}{H} \right) ``` - where $`h_0`$ is a reference height, $`\rho_0`$ is an air density at the reference height, and $`H`$ is a scale height. These parameters are referred from the table shown in References [2]. + where $h_0$ is a reference height, $\rho_0$ is an air density at the reference height, and $H$ is a scale height. These parameters are referred from the table shown in References [2]. + NRLMSISE00 1. Read Space weather table diff --git a/Specifications/Environment/Spec_CelestialRotation.md b/Specifications/Environment/Spec_CelestialRotation.md index 58ad4fbb..d385e4ed 100644 --- a/Specifications/Environment/Spec_CelestialRotation.md +++ b/Specifications/Environment/Spec_CelestialRotation.md @@ -11,11 +11,11 @@ 3. How to use - Make an instance of the `CelestialRotation` class in `CelestialInformation` class. - Select rotation mode in `SampleSimBase.ini` - - `Idle` : no motion ($`\mathrm{DCM_{ECItoECEF}}`$ is set as the unit matrix) + - `Idle` : no motion ($\mathrm{DCM_{ECItoECEF}}$ is set as the unit matrix) - If the rotation mode input is neither `Full` nor `Simple`, the `Idle` mode is set. - - `Simple` : axial rotation only ($`\mathrm{DCM_{ECItoECEF}} = \bf{R}`$) - - `Full` : Precession and Nutation are taken into account ($`\mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P}`$) - - $`\bf{R}`$, $`\bf{N}`$, $`\bf{P}`$ stand for the DCM of axial rotation, nutation, precession, respectively. + - `Simple` : axial rotation only ($\mathrm{DCM_{ECItoECEF}} = \bf{R}$) + - `Full` : Precession and Nutation are taken into account ($\mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P}$) + - $\bf{R}$, $\bf{N}$, $\bf{P}$ stand for the DCM of axial rotation, nutation, precession, respectively. ## 2. Explanation of Algorithm @@ -28,7 +28,7 @@ ```math \mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P} ``` - - where $`\bf{R}`$, $`\bf{N}`$, $`\bf{P}`$ stand for the DCM of axial rotation, nutation, precession, respectively. + - where $\bf{R}$, $\bf{N}$, $\bf{P}$ stand for the DCM of axial rotation, nutation, precession, respectively. 2. inputs and outputs - Input @@ -51,8 +51,8 @@ - JDJ2000 = 2451545.0 [day] - JC = 36525 [day/century] - - By using tTT, we get the DCM of precession ($`\bf{P}`$) and nutation ($`\bf{N}`$) with `Precession` and `Nutation` functions. - - $`\varepsilon`$，$`\Delta \varepsilon`$，$`\Delta \psi`$ are calculated in `Nutation` function. + - By using tTT, we get the DCM of precession ($\bf{P}$) and nutation ($\bf{N}$) with `Precession` and `Nutation` functions. + - $\varepsilon$，$\Delta \varepsilon$，$\Delta \psi$ are calculated in `Nutation` function. ```math \mathrm{E_q} = \Delta \psi \cos{(\varepsilon + \Delta \varepsilon)} \\ @@ -62,7 +62,7 @@ - GAST is calculated from julian date in `gstime` function in `src/Library/sgp4/sgp4unit.h`. - - By using GMST, We get the DCM of axial rotation ($`\bf{R}`$) with the `Rotation` function. The coordinate transformation from ECI to ECEF is calculated. + - By using GMST, We get the DCM of axial rotation ($\bf{R}$) with the `Rotation` function. The coordinate transformation from ECI to ECEF is calculated. ```math \mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P} ``` @@ -79,7 +79,7 @@ - Input - Greenwich Apparent Sidereal Time (GAST) - Output - - the DCM of axial rotation ($`\bf{R}`$) + - the DCM of axial rotation ($\bf{R}$) 3. algorithm ```math @@ -100,7 +100,7 @@ - Input - Julian century for terrestrial time (tTT) - Output - - the DCM of precession ($`\bf{P}`$) + - the DCM of precession ($\bf{P}$) 3. algorithm - Precession angles are calculated as follows. @@ -137,10 +137,10 @@ - Input - Julian century for terrestrial time (tTT) - Output - - Return: the DCM of precession ($`\bf{N}`$) - - $`\varepsilon`$: mean obliquity of the ecliptic - - $`\Delta \varepsilon`$: nutation in obliquity - - $`\Delta \psi`$: nutation in longitude + - Return: the DCM of precession ($\bf{N}$) + - $\varepsilon$: mean obliquity of the ecliptic + - $\Delta \varepsilon$: nutation in obliquity + - $\Delta \psi$: nutation in longitude 3. algorithm Delaunay angles are calculated as follows. @@ -157,9 +157,9 @@ - l' : mean anomaly of the sun - F : mean argument of latitude of the moon - D : mean elongation of the moon from the sun - - $`\Omega`$ : mean longitude of ascending node of the moon + - $\Omega$ : mean longitude of ascending node of the moon - $`\varepsilon`$ and $`\Delta \varepsilon`$ and $`\Delta \psi`$ are calculated as follows. + $\varepsilon$ and $\Delta \varepsilon$ and $\Delta \psi$ are calculated as follows. ```math \varepsilon = 23^\circ26'21".448 - 46".8150\mathrm{tTT} - 0".00059\mathrm{tTT}^2 + 0".001813\mathrm{tTT}^3 \\ @@ -167,7 +167,7 @@ \Delta \psi = -17.206\sin{\Omega} - 1.317\sin{2L'} + 0.207\sin{2\Omega} - 0.228\sin{2L} + 0.148\sin{l'}+0.071\sin{l}-0.052\sin{(2L'+l')} - 0.030\sin{(2L+l)}+0.022\sin{(2L'-l')} \, [\mathrm{arcsec}] \\ ``` - where $`L = F + \Omega`$，$`L' = L - D`$ + where $L = F + \Omega$，$L' = L - D$ ```math \bf{N} = @@ -191,10 +191,10 @@ ## 3. Results of verifications -1. $`\mathrm{DCM_{ECItoECEF}}`$ calculation in `Update` function +1. $\mathrm{DCM_{ECItoECEF}}$ calculation in `Update` function 1. overview - - The $`\mathrm{DCM_{ECItoECEF}}`$ calculation is compared with [Matlab's dcmeci2ecef function](https://jp.mathworks.com/help/aerotbx/ug/dcmeci2ecef.html#d123e38055) + - The $\mathrm{DCM_{ECItoECEF}}$ calculation is compared with [Matlab's dcmeci2ecef function](https://jp.mathworks.com/help/aerotbx/ug/dcmeci2ecef.html#d123e38055) 2. conditions for the verification 1. input value diff --git a/Specifications/Environment/Spec_MagEnvironment.md b/Specifications/Environment/Spec_MagEnvironment.md index cefd6097..61a6190a 100644 --- a/Specifications/Environment/Spec_MagEnvironment.md +++ b/Specifications/Environment/Spec_MagEnvironment.md @@ -33,8 +33,8 @@ ## 2. Explanation of Algorithm -+ IGRF calculates the magnetic field based on a spherical harmonic expansion of magnetic scalar potential $`V`$. - + The coefficients of the spherical harmonic expansion of $`V`$ is updated by [IAGA](https://www.ngdc.noaa.gov/IAGA/vmod/index.html). The latest version is the 13th generation, and S2E uses this version. ++ IGRF calculates the magnetic field based on a spherical harmonic expansion of magnetic scalar potential $V$. + + The coefficients of the spherical harmonic expansion of $V$ is updated by [IAGA](https://www.ngdc.noaa.gov/IAGA/vmod/index.html). The latest version is the 13th generation, and S2E uses this version. + Please refer [here](https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html#:~:text=The%20International%20Association%20of%20Geomagnetism,interior%2C%20its%20crust%20and%20its) for the details of IGRF. ## 3. Verification diff --git a/Specifications/Environment/Spec_SRPEnvironment.md b/Specifications/Environment/Spec_SRPEnvironment.md index 8cf88546..d3ea762d 100644 --- a/Specifications/Environment/Spec_SRPEnvironment.md +++ b/Specifications/Environment/Spec_SRPEnvironment.md @@ -19,10 +19,10 @@ - `CalcTruePressure`: Return solar pressure (N/m2) with eclipse effect for SRP disturbance calculation. - `CalcPowerDensity`: Return solar power density (W/m2) with eclipse effect for Electrical Power System calculation. - `GetPressure`: Return solar pressure (N/m2) without eclipse effect. - - `GetShadowFunction`: Return shadow function $`\nu`$. - - When the spacecraft is in umbra, $`\nu=0`$. - - When the spacecraft is in sunlight, $`\nu=1`$. - - When the spacecraft is in penumbra, $`0<\nu<1`$. + - `GetShadowFunction`: Return shadow function $\nu$. + - When the spacecraft is in umbra, $\nu=0$. + - When the spacecraft is in sunlight, $\nu=1$. + - When the spacecraft is in penumbra, $0<\nu<1$. @@ -34,14 +34,14 @@ 2. inputs and outputs - Constants - - Solar constant: $`P_{\odot} = 1366`$ W/m2 - - Speed of light: $`c = 299792458`$ m/s - - Astronomical Unit: $`AU = 149597870700`$ m + - Solar constant: $P_{\odot} = 1366$ W/m2 + - Speed of light: $c = 299792458$ m/s + - Astronomical Unit: $AU = 149597870700$ m - Input variables - - The sun position in the body-fixed frame of the spacecraft: $`\boldsymbol{r}_{\odot-sc}`$ m - - Unbold $`r_{\odot-sc}`$ is the norm of $`\boldsymbol{r}_{\odot-sc}`$ + - The sun position in the body-fixed frame of the spacecraft: $\boldsymbol{r}_{\odot-sc}$ m + - Unbold $r_{\odot-sc}$ is the norm of $\boldsymbol{r}_{\odot-sc}$ - Output - - Solar radiation pressure: $`P_{sc}`$ N/m2 + - Solar radiation pressure: $P_{sc}$ N/m2 3. algorithm ```math @@ -57,11 +57,11 @@ 2. inputs and outputs - Constants - - Radius of the earth: $`r_{\oplus}=6378137`$ m - - Radius of the sun: $`r_{\odot}=6.96\times10^{8}`$ m + - Radius of the earth: $r_{\oplus}=6378137$ m + - Radius of the sun: $r_{\odot}=6.96\times10^{8}$ m - Input variables - - The sun position in the body-fixed frame of the spacecraft: $`\boldsymbol{r}_{\odot-sc}`$ m - - The earth position in the body-fixed frame of the spacecraft: $`\boldsymbol{r}_{\oplus-sc}`$ m + - The sun position in the body-fixed frame of the spacecraft: $\boldsymbol{r}_{\odot-sc}$ m + - The earth position in the body-fixed frame of the spacecraft: $\boldsymbol{r}_{\oplus-sc}$ m - Output - none @@ -82,24 +82,24 @@ 2. inputs and outputs - Input - - The apparent radius of the Sun: $`A_{\odot}`$ - - The apparent radius of the Earth: $`A_{\oplus}`$ - - The apparent separation of the centers of the Sun and the Earth: $`\delta`$ - - The angle between the center of the Sun and the common chord: $`x`$ - - The length of the common chord of the apparent solar disk and apparent celestial disk: $`y`$ + - The apparent radius of the Sun: $A_{\odot}$ + - The apparent radius of the Earth: $A_{\oplus}$ + - The apparent separation of the centers of the Sun and the Earth: $\delta$ + - The angle between the center of the Sun and the common chord: $x$ + - The length of the common chord of the apparent solar disk and apparent celestial disk: $y$ - Output - - The shadow function: $`\nu`$ + - The shadow function: $\nu$ 3. algorithm - - If the occultation is total, then $`\nu=0`$. - - If the occultation is partial but maximum, then $`\nu=1-\left(\frac{A_{\oplus}}{A_{\odot}}\right)^2`$ - - If the occultation is partial, then $`\nu = 1-\frac{S}{\pi A^2_{\odot}}`$ + - If the occultation is total, then $\nu=0$. + - If the occultation is partial but maximum, then $\nu=1-\left(\frac{A_{\oplus}}{A_{\odot}}\right)^2$ + - If the occultation is partial, then $\nu = 1-\frac{S}{\pi A^2_{\odot}}$ - S is given by the following calculation. ```math S=A_{\odot}^2\arccos\left(\frac{x}{A_{\odot}}\right)+A_{\oplus}^2\arccos\left(\frac{\delta-x}{A_{\oplus}}\right)-\delta\cdot y ``` - - In other cases, since it means that no occultation takes place, then $`\nu=1`$. + - In other cases, since it means that no occultation takes place, then $\nu=1$. ## 3. Results of verifications From 00b771e5a80d4a9635520dbd5d23d8f7e17b74b6 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Tue, 15 Nov 2022 16:24:28 +0100 Subject: [PATCH 02/19] test fix newline in math block --- Specifications/Disturbance/Spec_GeoPotential.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 2bcfd143..5b3f3c6c 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -84,9 +84,9 @@ For unnormalized algorithms, see chapter 3.2.4 of [Satellite Orbits](https://www For normalized algorithm, we use following normalizing relation for Legendre polynomials, ```math -\bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m}\\ -\bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m}\\ -\bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m}\\ +\bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m} \\ +\bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m} \\ +\bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m} \\ ``` The recursion calculation of V and W can be changed to a normalized version as follows From 53d81f3da89a27ca6eed178810e3873748426a6e Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Tue, 15 Nov 2022 16:26:14 +0100 Subject: [PATCH 03/19] test fix newline in math block --- Specifications/Disturbance/Spec_GeoPotential.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 5b3f3c6c..7fd35268 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -84,9 +84,9 @@ For unnormalized algorithms, see chapter 3.2.4 of [Satellite Orbits](https://www For normalized algorithm, we use following normalizing relation for Legendre polynomials, ```math -\bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m} \\ -\bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m} \\ -\bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m} \\ +\bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m} \newline +\bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m} \newline +\bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m} \newline ``` The recursion calculation of V and W can be changed to a normalized version as follows From a42df9d42310f76abb4085eff34d6906d04806b4 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Tue, 15 Nov 2022 16:35:06 +0100 Subject: [PATCH 04/19] test fix newline in math block --- Specifications/Disturbance/Spec_GeoPotential.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 7fd35268..12efb6c4 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -84,9 +84,11 @@ For unnormalized algorithms, see chapter 3.2.4 of [Satellite Orbits](https://www For normalized algorithm, we use following normalizing relation for Legendre polynomials, ```math -\bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m} \newline -\bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m} \newline -\bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m} \newline +\begin{align*} + \bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m}\\ + \bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m}\\ + \bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m}\\ +\end{align*} ``` The recursion calculation of V and W can be changed to a normalized version as follows From 91d0fc1959107a16557831a679cc1970674cb21e Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:06:25 +0100 Subject: [PATCH 05/19] fix multiple line equation in geopotential --- .../Disturbance/Spec_GeoPotential.md | 90 ++++++++++++------- 1 file changed, 59 insertions(+), 31 deletions(-) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 12efb6c4..13842ae4 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -84,63 +84,81 @@ For unnormalized algorithms, see chapter 3.2.4 of [Satellite Orbits](https://www For normalized algorithm, we use following normalizing relation for Legendre polynomials, ```math -\begin{align*} +\begin{align} \bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m}\\ \bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m}\\ \bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m}\\ -\end{align*} +\end{align} ``` The recursion calculation of V and W can be changed to a normalized version as follows ```math -N_{m,m}\bar{V}_{m,m}=(2m-1)(\frac{xR_{e}}{r^2}N_{m-1,m-1}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}N_{m-1,m-1}\bar{W}_{m-1,m-1})\\ -\bar{V}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\\ -\bar{W}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\\ +\begin{align} + N_{m,m}\bar{V}_{m,m}=(2m-1)(\frac{xR_{e}}{r^2}N_{m-1,m-1}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}N_{m-1,m-1}\bar{W}_{m-1,m-1})\\ + \bar{V}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\\ + \bar{W}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\\ +\end{align} ``` ```math -N_{n,m}\bar{V}_{n,m}=\frac{2n-1}{n-m}(\frac{zR_{e}}{r^2}N_{n-1,m}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}N_{n-2,m}\bar{V}_{n-2,m})\\ -\bar{V}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{V}_{n-2,m}\\ -\bar{W}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{W}_{n-2,m}\\ +\begin{align} + N_{n,m}\bar{V}_{n,m}=\frac{2n-1}{n-m}(\frac{zR_{e}}{r^2}N_{n-1,m}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}N_{n-2,m}\bar{V}_{n-2,m})\\ + \bar{V}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{V}_{n-2,m}\\ + \bar{W}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{W}_{n-2,m}\\ +\end{align} ``` The recurrence relation of normalize function can be expressed as follows ```math -N_{0,0}=1\\ -N_{m,m}=(2m-1)\sqrt{\frac{2m}{2m+1}}N_{m-1,m-1}\quad(m\geq1) +\begin{align} + N_{0,0}=1\\ + N_{m,m}=(2m-1)\sqrt{\frac{2m}{2m+1}}N_{m-1,m-1}\quad(m\geq1) +\end{align} ``` ```math -N_{n,m}=\sqrt{\frac{2n-1}{2n+1}}\sqrt{\frac{n+m}{n-m}}N_{n-1,m}\quad(n\geq1,0\leq m \leq n) +\begin{align} + N_{n,m}=\sqrt{\frac{2n-1}{2n+1}}\sqrt{\frac{n+m}{n-m}}N_{n-1,m}\quad(n\geq1,0\leq m \leq n) +\end{align} ``` So, the divisions of the normalized functions are described as follows ```math -\frac{N_{0,0}}{N_{1,1}}=\sqrt{2m+1}=\sqrt{3}\\ -\frac{N_{m-1,m-1}}{N_{m,m}}=\frac{1}{2m-1}\sqrt{\frac{2m+1}{2m}}\quad(m\geq2)\\ -\frac{N_{n-1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}\quad(n\geq1,0\leq m \leq n)\\ -\frac{N_{n-2,m}}{N_{n,m}}=\frac{N_{n-2,m}}{N_{n-1,m}}\frac{N_{n-1,m}}{N_{n,m}}\\ -\frac{N_{n-2,m}}{N_{n,m}}=\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\frac{N_{n-1,m}}{N_{n,m}}\quad(n\geq2,0\leq m \leq n) +\begin{align} + \frac{N_{0,0}}{N_{1,1}}=\sqrt{2m+1}=\sqrt{3}\\ + \frac{N_{m-1,m-1}}{N_{m,m}}=\frac{1}{2m-1}\sqrt{\frac{2m+1}{2m}}\quad(m\geq2)\\ + \frac{N_{n-1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}\quad(n\geq1,0\leq m \leq n)\\ + \frac{N_{n-2,m}}{N_{n,m}}=\frac{N_{n-2,m}}{N_{n-1,m}}\frac{N_{n-1,m}}{N_{n,m}}\\ + \frac{N_{n-2,m}}{N_{n,m}}=\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\frac{N_{n-1,m}}{N_{n,m}}\quad(n\geq2,0\leq m \leq n) +\end{align} ``` The recurrence relations for normalized V and W are derived as follows ```math -\bar{V}_{0,0}=\frac{Re}{r}\\ -\bar{W}_{0,0}=0\\ -\bar{V}_{1,1}=\sqrt{3}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{0,0}-\frac{yR_e}{r^2}\bar{W}_{0,0})\\ +\begin{align} + \bar{V}_{0,0}=\frac{Re}{r}\\ + \bar{W}_{0,0}=0\\ + \bar{V}_{1,1}=\sqrt{3}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{0,0}-\frac{yR_e}{r^2}\bar{W}_{0,0})\\ +\end{align} ``` ```math -\bar{V}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\quad(m\geq2)\\ -\bar{W}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\quad(m\geq2)\\ +\begin{align} + \bar{V}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\quad(m\geq2)\\ + \bar{W}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\quad(m\geq2)\\ +\end{align} ``` ```math -\bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ -\bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ +\begin{align} + \bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ + \bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ +\end{align} ``` ```math -\bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{V}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ -\bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{W}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ +\begin{align} + \bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{V}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ + \bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{W}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ +\end{align} ``` @@ -160,26 +178,36 @@ For unnormalized algorithms, See chapter 3.2.5 of [Satellite Orbits](https://ww When we use the normalized coefficients $\bar{C}_{n,m}$ and $\bar{S}_{n,m}$ and $\bar{V}_{n,m}$ and $\bar{W}_{n,m}$ functions, the acceleration calculation is described like follows ```math -\ddot{x}_{n,m}=-\frac{GM}{Re^{2}}\bar{C}_{n,0}\bar{V}_{n+1,1} =-\frac{GM}{Re^{2}} C_{n,0}V_{n+1,1} \frac{N_{n,0}}{N_{n+1,1}}\quad(m=0) +\begin{align} + \ddot{x}_{n,m}=-\frac{GM}{Re^{2}}\bar{C}_{n,0}\bar{V}_{n+1,1} =-\frac{GM}{Re^{2}} C_{n,0}V_{n+1,1} \frac{N_{n,0}}{N_{n+1,1}}\quad(m=0) +\end{align} ``` The division of normalized function $\frac{N_{n,0}}{N_{n+1,1}}$ should be removed, so we have to multiply following correction factors into the equation. When $m=0$, following correction factors are useful for x and y acceleration ```math -\frac{N_{n+1,1}}{N_{n,0}}=\sqrt{\frac{1}{2}}\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+2)(n+1)} +\begin{align} + \frac{N_{n+1,1}}{N_{n,0}}=\sqrt{\frac{1}{2}}\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+2)(n+1)} +\end{align} ``` When $m=1$, following correction factors are useful for x and y acceleration ```math -\frac{N_{n+1,0}}{N_{n,1}}=\sqrt{2}\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{n(n+1)}}\quad(m=1)\\ +\begin{align} + \frac{N_{n+1,0}}{N_{n,1}}=\sqrt{2}\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{n(n+1)}}\quad(m=1)\\ +\end{align} ``` When $m>1$, following correction factors are useful for x and y acceleration ```math -\frac{N_{n+1,m+1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+m+1)(n+m+2)}\\ -\frac{N_{n+1,m-1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{(n-m+1)(n-m+2)}}\\ +\begin{align} + \frac{N_{n+1,m+1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+m+1)(n+m+2)}\\ + \frac{N_{n+1,m-1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{(n-m+1)(n-m+2)}}\\ +\end{align} ``` When $m>=0$, following correction factors are useful for z acceleration ```math -\frac{N_{n+1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{n+m+1}{n-m+1}} +\begin{align} + \frac{N_{n+1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{n+m+1}{n-m+1}} +\end{align} ``` 4. note - To accelerate the calculation, the double `for loop` of acceleration calculation and the recursion loop need to be integrated in future. From c05881a301127a347720dcfbc65a0a30d83473fe Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:11:59 +0100 Subject: [PATCH 06/19] fix equation in component --- Specifications/Component/AOCS/Spec_STT.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Specifications/Component/AOCS/Spec_STT.md b/Specifications/Component/AOCS/Spec_STT.md index 06c1d647..a67ba320 100644 --- a/Specifications/Component/AOCS/Spec_STT.md +++ b/Specifications/Component/AOCS/Spec_STT.md @@ -27,12 +27,12 @@ ```math - \theta_{ce} = \arctan{\left(\frac{|\boldsymbol{r_{se}}|}{R_e}\right)} - \\ - \theta_{cs} = \arccos{(\boldsymbol{r_{se}}*\boldsymbol{r_{sight}})} - \\ + \begin{align} + \theta_{ce} = \arctan{\left(\frac{|\boldsymbol{r_{se}}|}{R_e}\right)}\\ + \theta_{cs} = \arccos{(\boldsymbol{r_{se}}*\boldsymbol{r_{sight}})}\\ \theta_{es} = \theta_{ce} - \theta_{cs} \tag{1} + \end{align} ``` ![](./figs/stt_earth_judgement.png) From 150d900f120789220f012492ed45fbdca913b1b1 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:20:23 +0100 Subject: [PATCH 07/19] fix equation in disturbances --- .../Disturbance/Spec_GeoPotential.md | 64 +++++++++---------- .../Disturbance/Spec_SurfaceForce.md | 6 +- .../Disturbance/Spec_SurfaceForce_AirDrag.md | 32 ++++++---- .../Spec_SurfaceForce_SolarRadiation.md | 8 ++- .../Disturbance/Spec_ThirdBodyGravity.md | 2 +- 5 files changed, 60 insertions(+), 52 deletions(-) diff --git a/Specifications/Disturbance/Spec_GeoPotential.md b/Specifications/Disturbance/Spec_GeoPotential.md index 13842ae4..0f2e8575 100644 --- a/Specifications/Disturbance/Spec_GeoPotential.md +++ b/Specifications/Disturbance/Spec_GeoPotential.md @@ -85,9 +85,9 @@ For normalized algorithm, we use following normalizing relation for Legendre pol ```math \begin{align} - \bar{P}_{n,m}=\frac{1}{N_{n,m}}P_{n,m}\\ - \bar{V}_{n,m}=\frac{1}{N_{n,m}}V_{n,m}\\ - \bar{W}_{n,m}=\frac{1}{N_{n,m}}W_{n,m}\\ + \bar{P}_{n,m} &= \frac{1}{N_{n,m}}P_{n,m}\\ + \bar{V}_{n,m} &= \frac{1}{N_{n,m}}V_{n,m}\\ + \bar{W}_{n,m} &= \frac{1}{N_{n,m}}W_{n,m}\\ \end{align} ``` @@ -95,17 +95,17 @@ The recursion calculation of V and W can be changed to a normalized version as f ```math \begin{align} - N_{m,m}\bar{V}_{m,m}=(2m-1)(\frac{xR_{e}}{r^2}N_{m-1,m-1}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}N_{m-1,m-1}\bar{W}_{m-1,m-1})\\ - \bar{V}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\\ - \bar{W}_{m,m}=\frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\\ + N_{m,m}\bar{V}_{m,m} &= (2m-1)(\frac{xR_{e}}{r^2}N_{m-1,m-1}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}N_{m-1,m-1}\bar{W}_{m-1,m-1})\\ + \bar{V}_{m,m} &= \frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\\ + \bar{W}_{m,m} &= \frac{N_{m-1,m-1}}{N_{m,m}}(2m-1)(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\\ \end{align} ``` ```math \begin{align} - N_{n,m}\bar{V}_{n,m}=\frac{2n-1}{n-m}(\frac{zR_{e}}{r^2}N_{n-1,m}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}N_{n-2,m}\bar{V}_{n-2,m})\\ - \bar{V}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{V}_{n-2,m}\\ - \bar{W}_{n,m}=\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{W}_{n-2,m}\\ + N_{n,m}\bar{V}_{n,m} &= \frac{2n-1}{n-m}(\frac{zR_{e}}{r^2}N_{n-1,m}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}N_{n-2,m}\bar{V}_{n-2,m})\\ + \bar{V}_{n,m} &= \frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{V}_{n-2,m}\\ + \bar{W}_{n,m} &= \frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\frac{N_{n-1,m}}{N_{n,m}}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\frac{N_{n-2,m}}{N_{n,m}}\bar{W}_{n-2,m}\\ \end{align} ``` @@ -113,51 +113,46 @@ The recurrence relation of normalize function can be expressed as follows ```math \begin{align} - N_{0,0}=1\\ - N_{m,m}=(2m-1)\sqrt{\frac{2m}{2m+1}}N_{m-1,m-1}\quad(m\geq1) -\end{align} -``` - -```math -\begin{align} - N_{n,m}=\sqrt{\frac{2n-1}{2n+1}}\sqrt{\frac{n+m}{n-m}}N_{n-1,m}\quad(n\geq1,0\leq m \leq n) + N_{0,0} &= 1\\ + N_{m,m} &= (2m-1)\sqrt{\frac{2m}{2m+1}}N_{m-1,m-1}\quad(m\geq1)\\ + N_{n,m} &= \sqrt{\frac{2n-1}{2n+1}}\sqrt{\frac{n+m}{n-m}}N_{n-1,m}\quad(n\geq1,0\leq m \leq n) \end{align} ``` So, the divisions of the normalized functions are described as follows ```math \begin{align} - \frac{N_{0,0}}{N_{1,1}}=\sqrt{2m+1}=\sqrt{3}\\ - \frac{N_{m-1,m-1}}{N_{m,m}}=\frac{1}{2m-1}\sqrt{\frac{2m+1}{2m}}\quad(m\geq2)\\ - \frac{N_{n-1,m}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}\quad(n\geq1,0\leq m \leq n)\\ - \frac{N_{n-2,m}}{N_{n,m}}=\frac{N_{n-2,m}}{N_{n-1,m}}\frac{N_{n-1,m}}{N_{n,m}}\\ - \frac{N_{n-2,m}}{N_{n,m}}=\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\frac{N_{n-1,m}}{N_{n,m}}\quad(n\geq2,0\leq m \leq n) + \frac{N_{0,0}}{N_{1,1}} &= \sqrt{2m+1}=\sqrt{3}\\ + \frac{N_{m-1,m-1}}{N_{m,m}} &= \frac{1}{2m-1}\sqrt{\frac{2m+1}{2m}}\quad(m\geq2)\\ + \frac{N_{n-1,m}}{N_{n,m}} &= \sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}\quad(n\geq1,0\leq m \leq n)\\ + \frac{N_{n-2,m}}{N_{n,m}} &= \frac{N_{n-2,m}}{N_{n-1,m}}\frac{N_{n-1,m}}{N_{n,m}}\\ + \frac{N_{n-2,m}}{N_{n,m}} &= \sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\frac{N_{n-1,m}}{N_{n,m}}\quad(n\geq2,0\leq m \leq n) \end{align} ``` The recurrence relations for normalized V and W are derived as follows ```math \begin{align} - \bar{V}_{0,0}=\frac{Re}{r}\\ - \bar{W}_{0,0}=0\\ - \bar{V}_{1,1}=\sqrt{3}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{0,0}-\frac{yR_e}{r^2}\bar{W}_{0,0})\\ + \bar{V}_{0,0} &= \frac{Re}{r}\\ + \bar{W}_{0,0} &= 0\\ + \bar{V}_{1,1} &= \sqrt{3}(2m-1)(\frac{xR_{e}}{r^2}\bar{V}_{0,0}-\frac{yR_e}{r^2}\bar{W}_{0,0})\\ \end{align} ``` ```math \begin{align} - \bar{V}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\quad(m\geq2)\\ - \bar{W}_{m,m}=\sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\quad(m\geq2)\\ + \bar{V}_{m,m} &= \sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{V}_{m-1,m-1}-\frac{yR_e}{r^2}\bar{W}_{m-1,m-1})\quad(m\geq2)\\ + \bar{W}_{m,m} &= \sqrt{\frac{2m+1}{2m}}(\frac{xR_{e}}{r^2}\bar{W}_{m-1,m-1}+\frac{yR_e}{r^2}\bar{V}_{m-1,m-1})\quad(m\geq2)\\ \end{align} ``` ```math \begin{align} - \bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ - \bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ + \bar{V}_{n,m} &= \sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ + \bar{W}_{n,m} &= \sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m})\quad(n=1,0\leq m \leq n)\\ \end{align} ``` ```math \begin{align} - \bar{V}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{V}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ - \bar{W}_{n,m}=\sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{W}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ + \bar{V}_{n,m} &= \sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{V}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{V}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ + \bar{W}_{n,m} &= \sqrt{\frac{2n+1}{2n-1}}\sqrt{\frac{n-m}{n+m}}(\frac{2n-1}{n-m}\frac{zR_{e}}{r^2}\bar{W}_{n-1,m}-\frac{n+m-1}{n-m}\frac{R_e^2}{r^2}\sqrt{\frac{2n-1}{2n-3}}\sqrt{\frac{n-m-1}{n+m-1}}\bar{W}_{n-2,m})\quad(n\geq2,0\leq m \leq n)\\ \end{align} ``` @@ -179,7 +174,8 @@ For unnormalized algorithms, See chapter 3.2.5 of [Satellite Orbits](https://ww When we use the normalized coefficients $\bar{C}_{n,m}$ and $\bar{S}_{n,m}$ and $\bar{V}_{n,m}$ and $\bar{W}_{n,m}$ functions, the acceleration calculation is described like follows ```math \begin{align} - \ddot{x}_{n,m}=-\frac{GM}{Re^{2}}\bar{C}_{n,0}\bar{V}_{n+1,1} =-\frac{GM}{Re^{2}} C_{n,0}V_{n+1,1} \frac{N_{n,0}}{N_{n+1,1}}\quad(m=0) + \ddot{x}_{n,m} &= -\frac{GM}{Re^{2}}\bar{C}_{n,0}\bar{V}_{n+1,1}\\ + &= -\frac{GM}{Re^{2}} C_{n,0}V_{n+1,1} \frac{N_{n,0}}{N_{n+1,1}}\quad(m=0) \end{align} ``` The division of normalized function $\frac{N_{n,0}}{N_{n+1,1}}$ should be removed, so we have to multiply following correction factors into the equation. @@ -199,8 +195,8 @@ When $m=1$, following correction factors are useful for x and y acceleration When $m>1$, following correction factors are useful for x and y acceleration ```math \begin{align} - \frac{N_{n+1,m+1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+m+1)(n+m+2)}\\ - \frac{N_{n+1,m-1}}{N_{n,m}}=\sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{(n-m+1)(n-m+2)}}\\ + \frac{N_{n+1,m+1}}{N_{n,m}} &= \sqrt{\frac{2n+1}{2n+3}}\sqrt{(n+m+1)(n+m+2)}\\ + \frac{N_{n+1,m-1}}{N_{n,m}} &= \sqrt{\frac{2n+1}{2n+3}}\sqrt{\frac{1}{(n-m+1)(n-m+2)}}\\ \end{align} ``` When $m>=0$, following correction factors are useful for z acceleration diff --git a/Specifications/Disturbance/Spec_SurfaceForce.md b/Specifications/Disturbance/Spec_SurfaceForce.md index 0185372b..1306c1bb 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce.md +++ b/Specifications/Disturbance/Spec_SurfaceForce.md @@ -63,8 +63,10 @@ - $\boldsymbol{r}_{cg}$ is the position vector of the center of mass ```math - \boldsymbol{F}=-C_{n}\boldsymbol{n}+C_{t}\boldsymbol{t}\\ - \boldsymbol{T}=(\boldsymbol{r}_{s}-\boldsymbol{r}_{cg})\times\boldsymbol{F} + \begin{align} + \boldsymbol{F} &= -C_{n}\boldsymbol{n}+C_{t}\boldsymbol{t}\\ + \boldsymbol{T} &= (\boldsymbol{r}_{s}-\boldsymbol{r}_{cg})\times\boldsymbol{F} + \end{align} ``` - Detail of the $C_{n}$ and $C_{t}$ are defined by subclasses by using `CalcCoef` function diff --git a/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md b/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md index 8a894d5d..1198f32a 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md +++ b/Specifications/Disturbance/Spec_SurfaceForce_AirDrag.md @@ -1,4 +1,4 @@ -# Surface Force ~Air Drag~ +# Surface Force: Air Drag ## 1. Overview @@ -25,9 +25,11 @@ - `CalcCoef` calculates the normal and in-plane coefficients for `SurfaceForce` calculation. The air drag force acting on a surface is expressed as the following equation ```math - \boldsymbol{F}=-C_{n}\boldsymbol{n}+C_{t}\boldsymbol{t}\\ - C_{n}=\frac{1}{2}\rho A v^2 C_{n}^{\prime}\\ - C_{t}=\frac{1}{2}\rho A v^2 C_{t}^{\prime} + \begin{align} + \boldsymbol{F}=-C_{n}\boldsymbol{n}+C_{t}\boldsymbol{t}\\ + C_{n}=\frac{1}{2}\rho A v^2 C_{n}^{\prime}\\ + C_{t}=\frac{1}{2}\rho A v^2 C_{t}^{\prime} + \end{align} ``` - This function mainly calculates the common part of the coefficient calculation. $C_{n}^{\prime}$ and $C_{t}^{\prime}$ are calculated in `CalCnCt` function, and they will be used in this function. @@ -68,9 +70,11 @@ - $C_{n}^{\prime}$ and $C_{t}^{\prime}$ are calculated as following equations ```math - C_{n}^{\prime} = \frac{2-\sigma_{d}}{\sqrt{\pi}}\frac{\Pi(S_{n})}{S^{2}}+\frac{\sigma_{d}}{2}\frac{\chi(S_{n})}{S^{2}}\sqrt{\frac{T_{w}}{T_{m}}}\\ - C_{t}^{\prime} =\frac{\sigma_{d}}{\sqrt{\pi}}\frac{\chi(S_{n})}{S^{2}}S_{t} - ``` + \begin{align} + C_{n}^{\prime} &= \frac{2-\sigma_{d}}{\sqrt{\pi}}\frac{\Pi(S_{n})}{S^{2}}+\frac{\sigma_{d}}{2}\frac{\chi(S_{n})}{S^{2}}\sqrt{\frac{T_{w}}{T_{m}}}\\ + C_{t}^{\prime} &=\frac{\sigma_{d}}{\sqrt{\pi}}\frac{\chi(S_{n})}{S^{2}}S_{t} + \end{align} + ``` - $S, S_{n}, S_{t}$ are defined as follows - $k=1.38064852E-23$ is the Boltzmann constant @@ -78,17 +82,21 @@ - $\cos{\theta}$ and $\sin{\theta}$ are calculated in `SurfaceForce` base class. ```math - S=\sqrt{\frac{Mv^{2}}{2kT_{w}}}\\ - S_{n}=S\cos{\theta}\\ - S_{t}=S\sin{\theta}\\ + \begin{align} + S &= \sqrt{\frac{Mv^{2}}{2kT_{w}}}\\ + S_{n} &= S\cos{\theta}\\ + S_{t} &= S\sin{\theta}\\ + \end{align} ``` - $\Pi(x)$ and $\chi(x)$ are defined as follows - where `erf` is the [Gauss error function](https://en.wikipedia.org/wiki/Error_function). ```math - \Pi(x)=x e^{-x^{2}}+\sqrt{\pi}(x^2+0.5)(1+erf(x))\\ - \chi(x)=e^{-x^{2}}+\sqrt{\pi}x(1+erf(x)) + \begin{align} + \Pi(x) &= x e^{-x^{2}}+\sqrt{\pi}(x^2+0.5)(1+erf(x))\\ + \chi(x) &= e^{-x^{2}}+\sqrt{\pi}x(1+erf(x)) + \end{align} ``` 4. note - Please see the reference document for more information on detailed calculations. diff --git a/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md b/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md index c81ae525..aae349ea 100644 --- a/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md +++ b/Specifications/Disturbance/Spec_SurfaceForce_SolarRadiation.md @@ -1,4 +1,4 @@ -# Surface Force ~Solar Radiation Pressure~ +# Surface Force: Solar Radiation Pressure ## 1. Overview @@ -41,8 +41,10 @@ - $\theta$ is the angle between the normal vector and the sun vector. ```math - C_{n} = AP \left((1+\nu\mu)\cos^{2}{\theta}+\frac{2}{3}\nu(1-\mu)\cos{\theta} \right)\\ - C_{t} = AP(1-\nu\mu)\cos{\theta}\sin{\theta} + \begin{align} + C_{n} &= AP \left((1+\nu\mu)\cos^{2}{\theta}+\frac{2}{3}\nu(1-\mu)\cos{\theta} \right)\\ + C_{t} &= AP(1-\nu\mu)\cos{\theta}\sin{\theta} + \end{align} ``` 4. note diff --git a/Specifications/Disturbance/Spec_ThirdBodyGravity.md b/Specifications/Disturbance/Spec_ThirdBodyGravity.md index de1b403e..c4103af1 100644 --- a/Specifications/Disturbance/Spec_ThirdBodyGravity.md +++ b/Specifications/Disturbance/Spec_ThirdBodyGravity.md @@ -38,7 +38,7 @@ 3. algorithm - Definition of the variables - - \boldsymbol{r}$: the position of the spacecraft (written as `sat_pos_i` in the code) + - $\boldsymbol{r}$: the position of the spacecraft (written as `sat_pos_i` in the code) - $\boldsymbol{s}$: the position of the third body (written as `third_body_pos_i` in the code) - $\boldsymbol{s}-\boldsymbol{r}$: the position of the third body seen from the spacecraft (written as `third_body_pos_from_sc_i` in the code) From 6301d21b0611ad30055efaf74b83a843e4005e0e Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:37:37 +0100 Subject: [PATCH 08/19] fix equation in dynamics --- README.md | 6 +- .../Dynamics/Spec_AttitudeDynamics.md | 38 +++++------ Specifications/Dynamics/Spec_EnckeMethod.md | 8 ++- Specifications/Dynamics/Spec_KeplerOrbit.md | 64 +++++++++++-------- Specifications/Dynamics/Spec_RelativeOrbit.md | 4 +- Specifications/Dynamics/Spec_Rk4Orbit.md | 16 +++-- 6 files changed, 77 insertions(+), 59 deletions(-) diff --git a/README.md b/README.md index a1c79659..e4a912b6 100644 --- a/README.md +++ b/README.md @@ -79,7 +79,11 @@ 1. [Attitude Dynamics](./Specifications/Dynamics/Spec_AttitudeDynamics.md) 2. [Controlled Attitude](./Specifications/Dynamics/Spec_ControlledAttitude.md) 2. [Orbit](./Specifications/Dynamics/Spec_Orbit.md) - 1. [Relative Orbit](./Specifications/Dynamics/Spec_RelativeOrbit.md) + 1. [Kepler Orbit](./Specifications/Dynamics/Spec_KeplerOrbit.md) + 2. [RK4 Orbit Propagation](./Specifications/Dynamics/Spec_Rk4Orbit.md) + 3. [SGP4 Orbit Propagation with TLE](./Specifications/Dynamics/Spec_Sgp4Orbit.md) + 4. [ENCKE method](./Specifications/Dynamics/Spec_EnckeOrbit.md) + 5. [Relative Orbit](./Specifications/Dynamics/Spec_RelativeOrbit.md) 3. Thermal 4. Environment diff --git a/Specifications/Dynamics/Spec_AttitudeDynamics.md b/Specifications/Dynamics/Spec_AttitudeDynamics.md index a2ef38b4..eb3b20e9 100644 --- a/Specifications/Dynamics/Spec_AttitudeDynamics.md +++ b/Specifications/Dynamics/Spec_AttitudeDynamics.md @@ -61,38 +61,36 @@ If the differential equation (1) is given, the state quantity in $n+1$ step can be calculated as (2). ```math - \hat{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x},t) - \tag{1} - ``` - ```math - \boldsymbol{x_{n+1}} = \boldsymbol{x_{n}} + \cfrac{\Delta t}{6}(\boldsymbol{k_1}+2\boldsymbol{k_2} - +2\boldsymbol{k_3}+\boldsymbol{k_4}) - \tag{2} + \begin{align} + \hat{\boldsymbol{x}} &= \boldsymbol{f}(\boldsymbol{x},t)\\ + \boldsymbol{x_{n+1}} &= \boldsymbol{x_{n}} + \cfrac{\Delta t}{6}(\boldsymbol{k_1}+2\boldsymbol{k_2}+2\boldsymbol{k_3}+\boldsymbol{k_4}) + \end{align} ``` where $\Delta t$ is a time step, which meets the equation (3). ```math t_{n+1} = t_{n} + \Delta t - \tag{3} ``` $k_i \, (i=1,2,3,4)$, which has the same number of elements, can be calculated as the equations (4). ```math - \boldsymbol{k_{1}} = \boldsymbol{f}(\boldsymbol{x_n},t_n) \\ - \boldsymbol{k_{2}} = \boldsymbol{f}\left(\boldsymbol{x_n}+\frac{\Delta t}{2} \boldsymbol{k_1},t_n+\frac{\Delta t}{2} \right) \\ - \boldsymbol{k_{3}} = \boldsymbol{f}\left(\boldsymbol{x_n}+\frac{\Delta t}{2} \boldsymbol{k_2},t_n+\frac{\Delta t}{2} \right) \\ - \boldsymbol{k_{4}} = \boldsymbol{f}\left(\boldsymbol{x_n} + \Delta t \boldsymbol{k_3},t_n+\Delta t \right) - \tag{4} + \begin{align} + \boldsymbol{k_{1}} &= \boldsymbol{f}(\boldsymbol{x_n},t_n) \\ + \boldsymbol{k_{2}} &= \boldsymbol{f}\left(\boldsymbol{x_n}+\frac{\Delta t}{2} \boldsymbol{k_1},t_n+\frac{\Delta t}{2} \right) \\ + \boldsymbol{k_{3}} &= \boldsymbol{f}\left(\boldsymbol{x_n}+\frac{\Delta t}{2} \boldsymbol{k_2},t_n+\frac{\Delta t}{2} \right) \\ + \boldsymbol{k_{4}} &= \boldsymbol{f}\left(\boldsymbol{x_n} + \Delta t \boldsymbol{k_3},t_n+\Delta t \right) + \end{align} ``` In this attitude propagation, the quantity of state $\boldsymbol{x}$ consists of 7 elements, including `Quaternion_i2b` and angular velocity $\boldsymbol{\omega}_b$. ```math - \boldsymbol{\omega}_b = [{\omega}_{bx} \, {\omega}_{by} \, {\omega}_{bz}]^T \\ - \boldsymbol{q}_{i2b} = [q_x \, q_y \, q_z \, q_w]^T \\ - \boldsymbol{x} = [\boldsymbol{\omega}_b, \boldsymbol{q}_{i2b}]^T - \tag{5} + \begin{align} + \boldsymbol{\omega}_b &= [{\omega}_{bx} \, {\omega}_{by} \, {\omega}_{bz}]^T \\ + \boldsymbol{q}_{i2b} &= [q_x \, q_y \, q_z \, q_w]^T \\ + \boldsymbol{x} &= [\boldsymbol{\omega}_b, \boldsymbol{q}_{i2b}]^T + \end{align} ``` 4. note @@ -115,7 +113,6 @@ ```math \dot{\boldsymbol{\omega}}_b = \boldsymbol{I}_b^{-1}(\boldsymbol{T}_b - \boldsymbol{\omega}_b \times \boldsymbol{h}_b) - \tag{6} ``` where $\boldsymbol{\omega}_b$[rad/s] is angular velocity in the body-fixed coordinate, $\boldsymbol{I}_b$[kgm$^2$] is inertia tensor of the satellite, $\boldsymbol{T}_b$[Nm] is torque in the body-fixed coordinate, $\boldsymbol{h}_b$[Nms] is angular momentum of the satellite in the body-fixed coordinate. @@ -125,12 +122,11 @@ \dot{\boldsymbol{q}}_{i2b} = \cfrac{1}{2} \begin{bmatrix} 0 & {\omega}_{bz} & -{\omega}_{by} & {\omega}_{bx} \\ - \- {\omega}_{bz} & 0 & {\omega}_{bx} & {\omega}_{by} \\ + - {\omega}_{bz} & 0 & {\omega}_{bx} & {\omega}_{by} \\ {\omega}_{by} & -{\omega}_{bx} & 0 & {\omega}_{bz} \\ - \- {\omega}_{bx} & -{\omega}_{by} & -{\omega}_{bz} & 0 + - {\omega}_{bx} & -{\omega}_{by} & -{\omega}_{bz} & 0 \end{bmatrix} \boldsymbol{q}_{i2b} - \tag{7} ``` ## 3. Results of verifications diff --git a/Specifications/Dynamics/Spec_EnckeMethod.md b/Specifications/Dynamics/Spec_EnckeMethod.md index a6d988a2..f97eaa2b 100644 --- a/Specifications/Dynamics/Spec_EnckeMethod.md +++ b/Specifications/Dynamics/Spec_EnckeMethod.md @@ -60,9 +60,11 @@ 3. Propagate the difference - Propagate the following differential equation. At this moment, we use the fourth-order Runge-Kutta method as a propagator. ```math - \ddot{\boldsymbol{\delta}} = -\frac{\mu}{r_{ref}^3}(\boldsymbol{\delta}+f(q)\boldsymbol{r})+\boldsymbol{a}_d\\ - f(q) = q \frac{q^2 + 3q + 3}{(1+q)^{1.5} + 1}\\ - q = \frac{\boldsymbol{\delta}\cdot(\boldsymbol{\delta}-2\boldsymbol{r}_i)}{r_i} + \begin{align} + \ddot{\boldsymbol{\delta}} &= -\frac{\mu}{r_{ref}^3}(\boldsymbol{\delta}+f(q)\boldsymbol{r})+\boldsymbol{a}_d\\ + f(q) &= q \frac{q^2 + 3q + 3}{(1+q)^{1.5} + 1}\\ + q &= \frac{\boldsymbol{\delta}\cdot(\boldsymbol{\delta}-2\boldsymbol{r}_i)}{r_i} + \end{align} ``` ## 3. Results of verifications diff --git a/Specifications/Dynamics/Spec_KeplerOrbit.md b/Specifications/Dynamics/Spec_KeplerOrbit.md index 003689d5..f339c038 100644 --- a/Specifications/Dynamics/Spec_KeplerOrbit.md +++ b/Specifications/Dynamics/Spec_KeplerOrbit.md @@ -66,7 +66,7 @@ ```math R_z(\theta) = \begin{pmatrix} \cos{\theta} & \sin{\theta} & 0 \\ - \-\sin{\theta} & \cos{\theta} & 0 \\ + -\sin{\theta} & \cos{\theta} & 0 \\ 0 & 0 & 1 \end{pmatrix} ``` @@ -93,21 +93,25 @@ - Details are described in `KeplerOrbit::SolveKeplerFirstOrder` - Calculate two dimensional position $x^{\*}$, $y^{\*}$ and velocity $\dot{x}^{\*}$, $\dot{y}^{\*}$ in the orbital plane ```math - x^* = a(\cos{u}-e)\\ - y^* = a\sqrt{1-e^2}\sin{u}\\ - \dot{x}^* = -an\frac{\sin{u}}{1^e\cos{u}}\\ - \dot{y}^* = an\sqrt{1-e^2}\frac{\cos{u}}{1^e\cos{u}}\\ + \begin{align} + x^* &= a(\cos{u}-e)\\ + y^* &= a\sqrt{1-e^2}\sin{u}\\ + \dot{x}^* &= -an\frac{\sin{u}}{1^e\cos{u}}\\ + \dot{y}^* &= an\sqrt{1-e^2}\frac{\cos{u}}{1^e\cos{u}}\\ + \end{align} ``` - Convert to the inertial frame ```math - \boldsymbol{r}_{i} = R_{p2i}\begin{bmatrix} x^* \\ - y^* \\ - 0 - \end{bmatrix}\\ - \boldsymbol{v}_{i} = R_{p2i}\begin{bmatrix} \dot{x}^* \\ - \dot{y}^* \\ - 0 - \end{bmatrix}\\ + \begin{align} + \boldsymbol{r}_{i} &= R_{p2i}\begin{bmatrix} x^* \\ + y^* \\ + 0 + \end{bmatrix}\\ + \boldsymbol{v}_{i} &= R_{p2i}\begin{bmatrix} \dot{x}^* \\ + \dot{y}^* \\ + 0 + \end{bmatrix}\\ + \end{align} ``` 3. `KeplerOrbit::SolveKeplerFirstOrder` function @@ -168,19 +172,25 @@ ``` - $x_{p}, y_{p}$ : Position in the orbital plane ```math - x_{p} = r_{ix} \cos{\Omega} + r_{iy} \sin{\Omega};\\ - y_{p} = (-r_{ix} \sin{\Omega} + r_{iy} \cos{\Omega})\cos{i} + r_{iz}\sin{i}; + \begin{align} + x_{p} &= r_{ix} \cos{\Omega} + r_{iy} \sin{\Omega};\\ + y_{p} &= (-r_{ix} \sin{\Omega} + r_{iy} \cos{\Omega})\cos{i} + r_{iz}\sin{i} + \end{align} ``` - $\dot{x}_{p}, \dot{y}_{p}$ : Velocity in the orbital plane ```math - \dot{x}_{p} = v_{ix} \cos{\Omega} + v_{iy} \sin{\Omega};\\ - \dot{y}_{p} = (-v_{ix} \sin{\Omega} + v_{iy} \cos{\Omega})\cos{i} + v_{iz}\sin{i}; + \begin{align} + \dot{x}_{p} &= v_{ix} \cos{\Omega} + v_{iy} \sin{\Omega};\\ + \dot{y}_{p} &= (-v_{ix} \sin{\Omega} + v_{iy} \cos{\Omega})\cos{i} + v_{iz}\sin{i} + \end{align} ``` - $e$ : Eccentricity ```math - c_1 = \frac{h}{\mu}\dot{y}_p - \frac{x_p}{r}\\ - c_2 = -\frac{h}{\mu}\dot{x}_p - \frac{y_p}{r}\\ - e = \sqrt{c_1^2 + c_2^2} + \begin{align} + c_1 &= \frac{h}{\mu}\dot{y}_p - \frac{x_p}{r}\\ + c_2 &= -\frac{h}{\mu}\dot{x}_p - \frac{y_p}{r}\\ + e &= \sqrt{c_1^2 + c_2^2} + \end{align} ``` - $\omega$ : Argument of Perigee ```math @@ -188,13 +198,17 @@ ``` - $t_{epoch}$ : Epoch [Julian day] ```math - f = \tan^{-1}\left(\frac{y_p}{x_p}\right) - \omega\\ - u = \tan^{-1}\frac{\frac{r \sin{f}}{\sqrt{1-e^2}}}{r\cos{f} + ae}\\ + \begin{align} + f &= \tan^{-1}\left(\frac{y_p}{x_p}\right) - \omega\\ + u &= \tan^{-1}\frac{\frac{r \sin{f}}{\sqrt{1-e^2}}}{r\cos{f} + ae}\\ + \end{align} ``` ```math - n = \sqrt{\frac{\mu}{a^3}}\\ - dt = \frac{u - e\sin{u}}{n}\\ - t_{epoch} = t - \frac{dt}{24*60*60} + \begin{align} + n &= \sqrt{\frac{\mu}{a^3}}\\ + dt &= \frac{u - e\sin{u}}{n}\\ + t_{epoch} &= t - \frac{dt}{24*60*60} + \end{align} ``` ## 3. Results of verifications diff --git a/Specifications/Dynamics/Spec_RelativeOrbit.md b/Specifications/Dynamics/Spec_RelativeOrbit.md index 38eb5e84..b531c81d 100644 --- a/Specifications/Dynamics/Spec_RelativeOrbit.md +++ b/Specifications/Dynamics/Spec_RelativeOrbit.md @@ -30,8 +30,8 @@ - Set up the configuration of the `[ORBIT]` section in the `Sat.ini`. + Set `propagate_mode =2` to use the relative orbit propagation + Choose `update_method`. - * `update_method = 0` means the orbit is updated with the propagation of the relative dynamics equation($\dot{\boldsymbol{x}}=\boldsymbol{Ax}+\boldsymbol{Bu}$, i.e., Hill equation). - * `update_method = 1` means the orbit is updated with the STM($\boldsymbol{x}(t)=\boldsymbol{\Phi}(t,t_0)\boldsymbol{x}(t_0)$, i.e., Clohessy-Wiltshire solution). + * `update_method = 0` means the orbit is updated with the propagation of the relative dynamics equation( $\dot{\boldsymbol{x}}=\boldsymbol{Ax}+\boldsymbol{Bu}$ , i.e., Hill equation). + * `update_method = 1` means the orbit is updated with the STM( $\boldsymbol{x}(t)=\boldsymbol{\Phi}(t,t_0)\boldsymbol{x}(t_0)$ , i.e., Clohessy-Wiltshire solution). + When you choose `update_method = 0`, set `relative_dynamics_model_type`. + When you choose `update_method = 1`, set `stm_model_type`. + Set the initial relative position of the satellite in the LVLH frame. LVLH frame definition is: diff --git a/Specifications/Dynamics/Spec_Rk4Orbit.md b/Specifications/Dynamics/Spec_Rk4Orbit.md index 4bfdcc73..adcacd28 100644 --- a/Specifications/Dynamics/Spec_Rk4Orbit.md +++ b/Specifications/Dynamics/Spec_Rk4Orbit.md @@ -31,13 +31,15 @@ - The position and velocity of the satellite are updated by using RK4. As the input of RK4, the six-state variables are set. These state variables are the three-dimensional position [$x$, $y$ ,$z$] and three-dimensional velocity [$v_x$, $v_y$, $v_z$] at the inertial coordinate. Here, the inertial coordinate is decided by the `PlanetSelect.ini` - As the force which works to the satellite motion is the external acceleration [$a_x$,$a_y$,$a_z$] calculated from the disturbance class or thruster class and the gravity force from the center planet, which is defined in `PlanetSelect.ini`. As a summary, the orbit is calculated as the following equation. ```math - \dot{x} = v_x\\ - \dot{y} = v_y\\ - \dot{z} = v_z\\ - \dot{v}_x = a_x-\mu\frac{x}{r^3}\\ - \dot{v}_y = a_y-\mu\frac{y}{r^3}\\ - \dot{v}_z = a_z-\mu\frac{z}{r^3}\\ - r = \sqrt{x^2+y^2+z^2} + \begin{align} + \dot{x} &= v_x\\ + \dot{y} &= v_y\\ + \dot{z} &= v_z\\ + \dot{v}_x &= a_x-\mu\frac{x}{r^3}\\ + \dot{v}_y &= a_y-\mu\frac{y}{r^3}\\ + \dot{v}_z &= a_z-\mu\frac{z}{r^3}\\ + r &= \sqrt{x^2+y^2+z^2} + \end{align} ``` ## 3. Results of verifications From b7179cdd91dad0a480e16742dfc6519268dae53b Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:38:37 +0100 Subject: [PATCH 09/19] fix link --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index e4a912b6..238198cd 100644 --- a/README.md +++ b/README.md @@ -81,8 +81,8 @@ 2. [Orbit](./Specifications/Dynamics/Spec_Orbit.md) 1. [Kepler Orbit](./Specifications/Dynamics/Spec_KeplerOrbit.md) 2. [RK4 Orbit Propagation](./Specifications/Dynamics/Spec_Rk4Orbit.md) - 3. [SGP4 Orbit Propagation with TLE](./Specifications/Dynamics/Spec_Sgp4Orbit.md) - 4. [ENCKE method](./Specifications/Dynamics/Spec_EnckeOrbit.md) + 3. [SGP4 Orbit Propagation with TLE](./Specifications/Dynamics/Spec_Sgp4.md) + 4. [ENCKE method](./Specifications/Dynamics/Spec_EnckeMethod.md) 5. [Relative Orbit](./Specifications/Dynamics/Spec_RelativeOrbit.md) 3. Thermal From 52f45da3dbdf241949f51c5f69659124e33529e2 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:47:49 +0100 Subject: [PATCH 10/19] fix equation in environment --- .../Environment/Spec_CelestialRotation.md | 63 ++++++++++--------- .../Environment/Spec_SRPEnvironment.md | 8 ++- 2 files changed, 40 insertions(+), 31 deletions(-) diff --git a/Specifications/Environment/Spec_CelestialRotation.md b/Specifications/Environment/Spec_CelestialRotation.md index d385e4ed..b82156a4 100644 --- a/Specifications/Environment/Spec_CelestialRotation.md +++ b/Specifications/Environment/Spec_CelestialRotation.md @@ -11,10 +11,10 @@ 3. How to use - Make an instance of the `CelestialRotation` class in `CelestialInformation` class. - Select rotation mode in `SampleSimBase.ini` - - `Idle` : no motion ($\mathrm{DCM_{ECItoECEF}}$ is set as the unit matrix) + - `Idle` : no motion ( $\mathrm{DCM_{ECItoECEF}}$ is set as the unit matrix) - If the rotation mode input is neither `Full` nor `Simple`, the `Idle` mode is set. - - `Simple` : axial rotation only ($\mathrm{DCM_{ECItoECEF}} = \bf{R}$) - - `Full` : Precession and Nutation are taken into account ($\mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P}$) + - `Simple` : axial rotation only ( $\mathrm{DCM_{ECItoECEF}} = \bf{R}$ ) + - `Full` : Precession and Nutation are taken into account ( $\mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P}$ ) - $\bf{R}$, $\bf{N}$, $\bf{P}$ stand for the DCM of axial rotation, nutation, precession, respectively. ## 2. Explanation of Algorithm @@ -51,18 +51,20 @@ - JDJ2000 = 2451545.0 [day] - JC = 36525 [day/century] - - By using tTT, we get the DCM of precession ($\bf{P}$) and nutation ($\bf{N}$) with `Precession` and `Nutation` functions. - - $\varepsilon$，$\Delta \varepsilon$，$\Delta \psi$ are calculated in `Nutation` function. + - By using tTT, we get the DCM of precession ( $\bf{P}$ ) and nutation ( $\bf{N}$ ) with `Precession` and `Nutation` functions. + - $\varepsilon$, $\Delta \varepsilon$, $\Delta \psi$ are calculated in `Nutation` function. ```math - \mathrm{E_q} = \Delta \psi \cos{(\varepsilon + \Delta \varepsilon)} \\ - \mathrm{GAST} = \mathrm{GMST} + \mathrm{E_q} + \begin{align} + \mathrm{E_q} &= \Delta \psi \cos{(\varepsilon + \Delta \varepsilon)} \\ + \mathrm{GAST} &= \mathrm{GMST} + \mathrm{E_q} + \end{align} ``` - where GAST is Greenwich Apparent Sidereal Time, GMST is Greenwich Mean Sidereal Time - GAST is calculated from julian date in `gstime` function in `src/Library/sgp4/sgp4unit.h`. - - By using GMST, We get the DCM of axial rotation ($\bf{R}$) with the `Rotation` function. The coordinate transformation from ECI to ECEF is calculated. + - By using GMST, We get the DCM of axial rotation ( $\bf{R}$ ) with the `Rotation` function. The coordinate transformation from ECI to ECEF is calculated. ```math \mathrm{DCM_{ECItoECEF}} = \bf{R}\bf{N}\bf{P} ``` @@ -79,7 +81,7 @@ - Input - Greenwich Apparent Sidereal Time (GAST) - Output - - the DCM of axial rotation ($\bf{R}$) + - the DCM of axial rotation ( $\bf{R}$ ) 3. algorithm ```math @@ -100,20 +102,22 @@ - Input - Julian century for terrestrial time (tTT) - Output - - the DCM of precession ($\bf{P}$) + - the DCM of precession ( $\bf{P}$ ) 3. algorithm - Precession angles are calculated as follows. ```math - \zeta = 2306.2181" \mathrm{tTT} + 0.30188" \mathrm{tTT}^2 + 0.017998" \mathrm{tTT}^3 \\ - \theta = 2004.3109" \mathrm{tTT} - 0.42665" \mathrm{tTT}^2 - 0.041833" \mathrm{tTT}^3 \\ - z = 2306.2181" \mathrm{tTT} + 1.09468" \mathrm{tTT}^2 + 0.018203" \mathrm{tTT}^3 \\ + \begin{align} + \zeta &= 2306.2181" \mathrm{tTT} + 0.30188" \mathrm{tTT}^2 + 0.017998" \mathrm{tTT}^3 \\ + \theta &= 2004.3109" \mathrm{tTT} - 0.42665" \mathrm{tTT}^2 - 0.041833" \mathrm{tTT}^3 \\ + z &= 2306.2181" \mathrm{tTT} + 1.09468" \mathrm{tTT}^2 + 0.018203" \mathrm{tTT}^3 \\ + \end{align} ``` ```math \bf{P} = \begin{pmatrix} \cos{(-z)} & \sin{(-z)} & 0 \\ - \- \sin{(-z)} & \cos{(-z)} & 0 \\ + - \sin{(-z)} & \cos{(-z)} & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} @@ -123,7 +127,7 @@ \end{pmatrix} \begin{pmatrix} \cos{(-\zeta)} & \sin{(-\zeta)} & 0 \\ - \- \sin{(-\zeta)} & \cos{(-\zeta)} & 0 \\ + - \sin{(-\zeta)} & \cos{(-\zeta)} & 0 \\ 0 & 0 & 1 \end{pmatrix} ``` @@ -137,21 +141,22 @@ - Input - Julian century for terrestrial time (tTT) - Output - - Return: the DCM of precession ($\bf{N}$) + - Return: the DCM of precession ( $\bf{N}$ ) - $\varepsilon$: mean obliquity of the ecliptic - $\Delta \varepsilon$: nutation in obliquity - $\Delta \psi$: nutation in longitude 3. algorithm Delaunay angles are calculated as follows. - - ```math - l = 134.96340251^\circ + 1717915923.2178"\mathrm{tTT} + 31.8792"\mathrm{tTT}^2 + 0.051635"\mathrm{tTT}^3 - 0.00024470"\mathrm{tTT}^4 \\ - l' = 357.52910918^\circ + 129596581.0481"\mathrm{tTT} - 0.5532"\mathrm{tTT}^2 + 0.000136"\mathrm{tTT}^3 - 0.00001149"\mathrm{tTT}^4 \\ - F = 93.27209062^\circ + 1739527262.8478"\mathrm{tTT} - 12.7512"\mathrm{tTT}^2 - 0.001037"\mathrm{tTT}^3 + 0.00000417"\mathrm{tTT}^4 \\ - D = 297.85019547^\circ + 1602961601.2090"\mathrm{tTT} - 6.3706"\mathrm{tTT}^2+0.006593"\mathrm{tTT}^3 -0.00003169"\mathrm{tTT}^4 \\ - \Omega = 125.04455501^\circ - 6962890.5431"\mathrm{tTT} + 7.4722"\mathrm{tTT}^2+0.007702"\mathrm{tTT}^3-0.00005939"\mathrm{tTT}^4 \\ - ``` + ```math + \begin{align} + l &= 134.96340251^\circ + 1717915923.2178"\mathrm{tTT} + 31.8792"\mathrm{tTT}^2 + 0.051635"\mathrm{tTT}^3 - 0.00024470"\mathrm{tTT}^4 \\ + l' &= 357.52910918^\circ + 129596581.0481"\mathrm{tTT} - 0.5532"\mathrm{tTT}^2 + 0.000136"\mathrm{tTT}^3 - 0.00001149"\mathrm{tTT}^4 \\ + F &= 93.27209062^\circ + 1739527262.8478"\mathrm{tTT} - 12.7512"\mathrm{tTT}^2 - 0.001037"\mathrm{tTT}^3 + 0.00000417"\mathrm{tTT}^4 \\ + D &= 297.85019547^\circ + 1602961601.2090"\mathrm{tTT} - 6.3706"\mathrm{tTT}^2+0.006593"\mathrm{tTT}^3 -0.00003169"\mathrm{tTT}^4 \\ + \Omega &= 125.04455501^\circ - 6962890.5431"\mathrm{tTT} + 7.4722"\mathrm{tTT}^2+0.007702"\mathrm{tTT}^3-0.00005939"\mathrm{tTT}^4 \\ + \end{align} + ``` - l : mean anomaly of the moon - l' : mean anomaly of the sun @@ -162,9 +167,11 @@ $\varepsilon$ and $\Delta \varepsilon$ and $\Delta \psi$ are calculated as follows. ```math - \varepsilon = 23^\circ26'21".448 - 46".8150\mathrm{tTT} - 0".00059\mathrm{tTT}^2 + 0".001813\mathrm{tTT}^3 \\ - \Delta \epsilon = 9.205\cos{\Omega} + 0.573\cos{2L'} - 0.090\cos{2\Omega} + 0.098\cos{2L}+0.007\cos{l'} - 0.001\cos{l} + 0.022\cos{(2L'+l')} + 0.013\cos{(2L+l)}-0.010\cos({2L'-l')} \, [\mathrm{arcsec}] \\ - \Delta \psi = -17.206\sin{\Omega} - 1.317\sin{2L'} + 0.207\sin{2\Omega} - 0.228\sin{2L} + 0.148\sin{l'}+0.071\sin{l}-0.052\sin{(2L'+l')} - 0.030\sin{(2L+l)}+0.022\sin{(2L'-l')} \, [\mathrm{arcsec}] \\ + \begin{align} + \varepsilon &= 23^\circ26'21".448 - 46".8150\mathrm{tTT} - 0".00059\mathrm{tTT}^2 + 0".001813\mathrm{tTT}^3 \\ + \Delta \epsilon &= 9.205\cos{\Omega} + 0.573\cos{2L'} - 0.090\cos{2\Omega} + 0.098\cos{2L}+0.007\cos{l'} - 0.001\cos{l} + 0.022\cos{(2L'+l')} + 0.013\cos{(2L+l)}-0.010\cos({2L'-l')} \, [\mathrm{arcsec}] \\ + \Delta \psi &= -17.206\sin{\Omega} - 1.317\sin{2L'} + 0.207\sin{2\Omega} - 0.228\sin{2L} + 0.148\sin{l'}+0.071\sin{l}-0.052\sin{(2L'+l')} - 0.030\sin{(2L+l)}+0.022\sin{(2L'-l')} \, [\mathrm{arcsec}] \\ + \end{align} ``` where $L = F + \Omega$，$L' = L - D$ @@ -178,7 +185,7 @@ \end{pmatrix} \begin{pmatrix} \cos{(-\Delta \psi)} & \sin{(-\Delta \psi)} & 0 \\ - \- \sin{(-\Delta \psi)} & \cos{(-\Delta \psi)} & 0 \\ + - \sin{(-\Delta \psi)} & \cos{(-\Delta \psi)} & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} diff --git a/Specifications/Environment/Spec_SRPEnvironment.md b/Specifications/Environment/Spec_SRPEnvironment.md index d3ea762d..9eec066e 100644 --- a/Specifications/Environment/Spec_SRPEnvironment.md +++ b/Specifications/Environment/Spec_SRPEnvironment.md @@ -67,9 +67,11 @@ 3. algorithm ```math - A_{\odot} = \sin^{-1}\left(\frac{r_{\odot}}{r_{\odot-sc}}\right)\\ - A_{\oplus} = \sin^{-1}\left(\frac{r_{\oplus}}{r_{\oplus-sc}}\right)\\ - \delta = \cos^{-1}\left(\frac{r_{\odot-sc}}{r_{\oplus-sc}}\cdot \boldsymbol{r}_{\oplus-sc}\cdot(\boldsymbol{r}_{\odot-sc}-\boldsymbol{r}_{\oplus-sc})\right)\\ + \begin{align} + A_{\odot} &= \sin^{-1}\left(\frac{r_{\odot}}{r_{\odot-sc}}\right)\\ + A_{\oplus} &= \sin^{-1}\left(\frac{r_{\oplus}}{r_{\oplus-sc}}\right)\\ + \delta &= \cos^{-1}\left(\frac{r_{\odot-sc}}{r_{\oplus-sc}}\cdot \boldsymbol{r}_{\oplus-sc}\cdot(\boldsymbol{r}_{\odot-sc}-\boldsymbol{r}_{\oplus-sc})\right)\\ + \end{align} ``` 4. note From 171ffe134b6f3ecbf4464aa2f11b6e1be9ad0d38 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:59:02 +0100 Subject: [PATCH 11/19] fix document format --- General/DocumentFormat.md | 24 ++++++++---------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/General/DocumentFormat.md b/General/DocumentFormat.md index 58e31d8b..4aae1a78 100644 --- a/General/DocumentFormat.md +++ b/General/DocumentFormat.md @@ -3,10 +3,7 @@ ## 0. General rule - The file name should be `Spec_CamelCase.md` in the Specifications directory. - Please use the markdown format. -- The markdown files are converted to a book with [mdBook](https://rust-lang.github.io/mdBook/) and published in Github Pages. - - The equations are rendered with [MathJax](https://www.mathjax.org/) - - e.g., Please use `\boldsymbol` for bold characters instead of `\bm`. -- So, please follow the writing rules of mdBook and MathJax. +- Because GitHub started math description ([link](https://github.blog/2022-05-19-math-support-in-markdown/)), we need to describe equations suit with the rule of GitHub and [MathJax](https://www.mathjax.org/). ## 1. Overview 1. Functions @@ -30,19 +27,14 @@ 3. Algorithm - MathJax description - please use equations as - - \\[ + $$ \dot{\boldsymbol{x}}=f(\boldsymbol{x},t) - \\] - - you can also use inline equation as \\( x=y \\) - - - Following description is automatically converted to the above style, but do not use this for new files. - - ```math - \dot{\boldsymbol{x}}=f(\boldsymbol{x},t) - ``` - - inline equation: $`x=y`$ - + $$ + or + ```math + \dot{\boldsymbol{x}}=f(\boldsymbol{x},t) + ``` + - you can also use inline equation as $ x=y $ 4. Note From 420b25b12fcb8a039e7d4a3c14c447b1c0b2b1ad Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 09:59:42 +0100 Subject: [PATCH 12/19] fix small --- General/DocumentFormat.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/General/DocumentFormat.md b/General/DocumentFormat.md index 4aae1a78..d7d2a11c 100644 --- a/General/DocumentFormat.md +++ b/General/DocumentFormat.md @@ -3,7 +3,7 @@ ## 0. General rule - The file name should be `Spec_CamelCase.md` in the Specifications directory. - Please use the markdown format. -- Because GitHub started math description ([link](https://github.blog/2022-05-19-math-support-in-markdown/)), we need to describe equations suit with the rule of GitHub and [MathJax](https://www.mathjax.org/). +- Because GitHub started to support math description ([link](https://github.blog/2022-05-19-math-support-in-markdown/)), we need to describe equations suit with the rule of GitHub and [MathJax](https://www.mathjax.org/). ## 1. Overview 1. Functions From 5ae68698344704f2b2b50a1b44eac3e4d8a03820 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:01:30 +0100 Subject: [PATCH 13/19] remove Github pages description --- README.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 238198cd..679c28c8 100644 --- a/README.md +++ b/README.md @@ -13,14 +13,13 @@ - This version is only for s2e-core developers - feature/branch-name - Writing documents before merge with the `develop` -- You can read the documents in the [GitHub Pages of this repository](https://ut-issl.github.io/s2e-documents). - If you have any questions or comments for S2E, feel free to ask us on the [discussions page of s2e-core](https://github.com/ut-issl/s2e-core/discussions). ## Index 1. General Information - 1. [Coding Convention](./General/CodingConvention.md) - 2. [Format of Documents](./General/DocumentFormat.md) + 1. [Coding Convention of S2E](./General/CodingConvention.md) + 2. [Format of S2E Documents](./General/DocumentFormat.md) 3. Mandatory set up and How to Build the S2E 1. [How to build with Visual Studio](./General/HowToCompileWithVisualStudio.md) 2. [How to build with Ubuntu in Docker](./General/HowToCompileWithUbuntuInDocker.md) for **both Windows and Mac** users From 0a1d3bdaa37857d3cf4229532ed50f332b9037e8 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:11:40 +0100 Subject: [PATCH 14/19] fix small --- General/DocumentFormat.md | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/General/DocumentFormat.md b/General/DocumentFormat.md index d7d2a11c..bea48565 100644 --- a/General/DocumentFormat.md +++ b/General/DocumentFormat.md @@ -25,16 +25,14 @@ 2. Inputs and Outputs 3. Algorithm - - MathJax description - - please use equations as - $$ - \dot{\boldsymbol{x}}=f(\boldsymbol{x},t) - $$ + - Math description + - please use equations as + $$\dot{\boldsymbol{x}}=f(\boldsymbol{x},t)$$ or ```math \dot{\boldsymbol{x}}=f(\boldsymbol{x},t) ``` - - you can also use inline equation as $ x=y $ + - you can also use inline equation as $x=y$ 4. Note From 8c2b44acc2e18f8d19a5b52041d38d14fdd72aba Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:26:39 +0100 Subject: [PATCH 15/19] fix bug in components --- README.md | 2 +- Specifications/Component/AOCS/Spec_RWJitter.md | 4 ++-- .../Component/Propulsion/Spec_SimpleThruster.md | 8 ++------ 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 679c28c8..03dc18a8 100644 --- a/README.md +++ b/README.md @@ -57,7 +57,7 @@ 3. CDH 4. Logic 5. Mission - 1. [Telescope](./Specifications/Component/Mission/Spec_Telescope_en.md) + 1. [Telescope](./Specifications/Component/Mission/Spec_Telescope_en.md) ([Japanese version](./Specifications/Component/Mission/Spec_Telescope_ja.md)) 6. Power 1. [PCU](./Specifications/Component/Power/Spec_PCU.md) 7. Propulsion diff --git a/Specifications/Component/AOCS/Spec_RWJitter.md b/Specifications/Component/AOCS/Spec_RWJitter.md index d4461d6f..c76c7d36 100644 --- a/Specifications/Component/AOCS/Spec_RWJitter.md +++ b/Specifications/Component/AOCS/Spec_RWJitter.md @@ -57,7 +57,7 @@ u(t)=\sum_{i=1}^n C_i\Omega^2\sin(2\pi h_i\Omega t+\alpha_i) ``` - - where $u(t)$ is the disturbance force and torque in Newton (N) or Newton-meters (Nm), $n$ is the number of harmonics included in the model, $C_i$ is the amplitude of the $i$th harmonic in $\mathrm{N^2/Hz}$ (or $\mathrm{(Nm)^2/Hz}$), $\Omega$ is the wheel speed in Hz, $h_i$ is the $i$th harmonic number and $\alpha_i$ is a random phase (assumed to be uniform over $[0, 2\pi]$) [1]. + - where $u(t)$ is the disturbance force and torque in Newton (N) or Newton-meters (Nm), $n$ is the number of harmonics included in the model, $C_i$ is the amplitude of the $i$ th harmonic in $\mathrm{N^2/Hz}$ (or $\mathrm{(Nm)^2/Hz}$), $\Omega$ is the wheel speed in Hz, $h_i$ is the $i$ th harmonic number and $\alpha_i$ is a random phase (assumed to be uniform over $[0, 2\pi]$) [1]. - $\alpha_i$ is generated as a uniform random number in the constructor. - When users want to use a more precise model, set `considers_structural_resonance` to ENABLE in `RW.ini` and use a model that takes structural resonance inside the RW into account. + If structural resonances are not taken into account, the RW disturbance will be underestimated, but it is not a significant change in general. @@ -72,7 +72,7 @@ - output: + jitter force and torque with structural resonance in component frame 3. algorithm - - The transfer function from disturbance by harmonics of RW without resonance ($u(t)$) to disturbance with resonance ($y(t)$) is modeled as following equation: + - The transfer function from disturbance by harmonics of RW without resonance ( $u(t)$ ) to disturbance with resonance ( $y(t)$ ) is modeled as following equation: ```math G(s)=\frac{s^2+2\zeta\omega_ns+\omega_n^2}{s^2+2d\zeta\omega_ns+\omega_n^2} ``` diff --git a/Specifications/Component/Propulsion/Spec_SimpleThruster.md b/Specifications/Component/Propulsion/Spec_SimpleThruster.md index 4fcb26d2..600224d0 100644 --- a/Specifications/Component/Propulsion/Spec_SimpleThruster.md +++ b/Specifications/Component/Propulsion/Spec_SimpleThruster.md @@ -45,10 +45,9 @@ ```math F_{thrust} = \epsilon * F_{max} + n_{f} - \tag{1} ``` - where $F_{thrust} $ is thrust magnitude, $\epsilon$ is the duty of thruster, $F_{max}$ is the maximum thrust magnitude, and $n_{f}$ is the error of thrust magnitude. + where $F_{thrust}$ is thrust magnitude, $\epsilon$ is the duty of thruster, $F_{max}$ is the maximum thrust magnitude, and $n_{f}$ is the error of thrust magnitude. Thrust direction can be calculated as follows: @@ -58,7 +57,6 @@ ```math \boldsymbol{d}_{thrust} = \boldsymbol{q}(\boldsymbol{d}_{err},n_{d}) * \boldsymbol{d}_{true} - \tag{2} ``` where $\boldsymbol{d}_{true}$ is the thrust vector without errors, $n$ is the random angles to rotate the direction of error $\boldsymbol{d}_{err}$, $\boldsymbol{d}_{x}$ is the vector which is not equal to $\boldsymbol{d}_{true}$, $n_d$ is the directional error, and $\boldsymbol{d}_{thrust}$ is the thrust vector with errors. $\boldsymbol{q}(\boldsymbol{d},n)$ is the quaternion which has the rotation axis $\boldsymbol{d}$ and the rotation angle $n$. @@ -67,7 +65,6 @@ ```math \boldsymbol{F}_{thrust} = F_{thrust} * \boldsymbol{d}_{thrust} - \tag{3} ``` where $\boldsymbol{F}_{thrust}$ is thrust. @@ -87,10 +84,9 @@ ```math \boldsymbol{T}_{thrust} = (\boldsymbol{v}_{thruster}-\boldsymbol{v}_{SC}) \times \boldsymbol{F}_{thrust} - \tag{4} ``` - where $\boldsymbol{T}_{thrust} $ is torque by the thruster, $\boldsymbol{v}_{thruster}$ is thruster position and $\boldsymbol{v}_{SC}$ is the mass center of spacecraft. + where $\boldsymbol{T}_{thrust}$ is torque by the thruster, $\boldsymbol{v}_{thruster}$ is thruster position and $\boldsymbol{v}_{SC}$ is the mass center of spacecraft. ## 3. Results of verifications 1. Case1 From 002ad7085dd278d0b3a46f72adaab366f71639f0 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:41:44 +0100 Subject: [PATCH 16/19] fix bug in thruster --- .../Component/Propulsion/Spec_SimpleThruster.md | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/Specifications/Component/Propulsion/Spec_SimpleThruster.md b/Specifications/Component/Propulsion/Spec_SimpleThruster.md index 600224d0..dfa7e5f1 100644 --- a/Specifications/Component/Propulsion/Spec_SimpleThruster.md +++ b/Specifications/Component/Propulsion/Spec_SimpleThruster.md @@ -54,12 +54,16 @@ ```math \boldsymbol{d}_{err} = \boldsymbol{q}(\boldsymbol{d}_{true},n) * \boldsymbol{d}_x ``` - ```math \boldsymbol{d}_{thrust} = \boldsymbol{q}(\boldsymbol{d}_{err},n_{d}) * \boldsymbol{d}_{true} ``` - where $\boldsymbol{d}_{true}$ is the thrust vector without errors, $n$ is the random angles to rotate the direction of error $\boldsymbol{d}_{err}$, $\boldsymbol{d}_{x}$ is the vector which is not equal to $\boldsymbol{d}_{true}$, $n_d$ is the directional error, and $\boldsymbol{d}_{thrust}$ is the thrust vector with errors. $\boldsymbol{q}(\boldsymbol{d},n)$ is the quaternion which has the rotation axis $\boldsymbol{d}$ and the rotation angle $n$. + where + - $\boldsymbol{d}_{true}$ is the thrust vector without errors + - $n$ is the random angles to rotate the direction of error $\boldsymbol{d}_{err}$ + - $\boldsymbol{d}_{x}$ is the vector which is not equal to $\boldsymbol{d}_{true}$ + - $n_d$ is the directional error, and $\boldsymbol{d}_{thrust}$ is the thrust vector with errors + - $\boldsymbol{q}(\boldsymbol{d},n)$ is the quaternion which has the rotation axis $\boldsymbol{d}$ and the rotation angle $n$ . Thrust can be calculated as follows: @@ -86,7 +90,10 @@ \boldsymbol{T}_{thrust} = (\boldsymbol{v}_{thruster}-\boldsymbol{v}_{SC}) \times \boldsymbol{F}_{thrust} ``` - where $\boldsymbol{T}_{thrust}$ is torque by the thruster, $\boldsymbol{v}_{thruster}$ is thruster position and $\boldsymbol{v}_{SC}$ is the mass center of spacecraft. + where + - $\boldsymbol{T}_{thrust}$ is torque by the thruster + - $\boldsymbol{v}_{thruster}$ is thruster position + - $\boldsymbol{v}_{SC}$ is the mass center of spacecraft. ## 3. Results of verifications 1. Case1 From df768afa3ccb1d3f3edfb8f34666dc0872ae4d55 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:46:50 +0100 Subject: [PATCH 17/19] fix small --- Specifications/Component/Propulsion/Spec_SimpleThruster.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Specifications/Component/Propulsion/Spec_SimpleThruster.md b/Specifications/Component/Propulsion/Spec_SimpleThruster.md index dfa7e5f1..79fedfe7 100644 --- a/Specifications/Component/Propulsion/Spec_SimpleThruster.md +++ b/Specifications/Component/Propulsion/Spec_SimpleThruster.md @@ -62,7 +62,8 @@ - $\boldsymbol{d}_{true}$ is the thrust vector without errors - $n$ is the random angles to rotate the direction of error $\boldsymbol{d}_{err}$ - $\boldsymbol{d}_{x}$ is the vector which is not equal to $\boldsymbol{d}_{true}$ - - $n_d$ is the directional error, and $\boldsymbol{d}_{thrust}$ is the thrust vector with errors + - $n_d$ is the directional error + - $\boldsymbol{d}_{thrust}$ is the thrust vector with errors - $\boldsymbol{q}(\boldsymbol{d},n)$ is the quaternion which has the rotation axis $\boldsymbol{d}$ and the rotation angle $n$ . Thrust can be calculated as follows: From b3d1db018691e6de7173281340d1de7f2eef9797 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 16 Nov 2022 10:52:53 +0100 Subject: [PATCH 18/19] fix small --- Specifications/Disturbance/Spec_GGTorque.md | 2 -- Specifications/Disturbance/Spec_MagDisturbance.md | 12 ++++++------ 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/Specifications/Disturbance/Spec_GGTorque.md b/Specifications/Disturbance/Spec_GGTorque.md index 579ed59b..b25ba715 100644 --- a/Specifications/Disturbance/Spec_GGTorque.md +++ b/Specifications/Disturbance/Spec_GGTorque.md @@ -26,11 +26,9 @@ Gravity gradient torque is calculated by the following equations (1) and (2). When two arguments are given, (1) is applied. When three arguments are given, (2) is used. ```math \boldsymbol{T}_{GG} = \cfrac{3 \mu}{R_0^5} \boldsymbol{R}_0 \times (\boldsymbol{I}\cdot \boldsymbol{R}_o) - \tag{1} ``` ```math \boldsymbol{T}_{GG} = \cfrac{3 \mu}{R_0^3} \boldsymbol{u} \times (\boldsymbol{I}\cdot \boldsymbol{u}) - \tag{2} ``` where $\mu$ is the gravitational constant of the Earth, $R_0$ is the distance between the Earth center and the satellite, $\boldsymbol{R_0}$ is the vector from Earth center to the satellite, $\boldsymbol{I}$ is the inertia tensor of the satellite. diff --git a/Specifications/Disturbance/Spec_MagDisturbance.md b/Specifications/Disturbance/Spec_MagDisturbance.md index 2797f96e..0fea6b74 100644 --- a/Specifications/Disturbance/Spec_MagDisturbance.md +++ b/Specifications/Disturbance/Spec_MagDisturbance.md @@ -26,10 +26,9 @@ - This function performs disturbance calculation and torque output simultaneously. The magnetic disturbance torque is added to the disturbance torque when this function is called. This process is performed in each loop of the posture calculation. 2. algorithm - Magnetic disturbance torque is calculated by the following equations (1). + Magnetic disturbance torque is calculated by the following equations. ```math \boldsymbol{T}_{mag} = \boldsymbol{M} \times \boldsymbol{B} - \tag{1} ``` where $\boldsymbol{M}$ is the residual magnetic moment in the body-fixed frame, $\boldsymbol{B}$ is the magnetic field in the body-fixed frame. @@ -47,11 +46,12 @@ - (Vector<3>) RMM: the residual magnetic moment in the body-fixed frame 3. algorithm - The residual magnetic moment is calculated by the following equations (2). + The residual magnetic moment is calculated by the following equations. ```math - \boldsymbol{M}(t_{k+1}) = \boldsymbol{M}_0 + \boldsymbol{r}(t_{k+1}) + \boldsymbol{w}_1(t_{k+1}) \\ - \boldsymbol{r}(t_{k+1}) = \boldsymbol{r}(t_{k}) +\boldsymbol{w_2}(t_{k+1}) - \tag{2} + \begin{align} + \boldsymbol{M}(t_{k+1}) &= \boldsymbol{M}_0 + \boldsymbol{r}(t_{k+1}) + \boldsymbol{w}_1(t_{k+1}) \\ + \boldsymbol{r}(t_{k+1}) &= \boldsymbol{r}(t_{k}) +\boldsymbol{w_2}(t_{k+1}) + \end{align} ``` where $\boldsymbol{M}_0$ is the average residual magnetic moment in the body-fixed frame, $\boldsymbol{r}$ is the random walk of RMM, and $\boldsymbol{w}_i \sim N([0,0,0],\Sigma_i)$ is the white noise. From 6877714ec24c393ab4565100219b72037ddd84c1 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Fri, 18 Nov 2022 20:22:25 +0100 Subject: [PATCH 19/19] fix inline equation --- Specifications/Component/AOCS/Spec_RWJitter.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Specifications/Component/AOCS/Spec_RWJitter.md b/Specifications/Component/AOCS/Spec_RWJitter.md index c76c7d36..e0101497 100644 --- a/Specifications/Component/AOCS/Spec_RWJitter.md +++ b/Specifications/Component/AOCS/Spec_RWJitter.md @@ -22,9 +22,9 @@ 3. how to use - Set the harmonics coefficients in `radial_force_harmonics_coef.csv` and `radial_torque_harmonics_coef.csv` - - The first column is an array of the $h_i$($i$-th harmonic number). The second column is an array of the $C_i$(amplitude of the $i$-th harmonic). + - The first column is an array of the $h_i$( $i$-th harmonic number). The second column is an array of the $C_i$ (amplitude of the $i$-th harmonic). - Set parameters in `RW.ini` - - When only the static imbalance and dynamic imbalance(correspond to $C_i$ at $h_i$≒1) is known according to the spec sheet, edit the files as follows. + - When only the static imbalance and dynamic imbalance(correspond to $C_i$ at $h_i\ne1$) is known according to the spec sheet, edit the files as follows. + `radial_force_harmonics_coef.csv` * Set $h_1$(the line 1 of the first column) as $1.0$. * Set $C_1$(the line 1 of the second column) as the static imbalance on the spec sheet.