In MATLAB, Computing a new TLE/orbit following a delta-v impulse?

Question

I have already worked on this project for a week, and I still can not solve this problem by myself. How can I update data for the SGP4 OrbitPropagator in Matlab's satellite function?

The thing I want to do its read the ISS(International Space Station) TLE file, and simulate its orbit with SGP4 Propagator. Until here, I did well, but we know the orbit will lose altitude to drag with time, so I want to give it a delta-V of 1 m/s to maintain its orbit.

I have tried three different ways from MathWorks help, but none of them were successful.

Method 1 :

ISSSGP4 = satellite(sc1, tleFile, "Name", "ISSSGP4", "OrbitPropagator", "sgp4");

The tleFile is the source data which the OrbitPropagator needs, and then I tried to update the TLE file after I give it a delta-V.

[positionSGP4, velocitySGP4] = states(satSGP4, "CoordinateFrame", "ecef");

we can use states to get current ISS velocity, add the delta-V to it, and do a coordinate transformation, using ijk2keplerian to get the Keplerian elements which we need to use in the TLE.

[a, ecc, incl, RAAN, argp, nu, truelon, arglat, lonper] = ijk2keplerian(r_ijk, v_ijk);

Mean anomaly and mean motion I can calculate by myself. The problem is, TLE requires very very rigorous format so I can't generate the valid TLE file by just using fprintf.

fid = fopen(tle_filename, 'w');
fprintf(fid, '%s\n', satellite_name);
fprintf(fid, '1 25544U 98067A   23232.99987622  .00013743  00000-0  25255-3 0  9994\n');
fprintf(fid, '2 25544  %.4f %.4f %s  %.4f %.4f %f 1847', incl,RAAN,ecc,argp,M_e,n);
fclose(fid);

And, I know trying to generate a TLE file may be a stupid idea.

Method 2:

In the satellite documentation, we can see they provide another way to do the OrbitPropagator:

satellite(scenario, semimajoraxis, eccentricity, inclination, RAAN, argofperiapsis, trueanomaly)

which means you can add a Satellite object from Keplerian elements defined in the Geocentric Celestial Reference Frame (GCRF) to the satellite scenario.

They also said

If you specify the satellite using Keplerian elements, the OrbitPropagator value can take one of these options:

"two-body-keplerian"
"sgp4"
"sdp4"

means this way should still support the sgp4 OrbitPropagator which I am using. However, when I tried to doing by myself

ISSSGP42 = satellite(sc2, semiMajorAxis, eccentricity, inclination, rightAscensionOfAscendingNode, ...
argumentOfPeriapsis, trueAnomaly, "Name", "ISSSGP4", "OrbitPropagator", "sgp4");
[positionSGP42, velocitySGP42] = states(ISSSGP42, "CoordinateFrame", "ecef");
[positionSGP4_n2, velocitySGP4_n2] = convert_positions(positionSGP42, velocitySGP42);

The operation completed but the response looks wrong to me

The original looked like this Why is the answer so different?

Method 3: In the satellite documentation, they provide a third way to use the OrbitPropagator.

satellite(scenario, positiontable, velocitytable)

but it cannot use the SGP4, so how does it work?

I know this is asking a lot. Thank you for any help!

Answer 1

You said,

"I already work on this project for a week, and I still can not solve this problem by myself."

No one could do it that fast, unless their solution was wrong. Getting it right would require a vast investment of effort. As David Hammen said in PseudoTLE creation ,

Creating TLEs is a highly nontrivial undertaking, and will take months (maybe even years) of effort to produce the required algorithms. ... You might want to rethink what you want to do, and why.

Before getting into the gory details, I want to suggest an alternate data source.

For precision work, TLEs are useless. Never use them unless you can't get anything else. If your only object of interest is the ISS, then you're in good shape, because NASA publishes much better ephemeris than you can get from a TLE at https://spotthestation.nasa.gov/trajectory_data.cfm . If you're interested in radar altimetry science missions , or GPS and other navigation systems, then spend your time learning how to use RINEX data from CDDIS (Suppose I wanted to compare TLEs to actual LEO satellite positions, what data is available? From which may it be easiest to extract X, Y, Z, T points?), but be aware of TLE and RINEX gps differences .

Now, back to answering your question. You said,

"we can use states to get current ISS velocity, and do coordinate transformation, adding deltaV and using ijk2keplerian to get the keplerian element which we need to use in TLE."

No, not that way. You could easily perform that calculation, but it will give you the wrong answer, because TLEs must be calculated in a very different way. If you take a state vector at a single point in time, and compute a Keplerian ellipse from it, you get osculating elements, describing an ellipse that is tangent to the real orbit at only that point. The variable names in TLEs sound the same, but the numbers have different meanings. TLEs use mean elements, as defined in Kozai (1959) and Brouwer (1959). To convert osculating elements into those mean elements, you have to do what Walter (1967) describes, which is not easy, as in Mean to Osculating conversion for non-J2 averaged elements , Calculated Classical Orbital Elements of the ISS seem to differ from the actual ones and ISS (ZARYA) TLE to Latitude Longitude Conversion .

"But we know the orbit will going down with time, so I want to give it a deltaV:1 m/s to maintain its orbit."

One of the many strange facts about TLEs is that although they do include atmospheric drag in the model and the published parameters, "The drag effect on eccentricity is modeled in such a way that perigee height remains constant." (see Differences between SGP8 and the standard SGP4? Is it ever used in practice?)

"I can't generate the valid TLE file by just using fprintf."

Sure you can! This is the only part that's actually easy. You just have to give it the proper field lengths, by using numbers before the dots as well as after them.

fprintf(fid, '2 25544 %8.4f %8.4f %7.0f %8.4f %8.4f %11.8f 1847', incl, RAAN, ecc*1e7, argp, M_e, n)

Note the multiplication of the eccentricity to properly deal with the absent but implied 0 and dot. You seem to be missing the final checksum digit, but happily for you STK doesn't validate it. Formatting the TLE is easy. Figuring out what numbers to put into the format is the hard part.

"I know generate TLE file may be a stupid idea."

Many have tried, and failed, because it is extremely difficult to get right. I strongly advise against it. So do How to generate TLE file? and How can I plot a satellite's orbit in 3D from a TLE using Python and Skyfield?