NavList:

Frank, thank you for your comprehensive reply, that's exactly the kind of feedback I was looking for in order to give me pointers for future learnings. Calculating a fix instead of plotting it is the logical next step and it doesn't seem very complicated to do, either based on the Lemoine point (for triangles) or centroid (for triangles or shapes with more than 3 sides). Something to investigate with my Sunday morning coffee. Can you point me to a reference on the subject?

Regarding your question about my software background, I am familiar with software and Github is a place that makes it easy to preserve software artifacts and control changes to them. It may seem arcane to some but it's really quite simple: if people click on the link to the repository, the home page contains a detailed description of what the program does and how it works. But to answer your question, there are two main aspects to the program: astronomical calculations and sight reduction.

Astronomical calculations include the calculation of the position of celestial bodies and other associated astronomical values. The positions of bodies are derived from Chebyshev polynomials for their heliocentric longitude, latitude, and distance, and the coefficients of those polynomials are long Fourier series. I reused the formulas from Henning Umland's superb AstroCalc program. I merely converted them from Javascript to the HP Prime's PPL language. They are very accurate when compared to the Almanac and Umland claims they are good from 1800 to 2200. Based on that data, the program is able to plot a star chart using a stereographic projection, list the visible bodies at morning and evening nautical twilight, etc.

Sight reduction functions include the ability to calculate the sextant's index error using a measure of the sun's diameter, averaging multiple sights, calculating lines of positions based on observations, translating lines of position to compute a running fix, plotting those of lines of position, etc.

The whole thing runs on a cheap but robust HP Prime scientific calculator that can live in the chart table or sextant box. You can compute a fix with just a sextant, a watch and a calculator in just a few minutes. If someone has or buys an HP Prime calculator, they can easily download my program and install it on their calculator using the instructions provided in the link above and turn it into a cheap celestial navigation computer. 

To answer your specific question about refraction, the function starts at line 1598 in this file: https://github.com/placidanomaly/HPPrimeSightReduction/blob/main/src/Sight%20Reduction.txt . It is based on a paper by G.G. Bennett which I happen to have found on this web site.

More about myself: I am an engineer and I've been working in the satellite industry for the past few years, planning and designing satellite constellations. I am quite familiar with perturbation theory and geodetic calculations (in fact I'm wondering how the actual ellipsoid shape of the earth impacts the accuracy of the simplified spherical trigonometry used in celestial navigation?). I know ephemerides can be calculated by integrating equations from the JPL but that would not be possible on a calculator so I was very happy to find the polynomial formulas on celnav.de. I have included a screen shot of another proof-of-concept program I wrote to orchestrate the flow of traffic between user terminals and ground stations in real time across an OISL-equipped satellite constellation that is currently being built.

File: