Am trying to do a simple test of 2-stars fix math with Skyfield python lib and Celestial Programming site at https://www.celestialprogramming.com/snippets/twoStarFix.html
So I’m getting the data of imagined observation about a 2 stars from Hipparcos using Skyfield for some geographic position at a time t. Then I’m calculating back the geographic position of observer using the math from Celestial Programming site. And I’m getting a big error, eg.:

in data: 49N, 25E
result: 48.966578719047284, 25.25305554808433

Where is the problem? Is there an error in code somewhere? Thanks for help!

The code itself:

# Greg Miller ([email protected]) 2023 http://www.celestialprogramming.com/

from skyfield.api import Topos, N, E, W, wgs84, load, Star
from skyfield.data import hipparcos
from skyfield.positionlib import position_of_radec
from math import pi, sin, cos, acos, atan2, sqrt

rad = pi/180.0

# Angles returned in radians
def rectToPolar(x, y, z):
    r = sqrt(x * x + y * y + z * z);
    lon = atan2(y, x);
    lat = acos(z / r);
    lat = pi/2 - lat # Make sure lat is in range +/-90deg
    return [lat, lon, r]

# star[0]=GP Longitude, star[1]=Declination, star[2]=Zenith Distance
def sightReduce(star1, star2):
    # From "An Analytical Solution of the Two Star Sight Problem of Celestial Navigation." 
    # James A. Van Allen

    a1 = cos(star1[0])*cos(star1[1])
    a2 = cos(star2[0])*cos(star2[1])

    b1=sin(star1[0])*cos(star1[1])
    b2=sin(star2[0])*cos(star2[1])

    c1=sin(star1[1])
    c2=sin(star2[1])

    p1=cos(star1[2])
    p2=cos(star2[2])

    # Eq 3
    A=a1*b2 - a2*b1
    B=b2*c1 - b1*c2
    C=b2*p1 - b1*p2
    D=a1*c2 - a2*c1
    E=b1*c2 - b2*c1
    F=c2*p1 - c1*p2

    # Eq 4 substituded in to eq 5
    # x^2(1+A^2/B^2 + D^2/E^2) + x(-2AC/B^2 - 2FD/E^2) + C^2/B^2 + F^2/E^2 - 1

    # Quadratic solve
    a=(1+(A*A)/(B*B)+(D*D)/(E*E))
    b=((-2*A*C)/(B*B) - (2*F*D)/(E*E))
    c=(C*C)/(B*B) + (F*F)/(E*E) - 1

    x1=(-b+sqrt(b*b - 4*a*c))/(2*a)
    x2=(-b-sqrt(b*b - 4*a*c))/(2*a)

    # eq  4
    y1=(F-D*x1)/E
    y2=(F-D*x2)/E

    z1=(C-A*x1)/B
    z2=(C-A*x2)/B

    s1=rectToPolar(x1,y1,z1)
    s2=rectToPolar(x2,y2,z2)

    return [s1[0]/rad, s1[1]/rad, s2[0]/rad, s2[1]/rad]

with load.open(hipparcos.URL) as f:
    df = hipparcos.load_dataframe(f)
planets = load('de421.bsp') # Load the JPL ephemeris DE421 (covers 1900-2050).
earth = planets['earth']

origLat, origLng = 49, 25

(txtLat, txtLatS) = (-origLat, "S") if (origLat < 0) else (origLat, "N")
(txtLng, txtLngS) = (-origLng, "W") if (origLng < 0) else (origLng, "E")
te = Topos(str(txtLat)+" "+txtLatS, str(txtLng)+" "+txtLngS)# elevation_m=300

tx = 1710621420+27
t = load.timescale().utc(1970, 1, 1, 0, 0, tx)
print(t.utc_strftime())
body = [4427, 67301]

def getHip(body):
    star = Star.from_dataframe(df.loc[body])
    pps = (earth + te).at(t).observe(star).apparent()
    alt0, az0, distance = pps.altaz()
    pps = (earth).at(t).observe(star)
    ra0, dec0, distance = pps.radec('date')
    xppq = position_of_radec(ra0.radians*12/pi, dec0.radians*180/pi, t=t, center=399)
    xq2 = wgs84.geographic_position_of(xppq)
    return [xq2.longitude.radians, dec0.radians, pi/2-alt0.radians] # [0]=GP Longitude, [1]=Declination, [2]=Zenith Distance

print(origLat, origLng)
print(sightReduce(getHip(body[0]), getHip(body[1])))

Have tried this python code example and got a big error.