cartesian_to_geodetic_triaxial_2 Subroutine

    subroutine cartesian_to_geodetic_triaxial_2(a,b,c,r,eps,phi,lambda,h)

    implicit none

    real(wp),intent(in) :: a    !! ellipsoid radii `a >= b >= c`
    real(wp),intent(in) :: b    !! ellipsoid radii `a >= b >= c`
    real(wp),intent(in) :: c    !! ellipsoid radii `a >= b >= c`
    real(wp),dimension(3),intent(in) :: r !! Cartesian coordinates (x,y,z)
    real(wp),intent(in)  :: eps  !! convergence tolerance
    real(wp),intent(out) :: phi !! latitude (rad)
    real(wp),intent(out) :: lambda !! longitude (rad)
    real(wp),intent(out) :: h !! altitude

    integer,parameter :: maxiter = 20 !! maximum number of iterations

    integer :: i  !! iteration counter
    real(wp),dimension(3,3) :: AA
    real(wp),dimension(3) :: bvec, xvec
    real(wp) :: a2,b2,c2,x,y,z,ex2,ee2,e,f,g,xo,yo,zo,j11,j12,j21,j23,rmag,omee2
    logical :: success

    x = r(1)
    y = r(2)
    z = r(3)

    if (a<b .or. b<c) error stop 'error in cartesian_to_geodetic_triaxial_2: invalid a,b,c'
    call special_cases(a,b,c,x,y,z,phi,lambda,h,success)
    if (success) return

    rmag  = norm2(r)
    a2    = a*a
    b2    = b*b
    c2    = c*c
    ex2   = (a2-c2)/a2
    ee2   = (a2-b2)/a2
    omee2 = one-ee2
    E     = one/a2
    F     = one/b2
    G     = one/c2
    xo    = a*x/rmag
    yo    = b*y/rmag
    zo    = c*z/rmag

    do i = 1, maxiter

        j11 = F*yo-(yo-y)*E
        j12 = (xo-x)*F-E*xo
        j21 = G*zo-(zo-z)*E
        j23 = (xo-x)*G-E*xo

        ! solve the linear system:
        AA = reshape(-[j11,j21,two*E*xo,&
                       j12,zero,two*F*yo,&
                       zero,j23,two*G*zo], [3,3])
        bvec = [ (xo-x)*F*yo-(yo-y)*E*xo, &
                 (xo-x)*G*zo-(zo-z)*E*xo, &
                 E*xo**2+F*yo**2+G*zo**2-one ]
        call linear_solver(AA,bvec,xvec,success)
        if (.not. success) then
            write(*,*) 'error in cartesian_to_geodetic_triaxial_2: matrix is singular'
            phi    = zero
            lambda = zero
            h      = zero
            return
        end if
        xo = xo + xvec(1)
        yo = yo + xvec(2)
        zo = zo + xvec(3)

        if (maxval(abs(xvec))<eps) exit

    end do

    ! outputs:
    phi = atan(zo*omee2/(one-ex2)/sqrt(omee2**2*xo**2+yo**2))
    lambda = atan2(yo, omee2*xo)
    h = sign(one,z-zo)*sign(one,zo)*sqrt((x-xo)**2+(y-yo)**2+(z-zo)**2)

    contains

    subroutine linear_solver(a,b,x,success)

        !!  Solve the 3x3 system: `A * x = b`
        !!  Reference: https://caps.gsfc.nasa.gov/simpson/software/m33inv_f90.txt

        implicit none

        real(wp),dimension(3,3),intent(in) :: a
        real(wp),dimension(3),intent(in)   :: b
        real(wp),dimension(3),intent(out)  :: x
        logical,intent(out)                :: success

        real(wp) :: det !! determinant of a
        real(wp),dimension(3,3) :: adj !! adjoint of a
        real(wp),dimension(3,3) :: ainv !! inverse of a

        det =   a(1,1)*a(2,2)*a(3,3)  &
              - a(1,1)*a(2,3)*a(3,2)  &
              - a(1,2)*a(2,1)*a(3,3)  &
              + a(1,2)*a(2,3)*a(3,1)  &
              + a(1,3)*a(2,1)*a(3,2)  &
              - a(1,3)*a(2,2)*a(3,1)

        success = abs(det) > dblmin ! check for singularity

        if (success) then
            adj(:,1) = [a(2,2)*a(3,3)-a(2,3)*a(3,2),&
                        a(2,3)*a(3,1)-a(2,1)*a(3,3),&
                        a(2,1)*a(3,2)-a(2,2)*a(3,1)]
            adj(:,2) = [a(1,3)*a(3,2)-a(1,2)*a(3,3),&
                        a(1,1)*a(3,3)-a(1,3)*a(3,1),&
                        a(1,2)*a(3,1)-a(1,1)*a(3,2)]
            adj(:,3) = [a(1,2)*a(2,3)-a(1,3)*a(2,2),&
                        a(1,3)*a(2,1)-a(1,1)*a(2,3),&
                        a(1,1)*a(2,2)-a(1,2)*a(2,1)]
            ainv = adj/det
            x = matmul(ainv,b)
        else
            x = zero
        end if

        end subroutine linear_solver

    end subroutine cartesian_to_geodetic_triaxial_2