xyz2philambda Subroutine

private subroutine xyz2philambda(ax, ay, b, x, y, z, phi, lambda)

Determination of the geodetic latitude and longitude

Note

This one has a different stopping criterion than the reference.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: ax
real(kind=wp), intent(in) :: ay
real(kind=wp), intent(in) :: b
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: y
real(kind=wp), intent(in) :: z
real(kind=wp), intent(inout) :: phi
real(kind=wp), intent(inout) :: lambda

Calls

proc~~xyz2philambda~~CallsGraph proc~xyz2philambda geodesy_module::xyz2philambda proc~geodetic_to_cartesian_triaxial_2 geodesy_module::geodetic_to_cartesian_triaxial_2 proc~xyz2philambda->proc~geodetic_to_cartesian_triaxial_2

Called by

proc~~xyz2philambda~~CalledByGraph proc~xyz2philambda geodesy_module::xyz2philambda proc~cartesian_to_geodetic_triaxial geodesy_module::cartesian_to_geodetic_triaxial proc~cartesian_to_geodetic_triaxial->proc~xyz2philambda

Source Code

    subroutine xyz2philambda(ax, ay, b, x, y, z, phi, lambda)

     real(wp),intent(in) :: ax, ay, b, x, y, z
     real(wp),intent(inout) :: phi, lambda !! input: initial guess, output: refined values

     real(wp) :: ee2,ex2,Sphi,Cphi,Slambda,Clambda,&
                 Den,NN,onemee2,onemex2,dndphi,dxdphi,&
                 dydphi,dzdphi,dndlam,dxdlam,dydlam,dzdlam
     integer :: n
     real(wp),dimension(3,2) :: J
     real(wp),dimension(2,3) :: Jt !! transpose of J
     real(wp),dimension(3,1) :: dl
     real(wp),dimension(2,2) :: Nmat, Ninv
     real(wp),dimension(2,1) :: dx
     real(wp),dimension(3) :: r0

    ! real(wp) :: s0, SS0
    ! real(wp),dimension(3,1) :: UU
    ! real(wp),dimension(1,1) :: tmp

     integer,parameter :: maxiter = 100 !! maximum number of iterations
     real(wp),parameter :: stop_tol = ten * epsilon(one) !! stopping tol for corrections

     ee2 = (ax*ax - ay*ay)/(ax*ax) ! eqn. 5
     ex2 = (ax*ax - b*b)/(ax*ax)   !
     onemee2 = one - ee2
     onemex2 = one - ex2

     !s0 = zero
     do n = 1, maxiter
        !SS0 = s0

        ! Design Matrix J
        Sphi = sin(phi)
        Cphi = cos(phi)
        Slambda = sin(lambda)
        Clambda = cos(lambda)

        NN  = ax/sqrt(one-ex2*Sphi*Sphi-ee2*Cphi*Cphi*Slambda*Slambda) ! eqn. 4
        Den = two*(one-ex2*Sphi**2-ee2*Cphi**2*Slambda**2)**(three/two)
        dndphi = -ax*sin(two*phi)*(ex2 - ee2*Slambda**2) / Den
        dxdphi = (dndphi*Cphi - NN*Sphi) * Clambda
        dydphi = onemee2*(dndphi*Cphi - NN*Sphi) * Slambda
        dzdphi = onemex2*(dndphi*Sphi + NN*Cphi)
        dndlam = -ax*ee2*Cphi**2*sin(two*lambda) / Den
        dxdlam = (dndlam*Clambda - NN*Slambda)*Cphi
        dydlam = onemee2*(dndlam*Slambda + NN*Clambda)*Cphi
        dzdlam = onemex2*dndlam*Sphi
        J = reshape([dxdphi,dydphi,dzdphi,dxdlam,dydlam,dzdlam],[3,2])

        ! Vector dl
        call geodetic_to_cartesian_triaxial_2(ax,ay,b,phi,lambda,0.0_wp,r0) ! just use the main one with alt=0
        dl(:,1) = [x,y,z] - r0 ! eqn. 51

        ! Solution
        Jt      = transpose(J)
        Nmat    = matmul(Jt,J) ! eqn. 53
        Ninv    = (one / (Nmat(1,1)*Nmat(2,2) - Nmat(1,2)*Nmat(2,1))) * &
                  reshape([Nmat(2,2),-Nmat(2,1),-Nmat(1,2),Nmat(1,1)], [2,2]) ! eqn. 54
        dx      = matmul(Ninv, matmul(Jt,dl)) ! eqn. 52
        phi     = phi    + dx(1,1)   ! corrections. eqn. 55
        lambda  = lambda + dx(2,1)   !

        ! ! original:
        ! UU      = matmul(J,dx) - dl
        ! tmp     = sqrt(matmul(transpose(UU),UU))
        ! s0      = tmp(1,1)
        ! if (s0 == SS0) exit

        ! JW: I think this is a better stopping criterion:
        if (all(abs(dx) <= stop_tol)) exit

     end do

    end subroutine xyz2philambda