The ellipsoidal stuff can be fun, but in the end the problem is that the Earth isn't really ellipsoidal either, so pretending it is doesn't get you much beyond pretending it is a sphere.
For efficient calculation of very good approximate answers on an ellipsoid, you'll want to take a look at geographiclib, the motivation and capabilities of which are described in C.F.F. Karney, "Algorithms for Geodesics", Journal of Geodesy 87(1) 43-55, 2013. I use this as an alternative for great circle route lengths, but there's really no point in worrying very hard about calculating paths that are accurate on an ellipsoid to less than a centimeter, when other things you've omitted can cause errors of kilometers.
The thing you're hunting is an equation that will give you the horizon, the locus of points at which the line of sight from your satellite to the ground just grazes the surface. The problem is that all the effects which throw you off get bigger as you get closer to the edge, so looking for the exact edge is inherently difficult. If the surface is purely convex, like an ellipsoid, then once you've been tangent to the surface, that's as far as you can go. Even so, you still have the problem of picking exactly what surface "tangent" means, and how to relate it to coordinates, as in the difference between geocentric latitude and geodetic latitude.
If, on the other hand, the curvature is more complicated, it is possible to encounter a local tangent but have that ray continue to meet the surface coming back around for another pass, if the bulge is large enough. On Earth, this kind of thing is very common when considering terrain: as you look just over the top of one hill, what you see behind it is a bigger hill, and then an even bigger hill behind that one. At some point, since the idealized ellipsoid keeps curving away and there is a highest mountain atop it, this process stops, and you can settle for a definition along the lines of "the point on the surface whose line of sight from the spacecraft intersects the surface at the largest of all cone angles for this clock angle", but then your definition is driven primarily by the messy details of staring at mountains from the side.
One of the several kinds of variation I mean is encoded in the term "geoid", which means the imaginary undulating surface we use to define "mean sea level".
You have to be careful to ask whether the precise height measurement someone is quoting uses ellipsoidal height (often referenced to WGS84, due to the influence of GPS, but not necessarily!) or orthometric height, as a surveyor calculates, relative to mean sea level. Some definitions of geoid claim it's the actual surface of the ocean, but it's not, because it has had all the tides and currents removed. If you know about them, you can add them back in, but you have to know a lot --- including the sea surface temperature, because water expands very slightly as it warms, so the water that's warmer is also higher by several centimeters than the colder water on the other side of the current. The geoid is also continued theoretically underneath the continents, where it is defined as the height the ocean would have if the continents were removed in height but their gravity was still there, so it is substantially higher than it would be if the continents just weren't there at all. In fact, the same thing happens under the ocean, for real: rock is denser than water, so underwater mountain ranges have more gravitational pull than equivalent volumes of water, so the ocean is higher above ridges and lower above trenches. The ocean surface, in that way of thinking, is a low-pass-filtered version of the sea floor.
Terrain is then added on top of that, and you need some kind of consistent source for it. The U.S. Government standard is Digital Terrain Elevation Data (DTED), as provided by https://earth-info.nga.mil/ along with the earth gravity models and geomagnetic models, but there are naturally restrictions on who can get data of just what accuracy from them. Much of the data underlying DTED was obtained by the Space Shuttle's Radar Tomography Mission, which was eventually released through NASA JPL at the equivalent resolution of the DTED Level Two that NGA won't distribute without a government permission slip. Even if you have the full roughly 8 GB dataset for global DTED-2, that still only gives you one point every thirty meters or so, which means you can't know what the height between the posts actually is --- you have to model it to interpolate it, and then you have errors from doing that. Cubic B-splines are a popular method, but I strongly suspect that is mainly because they are particularly easy, not because they are particularly accurate.
One of the fun and annoying things that then starts to happen is even if you believe you have mastered physical geography, there enters human geography. Consider trying to figure out which height in a certain region is the best one to use for the height of an antenna in your simulation. You shouldn't use the average, you should use the maximum --- because humans don't scatter antennas randomly, they intentionally seek out the locally highest point, in order to give them a better view. They also don't just leave them lying on the ground: antennas are on towers, on top of buildings, or both at once, and those buildings and towers are built right along the peak of the ridge line your horizon follows through the region, messing up its contours with their own profiles. Buildings introduce a whole catalog of horrors. For example, I rarely need to shovel snow off my sidewalk, because the street runs north-south, and is sunny most of the day; but I always have to shovel my driveway, because it is hidden from the sun by the large house next door until the sun is very high in the sky, which no longer happens in the winter. Does that matter? It definitely affects my likelihood of getting good GPS signals, just by moving my car 20 feet, which might well be part of what you have in mind: over exactly what region can some simulated satellite system satisfy its customers?
Understanding what happens at the edge of the visible region is really hard.
reference-requesttag. $\endgroup$