!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.h"

!===================================================================================================================================
!>
!!# Module ** hmap_frenet **
!!
!! This map uses the Frenet frame of a given periodic input curve X0(zeta) along the curve parameter zeta in [0,2pi].
!! It uses the signed orthonormal Frenet-Serret frame (TNB frame) that can be computed from derivatives of X0  in zeta.
!! h:  X_0(\zeta) + q_1 \sigma N(\zeta) + q_2 \sigma B(\zeta)
!! with a sign switching function \sigma(\zeta), that accounts for points of zero curvature.
!! the tangent is T=X_0' / |X_0'|, the bi-normal is B= (X_0' x X_0'') / |X_0' x X_0''|, and the normal N= B X T
!! Derivatives use the Frenet-Serret formulas:
!!
!! dT/dl = k N
!! dN/dl = -kappa T + tau B
!! dB/dl = -tau N
!!
!! With  l(\zeta) being the arc-length, and l' = |X_0'|.
!! the curvature is kappa=  |X_0' x  X_0''| / (l')^3,
!! and the torsion tau= (X_0' x X_0'').X_0''' /  |X_0' x X_0''|^2
!!
!! As a first representation of the curve X0(\zeta), we choose zeta to be the geometric toroidal angle zeta=phi, such that
!!             R0(zeta)*cos(zeta)
!!  X0(zeta)=( R0(zeta)*sin(zeta)  )
!!             Z0(zeta)
!! and both R0,Z0 are represented as a real Fourier series with modes 0... n_max and number of Field periods Nfp
!! R0(zeta) = sum_{n=0}^{n_{max}} rc(n)*cos(n*Nfp*zeta) + rs(n)*sin(n*Nfp*zeta)
!===================================================================================================================================
MODULE MODgvec_hmap_frenet
! MODULES
USE MODgvec_Globals, ONLY:PI,TWOPI,CROSS,wp,Unit_stdOut,abort,MPIroot
USE MODgvec_c_hmap,    ONLY:c_hmap, c_hmap_auxvar
IMPLICIT NONE

PUBLIC

!> Store data that can be precomputed on a set ot zeta points
!> depends on hmap_frenet, but could be used for different point sets in zeta
!
TYPE,EXTENDS(c_hmap_auxvar) :: t_hmap_frenet_auxvar
  REAL(wp)  :: lp,kappa,tau,sigma,lp_p,kappa_p,tau_p
  REAL(wp),DIMENSION(3)::X0,T,N,B
END TYPE t_hmap_frenet_auxvar

TYPE,EXTENDS(c_hmap) :: t_hmap_frenet
  !---------------------------------------------------------------------------------------------------------------------------------
  LOGICAL  :: initialized=.FALSE.
  !---------------------------------------------------------------------------------------------------------------------------------
  ! parameters for hmap_frenet:
  !curve description
  !INTEGER             :: nfp  !already part of c_hmap. Is overwritten in init
  INTEGER              :: n_max=0  !! input maximum mode number (without nfp), 0...n_max,
  REAL(wp),ALLOCATABLE :: rc(:)  !! input cosine coefficients of R0 as array (0:n_max) of modes (0,1,...,n_max)*nfp
  REAL(wp),ALLOCATABLE :: rs(:)  !! input   sine coefficients of R0 as array (0:n_max) of modes (0,1,...,n_max)*nfp
  REAL(wp),ALLOCATABLE :: zc(:)  !! input cosine coefficients of Z0 as array (0:n_max) of modes (0,1,...,n_max)*nfp
  REAL(wp),ALLOCATABLE :: zs(:)  !! input   sine coefficients of Z0 as array (0:n_max) of modes (0,1,...,n_max)*nfp
  INTEGER,ALLOCATABLE  :: Xn(:)   !! array of mode numbers,  local variable =(0,1,...,n_max)*nfp
  LOGICAL              :: omnig=.FALSE.   !! omnigenity. True: sign change of frame at pi/nfp , False: no sign change
  !---------------------------------------------------------------------------------------------------------------------------------

  CONTAINS

  FINAL :: hmap_frenet_free
  PROCEDURE :: eval_all        => hmap_frenet_eval_all
  PROCEDURE :: eval            => hmap_frenet_eval
  PROCEDURE :: eval_aux        => hmap_frenet_eval_aux
  PROCEDURE :: eval_dxdq       => hmap_frenet_eval_dxdq
  PROCEDURE :: eval_dxdq_aux   => hmap_frenet_eval_dxdq_aux
  PROCEDURE :: eval_Jh         => hmap_frenet_eval_Jh
  PROCEDURE :: eval_Jh_aux     => hmap_frenet_eval_Jh_aux
  PROCEDURE :: eval_Jh_dq      => hmap_frenet_eval_Jh_dq
  PROCEDURE :: eval_Jh_dq_aux  => hmap_frenet_eval_Jh_dq_aux
  PROCEDURE :: eval_gij        => hmap_frenet_eval_gij
  PROCEDURE :: eval_gij_aux    => hmap_frenet_eval_gij_aux
  PROCEDURE :: eval_gij_dq     => hmap_frenet_eval_gij_dq
  PROCEDURE :: eval_gij_dq_aux => hmap_frenet_eval_gij_dq_aux
  PROCEDURE :: get_dx_dqi       => hmap_frenet_get_dx_dqi
  PROCEDURE :: get_dx_dqi_aux   => hmap_frenet_get_dx_dqi_aux
  PROCEDURE :: get_ddx_dqij     => hmap_frenet_get_ddx_dqij
  PROCEDURE :: get_ddx_dqij_aux => hmap_frenet_get_ddx_dqij_aux

  !---------------------------------------------------------------------------------------------------------------------------------
  ! procedures for hmap_frenet:
  PROCEDURE :: eval_X0       => hmap_frenet_eval_X0_fromRZ
  PROCEDURE :: sigma         => hmap_frenet_sigma
END TYPE t_hmap_frenet

!INITIALIZATION FUNCTION:
INTERFACE t_hmap_frenet
  MODULE PROCEDURE hmap_frenet_init,hmap_frenet_init_params
END INTERFACE t_hmap_frenet

INTERFACE t_hmap_frenet_auxvar
  MODULE PROCEDURE hmap_frenet_init_aux
END INTERFACE t_hmap_frenet_auxvar

LOGICAL :: test_called=.FALSE.

!===================================================================================================================================

CONTAINS


!===================================================================================================================================
!> initialize the type hmap_frenet, reading from parameterfile and then call init_params
!===================================================================================================================================
FUNCTION hmap_frenet_init() RESULT(sf)
  ! MODULES
  USE MODgvec_ReadInTools, ONLY: GETLOGICAL,GETINT, GETREALARRAY
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
    TYPE(t_hmap_frenet) :: sf !! self
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
    INTEGER :: nfp,n_max,nvisu
    REAL(wp),ALLOCATABLE :: Rc(:),Rs(:),Zc(:),Zs(:)
    LOGICAL :: omnig
  !===================================================================================================================================
    SWRITE(UNIT_stdOut,'(4X,A)')'INIT HMAP :: FRENET FRAME OF A CLOSED CURVE . GET PARAMETERS:'

    nfp   = GETINT("hmap_nfp")
    n_max = GETINT("hmap_n_max")
    nvisu = GETINT("hmap_nvisu",-1)

    ALLOCATE(Rc(0:n_max)) ; Rc=0.0_wp ; Rc=GETREALARRAY("hmap_rc",n_max+1,Rc)
    ALLOCATE(Rs(0:n_max)) ; Rs=0.0_wp ; Rs=GETREALARRAY("hmap_rs",n_max+1,Rs)
    ALLOCATE(Zc(0:n_max)) ; Zc=0.0_wp ; Zc=GETREALARRAY("hmap_zc",n_max+1,Zc)
    ALLOCATE(Zs(0:n_max)) ; Zs=0.0_wp ; Zs=GETREALARRAY("hmap_zs",n_max+1,Zs)

    omnig=GETLOGICAL("hmap_omnig",.FALSE.) !omnigenity


    sf=hmap_frenet_init_params(nfp,n_max,nvisu,Rc,Rs,Zc,Zs,omnig)
    DEALLOCATE(rc,rs,zc,zs)
  END FUNCTION hmap_frenet_init

!===================================================================================================================================
!> initialize the type hmap_frenet with number of elements
!!
!===================================================================================================================================
FUNCTION hmap_frenet_init_params(nfp,n_max,nvisu,Rc,Rs,Zc,Zs,omnig) RESULT(sf)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER, INTENT(IN) :: nfp         !! number of field periods
  INTEGER, INTENT(IN) :: n_max       !! maximum mode number of guiding curve
  INTEGER, INTENT(IN) :: nvisu       !! number of visualization points per field period (-1: no visualization)
  REAL(wp),INTENT(IN) :: Rc(0:n_max) !! R cos(-n*zeta) modes of guiding curve, 0..n_max
  REAL(wp),INTENT(IN) :: Rs(0:n_max) !! R sin(-n*zeta) modes of guiding curve, 0..n_max
  REAL(wp),INTENT(IN) :: Zc(0:n_max) !! Z cos(-n*zeta) modes of guiding curve, 0..n_max
  REAL(wp),INTENT(IN) :: Zs(0:n_max) !! Z sin(-n*zeta) modes of guiding curve, 0..n_max
  LOGICAL ,INTENT(IN) :: omnig       !! omnigeneity, gives sign function of Frenet frame. False: sigma=1,
                                     !! True: sigma=+1 for 0<=zeta<=pi/nfp, and -1 for pi/nfp<zeta<2pi
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  TYPE(t_hmap_frenet) :: sf !! self
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER :: n
!===================================================================================================================================
  SWRITE(UNIT_stdOut,'(4X,A)')'INIT HMAP :: FRENET FRAME OF A CLOSED CURVE ...'

  sf%nfp=nfp
  IF(sf%nfp.LE.0) &
     CALL abort(__STAMP__, &
          "hmap_frenet init: nfp > 0 not fulfilled!")

  sf%n_max=n_max
  ALLOCATE(sf%Xn(0:sf%n_max))
  DO n=0,sf%n_max
    sf%Xn(n)=n*sf%nfp
  END DO
  ALLOCATE(sf%rc(0:sf%n_max));sf%rc=Rc
  ALLOCATE(sf%rs(0:sf%n_max));sf%rs=Rs
  ALLOCATE(sf%zc(0:sf%n_max));sf%zc=Zc
  ALLOCATE(sf%zs(0:sf%n_max));sf%zs=Zs
  sf%omnig=omnig

  IF (.NOT.(sf%rc(0) > 0.0_wp)) THEN
     CALL abort(__STAMP__, &
          "hmap_frenet init: condition rc(n=0) > 0 not fulfilled!")
  END IF

  IF(MPIroot)THEN
    IF(nvisu.GT.0) CALL VisuFrenet(sf,nvisu*sf%nfp)

    CALL CheckZeroCurvature(sf)
  END IF

  sf%initialized=.TRUE.
  SWRITE(UNIT_stdOut,'(4X,A)')'...DONE.'
  IF(.NOT.test_called) CALL hmap_frenet_test(sf)

END FUNCTION hmap_frenet_init_params

!===================================================================================================================================
!> finalize the type hmap_frenet
!!
!===================================================================================================================================
SUBROUTINE hmap_frenet_free( sf )
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  TYPE(t_hmap_frenet), INTENT(INOUT) :: sf !! self
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
!===================================================================================================================================
  IF(.NOT.sf%initialized) RETURN
  DEALLOCATE(sf%rc)
  DEALLOCATE(sf%rs)
  DEALLOCATE(sf%zc)
  DEALLOCATE(sf%zs)

  sf%initialized=.FALSE.

END SUBROUTINE hmap_frenet_free


!===================================================================================================================================
!> Sample axis and check for zero (<1.e-12) curvature
!!
!===================================================================================================================================
SUBROUTINE checkZeroCurvature( sf)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER               :: iz,nz
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,B
  REAL(wp)              :: lp,absB
  REAL(wp),DIMENSION((sf%n_max+1)*8) :: zeta,kappa
  LOGICAL ,DIMENSION((sf%n_max+1)*8) :: checkzero
!===================================================================================================================================
  nz=(sf%n_max+1)*8
  DO iz=1,nz
    zeta(iz)=REAL(iz-1,wp)/REAL(nz,wp)*TWOPI/sf%nfp  !0...2pi/nfp without endpoint
    CALL sf%eval_X0(zeta(iz),X0,X0p,X0pp,X0ppp)
    lp=SQRT(SUM(X0p*X0p))
    B=CROSS(X0p,X0pp)
    absB=SQRT(SUM(B*B))
    kappa(iz)=absB/(lp**3)
  END DO !iz
  checkzero=(kappa.LT.1.0e-8)
  IF(ANY(checkzero))THEN
    IF(sf%omnig)THEN
      !omnig=True: kappa can only be zero once, at 0,pi/nfp,[2pi/nfp...]
      IF(.NOT.(checkzero(1).AND.checkzero(nz/2+1).AND.(COUNT(checkzero).EQ.2)))THEN
        DO iz=1,nz
          IF(checkzero(iz)) WRITE(UNIT_StdOut,'(A,E15.5)')'         ...curvature <1e-8 at zeta/(2pi/nfp)=',zeta(iz)*sf%nfp/TWOPI
        END DO
        CALL abort(__STAMP__, &
             "hmap_frenet checkZeroCurvature with omnig=True: found additional points with zero curvature")
      END IF
    ELSE
      DO iz=1,nz
        IF(checkzero(iz)) WRITE(UNIT_StdOut,'(A,E15.5)')'         ...curvature <1e-8 at zeta/(2pi/nfp)=',zeta(iz)*sf%nfp/TWOPI
      END DO
      CALL abort(__STAMP__, &
           "hmap_frenet checkZeroCurvature with omnig=False: found points with zero curvature")
    END IF
  END IF
END SUBROUTINE CheckZeroCurvature

!===================================================================================================================================
!> Write evaluation of the axis and signed frenet frame to file
!!
!===================================================================================================================================
SUBROUTINE VisuFrenet( sf ,nvisu)
! MODULES
USE MODgvec_Output_CSV,     ONLY: WriteDataToCSV
USE MODgvec_Output_vtk,     ONLY: WriteDataToVTK
USE MODgvec_Output_netcdf,     ONLY: WriteDataToNETCDF
USE MODgvec_Analyze_vars,     ONLY: outfileType
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  INTEGER             , INTENT(IN) :: nvisu     !!
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,T,N,B
  REAL(wp)              :: zeta,sigma,lp,absB,kappa,tau,eps
  INTEGER               :: iVar,ivisu,itest
  INTEGER,PARAMETER     :: nVars=26
  CHARACTER(LEN=20)     :: VarNames(1:nVars)
  REAL(wp)              :: values(1:nVars,1:nvisu*sf%nfp+1)
!===================================================================================================================================
  IF(nvisu.LE.0) RETURN
  iVar=0
  VarNames(ivar+1:iVar+3)=(/ "x", "y", "z"/);iVar=iVar+3
  VarNames(ivar+1:iVar+3)=(/"TX","TY","TZ"/);iVar=iVar+3
  VarNames(ivar+1:iVar+3)=(/"NX","NY","NZ"/);iVar=iVar+3
  VarNames(ivar+1:iVar+3)=(/"BX","BY","BZ"/);iVar=iVar+3
  VarNames(iVar+1       )="zeta_norm"       ;iVar=iVar+1
  VarNames(iVar+1       )="sigma_sign"      ;iVar=iVar+1
  VarNames(iVar+1       )="lprime"          ;iVar=iVar+1
  VarNames(iVar+1       )="kappa"           ;iVar=iVar+1
  VarNames(iVar+1       )="tau"             ;iVar=iVar+1
  VarNames(ivar+1:iVar+3)=(/ "X0pX", "X0pY", "X0pZ"/);iVar=iVar+3
  VarNames(ivar+1:iVar+3)=(/ "X0ppX", "X0ppY", "X0ppZ"/);iVar=iVar+3
  VarNames(ivar+1:iVar+3)=(/ "X0pppX", "X0pppY", "X0pppZ"/);iVar=iVar+3

!  values=0.
  DO ivisu=1,nvisu*sf%nfp+1
    eps=0.
    kappa=0.
    itest=0
    DO WHILE(kappa.LT.1.0e-6)! for tau being meaningful
      zeta=(REAL(ivisu-1,wp)+eps)/REAL(nvisu*sf%nfp,wp)*TWOPI
      CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
      lp=SQRT(SUM(X0p*X0p))
      T=X0p/lp
      B=CROSS(X0p,X0pp)
      absB=SQRT(SUM(B*B))
      kappa=absB/(lp**3)
      itest=itest+1
      eps=10**REAL(-16+itest,wp)
      IF(itest.EQ.15)THEN !meaningful kappa not found
        B=0.
        kappa=1.0e-6 !-12
        absB=1.
      END IF
    END DO
    tau=SUM(X0ppp*B)/(absB**2)
    B=B/absB
    N=CROSS(B,T)
    sigma=sf%sigma(zeta)
    iVar=0
    values(ivar+1:iVar+3,ivisu)=X0                ;iVar=iVar+3
    values(ivar+1:iVar+3,ivisu)=T                 ;iVar=iVar+3
    values(ivar+1:iVar+3,ivisu)=N*sigma           ;iVar=iVar+3
    values(ivar+1:iVar+3,ivisu)=B*sigma           ;iVar=iVar+3
    values(iVar+1       ,ivisu)=zeta*sf%nfp/TWOPI ;iVar=iVar+1
    values(iVar+1       ,ivisu)=sigma             ;iVar=iVar+1
    values(iVar+1       ,ivisu)=lp                ;iVar=iVar+1
    values(iVar+1       ,ivisu)=kappa             ;iVar=iVar+1
    values(iVar+1       ,ivisu)=tau               ;iVar=iVar+1
    values(ivar+1:iVar+3,ivisu)=X0p               ;iVar=iVar+3
    values(ivar+1:iVar+3,ivisu)=X0pp              ;iVar=iVar+3
    values(ivar+1:iVar+3,ivisu)=X0ppp             ;iVar=iVar+3
  END DO !ivisu
  IF((outfileType.EQ.1).OR.(outfileType.EQ.12))THEN
    CALL WriteDataToVTK(1,3,nVars-3,(/nvisu*sf%nfp/),1,VarNames(4:nVars),values(1:3,:),values(4:nVars,:),"visu_hmap_frenet.vtu")
  END IF
  IF((outfileType.EQ.2).OR.(outfileType.EQ.12))THEN
#if NETCDF
    CALL WriteDataToNETCDF(1,3,nVars-3,(/nvisu*sf%nfp/),(/"dim_zeta"/),VarNames(4:nVars),values(1:3,:),values(4:nVars,:), &
         "visu_hmap_frenet")
#else
    CALL WriteDataToCSV(VarNames(:) ,values, TRIM("out_visu_hmap_frenet.csv") ,append_in=.FALSE.)
#endif
  END IF
END SUBROUTINE VisuFrenet


!===================================================================================================================================
!> initialize the aux variable
!!
!===================================================================================================================================
FUNCTION hmap_frenet_init_aux( sf ,zeta,do_2nd_der) RESULT(xv)
! MODULES
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet),INTENT(IN) :: sf !! self
  REAL(wp)            ,INTENT(IN) :: zeta
  LOGICAL             ,INTENT(IN) :: do_2nd_der !! compute second derivative and store second derivative terms
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  TYPE(t_hmap_frenet_auxvar)      :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0p,X0pp,X0ppp,X0p4,Bloc
  REAL(wp)  :: absB,absB_p
