I am currently making up my mind about revisit analysis. I am searching for clever algorithms to implement this in Matlab. I want to define targets on Earth surface and get the revisit times for each of these targets individually.

The current idea is this:

  1. define x targets on earth surface, defined in geodetic coordinates 1.1. transform the targets in cartesian coordinates earth fixed

  2. propagate the satellite's orbit, and get the state vectors in cartesian coordinates 2.1 transform them into cartesian coordinates earth fixed

  3. calculate the connecting vector between each target and the satellite position vector. 3.1. calculate the angle between this connecting vector and the satellite position vector. 3.2. check if this angle < off-nadir angle -> if yes, the target is visibile from this satellite position.

So this is roughly the idea, which I also already implemented and it works. The problem is, that it is very costy in terms of performance. As it needs to iterate through all satellite position vectors, and for each of them it needs to iterate through all targets.

So if the orbit propgagation delivers 20e^3 orbital states and we have 1e^06 targets on Earth's surface, the algo has to iterate (calculate and compare the angles) 20e^3 * 1e^06 times. Which doesn't seem to be the way to go.

So I would love to here your inputs and ideas on this. There has to be a more clever way then just iterating through all possible combinations.

Thanks and Cheers, af_ab


2 Answers 2


Too long to post as a comment, not enough for a complete answer, but hopefully contributes a bit:

Okay I see what you mean.

Calculating the origin-satellite-target angle for $10^6$ locations on Earth for $2 \times 10^4$ points along an orbital trajectory does feel unnecessarily cumbersome. The satellite's off-nadir angle defines a cone, and the Earth is an ellipsoid, and the curve defined by the intersection of those two (the satellite's instantaneous "footprint") can be calculated analytically I believe.

See for example

and all the references therein.

It certainly can for a spherical Earth at least. So then for each state vector it's just a problem of finding which sites are in the footprint.

If the satellite is in LEO then that footprint is a very small fraction of the Earth's surface, so any kind of pre-sorting/filtering will reduce the number of cases where you need to call a lot of transcendentals and FLOPs.

But setting up all those rules cleverly and correctly will take some head scratching if you can't find a paper that explains how someone did it before, thus your question in the first place.

