CN108334683A - A kind of Analytic Calculation Method of spacecraft trajectory characteristic point - Google Patents
A kind of Analytic Calculation Method of spacecraft trajectory characteristic point Download PDFInfo
- Publication number
- CN108334683A CN108334683A CN201810069338.0A CN201810069338A CN108334683A CN 108334683 A CN108334683 A CN 108334683A CN 201810069338 A CN201810069338 A CN 201810069338A CN 108334683 A CN108334683 A CN 108334683A
- Authority
- CN
- China
- Prior art keywords
- parameter
- formula
- characteristic point
- calculated
- pseudo
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
Abstract
The present invention provides a kind of Analytic Calculation Method of spacecraft trajectory characteristic point, including：The orbit parameter of current point in time is converted, the first track six roots of sensation number is obtained；The flight time is calculated according to the trajectory characteristic point under aircraft parameter current and disome theory；In conjunction with the flight time, the parameter of pseudo-random numbers generation is calculated using the extrapolation algorithm based on intermediate orbit theory；The parameter of pseudo-random numbers generation is converted according to disome theory, obtains the second track six roots of sensation number；According to the parameter of pseudo-random numbers generation and the second track six roots of sensation number, true characteristic point parameter is calculated.Current location is obtained to the flight time of earth-satellite orbit characteristic point using disome theoretical calculation, obtains pseudo-random numbers generation parameter using the extrapolation algorithm based on intermediate orbit theory, then obtain true characteristic point parameter using disome theory with this as the starting point.Computational methods provided by the invention can greatly improve the calculating speed of trajectory characteristic point under the premise of ensuring precision compared to the prior art.
Description
Technical field
The present invention relates to spacecrafts.More particularly, to a kind of analytical Calculation of spacecraft trajectory characteristic point
Method.
Background technology
In the Track desigh of spacecraft, it is often necessary to calculate track near-earth according to the flight parameter at current time
Point, apogee will be fitted with trajectory characteristic points, the time of these trajectory characteristic points, position, speed such as equatorial plane intersection points as track
The input of suitable property inspection, therefore its calculating speed and precision have a significant impact to the efficiency of Track desigh.
Traditional computational methods generally use integral way, to meet computational accuracy, need that smaller material calculation is arranged,
This calculating speed when the flight time is shorter is fine, but when flight time is longer will cause to calculate overlong time and cannot be satisfied and want
It asks, it is especially more prominent in a large amount of Track desigh tasks.
Invention content
At least one in order to solve the above problem, the present invention provides a kind of parsing meter of spacecraft trajectory characteristic point
Calculation method, including：
The orbit parameter of current point in time is converted, the first track six roots of sensation number is obtained；
The flight time is calculated according to the trajectory characteristic point under aircraft parameter current and disome theory；
In conjunction with the flight time, the parameter of pseudo-random numbers generation is calculated using the extrapolation algorithm based on intermediate orbit theory；
The parameter of pseudo-random numbers generation is converted according to disome theory, obtains the second track six roots of sensation number；
According to the parameter of pseudo-random numbers generation and the second track six roots of sensation number, true characteristic point parameter is calculated.
Preferably, the orbit parameter of the current point in time includes the position and speed at current time.
Preferably, the first track six roots of sensation number include the first orbit inclination angle, the first right ascension of ascending node, the first eccentricity,
First semi-major axis, the first argument of perigee and the first mean anomaly.
Preferably, the second track six roots of sensation number include the second orbit inclination angle, the second right ascension of ascending node, the second eccentricity,
Second semi-major axis, the second argument of perigee and the second mean anomaly.
Preferably, first orbit inclination angle or the second orbit inclination angle are calculated using orbit inclination angle formula, and the track inclines
Angle formula is：
The right ascension of ascending node is calculated using right ascension of ascending node formula, and the right ascension of ascending node formula is：
The eccentricity is calculated using eccentricity formula, and the eccentricity formula is：
The semi-major axis is calculated using semi-major axis formula, and the semi-major axis formula is：
The argument of perigee is calculated using argument of perigee formula, and the argument of perigee formula is：
The mean anomaly is calculated using mean anomaly formula, and the mean anomaly formula is：
M0=E-esinE；
Wherein, hzFor the z-axis component of h, hx is the x-axis component of h, and hy is the y-axis component of h；The calculation formula of h is：
E is eccentric anomaly, and the calculation formula of E is：
Wherein, f is true anomaly, and the calculation formula of f is：F=μ-ω, μ are ascending node argument, and ω is argument of perigee, μ
Calculation formula be：
Wherein X is the x-axis component of position vector, and Y is the y-axis component of position vector, and Z is the z-axis component of position vector.
Preferably, the type of the trajectory characteristic point is the earth's core away from known type or mean anomaly known type.
Preferably, when the type of the trajectory characteristic point be the earth's core away from known type when, the calculation formula of the true anomaly
For：
Wherein, P is semi-latus rectum, and the calculation formula of P is：
Wherein, GM is Gravitational coefficient of the Earth；
The calculation formula of the flight time is：
Wherein, T is the orbital period, and the calculation formula of the T is：
Preferably, the step of calculating true characteristic point parameter includes：
It calculates from pseudo-random numbers generation to the residual non-uniformity of true characteristic point；And
Go out the parameter of true characteristic point according to disome theoretical calculation；
The calculation formula of the residual non-uniformity is identical as the calculation formula of flight time.
Preferably, the parametric solution model of the true characteristic point is：
Wherein,
R=a (1-emxcos Δ E+emysin Δ E)
Wherein, Δ E is eccentric anomaly increment, meets Kepler equations：
Beneficial effects of the present invention are as follows：
The present invention provides a kind of Analytic Calculation Method of spacecraft trajectory characteristic point, is obtained using disome theoretical calculation
Current location obtains pseudo-random numbers generation to the flight time of earth-satellite orbit characteristic point using the extrapolation algorithm based on intermediate orbit theory
Parameter, then with this as the starting point the parameter of true characteristic point is obtained using disome theory.Thus under the premise of ensureing computational accuracy, greatly
Width improves the calculating speed of characteristic point.
Description of the drawings
Specific embodiments of the present invention will be described in further detail below in conjunction with the accompanying drawings.
Fig. 1 shows the computational methods flow chart in one embodiment of the present invention.
Fig. 2 shows the computational accuracies of the computational methods in one embodiment of the present invention.
Fig. 3 shows the computational accuracy that the prior art uses integral algorithm to calculate.
Specific implementation mode
In order to illustrate more clearly of the present invention, the present invention is done further with reference to preferred embodiments and drawings
It is bright.Similar component is indicated with identical reference numeral in attached drawing.It will be appreciated by those skilled in the art that institute is specific below
The content of description is illustrative and be not restrictive, and should not be limited the scope of the invention with this.
The present invention provides a kind of Analytic Calculation Method of spacecraft trajectory characteristic point, as shown in Figure 1, including：
S1：The orbit parameter of current point in time is converted, the first track six roots of sensation number is obtained.
Input parameter is the lower position of WGS84 systems at current timeSpeedIt is converted under J2000 coordinate systems
PositionSpeedThen, it is track six roots of sensation number by the position at current time, rate conversion using disome theory.
Specifically, the part includes the following steps：
S11：Initial conversion formula is calculated using formula,
S12：Orbit inclination angle is calculated,
S13：Right ascension of ascending node is calculated,
S14：Eccentricity is calculated,
S15：Semi-major axis is calculated,
S16：Argument of perigee
S17：Mean anomaly is calculated,
M0=E-esinE.
Wherein, hzFor the z-axis component of h, hx is the x-axis component of h, and hy is the y-axis component of h；The calculation formula of h is：
E is eccentric anomaly, and the calculation formula of E is：
Wherein, f is true anomaly, and the calculation formula of f is：F=μ-ω, μ are ascending node argument, and ω is argument of perigee, μ
Calculation formula be：
Wherein X is the x-axis component of position vector, and Y is the y-axis component of position vector, and Z is the z-axis component of position vector.
S2：The flight time is calculated according to the trajectory characteristic point under aircraft parameter current and disome theory；
Trajectory characteristic point is divided into two classes, and one kind is the earth's core away from known type, and one kind is mean anomaly known type.For the former,
True anomaly need to first be calculated：
It obtains substituting into (9), (10) after true anomaly and obtains mean anomaly M1.
Flight time tfComputational methods：
S3：In conjunction with the flight time, pseudo-random numbers generation is calculated using the extrapolation algorithm based on intermediate orbit theory；
According to transfer time, position, the speed of aircraft are calculated using extrapolation algorithm, transfer time is to be based on disome
What theory obtained, by Perturbation Effect there are deviation, the parameter drift-out the extrapolated parameter of real features points, therefore referred to herein as
For " pseudo-random numbers generation ", calculate that step includes：
S31：By the vector x (t of current time J2000 coordinate systemi), be converted to ellipsoidal coordinates.
Specifically, converted as follows,
Wherein x (ti) by xi, yi, zi,WithComposition, and
D=(ri 2-c2)+δ(2zi+δ)
S32：Calculate first three constant of Jacobi (α1, α3, α2)。
S33：Calculate the quarternary quantic of F (ρ) and G (η).
F (ρ)=μ [c2p0(1-S0)+(ρ2+c2)(γ0ρ2+2ρ-p0)]
F (ρ)=μ γ1(γρ2+2ρ-p)(ρ2+2A1ρ+B1)
G (η)=μ [- p0(1-S0)+(1-η2)(p0+2δη+c2γ0η2)]
G (η)=μ S1p0(S+2Pη-η2)(1+P1η-Q1η2)
S34：Initialize six integral coefficient (R1, R2, R3, N1, N2, N3)。
Calculate R1, R2, R3, need AkAnd Wk, wherein A0=1, A is obtained by Factorization1And other, W0=V0=W is true
Perigee angle, W1=(W+eV1)/p, V1=sinW.
N2=D1[u+k1T2/2+k1T2/2+3k1 2T4/8+5k1 3T6/16]
Calculate N1, N2, N3, needC1k, C2kAnd Tk, wherein
T0=u
T1=1-cosu
Tk=[(k-1) Tk-2-cosusink-1U]/k, k=2 ..., 6
S35：Three constant (β after calculating Jacobi1, β2, β3)
Specifically, being calculated using following formula：
ti+β1=R1(ρi)+c2N1(ηi)
β2=-α2R2(ρi)+α2N2(ηi)
β3=φi+c2α3R3(ρi)-α3N3(ηi)
Using initial environment and initialized coefficient, it is arranged,
β1=-T time
β2=ω argument of perigees
β3=Ω right ascension of ascending node
S36：Use Jacobi constants (α1, α3, α2, β1, β2, β3) equation of motion is replaced, solve given time tfWhen ρf,
ηfAnd φf。
Specifically, being calculated using following formula：
tf+β1=R1(ρf)+c2N1(ηf)
β2=-α2R2(ρf)+α2N2(ηf)
β3=φi+c2α3R3(ρf)-α3N3(ηf)
First equation of motion is the general formula of Kepler's equations.Ellipsoidal coordinates are used as long as calculating, corresponding value is under
Formula calculates：
Wherein
S37：By time tfWhen ellipsoidal coordinates under position vector X (tf), be converted to the vector x (t of J2000 coordinate systemsf)
Wherein
S4：The parameter of pseudo-random numbers generation is converted according to disome theory, obtains the second track six roots of sensation number.
The computational methods of the second track six roots of sensation number are identical as the first track six roots of sensation number calculating method, and the present invention is no longer superfluous
It states.
It is walked according to third calculatedWithIt is converted, obtains the track six roots of sensation number of " pseudo-random numbers generation "：a1、
Ω1、i1、e1、ω1、M1。
S5：According to the parameter of pseudo-random numbers generation and the second track six roots of sensation number, true characteristic point parameter is calculated.
First, it according to the parameter of " pseudo-random numbers generation " and track six roots of sensation number, is calculated from " pseudo-random numbers generation " to " true feature
The residual non-uniformity of point ", the computational methods of residual non-uniformity are identical as the computational methods of aforementioned flight time, and the present invention is not
It repeats again.
But it is indefinite to should be noted that orbit perturbation influences, and causes " pseudo-random numbers generation " may be before " true characteristic point "
Face, it is also possible to which behind " true characteristic point ", therefore residual non-uniformity may be just, it is also possible to be negative.
Last basis " pseudo-random numbers generation " parameter and residual non-uniformity, the ginseng of " true characteristic point " is gone out according to disome theoretical calculation
Number.Solving model is as follows：
Wherein：
And：
R=a (1-emxcos Δ E+emysin Δ E)
Δ E is eccentric anomaly increment, meets Kepler equations：
The present invention provides a kind of Analytic Calculation Method of spacecraft trajectory characteristic point, is obtained using disome theoretical calculation
Current location obtains pseudo-random numbers generation to the flight time of earth-satellite orbit characteristic point using the extrapolation algorithm based on intermediate orbit theory
Parameter, then with this as the starting point the parameter of true characteristic point is obtained using disome theory.Thus under the premise of ensureing computational accuracy, greatly
Width improves the calculating speed of characteristic point.
In addition, the present invention also provides the advantageous effects that verification test proves the computational methods, as described below：
Below by taking certain spacecraft as an example, the current flight parameter of the spacecraft is shown in Table 1.
1 spacecraft flight parameter of table
| Serial number | Position X (m) | Position Y (m) | Position Z (m) | Speed X (m/s) | Speed Y (m/s) | Speed Z (m/s) |
| 1 | -2313413 | 3630020 | 4880225 | 1911.236 | -3074.557 | 4167.331 |
Integral algorithm (integration step 2s) is respectively adopted and this algorithm calculates trajectory characteristic point (by taking apogee as an example), with
The standard results that STK is calculated are compared, and precision is as shown in table 2.
2 result of calculation of table compares
According to table 2 as a result, this algorithm is suitable with the characteristic point of integral algorithm calculating error.
Characteristic point calculating has been carried out to 100 tracks of spacecraft, has been seen respectively shown in Fig. 2 and Fig. 3.It is calculated using integral
When method, with the increase of orbital flight time, the characteristic point calculating time of single track progressively increases to 8000 from 2000 microseconds
Microsecond.When using this algorithm, the characteristic point of single track calculates time about 65 microsecond to 80 microseconds, and it is winged with track to calculate the time
The row time is substantially unrelated.
From the above mentioned, it is seen that this method under the premise of ensureing computational accuracy, can greatly improve the calculating speed of characteristic point
Degree.
Obviously, the above embodiment of the present invention be only to clearly illustrate example of the present invention, and not be pair
The restriction of embodiments of the present invention may be used also on the basis of the above description for those of ordinary skill in the art
To make other variations or changes in different ways, all embodiments can not be exhaustive here, it is every to belong to this hair
Row of the obvious changes or variations that bright technical solution is extended out still in protection scope of the present invention.
Claims (10)
1. a kind of Analytic Calculation Method of spacecraft trajectory characteristic point, which is characterized in that including：
The orbit parameter of current point in time is converted, the first track six roots of sensation number is obtained；
The flight time is calculated according to the trajectory characteristic point under aircraft parameter current and disome theory；
In conjunction with the flight time, the parameter of pseudo-random numbers generation is calculated using the extrapolation algorithm based on intermediate orbit theory；
The parameter of pseudo-random numbers generation is converted according to disome theory, obtains the second track six roots of sensation number；
According to the parameter of pseudo-random numbers generation and the second track six roots of sensation number, true characteristic point parameter is calculated.
2. method according to claim 1, which is characterized in that the orbit parameter of the current point in time includes current time
Position and speed.
3. method according to claim 1, which is characterized in that the parameter of the pseudo-random numbers generation is to be calculated using extrapolation algorithm
It arrives.
4. method according to claim 2, which is characterized in that the first track six roots of sensation number includes the first orbit inclination angle, the
One right ascension of ascending node, the first eccentricity, the first semi-major axis, the first argument of perigee and the first mean anomaly.
5. method according to claim 2, which is characterized in that the second track six roots of sensation number include the second orbit inclination angle, second liter
Intersection point right ascension, the second eccentricity, the second semi-major axis, the second argument of perigee and the second mean anomaly.
6. the method according to claim 4 or 5, which is characterized in that first orbit inclination angle or the second orbit inclination angle use
Orbit inclination angle formula calculates, and the orbit inclination angle formula is：
The right ascension of ascending node is calculated using right ascension of ascending node formula, and the right ascension of ascending node formula is：
The eccentricity is calculated using eccentricity formula, and the eccentricity formula is：
The semi-major axis is calculated using semi-major axis formula, and the semi-major axis formula is：
The argument of perigee is calculated using argument of perigee formula, and the argument of perigee formula is：
The mean anomaly is calculated using mean anomaly formula, and the mean anomaly formula is：
M0=E-esinE；
Wherein, hzFor the z-axis component of h, hx is the x-axis component of h, and hy is the y-axis component of h；The calculation formula of h is：
E is eccentric anomaly, and the calculation formula of E is：
Wherein, f is true anomaly, and the calculation formula of f is：F=μ-ω, μ are ascending node argument, and ω is argument of perigee, the meter of μ
Calculating formula is：
Wherein X is the x-axis component of position vector, and Y is the y-axis component of position vector, and Z is the z-axis component of position vector.
7. method according to claim 1, which is characterized in that the type of the trajectory characteristic point is the earth's core away from known type or puts down
Anomaly known type.
8. method according to claim 7, which is characterized in that when the trajectory characteristic point type be the earth's core away from known type
When, the calculation formula of the true anomaly is：
Wherein, P is semi-latus rectum, and the calculation formula of P is：
Wherein, GM is Gravitational coefficient of the Earth；
The calculation formula of the flight time is：
Wherein, T is the orbital period, and the calculation formula of the T is：
9. according to claim 1 or any one of 7 the methods, which is characterized in that the step of calculating true characteristic point parameter wraps
It includes：
It calculates from pseudo-random numbers generation to the residual non-uniformity of true characteristic point；And
Go out the parameter of true characteristic point according to disome theoretical calculation；
The computational methods of the residual non-uniformity are identical as the computational methods of the flight time.
10. method according to claim 1, which is characterized in that the parametric solution model of the true characteristic point is：
Wherein,
R=a (1-emxcos Δ E+emysin Δ E)
Wherein, Δ E is eccentric anomaly increment, meets Kepler equations：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810069338.0A CN108334683B (en) | 2018-01-24 | 2018-01-24 | Analytic calculation method for spacecraft orbit feature points |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810069338.0A CN108334683B (en) | 2018-01-24 | 2018-01-24 | Analytic calculation method for spacecraft orbit feature points |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108334683A true CN108334683A (en) | 2018-07-27 |
| CN108334683B CN108334683B (en) | 2022-02-18 |
Family
ID=62925580
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810069338.0A Active CN108334683B (en) | 2018-01-24 | 2018-01-24 | Analytic calculation method for spacecraft orbit feature points |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108334683B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111024094A (en) * | 2019-12-23 | 2020-04-17 | 北京电子工程总体研究所 | Method for judging autonomous allowable derailment of aircraft |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101226062A (en) * | 2007-12-26 | 2008-07-23 | 北京控制工程研究所 | Method for calculating lunar orbit real-time in star |
| WO2011002534A1 (en) * | 2009-06-30 | 2011-01-06 | Georgetown Rail Equipment Company | Methods for gps milepost mapping |
| CN104484493A (en) * | 2014-10-29 | 2015-04-01 | 中国人民解放军63920部队 | Design method for airship returning back to predetermined fall point recursive orbit |
-
2018
- 2018-01-24 CN CN201810069338.0A patent/CN108334683B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101226062A (en) * | 2007-12-26 | 2008-07-23 | 北京控制工程研究所 | Method for calculating lunar orbit real-time in star |
| WO2011002534A1 (en) * | 2009-06-30 | 2011-01-06 | Georgetown Rail Equipment Company | Methods for gps milepost mapping |
| CN104484493A (en) * | 2014-10-29 | 2015-04-01 | 中国人民解放军63920部队 | Design method for airship returning back to predetermined fall point recursive orbit |
Non-Patent Citations (3)
| Title |
|---|
| MOHAMMED CHESSAB MAHDI: "ORBIT DESIGN AND SIMULATION FOR KUFASAT NANOSATELLITE", 《ARTIFICIAL SATELLITES》 * |
| SHEN HONGXIN等: "Point return orbit design and characteristics analysis for manned", 《SCIENCE CHINA》 * |
| 李丹等: "基于轨道根数的低轨卫星轨道预测算法", 《光学精密工程》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111024094A (en) * | 2019-12-23 | 2020-04-17 | 北京电子工程总体研究所 | Method for judging autonomous allowable derailment of aircraft |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108334683B (en) | 2022-02-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yagasaki | Sun-perturbed Earth-to-Moon transfers with low energy and moderate flight time | |
| Rasotto et al. | Differential algebra space toolbox for nonlinear uncertainty propagation in space dynamics | |
| CN104076819B (en) | The boundary control method of satellite bounded accompanying flying under a kind of round reference orbit | |
| Zhao et al. | Progress in reentry trajectory planning for hypersonic vehicle | |
| CN104657559B (en) | The ground moon free Entry trajectory design method based on cylindrical type speed parameter section | |
| WO2012031398A1 (en) | Numerical method for simulating subsonic flows based on euler equations in lagrangian formulation | |
| Li et al. | Spacecraft close-range trajectory planning via convex optimization and multi-resolution technique | |
| CN107423511A (en) | Meet to immerse border implicit iterative solving method without sliding boundary condition and the condition of continuity | |
| Dang | New state transition matrix for relative motion on an arbitrary Keplerian orbit | |
| CN111259325A (en) | Improved level set method based on local curvature adaptive correction | |
| CN108334683A (en) | A kind of Analytic Calculation Method of spacecraft trajectory characteristic point | |
| Oguri et al. | Stochastic sequential convex programming for robust low-thrust trajectory design under uncertainty | |
| CN113093246B (en) | Ground multi-target point imaging rapid judging and task parameter calculating method | |
| Caughey et al. | Accelerated iterative calculation of transonic nacelle flowfields | |
| Dugaro et al. | Physical properties of terrestrial planets and water delivery in the habitable zone using N-body simulations with fragmentation | |
| CN110046439B (en) | Trajectory deviation analysis forecasting algorithm considering high-order disturbance gravitation influence | |
| Greco et al. | An intrusive polynomial algebra multiple shooting approach to the solution of optimal control problems | |
| Huo et al. | Accurate approximation of in-ecliptic trajectories for E-sail with constant pitch angle | |
| Kasahara et al. | Numerical studies of frontal motion in the atmosphere‐I | |
| CN113093776B (en) | Off-orbit parameter determination method and device for spacecraft | |
| Azimov et al. | Integrated targeting and guidance for powered planetary descent | |
| Prevereaud et al. | Debris aerodynamic interactions during uncontrolled atmospheric reentry | |
| Pontani et al. | Neighboring optimal guidance and attitude control of low-thrust earth orbit transfers | |
| Bergsma et al. | Application of Taylor-series integration to reentry problems with wind | |
| Jiang | Robust optimization of Mars entry trajectory under uncertainty |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |