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?
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 :)
– af_ab Mar 03 '22 at 08:44a = [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.
– af_ab Mar 03 '22 at 09:31