Another angle of attack would be to consider that between one state vector and the next (assuming they're fairly closely space) most points in the foot print stay in the footprint and only the ones near the extreme edge have to be recalculated. But again a lot of head scratching to set that up cleverly and correctly.

Where this problem has been addressed ad nauseam is in graphics and rendering; how a point light source illuminates an ellipsoid or how a camera sees it. Decades of work by hundreds of applied mathematicians have contributed to computer graphics software and journals trying to find the fastest way to solve exactly these kinds of massively parallel geometrical calculations.

In those cases folks have already done the head scratching and found all the "most cleverly and correctly" implemented solutions for pre-sorting in order to minimize the FLOPs and transcendentals. However not all solutions will be purely rigorous and absolutely accurate so they would need to be examined.

My guess is that for each of your $10^6$ locations on Earth and $2 \times 10^4$ points along an orbital trajectory you store both a position and its normal.

Computer graphics folks always tend to have the normals of everything handy and ready to go.

Loop through the state vectors one by one.

For each, do a single calculation of the minimum dot product of the two normals that could potentially put locations on Earth within the satellite's footprint. Hopefully your satellite is in LEO and this number will be like 0.95 for example.

Then only for these locations do the full-blown calculation based on the ellipsoidal Earth and whatever realistic field of view your satellite has (it might not be pointing exactly at the Earth's geocenter at all times, and might not be perfectly round).

I don't think this is as clever or as time-saving as you are hoping for, but it's a start?

  • 1
    $\begingroup$ "For each, do a single calculation of the minimum dot product of the two normals that could potentially put locations on Earth within the satellite's footprint. Hopefully your satellite is in LEO and this number will be like 0.95 for example." $\endgroup$
    – af_ab
    Commented Mar 3, 2022 at 8:38
  • 1
    $\begingroup$ What do you mean with "minimum dot product"? Shouldn't it be the maximum, so the case when they are parallel to each other? So if they are parallel (or nearly) to each other target and satellite are nearly in line, so probably inside the FOV. From this I could start the actual FOV calculation and visibility check. But of course this way I still need to calculate all dot products of each target to each satellite position. Which can become a lot :) $\endgroup$
    – af_ab
    Commented Mar 3, 2022 at 8:44
  • 1
    $\begingroup$ @af_ab My focus is on throwing away as much of the Earth as possible as quickly as possible to decrease the total computation (FLOPS and transcendentals). You first calculate the "minimum dot product" that "could potentially put locations on Earth within the satellite's footprint" and throw away everything below that minimum. It might be 0.99 or 0.95 for example. It doesn't mean any site above this minimum is definitely in the FOV, but it might be so you can't throw it away. $\endgroup$
    – uhoh
    Commented Mar 3, 2022 at 8:51
  • 1
    $\begingroup$ Sure I got the idea. I like it. I just tested it briefly (with matlab build in functions). a = [1 2 3]; % satellite state vec b = randi([-6000000,6000000],1e07,3); % target state vecs tic a_norm = normr(a); b_norm = normr(b); dot = a_norm*b_norm'; toc --------------------------------- So 1e07 targets need 0.826 secs to be normalized and the dot product to a sattelite state vector to be computed. So for 20k satellite states this would in total take around 4.6 hours. Without the actual FOV and visibility check calculations. $\endgroup$
    – af_ab
    Commented Mar 3, 2022 at 9:31
  • 1
    $\begingroup$ By testing the calculation of these "core" steps of the algorithm, I can estimate if it's worth it to continue in this direction. I have to see wether it's faster to use non-matalb build in functions. $\endgroup$
    – af_ab
    Commented Mar 3, 2022 at 9:35


In this answer, I'll try to address methods you may employ to reduce the computational burden, but some will lead to an approximate solution.

For an exact solution, you may wish to check out this project I wrote for a friend a few months ago: https://gitlab.com/nyx-space/showcase/coverage . This code will propagate each of the 250 spacecraft for one day (with J2 dynamics), generate an ephemeris for each, and search for the visibility of each spacecraft for each of the 41,000 cities for every 30 seconds of that ephemeris. The code is highly parallelized which explains its blazing fast speed: all 250 spacecraft are propagated in 5.5 seconds (in total, not 5.5s per spacecraft). The visibility computation is what's long.


With such a large problem, you might want to use the principle of "a cache", which is common in web systems under heavy load. The idea is to save some computation for later to avoid having to compute the same thing again. The trade off is the memory usage to store these calculations, and with it, the likelihood that this exact calculation will be used again.

One approach consists in storing the results of some calculation in a "key-value store". This will lead to O(1) access for the calculations. For example, you could use the (X, Y, Z, Epoch) position of a spacecraft in the ECI frame as a key and the value could be the list of locations visible at that time. This would work if it's likely for two spacecraft to pass near the same location in inertial space near the same time: you might do an approximation on the position and epoch in the key, e.g. round the position to the nearest 5 km and the epoch to the nearest 4 minutes (about one degree of rotation of the Earth if my math is right). These results will be approximate, but possibly good enough for the first few analyses.

Divide and conquer

If you have a network Matlab license, you should be able to run a script on several machines at once. An approach that I commonly use for large scale Monte Carlo analyses is to dispatch the work on several computing nodes (basically remote servers that run the same code) so that each of the computers operates on their own dataset. For example, node #1 might get the first 10,000 locations, node #2 locations numbered 10,001 to 20,000 and so on. Then, you should have divided the computational time by as many nodes as you have. The burden is then to correctly aggregate all of the results.

Dynamics approximations

What dynamics are the spacecraft subjected to? Are all of these spacecraft in near circular orbits and subjected only to the Earth gravity? If so, you may skip the orbital propagation all together and use the mean anomaly of each vehicle to compute its new approximate position in the future.


In computer systems, one must always trade off accuracy for speed. In very large databases, you can compute the exact number of items by counting them all one by one (which is long); or, you can continuously store all of the items in buckets (say one bucket as 10,000 items) and when the count is requested, simply count the number of buckets which will return you the number of items rounded to the nearest 10,000 items.

The same approach should be used for such a large problem as yours.

  1. Do you need extreme precision at the cost being able to run this simulations just once in the next week? Or ...
  2. Are you still trading the orbital planes and need to run many such simulations, in which case you might want to sacrifice the precision until you can eliminate most of the orbital plane permutations.

I hope this helps.

  • 1
    $\begingroup$ The orbit propagation is absolutely not the problem. It's a matter of seconds. It really is about the geometric visibility conditions to be checked. Before I think about optimizing the computation methods (like using parallel computing etc.) I would like to think about the initial logic to apply. And optimize this. $\endgroup$
    – af_ab
    Commented Mar 3, 2022 at 8:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.