Conversion ECEF XYZ to LLH (LAT/LONG/HEIGHT) and translation back - not accurate / possible error in IronPython script

Question

I've modeled a 3D earth with gridpoints, as below:

The points are represented in 3D space as XYZ coordinates. I then convert XYZ to Lat/Long/Height(elevation) based on the script I took from here: JSFiddle

For some reason I got really strange results when trying to find XYZ of LLH not from my set, so I tried to verify the initial script by converting XYZ to LLH and then the same LLH back to XYZ to see if I get the same coordinate.

Instead, the resulting coordinate is some XYZ on earth, unrelated to the original XYZ position.

XYZ to LLH script:
Source: JSFiddle

def xyzllh(x,y,z):
    """        xyz vector  to  lat,lon,height
         output:
    
              llhvec[3] with components
    
              flat      geodetic latitude in deg
              flon      longitude in deg
              altkm     altitude in km
    """
    dtr =  math.pi/180.0
    rrnrm = [0.0] * 3
    llhvec = [0.0] * 3
    
    geodGBL()
    
    esq = EARTH_Esq

    rp = math.sqrt( x*x + y*y + z*z )
    flatgc = math.asin( z / rp )/dtr
    testval= abs(x) + abs(y)
    
    if ( testval < 1.0e-10):
        flon = 0.0
    else:
        flon = math.atan2( y,x )/dtr 
    if (flon < 0.0 ):
        flon = flon + 360.0
    
    p = math.sqrt( x*x + y*y )
    
    # on pole special case
    
    if ( p < 1.0e-10 ):
        flat = 90.0
        if ( z < 0.0 ):
            flat = -90.0
        altkm = rp - rearth(flat)
        llhvec[0]  = flat
        llhvec[1]  = flon
        llhvec[2]  = altkm
        return  llhvec

    # first iteration, use flatgc to get altitude 
    # and alt needed to convert gc to gd lat.
    
    rnow = rearth(flatgc)
    altkm = rp - rnow
    flat = gc2gd(flatgc,altkm)
    
    rrnrm = radcur(flat)
    rn = rrnrm[1]
    
    for x in range(5):
        slat = math.sin(dtr*flat)
        tangd = ( z + rn*esq*slat ) / p
        flatn = math.atan(tangd)/dtr
    
        dlat = flatn - flat
        flat = flatn
        clat = math.cos( dtr*flat )
    
        rrnrm = radcur(flat)
        rn = rrnrm[1]
    
        altkm = (p/clat) - rn
    
        if ( abs(dlat) < 1.0e-12 ):
            break

    llhvec[0]  = flat
    llhvec[1]  = flon
    llhvec[2]  = altkm
    return llhvec 
# globals
EARTH_A = 0
EARTH_B = 0
EARTH_F = 0
EARTH_Ecc = 0
EARTH_Esq = 0

# starting function do_llhxyz()
CallCount = 0
llh  = [0.0] * 3
dtr = math.pi/180

CallCount = CallCount + 1

sans = " \n"

llh = xyzllh(x,y,z)

latitude = llh[0]
longitude= llh[1]
hkm = llh[2]
height = 1000.0 * hkm

latitude = fformat(latitude,5)
longitude = fformat(longitude,5)
height = fformat(height,1)

sans = sans +"Latitude,Longitude, Height (ellipsoidal) from ECEF\n"
sans = sans + "\n"
sans = sans +"Latitude  : " + str(latitude)  + "   deg N\n"
sans = sans +"Longitude : " + str(longitude - 180) + "   deg E\n"
sans = sans +"Height    : " + str(height)    + "   m\n"

lats = []
longs = []
heights = []
lats.append(str(latitude))
longs.append(str(longitude - 180))
heights.append(str(height))

And this is the LLH to XYZ script:
Source: www.mathworks.com

a = 6378137
t = 8.1819190842622e-2

# (prime vertical radius of curvature)
N = a / math.sqrt(1 - (t*t) * (math.sin(lat)*math.sin(lat)))

x = []
y = []
z = []

# results:
x.append( ((N+height) * math.cos(lat) * math.cos(long))/1000 )
y.append( ((N+height) * math.cos(lat) * math.sin(long))/1000 )
z.append( (((1-t*t) * N + height) * math.sin(lat))/1000 )

Anyone know what I'm doing wrong here? Thanks!