def osgb36_to_wgs84(easting, northing)
  # Ellipsoid parameters
  a_airy = 6377563.396      # Semi-major axis (Airy 1830)
  b_airy = 6356256.909      # Semi-minor axis (Airy 1830)
  a_wgs84 = 6378137.0       # Semi-major axis (WGS84)
  b_wgs84 = 6356752.3142    # Semi-minor axis (WGS84)

  # OSGB36 Helmert transformation parameters to WGS84
  tx = 446.448       # Translation along X-axis (meters)
  ty = -125.157      # Translation along Y-axis (meters)
  tz = 542.06        # Translation along Z-axis (meters)
  rx = -0.1502 / 3600 * Math::PI / 180  # Rotation about X-axis (radians)
  ry = -0.2470 / 3600 * Math::PI / 180  # Rotation about Y-axis (radians)
  rz = -0.8421 / 3600 * Math::PI / 180  # Rotation about Z-axis (radians)
  s = 20.4894 * 10**-6 # Scale factor

  # Convert eastings and northings to latitude and longitude on OSGB36
  phi, lambda = en_to_lat_lon(easting, northing, a_airy, b_airy)

  # Convert to cartesian coordinates
  x_osgb36, y_osgb36, z_osgb36 = lat_lon_to_cartesian(phi, lambda, a_airy, b_airy)

  # Apply Helmert transformation to convert to WGS84 cartesian coordinates
  x_wgs84 = tx + (1 + s) * x_osgb36 + (-rz) * y_osgb36 + (ry) * z_osgb36
  y_wgs84 = ty + (rz) * x_osgb36 + (1 + s) * y_osgb36 + (-rx) * z_osgb36
  z_wgs84 = tz + (-ry) * x_osgb36 + (rx) * y_osgb36 + (1 + s) * z_osgb36

  # Convert back to latitude and longitude on WGS84
  cartesian_to_lat_lon(x_wgs84, y_wgs84, z_wgs84, a_wgs84, b_wgs84)
end

# Convert easting and northing to latitude and longitude
# based on the transverse Mercator projection
# Returns latitude (phi) and longitude (lambda)
def en_to_lat_lon(easting, northing, a, b)
  e0 = 400000       # Easting of false origin
  n0 = -100000      # Northing of false origin
  f0 = 0.9996012717 # Scale factor on the central meridian
  phi0 = 49 * Math::PI / 180  # Latitude of false origin
  lambda0 = -2 * Math::PI / 180 # Longitude of false origin

  e2 = (a**2 - b**2) / a**2 # Eccentricity squared
  n = (a - b) / (a + b)

  phi = phi0
  m = 0

  while northing - n0 - m >= 0.001
    phi = (northing - n0 - m) / (a * f0) + phi
    m = meridional_arc(phi, phi0, a, b, n)
  end

  v = a * f0 * (1 - e2 * Math.sin(phi)**2)**-0.5
  rho = a * f0 * (1 - e2) * (1 - e2 * Math.sin(phi)**2)**-1.5
  eta2 = v / rho - 1

  vii = Math.tan(phi) / (2 * rho * v)
  viii = Math.tan(phi) / (24 * rho * v**3) * (5 + 3 * Math.tan(phi)**2 + eta2 - 9 * eta2 * Math.tan(phi)**2)
  ix = Math.tan(phi) / (720 * rho * v**5) * (61 + 90 * Math.tan(phi)**2 + 45 * Math.tan(phi)**4)

  phi -= (easting - e0)**2 * vii - (easting - e0)**4 * viii + (easting - e0)**6 * ix

  x = 1 / Math.cos(phi) / v
  xi = 1 / Math.cos(phi) / (6 * v**3) * (v / rho + 2 * Math.tan(phi)**2)
  xii = 1 / Math.cos(phi) / (120 * v**5) * (5 + 28 * Math.tan(phi)**2 + 24 * Math.tan(phi)**4)

  lambda = lambda0 + (easting - e0) * x - (easting - e0)**3 * xi + (easting - e0)**5 * xii

  [phi, lambda]
end

# Calculate meridional arc
def meridional_arc(phi, phi0, a, b, n)
  n2 = n**2
  n3 = n**3
  b1 = (1 + n + 5 / 4 * n2 + 5 / 4 * n3) * (phi - phi0)
  b2 = (3 * n + 3 * n2 + 21 / 8 * n3) * Math.sin(phi - phi0) * Math.cos(phi + phi0)
  b3 = (15 / 8 * n2 + 15 / 8 * n3) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0))
  b4 = 35 / 24 * n3 * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0))

  b * a * (b1 - b2 + b3 - b4)
end

# Convert latitude and longitude to cartesian coordinates
def lat_lon_to_cartesian(phi, lambda, a, b)
  e2 = (a**2 - b**2) / a**2
  nu = a / Math.sqrt(1 - e2 * Math.sin(phi)**2)

  x = nu * Math.cos(phi) * Math.cos(lambda)
  y = nu * Math.cos(phi) * Math.sin(lambda)
  z = (nu * (1 - e2)) * Math.sin(phi)

  [x, y, z]
end

# Convert cartesian coordinates to latitude and longitude
def cartesian_to_lat_lon(x, y, z, a, b)
  e2 = (a**2 - b**2) / a**2
  p = Math.sqrt(x**2 + y**2)
  phi = Math.atan2(z, p * (1 - e2))

  phi_prev = 0
  while (phi - phi_prev).abs > 1e-12
    nu = a / Math.sqrt(1 - e2 * Math.sin(phi)**2)
    phi_prev = phi
    phi = Math.atan2(z + e2 * nu * Math.sin(phi), p)
  end

  lambda = Math.atan2(y, x)

  [phi * 180 / Math::PI, lambda * 180 / Math::PI]
end