!===================================================================================================================================
  xv%do_2nd_der=do_2nd_der
  xv%zeta = zeta
  IF(xv%do_2nd_der)THEN
    CALL sf%eval_X0(zeta, xv%X0, X0p, X0pp, X0ppp,X0p4=X0p4)
  ELSE
    CALL sf%eval_X0(zeta, xv%X0, X0p, X0pp, X0ppp)
  END IF

  xv%lp = SQRT(SUM(X0p*X0p))
  xv%T = X0p / xv%lp
  Bloc = CROSS(X0p, X0pp)
  absB = SQRT(SUM(Bloc*Bloc))
  xv%kappa = absB / (xv%lp**3)
  IF(xv%kappa.LT.1.0e-8) &
      CALL abort(__STAMP__, &
          "hmap_frenet cannot evaluate frame at curvature < 1e-8 !", RealInfo=zeta)
  xv%sigma = sf%sigma(zeta)
  xv%tau = SUM(X0ppp*Bloc) / (absB**2)
  xv%B = Bloc / absB
  xv%N = CROSS(xv%B, xv%T)

  IF(xv%do_2nd_der)THEN
    xv%lp_p = SUM(X0pp*X0p) / xv%lp
    absB_p = SUM(Bloc*CROSS(X0p, X0ppp)) / absB
    xv%kappa_p = (absB_p*xv%lp -3*absB * xv%lp_p) / (xv%lp**4)
    xv%tau_p   = (SUM(X0p4*Bloc)*absB -2*SUM(X0ppp*Bloc)*absB_p) / (absB**3)
  END IF

  END FUNCTION hmap_frenet_init_aux


!===================================================================================================================================
!> evaluate all metrics necessary for optimizer
!!
!===================================================================================================================================
SUBROUTINE hmap_frenet_eval_all(sf,ndims,dim_zeta,xv,&
                                q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
                                Jh,    g_tt,    g_tz,    g_zz,&
                                Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1, &
                                Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2, &
                                g_t1,g_t2,g_z1,g_z2,Gh11,Gh22  )
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN)   :: sf
  INTEGER             , INTENT(IN)   :: ndims(3)    !! 3D dimensions of input arrays
  INTEGER             , INTENT(IN)   :: dim_zeta    !! which dimension is zeta dependent
  CLASS(c_hmap_auxvar), INTENT(IN)   :: xv(ndims(dim_zeta))  !! zeta point positions
  REAL(wp),DIMENSION(ndims(1),ndims(2),ndims(3)),INTENT(IN) :: q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp),DIMENSION(ndims(1),ndims(2),ndims(3)),INTENT(OUT):: Jh,g_tt    ,g_tz    ,g_zz    , &
                                                               Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1, &
                                                               Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2, &
                                                               g_t1,g_t2,g_z1,g_z2,Gh11,Gh22
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER :: i,j,k
  !===================================================================================================================================
  SELECT TYPE(xv)
  TYPE IS(t_hmap_frenet_auxvar)
    SELECT CASE(dim_zeta)
    CASE(1)
      !$OMP PARALLEL DO COLLAPSE(3) SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i,j,k)
      DO k=1,ndims(3); DO j=1,ndims(2); DO i=1,ndims(1)
        CALL hmap_frenet_eval_all_e(xv(i), &
                 q1(i,j,k),q2(i,j,k),dX1_dt(i,j,k),dX2_dt(i,j,k),dX1_dz(i,j,k),dX2_dz(i,j,k), &
                 Jh(i,j,k)    ,g_tt(i,j,k)    ,g_tz(i,j,k)    ,g_zz(i,j,k), &
                 Jh_dq1(i,j,k),g_tt_dq1(i,j,k),g_tz_dq1(i,j,k),g_zz_dq1(i,j,k), &
                 Jh_dq2(i,j,k),g_tt_dq2(i,j,k),g_tz_dq2(i,j,k),g_zz_dq2(i,j,k), &
                 g_t1(i,j,k),g_t2(i,j,k),g_z1(i,j,k),g_z2(i,j,k),Gh11(i,j,k),Gh22(i,j,k) )
      END DO; END DO; END DO
      !$OMP END PARALLEL DO
    CASE(2)
       !$OMP PARALLEL DO COLLAPSE(3) SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i,j,k)
      DO k=1,ndims(3); DO j=1,ndims(2); DO i=1,ndims(1)
        CALL hmap_frenet_eval_all_e(xv(j), &
                 q1(i,j,k),q2(i,j,k),dX1_dt(i,j,k),dX2_dt(i,j,k),dX1_dz(i,j,k),dX2_dz(i,j,k), &
                 Jh(i,j,k)    ,g_tt(i,j,k)    ,g_tz(i,j,k)    ,g_zz(i,j,k), &
                 Jh_dq1(i,j,k),g_tt_dq1(i,j,k),g_tz_dq1(i,j,k),g_zz_dq1(i,j,k), &
                 Jh_dq2(i,j,k),g_tt_dq2(i,j,k),g_tz_dq2(i,j,k),g_zz_dq2(i,j,k), &
                 g_t1(i,j,k),g_t2(i,j,k),g_z1(i,j,k),g_z2(i,j,k),Gh11(i,j,k),Gh22(i,j,k) )
      END DO; END DO; END DO
       !$OMP END PARALLEL DO
    CASE(3)
       !$OMP PARALLEL DO COLLAPSE(3) SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i,j,k)
      DO k=1,ndims(3); DO j=1,ndims(2); DO i=1,ndims(1)
        CALL hmap_frenet_eval_all_e(xv(k), &
                 q1(i,j,k),q2(i,j,k),dX1_dt(i,j,k),dX2_dt(i,j,k),dX1_dz(i,j,k),dX2_dz(i,j,k), &
                 Jh(i,j,k)    ,g_tt(i,j,k)    ,g_tz(i,j,k)    ,g_zz(i,j,k), &
                 Jh_dq1(i,j,k),g_tt_dq1(i,j,k),g_tz_dq1(i,j,k),g_zz_dq1(i,j,k), &
                 Jh_dq2(i,j,k),g_tt_dq2(i,j,k),g_tz_dq2(i,j,k),g_zz_dq2(i,j,k), &
                 g_t1(i,j,k),g_t2(i,j,k),g_z1(i,j,k),g_z2(i,j,k),Gh11(i,j,k),Gh22(i,j,k) )
      END DO; END DO; END DO
       !$OMP END PARALLEL DO
    END SELECT
  END SELECT !type(xv)

END SUBROUTINE hmap_frenet_eval_all

!===================================================================================================================================
!> evaluate all quantities at one given point (elemental)
!!
!===================================================================================================================================
PURE SUBROUTINE hmap_frenet_eval_all_e(xv,q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
                                       Jh,    g_tt,    g_tz,    g_zz,     &
                                       Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1, &
                                       Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2, &
                                       g_t1,g_t2,g_z1,g_z2,Gh11,Gh22  )
! MODULES
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  TYPE(t_hmap_frenet_auxvar),INTENT(IN) :: xv    !! precomputed auxiliary variables
  REAL(wp),INTENT(IN)  :: q1,q2       !! solution variables q1,q2
  REAL(wp),INTENT(IN)  :: dX1_dt,dX2_dt  !! theta derivative of solution variables q1,q2
  REAL(wp),INTENT(IN)  :: dX1_dz,dX2_dz  !!  zeta derivative of solution variables q1,q2
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp),INTENT(OUT) :: Jh,g_tt,g_tz,g_zz              !! Jac,1/Jac,g_{ab} with a=theta/zeta b=theta/zeta
  REAL(wp),INTENT(OUT) :: Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1  !! and their variation vs q1
  REAL(wp),INTENT(OUT) :: Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2  !! and their variation vs q2
  REAL(wp),INTENT(OUT) :: g_t1,g_t2,g_z1,g_z2,Gh11,Gh22  !! dq^{i}/dtheta*G^{i1}, dq^{i}/dtheta*G^{i2}, and dq^{i}/dzeta*G^{i1}, dq^{i}/dzeta*G^{i2} and G^{11},G^{22}
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp) :: Gh31,Gh32,Gh33
!===================================================================================================================================
  ASSOCIATE(  lp=>xv%lp, tau=>xv%tau, sigma=>xv%sigma, kappa=>xv%kappa)
  Gh11=1.0_wp
  !Gh21=0.0_wp
  Gh22=1.0_wp
  Gh31 =-lp*tau*q2
  Gh32 = lp*tau*q1

  !Jh=lp*(1.0_wp-sigma*kappa*q1)
  Jh_dq1=-lp*sigma*kappa
  Jh_dq2=0.0_wp

  Jh=lp+Jh_dq1*q1
  ! Gh33 = (lp**2)*((1.0_wp-sigma*kappa*q1)**2+tau**2*(q1**2+q2**2))
  Gh33 = Jh*Jh + Gh31*Gh31 + Gh32*Gh32

  g_t1 = dX1_dt
  g_t2 = dX2_dt
  g_z1 = dX1_dz + Gh31
  g_z2 = dX2_dz + Gh32

  g_tt =   dX1_dt *  g_t1         +  dX2_dt *  g_t2
  g_tz =   dX1_dt *  g_z1         +  dX2_dt *  g_z2
  g_zz =   dX1_dz * (g_z1 + Gh31) +  dX2_dz * (g_z2 + Gh32)  + Gh33

  !Gh11/dq1 =0 Gh12/dq1 =0 Gh13/dq1 = 0
  !            Gh22/dq1 =0 Gh23/dq1 = lp*tau
  !                        Gh33/dq1 = 2*(lp**2)*((1.0_wp-sigma*kappa*q1)*(-sigma*kappa)+tau**2 *(q1))
  !Gh11/dq2 =0 Gh12/dq2 =0 Gh13/dq2 = -lp*tau
  !            Gh22/dq2 =0 Gh23/dq2 = 0
  !                        Gh33/dq2 = 2*(lp*tau)**2*(q2)
  ! => g_t1 /dq1 =0, g_t1/dq2 =0, g_t2/dq1 =0, g_t2/dq2 =0
  ! => g_z1 /dq1 = Gh31/dq1, g_z1/dq2 =Gh31/dq2, g_z2/dq1 =Gh32/dq1, g_z2/dq2 =Gh32/dq2
  g_tt_dq1 = 0.0_wp
  g_tt_dq2 = 0.0_wp

  g_tz_dq1 =  lp*tau*dX2_dt
  g_tz_dq2 = -lp*tau*dX1_dt

  g_zz_dq1 =  2.0_wp*(lp*tau*(dX2_dz + Gh32)+Jh*Jh_dq1)
  g_zz_dq2 = -2.0_wp*lp*tau*(dX1_dz + Gh31)
  END ASSOCIATE
END SUBROUTINE hmap_frenet_eval_all_e


!===================================================================================================================================
!> sign function depending on zeta,
!! if omnig=False, sigma=1
!! if omnig=True, sigma=+1 for 0<=zeta<=pi/nfp, and -1 for pi/nfp<zeta<2pi
!!
!===================================================================================================================================
FUNCTION hmap_frenet_sigma(sf,zeta) RESULT(sigma)
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: zeta
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: sigma
!===================================================================================================================================
  sigma=MERGE(SIGN(1.0_wp,SIN(sf%nfp*zeta)),1.0_wp,sf%omnig)
END FUNCTION hmap_frenet_sigma

!===================================================================================================================================
!> evaluate the mapping h (q1,q2,zeta) -> (x,y,z)
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval( sf ,q_in) RESULT(x_out)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: x_out(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,T,N,B
  REAL(wp)          :: lp,absB,kappa,sigma,Jh
!===================================================================================================================================
  ! q(:) = (q1,q2,zeta) are the variables in the domain of the map
  ! X(:) = (x,y,z) are the variables in the range of the map
  !
  !  |x |
  !  |y |=  X0(zeta) + sigma*(N(zeta)*q1 + B(zeta)*q2)
  !  |z |

  ASSOCIATE(q1=>q_in(1),q2=>q_in(2),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  T=X0p/lp
  B=CROSS(X0p,X0pp)
  absB=SQRT(SUM(B*B))
  kappa=absB/(lp**3)
  IF(kappa.LT.1.0e-8) &
      CALL abort(__STAMP__, &
           "hmap_frenet cannot evaluate frame at curvature < 1e-8 !",RealInfo=zeta*sf%nfp/TWOPI)
  sigma=sf%sigma(zeta)
  Jh=lp*(1.0_wp-sigma*kappa*q1)
  IF(Jh.LT.1.0e-12) &
      CALL abort(__STAMP__, &
           "hmap_frenet, evaluation outside curvature radius (sigma*q1 >= 1./(kappa))",RealInfo=zeta*sf%nfp/TWOPI)
  B=B/absB
  N=CROSS(B,T)
  x_out=X0 +sigma*(q1*N + q2*B)
  END ASSOCIATE
END FUNCTION hmap_frenet_eval


!===================================================================================================================================
!> evaluate the mapping h (q1,q2,zeta) -> (x,y,z)
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_aux( sf ,q1,q2,xv) RESULT(x_out)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: x_out(3)
!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  x_out=xv%X0 +xv%sigma*(q1*xv%N + q2*xv%B)
  END SELECT !type(xv)
END FUNCTION hmap_frenet_eval_aux


!===================================================================================================================================
!> evaluate total derivative of the mapping  sum k=1,3 (dx(1:3)/dq^k) q_vec^k,
!! where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_dxdq( sf ,q_in,q_vec) RESULT(dxdq_qvec)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
  REAL(wp)            , INTENT(IN) :: q_vec(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                        :: dxdq_qvec(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,T,N,B
  REAL(wp)          :: lp,absB,kappa,tau,sigma,Jh
!===================================================================================================================================
  !  |x |
  !  |y |=  X0(zeta) + sigma*(N(zeta)*q1 + B(zeta)*q2)
  !  |z |
  !  dh/dq1 =sigma*N , dh/dq2=sigma*B
  !  dh/dq3 = l' [(1-sigma*kappa*q1)T + sigma*tau*(B*q1-N*q2) ]
  ASSOCIATE(q1=>q_in(1),q2=>q_in(2),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  T=X0p/lp
  B=CROSS(X0p,X0pp)
  absB=SQRT(SUM(B*B))
  kappa=absB/(lp**3)
  IF(kappa.LT.1.0e-8) &
      CALL abort(__STAMP__, &
           "hmap_frenet cannot evaluate frame at curvature < 1e-8 !",RealInfo=zeta*sf%nfp/TWOPI)

  sigma=sf%sigma(zeta)
  Jh=lp*(1.0_wp-sigma*kappa*q1)
  IF(Jh.LT.1.0e-12) &
      CALL abort(__STAMP__, &
           "hmap_frenet, evaluation outside curvature radius (sigma*q1 >= 1/(kappa))",RealInfo=zeta*sf%nfp/TWOPI)

  tau=SUM(X0ppp*B)/(absB**2)
  B=B/absB
  N=CROSS(B,T)
  dxdq_qvec(1:3)= sigma*(N*q_vec(1)+B*q_vec(2))+(Jh*T +sigma*lp*tau*(B*q1-N*q2))*q_vec(3)

  END ASSOCIATE !zeta
END FUNCTION hmap_frenet_eval_dxdq


!===================================================================================================================================
!> evaluate total derivative of the mapping  sum k=1,3 (dx(1:3)/dq^k) q_vec^k,
!! where dx(1:3)/dq^k, k=1,2,3 is evaluated at q_in=(X^1,X^2,zeta) ,
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_dxdq_aux( sf ,q1,q2,q1_vec,q2_vec,q3_vec,xv) RESULT(dxdq_qvec)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  REAL(wp)            , INTENT(IN) :: q1_vec,q2_vec,q3_vec
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                        :: dxdq_qvec(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp)          :: Jh
!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  Jh=xv%lp*(1.0_wp-xv%sigma*xv%kappa*q1)
  dxdq_qvec(1:3)= xv%sigma*(xv%N*q1_vec+xv%B*q2_vec) &
                  +(Jh*xv%T +xv%sigma*xv%lp*xv%tau*(xv%B*q1-xv%N*q2))*q3_vec
  END SELECT !type(xv)
END FUNCTION hmap_frenet_eval_dxdq_aux

!===============================================================================================================================
!> evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),
!!
!===============================================================================================================================
SUBROUTINE hmap_frenet_get_dx_dqi( sf ,q_in,dx_dq1,dx_dq2,dx_dq3)
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
  REAL(wp),DIMENSION(3),INTENT(OUT) :: dx_dq1,dx_dq2,dx_dq3
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,T,N,B
  REAL(wp)          :: lp,absB,kappa,tau,sigma
  !===================================================================================================================================
  !  |x |
  !  |y |=  X0(zeta) + sigma*(N(zeta)*q1 + B(zeta)*q2)
  !  |z |
  !  dh/dq1 =sigma*N , dh/dq2=sigma*B
  !  dh/dq3 = l' [(1-sigma*kappa*q1)T + sigma*tau*(B*q1-N*q2) ]
  ASSOCIATE(q1=>q_in(1),q2=>q_in(2),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  T=X0p/lp
  B=CROSS(X0p,X0pp)
  absB=SQRT(SUM(B*B))
  kappa=absB/(lp**3)
  sigma=sf%sigma(zeta)
  tau=SUM(X0ppp*B)/(absB**2)
  B=B/absB
  N=CROSS(B,T)
  dx_dq1(1:3)= sigma*N
  dx_dq2(1:3)= sigma*B
  dx_dq3(1:3)=lp*((1.0_wp-sigma*kappa*q1)*T +sigma*tau*(B*q1-N*q2))
  END ASSOCIATE !zeta
END SUBROUTINE hmap_frenet_get_dx_dqi

!===============================================================================================================================
!> evaluate all first derivatives dx(1:3)/dq^i, i=1,2,3 , at q_in=(X^1,X^2,zeta),
!!
!===============================================================================================================================
SUBROUTINE hmap_frenet_get_dx_dqi_aux( sf ,q1,q2,xv,dx_dq1,dx_dq2,dx_dq3)
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
  REAL(wp),DIMENSION(3),INTENT(OUT) :: dx_dq1,dx_dq2,dx_dq3
  !=================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  dx_dq1(1:3)= xv%sigma*xv%N
  dx_dq2(1:3)= xv%sigma*xv%B
  dx_dq3(1:3)=xv%lp*((1.0_wp-xv%sigma*xv%kappa*q1)*xv%T +xv%sigma*xv%tau*(xv%B*q1-xv%N*q2))
  END SELECT !type(xv)
END SUBROUTINE hmap_frenet_get_dx_dqi_aux

!=================================================================================================================================
!> evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),
!!
!===============================================================================================================================
SUBROUTINE hmap_frenet_get_ddx_dqij( sf ,q_in,ddx_dq11,ddx_dq12,ddx_dq13,ddx_dq22,ddx_dq23,ddx_dq33)
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
  REAL(wp),DIMENSION(3),INTENT(OUT) :: ddx_dq11,ddx_dq12,ddx_dq13,ddx_dq22,ddx_dq23,ddx_dq33
  !-----------------------------------------------------------------------------------------------------------------------------------
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,X0p4,Bloc,T,N,B
  REAL(wp)          :: lp,absB,kappa,tau,sigma
  REAL(wp)          :: lp_p,absB_p,kappa_p,tau_p
!===================================================================================================================================
  !  |x |
  !  |y |=  X0(zeta) + sigma*(N(zeta)*q1 + B(zeta)*q2)
  !  |z |
  ASSOCIATE(q1=>q_in(1),q2=>q_in(2),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp,X0p4=X0p4)
  lp=SQRT(SUM(X0p*X0p))
  T=X0p/lp
  Bloc=CROSS(X0p,X0pp)
  absB=SQRT(SUM(Bloc*Bloc))
  kappa=absB/(lp**3)
  sigma=sf%sigma(zeta)
  tau  = SUM(X0ppp*Bloc)/(absB**2)
  lp_p = SUM(X0pp*X0p) / lp
  absB_p = SUM(Bloc*CROSS(X0p, X0ppp)) / absB
  kappa_p = (absB_p*lp -3*absB * lp_p) / (lp**4)
  tau_p   = (SUM(X0p4*Bloc)*absB -2*SUM(X0ppp*Bloc)*absB_p) / (absB**3)
  B=Bloc/absB
  N=CROSS(B,T)
  ddx_dq11=0.0_wp
  ddx_dq12=0.0_wp
  ddx_dq13=sigma*lp*(-kappa*T +tau*B)

  ddx_dq22=0.0_wp
  ddx_dq23=-sigma*lp*tau*N
  ddx_dq33(1:3)= lp_p*((1.0_wp-sigma*kappa*q1)*T +sigma*tau*(B*q1-N*q2)) &
                 +lp*sigma*( -kappa_p*q1*T +tau_p*(B*q1-N*q2)                     &
                            +lp*( (1.0_wp-sigma*kappa*q1)*(kappa*N)             &
                                 +sigma*tau*((-tau*N)*q1-((-kappa*T +tau*B))*q2)) )
  END ASSOCIATE
END SUBROUTINE hmap_frenet_get_ddx_dqij

!=================================================================================================================================
!> evaluate all second derivatives d^2x(1:3)/(dq^i dq^j), i,j=1,2,3 is evaluated at q_in=(X^1,X^2,zeta),
!!
!===============================================================================================================================
SUBROUTINE hmap_frenet_get_ddx_dqij_aux( sf ,q1,q2,xv,ddx_dq11,ddx_dq12,ddx_dq13,ddx_dq22,ddx_dq23,ddx_dq33)
  IMPLICIT NONE
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
  REAL(wp),DIMENSION(3),INTENT(OUT) :: ddx_dq11,ddx_dq12,ddx_dq13,ddx_dq22,ddx_dq23,ddx_dq33
  !===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  ddx_dq11=0.0_wp
  ddx_dq12=0.0_wp
  ddx_dq13=xv%sigma*xv%lp*(-xv%kappa*xv%T +xv%tau*xv%B)
  ddx_dq22=0.0_wp
  ddx_dq23=-xv%sigma*xv%lp*xv%tau*xv%N
  ddx_dq33(1:3)= xv%lp_p*((1.0_wp-xv%sigma*xv%kappa*q1)*xv%T +xv%sigma*xv%tau*(xv%B*q1-xv%N*q2)) &
                 +xv%lp*xv%sigma*( -xv%kappa_p*q1*xv%T +xv%tau_p*(xv%B*q1-xv%N*q2)                             &
                                  +xv%lp*( (1.0_wp-xv%sigma*xv%kappa*q1)*(xv%kappa*xv%N)                      &
                                          +xv%sigma*xv%tau*((-xv%tau*xv%N)*q1-((-xv%kappa*xv%T +xv%tau*xv%B))*q2)) )
  END SELECT !type(xv)
END SUBROUTINE hmap_frenet_get_ddx_dqij_aux

!===================================================================================================================================
!> evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_Jh( sf ,q_in) RESULT(Jh)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: Jh
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,B
  REAL(wp)          :: lp,absB,kappa,sigma
!===================================================================================================================================
  ASSOCIATE(q1=>q_in(1),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  B=CROSS(X0p,X0pp)
  absB=SQRT(SUM(B*B))
  kappa=absB/(lp**3)
  IF(kappa.LT.1.0e-8) &
      CALL abort(__STAMP__, &
           "hmap_frenet cannot evaluate frame at curvature < 1e-8 !",RealInfo=zeta*sf%nfp/TWOPI)
  sigma=sf%sigma(zeta)

  Jh=lp*(1.0_wp-sigma*kappa*q1)
  IF(Jh .LT. 1.0e-8) &
      CALL abort(__STAMP__, &
           "hmap_frenet, evaluation outside curvature radius, Jh<0",RealInfo=zeta*sf%nfp/TWOPI)

  END ASSOCIATE !zeta
END FUNCTION hmap_frenet_eval_Jh

!===================================================================================================================================
!> evaluate Jacobian of mapping h: J_h=sqrt(det(G)) at q=(q^1,q^2,zeta)
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_Jh_aux( sf ,q1,q2,xv) RESULT(Jh)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: Jh
!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  Jh=xv%lp*(1.0_wp-xv%sigma*xv%kappa*q1)
  END SELECT !type(xv)
END FUNCTION hmap_frenet_eval_Jh_aux

!===================================================================================================================================
!> evaluate derivative of Jacobian of mapping h: sum_k q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_Jh_dq( sf ,q_in,q_vec) RESULT(Jh_dq)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q_in(3)
  REAL(wp)            , INTENT(IN) :: q_vec(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: Jh_dq
!-----------------------------------------------------------------------------------------------------------------------------------
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,Bloc
  REAL(wp)          :: lp,absB,kappa,sigma
  REAL(wp)          :: lp_p,absB_p,kappa_p
!===================================================================================================================================
  !  |x |
  !  |y |=  X0(zeta) + sigma*(N(zeta)*q1 + B(zeta)*q2)
  !  |z |
  ASSOCIATE(q1=>q_in(1),q2=>q_in(2),zeta=>q_in(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  Bloc=CROSS(X0p,X0pp)
  absB=SQRT(SUM(Bloc*Bloc))
  kappa=absB/(lp**3)
  sigma=sf%sigma(zeta)
  lp_p = SUM(X0pp*X0p) / lp
  absB_p = SUM(Bloc*CROSS(X0p, X0ppp)) / absB
  kappa_p = (absB_p*lp -3*absB * lp_p) / (lp**4)
  !tau_p   = (SUM(X0p4*Bloc)*absB -2*SUM(X0ppp*Bloc)*absB_p) / (absB**3)

  Jh_dq=-lp*sigma*kappa*q_vec(1) + (lp_p-sigma*q1*(lp_p*kappa+lp*kappa_p))*q_vec(3) !dsigma/dzeta is a dirac at kappa=0, so it is not evaluated
  END ASSOCIATE !zeta
END FUNCTION hmap_frenet_eval_Jh_dq

!===================================================================================================================================
!> evaluate derivative of Jacobian of mapping h: sum_k q_vec^k * dJ_h/dq^k, k=1,2,3 at q=(q^1,q^2,zeta)
!!
!! NOTE: needs auxvar with do_2nd_der=.TRUE.!! not checked for performance reasons.
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_Jh_dq_aux( sf ,q1,q2,q1_vec,q2_vec,q3_vec,xv) RESULT(Jh_dq)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: q1,q2
  REAL(wp)            , INTENT(IN) :: q1_vec,q2_vec,q3_vec
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: Jh_dq
!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)

  Jh_dq=-xv%lp*xv%sigma*xv%kappa*q1_vec + (xv%lp_p-xv%sigma*q1*(xv%lp_p*xv%kappa+xv%lp*xv%kappa_p))*q3_vec
  END SELECT !type(xv)
END FUNCTION hmap_frenet_eval_Jh_dq_aux


!===================================================================================================================================
!>  evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,,
!! where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and
!! dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_gij( sf ,qL_in,q_G,qR_in) RESULT(g_ab)
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: qL_in(3)
  REAL(wp)            , INTENT(IN) :: q_G(3)
  REAL(wp)            , INTENT(IN) :: qR_in(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                        :: g_ab
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,Bloc
  REAL(wp)              :: lp,absB,kappa,tau,sigma
  REAL(wp)              :: Ga, Gb, Gc
!===================================================================================================================================
  ! A = -q2*l' * tau
  ! B =  q1*l' * tau
  ! C = Jh^2 + (l'*tau)^2(q1^2+q2^2)
  !                       |q1  |   |1   0   Ga|        |q1  |
  !q_i G_ij q_j = (dalpha |q2  | ) |0   1   Gb| (dbeta |q2  | )
  !                       |q3  |   |Ga  Gb  Gc|        |q3  |
  ASSOCIATE(q1=>q_G(1),q2=>q_G(2),zeta=>q_G(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp)
  lp=SQRT(SUM(X0p*X0p))
  Bloc=CROSS(X0p,X0pp)
  absB=SQRT(SUM(Bloc*Bloc))
  kappa=absB/(lp**3)
  sigma=sf%sigma(zeta)
  tau=SUM(X0ppp*Bloc)/(absB**2)

  Ga = -lp*tau*q2
  Gb =  lp*tau*q1
  Gc = (lp**2)*((1.0_wp-sigma*kappa*q1)**2+tau**2*(q1**2+q2**2))
  g_ab=      qL_in(1)*qR_in(1) &
            +qL_in(2)*qR_in(2) &
       + Gc* qL_in(3)*qR_in(3) &
       + Ga*(qL_in(1)*qR_in(3)+qL_in(3)*qR_in(1)) &
       + Gb*(qL_in(2)*qR_in(3)+qL_in(3)*qR_in(2))
  END ASSOCIATE
END FUNCTION hmap_frenet_eval_gij

!===================================================================================================================================
!>  evaluate sum_ij (qL_i (G_ij(q_G)) qR_j) ,,
!! where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and
!! dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_gij_aux( sf ,qL1,qL2,qL3,q1,q2,qR1,qR2,qR3,xv) RESULT(g_ab)
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: qL1,qL2,qL3
  REAL(wp)            , INTENT(IN) :: q1,q2
  REAL(wp)            , INTENT(IN) :: qR1,qR2,qR3
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                        :: g_ab
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp)              :: Ga, Gb, Gc
!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  Ga = -xv%lp*xv%tau*q2
  Gb =  xv%lp*xv%tau*q1
  !Gc = (xv%lp**2)*((1.0_wp-xv%sigma*xv%kappa*q1)**2+xv%tau**2*(q1**2+q2**2))
  Gc = (xv%lp*(1.0_wp-xv%sigma*xv%kappa*q1))**2+ Ga*Ga+ Gb*Gb
  g_ab=      qL1*qR1 &
            +qL2*qR2 &
       + Gc* qL3*qR3 &
       + Ga*(qL1*qR3+qL3*qR1) &
       + Gb*(qL2*qR3+qL3*qR2)
  END SELECT !type(xv)
END FUNCTION hmap_frenet_eval_gij_aux

!===================================================================================================================================
!>  evaluate sum_k=1,3 q_vec^k * sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2
!! where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and
!! where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and
!! dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_gij_dq( sf ,qL_in,q_G,qR_in,q_vec) RESULT(g_ab_dq)
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: qL_in(3)
  REAL(wp)            , INTENT(IN) :: q_G(3)
  REAL(wp)            , INTENT(IN) :: qR_in(3)
  REAL(wp)            , INTENT(IN) :: q_vec(3)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: g_ab_dq
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(3) :: X0,X0p,X0pp,X0ppp,X0p4,Bloc
  REAL(wp)              :: lp,absB,kappa,tau,sigma
  REAL(wp)              :: lp_p,absB_p,kappa_p,tau_p
!===================================================================================================================================
  !                       |q1  |   |0  0        0           |        |q1  |
  !q_i G_ij q_j = (dalpha |q2  | ) |0  0      l'*tau        | (dbeta |q2  | )
  !                       |q3  |   |0  l'*tau  dG33/dq1     |        |q3  |
  ASSOCIATE(q1=>q_G(1),q2=>q_G(2),zeta=>q_G(3))
  CALL sf%eval_X0(zeta,X0,X0p,X0pp,X0ppp,X0p4=X0p4)
  lp=SQRT(SUM(X0p*X0p))
  Bloc=CROSS(X0p,X0pp)
  absB=SQRT(SUM(Bloc*Bloc))
  kappa=absB/(lp**3)
  sigma=sf%sigma(zeta)
  tau=SUM(X0ppp*Bloc)/(absB**2)
  lp_p = SUM(X0pp*X0p) / lp
  absB_p = SUM(Bloc*CROSS(X0p, X0ppp)) / absB
  kappa_p = (absB_p*lp -3*absB * lp_p) / (lp**4)
  tau_p   = (SUM(X0p4*Bloc)*absB -2*SUM(X0ppp*Bloc)*absB_p) / (absB**3)

  !G13 = -lp*tau*q2
  !G23 =  lp*tau*q1
  !G33 = (lp**2)*((1.0_wp-sigma*kappa*q1)**2+tau**2*(q1**2+q2**2))

  g_ab_dq =-(lp*tau*q_vec(2)+(lp_p*tau+lp*tau_p)*q2*q_vec(3))*(qL_in(1)*qR_in(3)+ qL_in(3)*qR_in(1)) &
           +(lp*tau*q_vec(1)+(lp_p*tau+lp*tau_p)*q1*q_vec(3))*(qL_in(2)*qR_in(3)+ qL_in(3)*qR_in(2)) &
           +2.0_wp*(qL_in(3)*qR_in(3))*( (lp**2)*( q_vec(1)*((tau**2+kappa**2)*q1-sigma*kappa)       &
                                                  +q_vec(2)*  tau**2*q2                        )     &
                                        +q_vec(3)*( lp_p*lp* ( (1.0_wp-sigma*kappa*q1)**2+tau**2*(q1**2+q2**2)) &
                                                   +(lp**2)* ( (1.0_wp-sigma*kappa*q1)*(-sigma*kappa_p*q1)      &
                                                              +tau*tau_p*(q1**2+q2**2) ) ) )
  END ASSOCIATE
END FUNCTION hmap_frenet_eval_gij_dq

!===================================================================================================================================
!>  evaluate sum_ij (qL_i d/dq^k(G_ij(q_G)) qR_j) , k=1,2
!! where qL=(dX^1/dalpha,dX^2/dalpha [,dzeta/dalpha]) and qR=(dX^1/dbeta,dX^2/dbeta [,dzeta/dbeta]) and
!! where qL=(dX^1/dalpha,dX^2/dalpha ,dzeta/dalpha) and qR=(dX^1/dbeta,dX^2/dbeta ,dzeta/dbeta) and
!! dzeta_dalpha then known to be either 0 of ds and dtheta and 1 for dzeta
!!
!! NOTE: needs auxvar with do_2nd_der=.TRUE.!! not checked for performance reasons.
!!
!===================================================================================================================================
FUNCTION hmap_frenet_eval_gij_dq_aux( sf ,qL1,qL2,qL3,q1,q2,qR1,qR2,qR3,q1_vec,q2_vec,q3_vec,xv) RESULT(g_ab_dq)
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN) :: sf
  REAL(wp)            , INTENT(IN) :: qL1,qL2,qL3
  REAL(wp)            , INTENT(IN) :: q1,q2
  REAL(wp)            , INTENT(IN) :: qR1,qR2,qR3
  REAL(wp)            , INTENT(IN) :: q1_vec,q2_vec,q3_vec
  CLASS(c_hmap_auxvar), INTENT(IN) :: xv
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                         :: g_ab_dq

!===================================================================================================================================
  SELECT TYPE(xv); TYPE IS(t_hmap_frenet_auxvar)
  g_ab_dq =-(xv%lp*xv%tau*q2_vec+(xv%lp_p*xv%tau+xv%lp*xv%tau_p)*q2*q3_vec)*(qL1*qR3+ qL3*qR1) &
           +(xv%lp*xv%tau*q1_vec+(xv%lp_p*xv%tau+xv%lp*xv%tau_p)*q1*q3_vec)*(qL2*qR3+ qL3*qR2) &
          +2.0_wp*(qL3*qR3)*( (xv%lp**2)*( q1_vec*((xv%tau**2+xv%kappa**2)*q1-xv%sigma*xv%kappa)           &
                                          +q2_vec*  xv%tau**2*q2                        )                  &
                             +q3_vec*( xv%lp_p*xv%lp* ( (1.0_wp-xv%sigma*xv%kappa*q1)**2+xv%tau**2*(q1**2+q2**2))   &
                                      +(xv%lp**2)*    ( (1.0_wp-xv%sigma*xv%kappa*q1)*(-xv%sigma*xv%kappa_p*q1)  &
                                                       +xv%tau*xv%tau_p*(q1**2+q2**2) ) ) )
  END SELECT ! TYPE(xv)
END FUNCTION hmap_frenet_eval_gij_dq_aux


!===================================================================================================================================
!> evaluate curve X0(zeta), position and first three derivatives, from given R0,Z0 Fourier
!!
!===================================================================================================================================
PURE SUBROUTINE hmap_frenet_eval_X0_fromRZ( sf,zeta,X0,X0p,X0pp,X0ppp,X0p4)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(IN ) :: sf
  REAL(wp)            , INTENT(IN ) :: zeta       !! position along closed curve parametrized in [0,2pi]
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)            , INTENT(OUT) :: X0(1:3)      !! curve position in cartesian coordinates
  REAL(wp)            , INTENT(OUT) :: X0p(1:3)     !! 1st derivative in zeta
  REAL(wp)            , INTENT(OUT) :: X0pp(1:3)    !! 2nd derivative in zeta
  REAL(wp)            , INTENT(OUT) :: X0ppp(1:3)   !! 3rd derivative in zeta
  REAL(wp)            , INTENT(OUT),OPTIONAL :: X0p4(1:3)  !! 4th derivative in zeta
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp) :: R0,R0p,R0pp,R0ppp,R0p4,Z0p4
  REAL(wp) :: coszeta,sinzeta
!===================================================================================================================================
  CALL eval_fourier1d(sf%n_max,sf%Xn,sf%rc,sf%rs,zeta,R0,R0p,R0pp,R0ppp,R0p4)
  CALL eval_fourier1d(sf%n_max,sf%Xn,sf%zc,sf%zs,zeta,X0(3),X0p(3),X0pp(3),X0ppp(3),Z0p4) !=Z0,Z0p,Z0pp,Z0ppp
  coszeta=COS(zeta)
  sinzeta=SIN(zeta)
  ASSOCIATE(x   =>X0(1)   ,y   =>X0(2)   , &
            xp  =>X0p(1)  ,yp  =>X0p(2)  , &
            xpp =>X0pp(1) ,ypp =>X0pp(2) , &
            xppp=>X0ppp(1),yppp=>X0ppp(2))
    !! angle zeta=geometric toroidal angle phi=atan(y/x)
    x=R0*coszeta
    y=R0*sinzeta

    xp = R0p*coszeta  - R0*sinzeta
    yp = R0p*sinzeta  + R0*coszeta
    !xp  = R0p*coszeta  -y
    !yp  = R0p*sinzeta  +x

    xpp = R0pp*coszeta - 2*R0p*sinzeta - R0*coszeta
    ypp = R0pp*sinzeta + 2*R0p*coszeta - R0*sinzeta
    !xpp  = R0pp*coszeta -2.0_wp*yp + x
    !ypp  = R0pp*sinzeta +2.0_wp*xp + y

    xppp = R0ppp*coszeta - 3*R0pp*sinzeta - 3*R0p*coszeta + R0*sinzeta
    yppp = R0ppp*sinzeta + 3*R0pp*coszeta - 3*R0p*sinzeta - R0*coszeta
    !xppp  = R0ppp*coszeta +3.0_wp*(xp-ypp) + y
    !yppp  = R0ppp*sinzeta +3.0_wp*(yp+xpp) + x
  IF(PRESENT(X0p4))THEN
    X0p4(1)  = R0p4*coszeta - 4*R0ppp*sinzeta - 6*R0pp*coszeta  + 4*R0p*sinzeta + R0*coszeta
    X0p4(2)  = R0p4*sinzeta + 4*R0ppp*coszeta - 6*R0pp*sinzeta  - 4*R0p*coszeta + R0*sinzeta
    X0p4(3)  = Z0p4
  END IF
  END ASSOCIATE !x,y,xp,yp,...

END SUBROUTINE hmap_frenet_eval_X0_fromRZ


!===================================================================================================================================
!> evaluate 1d fourier series from given cos/sin coefficients and mode numbers xn
!! SUM(xc(0:n_max)*COS(xn(0:n_max)*zeta)+xs(0:n_max)*SIN(xn(0:n_max)*zeta)
!! evaluate all derivatives 1,2,3 alongside
!!
!===================================================================================================================================
PURE SUBROUTINE eval_fourier1d(n_max,xn,xc,xs,zeta,x,xp,xpp,xppp,xp4)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER  , INTENT(IN ) :: n_max        !! number of modes is n_max+1  (0...n_max)
  INTEGER  , INTENT(IN ) :: xn(0:n_max)  !! array of mode numbers
  REAL(wp) , INTENT(IN ) :: xc(0:n_max)  !! cosine coefficients
  REAL(wp) , INTENT(IN ) :: xs(0:n_max)  !!   sine coefficients
  REAL(wp) , INTENT(IN ) :: zeta         !! angular position [0,2pi]
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp) , INTENT(OUT) :: x      !! value at zeta
  REAL(wp) , INTENT(OUT) :: xp     !! 1st derivative in zeta
  REAL(wp) , INTENT(OUT) :: xpp    !! 2nd derivative in zeta
  REAL(wp) , INTENT(OUT) :: xppp   !! 3rd derivative in zeta
  REAL(wp) , INTENT(OUT) :: xp4    !! 4th derivative in zeta
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  REAL(wp),DIMENSION(0:n_max) :: cos_nzeta,sin_nzeta,xtmp,xptmp
!===================================================================================================================================
  cos_nzeta=COS(REAL(xn,wp)*zeta)
  sin_nzeta=SIN(REAL(xn,wp)*zeta)
  xtmp = xc*cos_nzeta+xs*sin_nzeta
  xptmp= REAL(xn,wp)*(-xc*sin_nzeta+xs*cos_nzeta)
  x    = SUM(xtmp)
  xp   = SUM(xptmp)
  xpp  = SUM(REAL(-xn*xn,wp)*xtmp)
  xppp = SUM(REAL(-xn*xn,wp)*xptmp)
  xp4  = SUM(REAL(xn**4,wp)*xtmp)

END SUBROUTINE eval_fourier1d


!===================================================================================================================================
!> test hmap_frenet - evaluation of the map
!!
!===================================================================================================================================
SUBROUTINE hmap_frenet_test( sf )
USE MODgvec_GLobals, ONLY: UNIT_stdOut,testdbg,testlevel,nfailedMsg,nTestCalled,testUnit
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_hmap_frenet), INTENT(INOUT) :: sf  !!self
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER            :: iTest,idir,jdir,kdir,qdir,izeta,i,j,k,ndims(1:3),ijk(3)
  INTEGER,PARAMETER  :: nzeta=5
  INTEGER,PARAMETER  :: ns=2
  INTEGER,PARAMETER  :: nthet=3
  REAL(wp)           :: refreal,checkreal,x(3),q_in(3),q_test(3,3),x_eps(3),dxdq(3),gij,gij_eps,Jh_0,Jh_eps
  REAL(wp),ALLOCATABLE :: zeta(:)
  REAL(wp)           :: qloc(3),q_thet(3),q_zeta(3)
  REAL(wp)           :: dxdq_eps(3),dx_dqi(3,3),ddx_dqij(3,3,3)
  REAL(wp),ALLOCATABLE,DIMENSION(:,:,:) :: q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
                                     Jh,g_tt,    g_tz,    g_zz,     &
                                     Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1, &
                                     Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2, &
                                     g_t1,g_t2,g_z1,g_z2,Gh11,Gh22
  REAL(wp),PARAMETER :: realtol=1.0E-11_wp
  REAL(wp),PARAMETER :: epsFD=1.0e-8
  REAL(wp),PARAMETER :: realtolFD=1.0e-4
  CHARACTER(LEN=10)  :: fail
  REAL(wp)           :: R0, Z0
  TYPE(t_hmap_frenet_auxvar),ALLOCATABLE :: xv(:)
!===================================================================================================================================
  test_called=.TRUE. ! to prevent infinite loop in this routine
  IF(testlevel.LE.0) RETURN
  IF(testdbg) THEN
     Fail=" DEBUG  !!"
  ELSE
     Fail=" FAILED !!"
  END IF
  nTestCalled=nTestCalled+1
  SWRITE(UNIT_stdOut,'(A,I4,A)')'>>>>>>>>> RUN hmap_frenet TEST ID',nTestCalled,'    >>>>>>>>>'
  IF(testlevel.GE.1)THEN

    !evaluate on the axis q1=q2=0
    iTest=101 ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    q_in=(/0.0_wp, 0.0_wp, 0.335_wp*PI/)
    R0 = SUM(sf%rc(:)*COS(sf%Xn(:)*q_in(3)) + sf%rs(:)*SIN(sf%Xn(:)*q_in(3)))
    Z0 = SUM(sf%zc(:)*COS(sf%Xn(:)*q_in(3)) + sf%zs(:)*SIN(sf%Xn(:)*q_in(3)))
    x = sf%eval(q_in )
    checkreal=SUM((x-(/R0*COS(q_in(3)),R0*SIN(q_in(3)),Z0/))**2)
    refreal = 0.0_wp

    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
       nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
            '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
       nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3))') &
     '\n =>  should be ', refreal,' : |y-eval_map(x)|^2= ', checkreal
    END IF !TEST

    q_test(1,:)=(/1.0_wp, 0.0_wp, 0.0_wp/)
    q_test(2,:)=(/0.0_wp, 1.0_wp, 0.0_wp/)
    q_test(3,:)=(/0.0_wp, 0.0_wp, 1.0_wp/)
    DO qdir=1,3
      !check dx/dq^i with FD
      iTest=101+qdir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
      q_in=(/0.0_wp, 0.0_wp, 0.335_wp*PI/)
      x = sf%eval(q_in )
      x_eps = sf%eval(q_in+epsFD*q_test(qdir,:))
      dxdq = sf%eval_dxdq(q_in,q_test(qdir,:))
      checkreal=SQRT(SUM((dxdq - (x_eps-x)/epsFD)**2)/SUM(x*x))
      refreal = 0.0_wp

      IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtolFD))) THEN
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
              '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),(A,I3))') &
       '\n =>  should be <',realtolFD,' : |dxdqFD-eval_dxdq|= ', checkreal,", dq=",qdir
      END IF !TEST
    END DO

    !! TEST G_ij
    DO idir=1,3; DO jdir=idir,3
      iTest=iTest+1 ; IF(testdbg)WRITE(*,*)'iTest=',iTest
      checkreal= SUM(sf%eval_dxdq(q_in,q_test(idir,:))*sf%eval_dxdq(q_in,q_test(jdir,:)))
      refreal  =sf%eval_gij(q_test(idir,:),q_in,q_test(jdir,:))
      IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
              '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),2(A,I3))') &
       '\n =>  should be ', refreal,' : sum|Gij-eval_gij|= ', checkreal,', i=',idir,', j=',jdir
      END IF !TEST
    END DO; END DO
    !! TEST dJh_dqk with FD
    DO qdir=1,3
      iTest=iTest+1; IF(testdbg)WRITE(*,*)'iTest=',iTest
      Jh_0    = sf%eval_Jh(   q_in                     )
      Jh_eps  = sf%eval_Jh(   q_in+epsFD*q_test(qdir,:))
      refreal = sf%eval_Jh_dq(q_in                     ,q_test(qdir,:))
      checkreal=(Jh_eps-Jh_0)/epsFD-refreal
      refreal=0.0_wp
      IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtolFD))) THEN
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
              '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),(A,I3))') &
       '\n =>  should be <', realtolFD,' : |dJh_dqFD-eval_Jh_dq|= ', checkreal,', dq=',qdir
      END IF !TEST
    END DO !qdir
    !! TEST dG_ij_dqk with FD
    DO qdir=1,3
    DO idir=1,3; DO jdir=1,3
      iTest=iTest+1; IF(testdbg)WRITE(*,*)'iTest=',iTest
      gij     = sf%eval_gij(   q_test(idir,:),q_in                     ,q_test(jdir,:))
      gij_eps = sf%eval_gij(   q_test(idir,:),q_in+epsFD*q_test(qdir,:),q_test(jdir,:))
      refreal = sf%eval_gij_dq(q_test(idir,:),q_in                     ,q_test(jdir,:),q_test(qdir,:))
      checkreal=(gij_eps-gij)/epsFD-refreal
      refreal=0.0_wp
      IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtolFD))) THEN
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
              '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),3(A,I3))') &
       '\n =>  should be <', realtolFD,' : |dGij_dqFD-eval_gij_dq|= ', checkreal,', i=',idir,', j=',jdir,', dq=',qdir
      END IF !TEST
    END DO; END DO
    END DO !qdir

    CALL sf%get_dx_dqi(q_in,dx_dqi(1,:),dx_dqi(2,:),dx_dqi(3,:))
    DO qdir=1,3
      !check dx/dq^i with FD
      iTest=10+qdir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
      dxdq = sf%eval_dxdq(q_in,q_test(qdir,:))
      checkreal=SUM(ABS(dxdq-dx_dqi(qdir,:)))
      refreal=0.0_wp
      IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
              '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
         nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),(A,I3))') &
       '\n =>  should be ', refreal,' : |dxdq-eval_dxdq|= ', checkreal,", dq=",qdir
      END IF !TEST
    END DO !qdir
    CALL sf%get_ddx_dqij(q_in,ddx_dqij(1,1,:),ddx_dqij(1,2,:),ddx_dqij(1,3,:), &
                              ddx_dqij(2,2,:),ddx_dqij(2,3,:),ddx_dqij(3,3,:))
    ddx_dqij(2,1,:)=ddx_dqij(1,2,:)
    ddx_dqij(3,1,:)=ddx_dqij(1,3,:)
    ddx_dqij(3,2,:)=ddx_dqij(2,3,:)
    DO qdir=1,3
      dxdq = sf%eval_dxdq(q_in,q_test(qdir,:))
      DO idir=1,3
        iTest=20+idir+3*(qdir-1) ; IF(testdbg)WRITE(*,*)'iTest=',iTest
        dxdq_eps = sf%eval_dxdq(q_in+epsFD*q_test(idir,:),q_test(qdir,:))
        checkreal=SUM(ABS( (dxdq_eps-dxdq)/epsFD-ddx_dqij(qdir,idir,:)))
        refreal=0.0_wp
        IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtolFD))) THEN
           nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
                '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
           nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),2(A,I3))') &
         '\n =>  should be <', realtolFD,' : |ddx_dqijFD-eval_ddx_dqij|= ', checkreal,", dqi=",qdir,", dqj=",idir
        END IF !TEST
      END DO !idir
    END DO !qdir

 END IF !testlevel >=1
 IF (testlevel .GE. 2) THEN
  DO idir=1,3
    SELECT CASE(idir)
    CASE(1)
      jdir=2; kdir=3
    CASE(2)
      jdir=1; kdir=3
    CASE(3)
      jdir=1; kdir=2
    END SELECT
    ndims(idir)=nzeta+idir
    ndims(jdir)=ns
    ndims(kdir)=nthet
    ALLOCATE(zeta(ndims(idir)),xv(ndims(idir)))
    DO izeta=1,ndims(idir)
      zeta(izeta)=0.333_wp+REAL(izeta-1,wp)/REAL(ndims(idir)-1,wp)*0.221_wp
      xv(izeta)=hmap_frenet_init_aux(sf,zeta(izeta),.TRUE.)
    END DO
    ALLOCATE(q1(ndims(1),ndims(2),ndims(3)))
    ALLOCATE(q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz,Jh,g_tt,g_tz,g_zz,Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1,Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2,g_t1,g_t2,g_z1,g_z2,Gh11,Gh22, &
             mold=q1)
    !assign somewhat randomly
    DO k=1,ndims(3); DO j=1,ndims(2); DO i=1,ndims(1)
      q1(i,j,k) = 0.11_wp -0.22_wp *REAL((i+j)*k,wp)/REAL((ndims(idir)+ndims(jdir))*ndims(kdir),wp)
      q2(i,j,k) = 0.15_wp -0.231_wp*REAL((i+k)*j,wp)/REAL((ndims(idir)+ndims(kdir))*ndims(jdir),wp)
      dX1_dt(i,j,k)=-0.1_wp  +0.211_wp*REAL((i+2*j)*k,wp)/REAL((ndims(idir)+2*ndims(jdir))*ndims(kdir),wp)
      dX2_dt(i,j,k)= 0.231_wp-0.116_wp*REAL((2*i+k)*j,wp)/REAL((2*ndims(idir)+ndims(kdir))*ndims(jdir),wp)
      dX1_dz(i,j,k)=-0.024_wp+0.013_wp*REAL((3*i+2*j)*k,wp)/REAL((3*ndims(idir)+2*ndims(jdir))*ndims(kdir),wp)
      dX2_dz(i,j,k)=-0.06_wp +0.031_wp*REAL((2*k+3*k)*i,wp)/REAL((2*ndims(kdir)+3*ndims(kdir))*ndims(idir),wp)
    END DO; END DO; END DO
    CALL sf%eval_all(ndims,idir,xv,q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
         Jh,g_tt,g_tz,g_zz,&
         Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1,&
         Jh_dq2,g_tt_dq2,g_tz_dq2,g_zz_dq2,&
         g_t1,g_t2,g_z1,g_z2,Gh11,Gh22)
    DO k=1,ndims(3); DO j=1,ndims(2); DO i=1,ndims(1)
      ijk=(/i,j,k/)
      izeta=ijk(idir)
      qloc=(/q1(i,j,k),q2(i,j,k),zeta(izeta)/)
      q_thet=(/dX1_dt(i,j,k),dX2_dt(i,j,k),0.0_wp/)
      q_zeta=(/dX1_dz(i,j,k),dX2_dz(i,j,k),1.0_wp/)
      Jh(i,j,k)       =Jh(i,j,k)       - sf%eval_Jh(qloc)
      g_tt(i,j,k)     =g_tt(i,j,k)     - sf%eval_gij(q_thet,qloc,q_thet)
      g_tz(i,j,k)     =g_tz(i,j,k)     - sf%eval_gij(q_thet,qloc,q_zeta)
      g_zz(i,j,k)     =g_zz(i,j,k)     - sf%eval_gij(q_zeta,qloc,q_zeta)
      Jh_dq1(i,j,k)   =Jh_dq1(i,j,k)   - sf%eval_Jh_dq(qloc,(/1.0_wp,0.0_wp,0.0_wp/))
      Jh_dq2(i,j,k)   =Jh_dq2(i,j,k)   - sf%eval_Jh_dq(qloc,(/0.0_wp,1.0_wp,0.0_wp/))
      g_tt_dq1(i,j,k) =g_tt_dq1(i,j,k) - sf%eval_gij_dq(q_thet,qloc,q_thet,(/1.0_wp,0.0_wp,0.0_wp/))
      g_tt_dq2(i,j,k) =g_tt_dq2(i,j,k) - sf%eval_gij_dq(q_thet,qloc,q_thet,(/0.0_wp,1.0_wp,0.0_wp/))
      g_tz_dq1(i,j,k) =g_tz_dq1(i,j,k) - sf%eval_gij_dq(q_thet,qloc,q_zeta,(/1.0_wp,0.0_wp,0.0_wp/))
      g_tz_dq2(i,j,k) =g_tz_dq2(i,j,k) - sf%eval_gij_dq(q_thet,qloc,q_zeta,(/0.0_wp,1.0_wp,0.0_wp/))
      g_zz_dq1(i,j,k) =g_zz_dq1(i,j,k) - sf%eval_gij_dq(q_zeta,qloc,q_zeta,(/1.0_wp,0.0_wp,0.0_wp/))
      g_zz_dq2(i,j,k) =g_zz_dq2(i,j,k) - sf%eval_gij_dq(q_zeta,qloc,q_zeta,(/0.0_wp,1.0_wp,0.0_wp/))
      g_t1(i,j,k)     =g_t1(i,j,k)     - sf%eval_gij(q_thet,qloc,(/1.0_wp,0.0_wp,0.0_wp/))
      g_t2(i,j,k)     =g_t2(i,j,k)     - sf%eval_gij(q_thet,qloc,(/0.0_wp,1.0_wp,0.0_wp/))
      g_z1(i,j,k)     =g_z1(i,j,k)     - sf%eval_gij(q_zeta,qloc,(/1.0_wp,0.0_wp,0.0_wp/))
      g_z2(i,j,k)     =g_z2(i,j,k)     - sf%eval_gij(q_zeta,qloc,(/0.0_wp,1.0_wp,0.0_wp/))
      Gh11(i,j,k)     =Gh11(i,j,k)     - sf%eval_gij((/1.0_wp,0.0_wp,0.0_wp/),qloc,(/1.0_wp,0.0_wp,0.0_wp/))
      Gh22(i,j,k)     =Gh22(i,j,k)     - sf%eval_gij((/0.0_wp,1.0_wp,0.0_wp/),qloc,(/0.0_wp,1.0_wp,0.0_wp/))
    END DO; END DO; END DO

    iTest=201+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(Jh))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|Jh_all-eval_Jh(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=202+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tt))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tt_all-eval_g_tt(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=203+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tz))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tz_all-eval_g_tz(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=203+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_zz))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_zz_all-eval_g_zz(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=204+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(Jh_dq1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|Jh_dq1_all-eval_Jh_dq1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=205+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(Jh_dq2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|Jh_dq2_all-eval_Jh_dq2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=206+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tt_dq1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tt_dq1_all-eval_g_tt_dq1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=207+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tz_dq1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tz_dq1_all-eval_g_tz_dq1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=208+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_zz_dq1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_zz_dq1_all-eval_g_zz_dq1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=209+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tt_dq2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tt_dq2_all-eval_g_tt_dq2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=210+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_tz_dq2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_tz_dq2_all-eval_g_tz_dq2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=211+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_zz_dq2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_zz_dq2_all-eval_g_zz_dq2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=212+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_t1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_t1_all-eval_g_t1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=213+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_t2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_t2_all-eval_g_t2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=214+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_z1))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_z1_all-eval_g_z1(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=215+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(g_z2))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|g_z2_all-eval_g_z2(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=216+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(Gh11))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|Gh11_all-eval_Gh11(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    iTest=217+20*idir ; IF(testdbg)WRITE(*,*)'iTest=',iTest
    checkreal=SUM(ABS(Gh22))/REAL(PRODUCT(ndims),wp)
    refreal=0.0_wp
    IF(testdbg.OR.(.NOT.( ABS(checkreal-refreal).LT. realtol))) THEN
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(A,2(I4,A))') &
           '\n!! hmap_frenet TEST ID',nTestCalled ,': TEST ',iTest,Fail
      nfailedMsg=nfailedMsg+1 ; WRITE(testUnit,'(2(A,E11.3),A,I4)') &
    '\n =>  should be ', refreal,' : |sum(|Gh22_all-eval_Gh22(xall)|)|= ', checkreal, " ,idir=",idir
    END IF

    DEALLOCATE(zeta,q1,q2,dX1_dt,dX2_dt,dX1_dz,dX2_dz, &
               Jh,g_tt,g_tz,g_zz,Jh_dq1,g_tt_dq1,g_tz_dq1,g_zz_dq1,Jh_dq2,&
               g_tt_dq2,g_tz_dq2,g_zz_dq2,g_t1,g_t2,g_z1,g_z2,Gh11,Gh22)
    DEALLOCATE(xv)
  END DO !idir
END IF

 test_called=.FALSE. ! to prevent infinite loop in this routine


END SUBROUTINE hmap_frenet_test

END MODULE MODgvec_hmap_frenet
