!> parameters for the ama hp adaptation
module ama_hp_DWR_params
  real :: dwr_ama_err_max
  real :: dwr_ama_err_min
  real :: dwr_ama_err_zero_level
  real :: dwr_ama_err_aver
  real :: dwr_ama_err_total
  real :: dwr_ama_err_tol_K
  real :: dwr_ama_err_total_reduced
  real :: dwr_ama_err_total_old
  real :: dwr_ama_target_tol
  integer :: dwr_ama_ncount
  real :: dwr_ama_min_area
  real, parameter, private :: dwr_ama_r_max = 0.1  ! the maximal decrease of area
  real, parameter, private :: dwr_ama_c_max = 2.5  ! the maximal increase of area
  !real, parameter, private :: dwr_ama_r_max = 0.5  ! the maximal decrease of area
  !real, parameter, private :: dwr_ama_c_max = 1.5  ! the maximal increase of area
  integer, dimension(:), allocatable :: dwr_ama_iest
  real, dimension(:), allocatable :: dwr_ama_est

  real :: dwr_ama_equi_hp_sum

contains
  function Get_dwr_ama_r_max( ) result (r_max)
    r_max = dwr_ama_r_max
  end function Get_dwr_ama_r_max

  function Get_dwr_ama_c_max( ) result (c_max)
    c_max = dwr_ama_c_max
  end function Get_dwr_ama_c_max

    !> refinement or recoarsening factor by [Balan, Woopwn, May: CaF 2016]
  function RefineFactor(eta) result (ratio)
    real, intent(in) :: eta
    real :: eta_c, r_max, c_max, eta_min, eta_max
    real :: ratio, ratioRWTH, ratioMFF
    real :: xi, pi

    pi = 2 * asin(1.)

    eta_max = dwr_ama_err_max
    eta_min = dwr_ama_err_min
    eta_c   = dwr_ama_err_tol_K
    r_max = dwr_ama_r_max
    c_max =  dwr_ama_c_max

    eta_min = max(eta_min, 1E-10 * eta_c)

    if(eta <= dwr_ama_err_zero_level) then  ! small error estimate
       ratio = c_max

    elseif(eta >= eta_max) then  ! too large error estimate
       ratio = r_max

    elseif(eta >=eta_c) then
       xi = (log(eta) - log(eta_c) ) / (log(eta_max) - log(eta_c) )
       !ratio = 1./((r_max - 1) * xi* xi + 1)

       ratioRWTH = ((r_max - 1) * xi* xi + 1)   ! RWTH
       ratioMFF = ( 1.- r_max)/2. * ( cos(pi * xi) + 1) + r_max

       !ratio = ratioRWTH
       ratio = ratioMFF
       !print*, '#E#E#:', xi, ratio
    else

       xi = (log(eta) - log(eta_c) ) / (log(eta_min) - log(eta_c) )

       ratioRWTH = (c_max - 1) * xi* xi + 1     ! RWTH
       ratioMFF = (c_max -1) / 2. * ( cos(pi *(xi + 1)) + 1) + 1

       !ratio = ratioRWTH
       ratio = ratioMFF
    end if

    !write(*, *) eta, ratioRWTH, ratioMFF, eta_c

  end function RefineFactor


end module ama_hp_DWR_params


!> subroutines for hp method based on DWR error estimates
module ama_hp_DWR

  use main_data  ! contains type(mesh) :: grid for computation
  use problem_oper
  use euler_problem
  use estimates
  use lapack_oper
  use plot_geom
  use eval_sol
  use geometry
  use marking
  use ama_L2interpol
  use regularity
  use pm_fluxes
  use ama_hp_interpol
  use ama_hp_DWR_params

  implicit none

  public:: ComputeDWRAnisotropicMetric
  public:: ComputeDWRAnisotropicMetric_Elem
  public:: Set_DWR_element_size
  public:: Set_DWR_metric
  public:: Set_DWR_anisotropy
  public:: SeekMinimumIteratively
  public:: SeekMinimumGlobally
  public:: FindAnisotropy
  public:: Draw_deriv
  public:: Draw_anis
  public:: Eval_Dir_Derivative
  public:: Eval_products
  public:: Set_Coeffs_Function
  public:: Set_Coeffs_Gradient
  public:: Find_Params_init
  public:: Anis_estim_tot
  public:: Anis_estim
  public:: IntegrateCircle
  public:: Eval_High_Order_Deriv_Elem

contains

  !> evaluate the new metric based on the anisotropic DWR error estimates
  subroutine ComputeDWRAnisotropicMetric( DWR, grid )
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    class( element ), pointer :: elem
    real, dimension(:, :), allocatable :: anisotropy
    real :: ratio_min
    integer :: i
    integer :: nelem

    call state%cpuTime%startAdaptTime()

    !!!call TEST_ANIS_ESTIM( )

    ! evaluation of the parameters for the ama_hp_adaptation
    call Eval_AMA_paramets( )

    !print*,'####', state%space%adapt%adapt_space
    !stop

    ! setting of the initial tolerance
    if( state%space%adapt%adapt_level == 0) then
         state%space%adapt%tol_init = sqrt(  state%estim( dwr_aver, 1) )
         !WS
         state%space%adapt%tol_min = state%space%adapt%tol_init * 1E-1 
         !state%space%adapt%tol_min = state%space%adapt%tol_init * 1E-2
         print*, "FR_ANI: what should be in state%space%adapt%tol_min?"
         !state%space%adapt%tol_min = state%space%adapt%tol_init * 1E-3

    endif


    !the target tolerance for adaptation
    !state%space%adapt%tol_min = 0.5 * sqrt( abs( state%estim( DWR%eta_index, 1 ) ))
    !!state%space%adapt%tol_min = sqrt( abs( state%estim( DWR%eta_index, 1 ) ))
    !state%space%adapt%tol_min = 1E-2 * 0.5**(state%space%adapt%adapt_level)
    !state%space%adapt%tol_min = max(state%space%adapt%tol_min, 0.2 * state%space%adapt%tol_max)

    ratio_min = 0.1
    if(state%space%adapt%tol_min > state%space%adapt%tol_max) then
       !WS
       state%space%adapt%tol_min = state%space%adapt%tol_min * 0.25
       !state%space%adapt%tol_min = state%space%adapt%tol_min * 0.5

    else
       !WS
       state%space%adapt%tol_min=  state%space%adapt%tol_min * 0.25 ! 0.75
       !state%space%adapt%tol_min=  state%space%adapt%tol_min * 0.75
       !state%space%adapt%tol_min=  state%space%adapt%tol_min * 0.5
       !state%space%adapt%tol_min=  state%space%adapt%tol_max * 0.1
    endif
    ! miniml possible tolerance
    !state%space%adapt%tol_min=  max( state%space%adapt%tol_min,  state%space%adapt%tol_max * ratio_min)


    !print*,'state%space%adapt%tol_min = ', state%space%adapt%tol_min
    
    nelem = grid%nelem

    ! array for the element ansotropy
    allocate( anisotropy(1:nelem, 1:5), source = 0.0 )

    !> computing of parameters for DWR_AMA adaptivity (minimal and maximal value)
    call ComputeDWR_adapt_params(  grid )



    ! computations of the metric for each elements
    do i=1, nelem
       elem => grid%elem(i)

       call ComputeDWRAnisotropicMetric_Elem( elem, anisotropy(i, 1:5) )
       !print*, "anis::", i, anisotropy(i, 1:5)
    enddo

!    print*, "anisotropy::", norm2(anisotropy(:, 1)), &
!    	norm2(anisotropy(:, 2)), norm2(anisotropy(:, 3))
    !call Elem_sizes( )



    ! output of the error estimates in file DWR_AMA_estim
    call DWR%writeDWRerrorsFile(nelem, state%nsize )


    !call Draw_metric( 'metricA' , state%space%adapt%adapt_level)

    ! smoothing of the metric
    !call SmoothMetric( )

    ! not not USE the folowing, metric not yet associated!!
    !call Draw_metric( 'metricS' , state%space%adapt%adapt_level)


    ! direct smoothing of the anisotropy, i.e. lambda and sigma, SETTING of elem%rgabc
    call SmoothAnisotropy(grid, nelem, anisotropy(1:nelem, 1:4) )
!    print*, "anisotropy after smooth::", norm2(anisotropy(:, 1)), &
!    	norm2(anisotropy(:, 2)), norm2(anisotropy(:, 3))

    !call Draw_metric( 'metricX' , state%space%adapt%adapt_level)

    ! stop "in subroutine computeDWRAnisotropicMetric (2) "

    ! smoothing of the metric
    !call SmoothMetric( )
    !call Draw_metric( 'metricS' , state%space%adapt%adapt_level)

    call state%cpuTime%addAdaptTime( )

    deallocate(anisotropy)

  end subroutine ComputeDWRAnisotropicMetric

  !> evaluate the parameters for DWR_AMA adaptation
  subroutine ComputeDWR_adapt_params(grid )
    class( mesh ), intent(inout) :: grid
    class(element), pointer :: elem
    integer :: nelem  ! number of mesh elements
    integer :: i, i1, i2, i3, i4, imax
    real , dimension(:), allocatable :: est
    integer , dimension(:), allocatable :: iest

    nelem = grid%nelem

    dwr_ama_err_max = maxval( grid%elem(:)%estim_loc )
    dwr_ama_err_min = minval( grid%elem(:)%estim_loc )
    dwr_ama_err_total = sum(grid%elem(:)%estim_loc )

    dwr_ama_err_tol_K =  state%space%adapt%tol_min / nelem

    dwr_ama_err_zero_level =  dwr_ama_err_tol_K* 1E-07 !-15

    ! ordering of the elements according the values of estimator
    allocate(iest(1:grid%nelem), est(1:grid%nelem) )
    do i=1,grid%nelem
       iest(i) = i
       est(i) = grid%elem(i)%estim_loc
    enddo

    call order_estims(grid%nelem, iest, est)

    imax = int(0.02 * nelem) + 1  ! max value corresponds to the n% of the elements

    !write(*,'(a10, 2i5, 30es14.6)') 'maxim', imax, iest(1),dwr_ama_err_max, est(1), grid%elem(iest(1))%estim_loc
    !write(*,'(a10, 2i5, 30es14.6)') 'maxim', nelem,iest(nelem),dwr_ama_err_min, est(nelem), grid%elem(iest(nelem))%estim_loc
    !write(*,'(a10, 2i5, 30es14.6)') 'maxim', imax, iest(imax), est(imax), grid%elem(iest(imax))%estim_loc

    !do i=1, imax
    !   print*,'d3d3d3',i,iest(i), est(i), dwr_ama_err_max
    !enddo


    !stop

    !dwr_ama_err_max =  sum(est(1:imax)) / imax

    deallocate(iest, est)

    !    ! optimal error for one element
    !    dwr_ama_err_tol_K = state%space%adapt%tol_min / nelem
    ! else
    !    dwr_ama_err_max = max(dwr_ama_err_max,  elem%estim_loc)
    !    dwr_ama_err_min = min(dwr_ama_err_min, elem%estim_loc )
    !    dwr_ama_err_total = dwr_ama_err_total + elem%estim_loc
    ! endif




    ! ! last element
    ! if(elem%i == nelem) then
    !print*,'estim(dwrEtaI_aver =', sqrt(state%estim(dwrEtaI_aver,1:ndim)), DWR%eta_index, dwrEtaI_aver

    !print*,'dwr_ama_err_total = ', dwr_ama_err_total
    !print*,'dwr_ama_err_max   = ', dwr_ama_err_max
    !print*,'dwr_ama_err_min   = ', dwr_ama_err_min
    !print*,'dwr_ama_err_tol_K = ', dwr_ama_err_tol_K
    !print*,'dwr_ama_err_zero  = ', dwr_ama_err_zero_level
    !if( dwr_ama_err_min < dwr_ama_err_zero_level ) then
    !   dwr_ama_err_min = dwr_ama_err_zero_level
    !   print*, "ZERO LEVEL ACHIEVED"
    !endif

    i1 = 0. ; i2 = 0. ; i3 = 0. ; i4 = 0.

    ! do i=1, nelem
    !    elem => grid%elem(i)
    !    if(elem%estim_loc < dwr_ama_err_min) then
    !       i1 = i1 + 1
    !    elseif(elem%estim_loc < dwr_ama_err_tol_K) then
    !       i2 = i2 + 1
    !    elseif(elem%estim_loc < dwr_ama_err_max) then
    !       i3 = i3 + 1
    !    else
    !       i4 = i4 + 1
    !    endif

    ! enddo
    ! !print*,'dwr_ama_err COUNT: ', i1, i2, i3, i4

    ! endif

  end subroutine ComputeDWR_Adapt_Params

  subroutine Elem_sizes( )
    real :: eta_c, eta, r_max, c_max, ratio, eta_min, eta_max, q
    integer :: i, n

    !c_max = 2.5
    !r_max = 0.5

    eta_c = 1E-5  !1E-7

    eta = 1E-10 !1E-12 !1E-10
    n = 80  !110 !80

    q = 10.**(1./8)


    eta_min = eta
    eta_max = eta_min* q**n

    dwr_ama_err_min = eta_min
    dwr_ama_err_max = eta_max
    dwr_ama_err_zero_level = eta_min
    dwr_ama_err_tol_K = eta_c
    do i= 0, n
       ratio = RefineFactor(eta)
       write(91, *) eta, ratio

       eta = eta * q
    enddo

    stop "y3h3i3dddee"

  end subroutine Elem_sizes


  !> evaluate the new metric based on the anisotropic DWR error estimates for one element
  subroutine ComputeDWRAnisotropicMetric_Elem( elem, anisotropy )
    class(element), target, intent(inout) :: elem
    real, dimension(1:5), intent(inout) :: anisotropy  ! anisotropy of the element
    !real, dimension(,:,:), allocatable :: Deriv  ! scaled directional derivative
    real :: lambda, sigma, phi, lambda_K, sigma_K, phi_K, eta
    integer :: i, j, k, degP

    !write(*,'(a10, 2i5, es12.4, a2, 30es12.4)') &
    !     ' EST:', elem%i, elem%dof, elem%eta( dwr_aver  ,:),'|', &
    !     elem%eta(dwrEtaKV:dwrEtaKD, :),  elem%eta(dwrEtaKV_dual:dwrEtaKD_dual, :)

    !print*
    !do j=0, 2
    !   write(*,'(a10, 2i5, 30es12.4)') ' primal:', elem%i, elem%dof, elem%wST_LS( 1, 1:, j )
    !enddo
    !print*
    !do j=0, 2
    !   write(*,'(a10, 2i5, 30es12.4)') ' dual:', elem%i, elem%dof, elem%zST_LS( 1, 1:, j )
    !enddo


    ! evaluation of the derivatives
    call Eval_High_Order_Deriv_Elem(elem)

    ! setting of the optimal element size from the error estimate and the tolerance
    call Set_DWR_element_size(elem, lambda, lambda_K, sigma_K, phi_K)

    ! setting of the anisotropy
    if(state%space%adapt%adapt_space == 'AMAh') then
       ! only h-variant
       j = 1
       degP = j ! polynomial degree: p_K + degP
       !do k=1,ndim
       !   write(*,'(a10, 2i5, 30es12.4)') 'Der primal:', elem%i, elem%dof, elem%wSS(k, j, : )
       !write(*,'(a10, 2i5, 30es12.4)') 'Der dual  :', elem%i, elem%dof, elem%wSS(ndim+k, j, : )
       !enddo
       

       ! setting of the optimal anisotropy
       call Set_DWR_anisotropy(elem, j, lambda_K, sigma_K, phi_K, lambda, sigma, phi, eta)

    else if(state%space%adapt%adapt_space == 'AMAhp') then
       ! full hp-variant

       ! setting of the optimal anisotropy
       call Set_DWR_hp_anisotropy(elem, lambda_K, sigma_K, phi_K, lambda, sigma, phi, degP)
    endif

    ! isotropic variant
    !sigma = 1.
    !phi = 0.

    anisotropy(1) = lambda
    anisotropy(2) = sigma
    anisotropy(3) = phi
    anisotropy(4) = degP

    ! NOT used at this place, it is called in SmoothAnisotropy
    ! setting of the metric
    !!call Set_DWR_metric(elem, j, lambda, sigma, phi)

    !write(*,'(a8, i5, 3es12.4, a2, 3es12.4)') 'metric', elem%i, lambda, sigma, phi, '|', elem%rgabc(1:3)

  end subroutine ComputeDWRAnisotropicMetric_Elem


  !> setting of the optimal element size from the error estimate and the tolerance
  subroutine Set_DWR_element_size(elem, lambda, lambda_K, sigma_K, phi_K)
    use  ama_hp_DWR_params
    class(element), target, intent(inout) :: elem
   !! integer, intent(inout) :: ideg  !? anisotropy driven by polynomial function of degree elem%deg+ideg
    real, intent(inout)  :: lambda  ! computed new size of the element
    real, intent(inout)  :: lambda_K, sigma_K, phi_K  ! computed original size and shape
    real, dimension(:,:), allocatable :: F_K, A, U, VT
    real, dimension(:), allocatable :: S, xi, Fxi
    real :: rmultiple, tol_K, tol_max, ratio, ratio_OLD
    real :: ratio_limit_min   ! maximal decrease of the element area
    real :: ratio_limit_max   ! maximal increase of the element area
    real :: lambda_max   ! maximal possible element size (corresponds to the domain ares)
    real :: lambda_min   ! minimal possible element size
    real :: pi, scale
    integer :: i, j, N, nn, ifile

    ! setting from module ama_hp_DWR_params
    ratio_limit_min = Get_dwr_ama_r_max( )   ! 0.5
    ratio_limit_max = Get_dwr_ama_c_max( )   ! 2.5

    if(state%space%adapt%adapt_space == 'AMAhp') then
       !ratio_limit_max = 5.
    endif

    !tol_max = state%space%adapt%tol_max
    tol_max = state%space%adapt%tol_min  !! * 0.9 !* 0.1 !0.5 !25

    !tolerance for one element
    tol_K = tol_max / grid%nelem  ! NOT used IN THE NEW VERSION: function RefineFactor

    !tolerance for one element - hp variant
    !tol_K = tol_max  * (elem%deg + 2) /  equi_hp_sum

    !if(elem%i < 5) print*, 'tol_K =;;', tol_K

    ! scaling factor = the  area of the equilateral reference triangle
    scale = 0.75 * sqrt(3.)

    ! mapping from hat{K} to K
    allocate(F_K(1:2, 1:3), xi(1:2), Fxi(1:2) )
    F_K(1:2, 1) = (grid%x(elem%face(idx, 3), 1:2) - grid%x(elem%face(idx,2), 1:2))/sqrt(3.)
    F_K(1:2, 2) = grid%x(elem%face( idx, 1) ,1:2) - elem%xc(1:2)
    F_K(1:2, 3) = elem%xc(1:2)

    pi =  asin(1.0)*2.0
    N = 100

    !do i=0, N
    !   xi(1) = cos(1.*i/N * 2*pi)
    !   xi(2) = sin(1.*i/N * 2*pi)
    !
    !   write(65,*) F_K(1, 1) * xi(1) + F_K(1, 2) * xi(2)  + F_K(1, 3), &
    !        F_K(2, 1) * xi(1) + F_K(2, 2) * xi(2)  + F_K(2, 3)
    !enddo
    !write(65,'(x)')


    ! SVD decomposition
    nn = 2
    allocate( A(1:nn, 1:nn), S(1:nn), U(1:nn, 1:nn), VT(1:nn, 1:nn) )
    A(1:2, 1:2) =  F_K(1:2, 1:2)


    !do i=1,nn
    !   write(*,'(a6, i5, 30es12.4)') 'FK:', i, F_K(i, 1:3)
    !enddo

    !print*,'___________________', elem%i, elem%xc(:), &
    !     ',  det = ', A(1,1)*A(2,2) - A(1,2)*A(2,1)
    !do i=1,nn
    !   write(*,'(a6, i5, 30es12.4)') 'A:', i, A(i, 1:nn)
    !enddo

    call SVD_matrix(nn, A(1:nn, 1:nn), S(1:nn), U(1:nn, 1:nn), VT(1:nn, 1:nn) )

    !print*
    !do i=1,nn
    !   write(*,'(a6, i5, 30es12.4)') 'U, S, V, :', i, U(i, 1:nn), S(i), VT(i, 1:nn)
    !enddo
    !print*

    ! verification of SVD
    A = 0.
    A(1,1) = S(1)
    A(2,2) = S(2)
    A(1:2, 1:2) = matmul( U(1:2, 1:2), matmul(A(1:2, 1:2), VT(1:2, 1:2) ) )

!    do i=1,nn
!       write(*,'(a6, i5, 30es12.4)') 'A:', i, A(i, 1:nn)
!    enddo
    !print*,'###################################################'


    !do i=0, N
    !   xi(1) = cos(1.*i/N * 2*pi) * S(1)
    !   xi(2) = sin(1.*i/N * 2*pi) * S(2)
    !   Fxi(1:2) = matmul( U(1:2, 1:2), xi(1:2) ) + elem%xc(1:2)
    !   write(66,*) Fxi(1:2)
    !enddo
    !write(66,'(x)')

    ! do i = 0, 3
    !    j = mod(i, 3) + 1

    !    select case (j)
    !    case(1)
    !       xi(1:2) = (/ 0., 1. /)
    !    case(2)
    !       xi(1:2) = (/ -sqrt(3.)/2, -0.5 /)
    !    case(3)
    !       xi(1:2) = (/ sqrt(3.)/2, -0.5 /)
    !    end select

    !    !write(70,*) xi(1:2)


    !    !write(70,*)  F_K(1, 1) * xi(1) + F_K(1, 2) * xi(2)  + F_K(1, 3), &
    !    !     F_K(2, 1) * xi(1) + F_K(2, 2) * xi(2)  + F_K(2, 3)
    !    !write(70, '(x)')

    !    xi(1:2) = matmul(VT(1:2, 1:2), xi(1:2) )

    !    xi(1) = xi(1) * S(1)
    !    xi(2) = xi(2) * S(2)

    !    Fxi(1:2) = matmul( U(1:2, 1:2), xi(1:2) ) + elem%xc(1:2)
    !    write(67,*) Fxi(1:2)

    !    !write(69,*) xi(1:2)  + elem%xc(1:2)
    !    !write(69,*) Fxi(1:2)
    !    !write(69, '(x)')

    ! enddo
    ! write(67, '(x)')
    ! write(69, '(x)')


    ! original values of the anisotropy of the actual element (before adaptation)
    lambda_K = sqrt( S(1) * S(2) )
    sigma_K = sqrt(  S(1) / S(2) )
    if(U(1,1) > 1) U(1,1) = 1
    if(U(1,1) < -1) U(1,1) = -1
    phi_K = acos(U(1,1) )


    !write(*,'(a8, i5, 4es12.4, a2, 10es12.4)') &
    !     ' ANIS:', elem%i, sqrt(elem%area /scale), lambda_K, sigma_K, phi_K !, '|',  &

    !stop '3t733uhh: ama-hp_DWR.f90'
!    if(ndim > 1 .and. elem%i == 1) &
!      print*, "Subroutine Set_DWR_element_size - does it work for Euler???"


    ! NEW value of the element size
    if(elem%estim_loc <= 0.) then
       ratio = ratio_limit_max  ! case when error is 0.

    else
       ratio = tol_K / elem%estim_loc   !!!eta( dwr_aver ,1) )

       !!ifile = 500 + 10 * state%space%adapt%adapt_level

       if(ratio < 1.) then
          ratio = ratio**(1./( 2. * elem%deg + 0 )) !+ 2) ) ! VD: 18/03/01  IS IT OK??
          !!write(ifile+1, *) elem%xc(1:2), ratio, elem%deg
       else
          ratio = ratio**(1./( 2. * elem%deg + 0) )
          !!write(ifile+2, *) elem%xc(1:2), ratio, elem%deg
       endif

    !ratio = 1.

!       if(ratio > 1) then
!          write(500+10*state%space%adapt%adapt_level+1, *) elem%xc(:), ratio,  tol_K / max(1E-25, elem%estim_loc)
!       else
!          write(500+10*state%space%adapt%adapt_level+2, *) elem%xc(:), ratio,  tol_K / max(1E-25, elem%estim_loc)
!       endif

       ratio = max( ratio_limit_min, min ( ratio_limit_max,  ratio) )

       ! NEW APPROACH, ESCO 2018
       ratio_old = ratio
       ratio =  RefineFactor(elem%estim_loc)

       ! write(100 + state%space%adapt%adapt_level, *) &
    endif

!    if(ratio > 1) then
!       write(500+10*state%space%adapt%adapt_level+3, *) elem%xc(:), ratio,  tol_K / max(1E-25, elem%estim_loc)
!    else
!       write(500+10*state%space%adapt%adapt_level+4, *) elem%xc(:), ratio,  tol_K / max(1E-25, elem%estim_loc)
!    endif


    ! new lambda, ratio-times has to be increases/descreased
    lambda = lambda_K * ratio

    ! maximal element size (at most 10% of the area of the computatinal domain)
    lambda_max = sqrt( 0.1  * state%space%domain_volume / scale )

    ! minimal element size (at least 1E-14 of the area of the computatinal domain)
    lambda_min = sqrt( 1E-14  * state%space%domain_volume / scale )

    !if( lambda < lambda_min) write(999,*) elem%xc(:), lambda, lambda_min

    lambda = max(lambda_min, min(lambda_max, lambda) )

    !if(lambda**2 * scale > 0.1 * state%space%domain_volume) &
    !     write(*,'(a8, 4es12.4, a2, 10es12.4)') &
    !     ' AREAS:', elem%area , S(1)*S(2)*scale, lambda_K,  sqrt(S(1)*S(2)), '|',  &
    !     lambda, lambda**2 * scale, state%space%domain_volume
    !!ratio, tol_K / elem%eta( dwr_aver ,1), (tol_K / elem%eta( dwr_aver ,1))**(elem%deg/4.)
    !!     ratio, lambda_old, lambda


    deallocate(A, S, U, VT, xi, Fxi, F_K)
    !print*,'_________________________________________________________________________676___'

  end subroutine Set_DWR_element_size


  !> setting of the optimal anisotropy based on the DWR error estimates
  subroutine Set_DWR_metric(elem, ideg, lambda,  sigma, phi, degR)
    class(element), target, intent(inout) :: elem
    integer, intent(in) :: ideg  ! anisotropy driven by polynomial function of degree elem%deg+ideg
    real, intent(in)  :: lambda      ! set size of the element
    real, intent(in) :: sigma        ! optimal aspect ratio
    real, intent(in) :: phi          ! optimal orientation
    real, intent(in) :: degR          ! optimal orientation
    real :: lam_max, lam_min
    integer :: imt, degP
    real :: sigma_set, pi2, phiP
    real :: h_min, h_max, l_max, l_min
    !real :: sigma_limit_max = 50. ! 50. correspondance with AMA%pos: sigma <= 1/(1.5*pos)
    real :: sigma_limit_max = 500. ! 50. correspondance with AMA%pos: sigma <= 1/(1.5*pos)

    sigma_set = min(sigma, sigma_limit_max )
    !sigma_set = 1.

    pi2 =  asin(1.0)
    phiP = phi + pi2

    lam_max =  ( 1./lambda * sigma_set)**2
    lam_min =  ( 1./lambda / sigma_set)**2

    ! maximal and minimal admissible length
    h_max =  min( state%space%diam / 4., 4.)
    h_min = h_max * 1E-8

    l_max = 1. / h_min**2
    l_min = 1. / h_max**2


    ! if( abs(elem%xc(1) -0.6) < 0.04 .and. abs(elem%xc(2) -0.05) < 0.05 &
    !      .and. state%space%adapt%adapt_level >= 4) then
    !    write(*,'(a10, 8(2es12.4,a2))') 'lams, sims:', &
    !         lambda, sigma_set, '|',lambda * sigma_set, lambda / sigma_set, '|',&
    !         lam_max, lam_min, '|',l_max, l_min, '|'!,1./sqrt(lam_max), 1./sqrt(lam_min)
    ! endif

    !if( lam_max > l_max) write(999,*) elem%xc(:), lam_max, l_max

    ! limitation of the metric eigenvalues with respect h_max, h_min
    lam_max = max(l_min, min (l_max, lam_max) )
    lam_min = max(l_min, min (l_max, lam_min) )

    ! if( abs(elem%xc(1) -0.6) < 0.04 .and. abs(elem%xc(2) -0.05) < 0.05 &
    !      .and. state%space%adapt%adapt_level >= 4) then
    !    write(*,'(a10, 8(2es12.4,a2))') 'lams, sims:', &
    !         lambda, sigma_set, '|',1./sqrt(lam_max), 1./sqrt(lam_min), '|',&
    !         lam_max, lam_min, '|', l_max, l_min
    ! endif

    ! setting of the metric
    elem%rgabc(1) =  lam_max * cos(phiP)**2 + lam_min * sin(phiP)**2
    elem%rgabc(2) = (lam_max - lam_min) * cos(phiP) * sin(phiP)
    elem%rgabc(3) =  lam_max * sin(phiP)**2 + lam_min * cos(phiP)**2


    degP = int(degR+0.5)
    elem%psplit = degP - 1
    elem%ama_p = elem%psplit

    !if(degP /= 1) write(*,'(a8, 2i5, 2es12.4, 4i5)') &
    !     'i_min:', elem%i, degP, degR, elem%ama_p, elem%deg, elem%deg + elem%psplit

    !imt = 800 + state%space%adapt%adapt_level
    !call DrawEllips(imt, elem%rgabc(1:3), elem%xc(1:2) )

  end subroutine Set_DWR_metric


  !> setting of the optimal hp-anisotropy based on the DWR error estimates
  !> testing \f$ p_K,\ p_K + 1,  \ p_K - 1 \f$
  subroutine Set_DWR_hp_anisotropy(elem,  lambda_K, sigma_K, phi_K, &
       lambda,  sigma_opt, phi_opt, deg_opt)
    class(element), target, intent(inout) :: elem
    real, intent(inout) ::  lambda_K, sigma_K, phi_K  ! original element anisotropy
    real, intent(inout) :: lambda    ! prescribed size of the element, can be modified
    real, intent(out) :: sigma_opt   ! optimal aspect ratio
    real, intent(out) :: phi_opt     ! optimal orientation
    integer, intent(out) :: deg_opt  ! optimal polynomial degree (elemdeg + deg_opt - 1)
    real, dimension(:, :), allocatable :: params
    real :: lam_in, sigma, phi, eta, eta_min
    integer :: ideg, DoF, DoFP, i_min
    logical :: iprint

    iprint = .false.
    !if(  (abs( elem%xc(1) - 0.1875) < 1./16 .and. abs( elem%xc(2) - 0.1875) < 1./16) .or. &
    !     (abs( elem%xc(1) - 0.8125) < 1./16 .and. abs( elem%xc(2) - 0.8125) < 1./16)) then
    !if( elem%i == 10 .or. elem%i == 236) then
    !   iprint = .true.
    !endif


    allocate(params(0:2, 1:5), source = 1.0  )

    DoF = (elem%deg +1 )*(elem%deg + 2)

    eta_min = 1E+20
    do ideg = 0, 2   ! testing p_K-1, p_K, p_K+1

         ! minimal & maximal degree for the HO derivatives
       if(elem%deg + ideg >  1 .and. elem%deg + ideg <= MaxDegreeImplemented) then

          DoFP = (elem%deg + ideg )*(elem%deg + ideg + 1)

          ! input lambda computed from lambda_K corresponding to p_K
          lam_in = lambda_K * sqrt(1. * DoFP / DoF)

          ! setting of the optimal anisotropy for p_K + 1 + ideg
          call Set_DWR_anisotropy(elem, ideg, lam_in, sigma_K, phi_K, &
               lambda, sigma, phi, eta, iprint)

          params(ideg, 1 ) = lam_in
          params(ideg, 2 ) = sigma
          params(ideg, 3 ) = phi
          params(ideg, 4 ) = eta

          !!!eta = eta *( 1 +  (2.-ideg)/10)  ! preference for higher degrees

          if(eta < eta_min) then
             eta_min = eta
             i_min = ideg
          endif

          if(iprint) then

             write(*,'(a8, i5, 8es12.4)') 'estL', ideg, lambda, lam_in, params(ideg, 1:4)
             if(ideg == 2) print*,'___________', elem%i, elem%xc(:)
             if(ideg == 2) print*
          endif

       endif
    enddo

    ! ONLY A SMALL DECREASE OF THE INTERPOLATION ERROR
    !if(eta_min * 1.5 > params(1, 4) ) i_min = 1

    if(elem%deg == 1) i_min = max(i_min, 1)   ! minimal allowed deg
    if(elem%deg + 1 == MaxDegreeImplemented) i_min = min(i_min, 1)   ! maximal allowed deg

    ! setting of optimal hp-anisotropy
    lambda    = lambda * params(i_min, 1) / lambda_K  ! modifying the target size
    sigma_opt = params(i_min, 2)
    phi_opt   = params(i_min, 3)
    deg_opt   = i_min

    !if(i_min /= 1) write(*,'(a8, 5i5)') 'i_min:', elem%i, i_min, elem%deg

    !write(100+state%space%adapt%adapt_level*10+i_min ,'(a8, 2es12.4, i5, 8es12.4)') &
    !     'est_opt', elem%xc(:), i_min, params(i_min, 4), params(0:2, 4)/ params(i_min, 4)
    !lambda, lambda, params(i_min, 1:4)
    !print*
    !if(elem%i <= 3) print*,"Verify lambda, lambda_opt"

  end subroutine Set_DWR_hp_anisotropy

    !> setting of the optimal anisotropy based on the DWR error estimates
  subroutine Set_DWR_anisotropy(elem, ideg,  lambda_K, sigma_K, phi_K, &
       lambda,  sigma_opt, phi_opt, eta_min, iprint)
    class(element), target, intent(inout) :: elem
    integer, intent(in) :: ideg  ! anisotropy driven by polynomial function of degree elem%deg+ideg
    real, intent(inout) ::  lambda_K, sigma_K, phi_K  ! original element anisotropy
    real, intent(in)  :: lambda      ! prescribed size of the element (not used?)
    real, intent(out) :: sigma_opt   ! optimal aspect ratio
    real, intent(out) :: phi_opt     ! optimal orientation
    real, intent(inout) :: eta_min   ! value of the minimized error estimate
    logical, intent(in), optional :: iprint
    real, dimension(:), allocatable :: ai, mi, Re_1, wR, yi
    real, dimension(:,:), allocatable :: mA, xi, Fxi, etasR, etasA, Du, Dz
    real, dimension(:,:,:), allocatable :: anis
    logical :: upper   ! correct upper bound (with expand of bound)
    real :: eta, pi, wi(1:ndim), Dwi(1:ndim, 1:2), val, val1, val2
    integer :: deg, degP, degQ, ieta, optim, dof1, dof2
    integer :: Qdof, i, j,l, kk
    real :: sigma_increase_max = 500. ! 500 == NO limitation
    logical :: diff_term

    if (state%model%Re > 0.0) then
      diff_term = .true.
    else
      diff_term = .false.
    end if

    upper = .true.
    !upper = .false.

    pi =  asin(1.0)*2.0

    optim = 0
    ! optim = -10  ==> output print
    !if( abs(elem%xc(1) - 1.) < 0.02 .and. abs(elem%xc(2) - 1.) < 0.02 ) optim = -10


    degP = elem%deg + ideg    ! degree of the interpolation error function
    deg = degP - 1
    degQ = 2*deg            ! degree of the square of the gradient of the interpolation error function

    ! array for the anisotropy
    !allocate(anis(1:4, 1:5) )
    ! 1:4 1 - w, 2 - grad(w), 3 - z, 4 - grad(z)
    ! 1:3 - A_w, \phi, \sigma,
    allocate(anis(1:ndim, 1:4, 1:3), source = 0.0 )

    allocate(mA( 1:2, 1:2), source = 0.0 ) ! diffusion matrix

    if(state%modelName == 'scalar' .and. diff_term) then
       allocate(Re_1(1:iRe), source = 0.0 )
       if(state%model%Re > 0.) Re_1(1) = 1./state%model%Re
       Re_1(2:iRe) =  elem%xi(0, 1, 2+1: 2+iRe-1)

       ! wi = elem%wST_LS(1, 1, 1) / sqrt(2) ! P_0 projection on element
       call Eval_aver_w_Elem(elem, wi(1:ndim))
       call Eval_aver_Dw_Elem(elem, Dwi(1:ndim, 1:2) )

       call Set_K_sk_scalar(ndim, nbDim, iRe, 1, wi(1:ndim), Dwi(1:ndim, 1:nbDim), &
            Re_1(1:iRe), &
            mA(1:nbDim, 1:nbDim), elem%xc(1:nbDim) )

       deallocate(Re_1)
       !write(*,'(a5, 5es12.4)') 'w,Dw', wi, Dwi, Re_1(1)
       !write(*,'(a5, 2es12.4)') 'mA=', mA(1,:)
       !write(*,'(a5, 2es12.4)') 'mA=', mA(2,:)
       !print*,'___________________________'
    ! EULER
    else if (ndim == 4 .and. state%model%Re == 0.0 .and. state%modelName == 'NSe') then
        ! we do not need mA and Dwi since there is no diffusion
    else
       stop "NOT implemented in subroutine Set_DWR_anisotropy, file ama-hp_interpol"
    endif

    allocate(Du( 1:ndim, 0:degP), source = 0.0 )  ! derivatives of order deg
    allocate(Dz( 1:ndim, 0:degP), source = 0.0 )  ! derivatives of order deg


    ! wSS(1:2*ndim,0:2,0:elem%deg+2)
    if (size(elem%wSS,1) /= 2*ndim) then
      stop "Problem in Set_DWR_anisotropy, wrong size of elem%wSS!"
    end if

    !print*, "For Euler we have to do it differently - eta(1:ndim)*weight(1:ndim)"
!    Du(0:degP) =  elem%wSS(1, ideg, 0:degP )
!    Dz(0:degP) =  elem%wSS(2, ideg, 0:degP )
    ! the derivatives of the highest order
    Du(1:ndim, 0:degP) =  elem%wSS(1:ndim, ideg, 0:degP )
    Dz(1:ndim, 0:degP) =  elem%wSS(ndim+1:2*ndim, ideg, 0:degP )

!    if (elem%i == 1) then
!      print*, "Du =", norm2(Du(1:ndim, 0:degP) )
!      print*, "Dz =", norm2(Dz(1:ndim, 0:degP) )
!    end if

    allocate( ai(0:2*degP), source = 0.0 )  ! coefficients

    ! function
    do kk = 1, ndim

      !! PRIMAL SOLUTION
      ! set \alpha_l coefficients in ESCO article under (28)
      call Set_Coeffs_Function(degP, Du(kk,0:degP), ai(0:degP) )
      !call Draw_deriv(ifile, degP, ai(0:degP) )
      call FindAnisotropy(degP, ai(0:degP), anis(kk,1, 1:3), upper )
      !  gradient
      if (diff_term) then
        call Set_Coeffs_Gradient(degP, degQ, Du(kk, 0:degP), mA(1:2, 1:2), ai(0:degQ) )
        !call Draw_deriv(ifileD, degQ, ai(0:degQ) )
        call FindAnisotropy(degQ, ai(0:degQ), anis(kk, 2, 1:3), upper )
        !call Draw_anis(ifileD2, degQ, anis(2, 1:3) )
      end if

      ! DUAL SOLUTION
      ! function
      call Set_Coeffs_Function(degP, Dz(kk, 0:degP), ai(0:degP) )
      call FindAnisotropy(degP, ai(0:degP), anis(kk, 3, 1:3), upper )
      ! its gradient
      if (diff_term) then
        call Set_Coeffs_Gradient(degP, degQ, Dz(kk,0:degP), mA(1:2, 1:2), ai(0:degQ) )
        call FindAnisotropy(degQ, ai(0:degQ), anis(kk, 4, 1:3), upper )
      end if

    end do

    !if(present(iprint) .and. iprint) then
    !   write(*,'(a12, 3es16.8)') 'ANIS of u:',    anis(1, 1:3)
    !   write(*,'(a12, 3es16.8)') 'ANIS of A Du:', anis(2, 1:3)
    !   write(*,'(a12, 3es16.8)') 'ANIS of z:',    anis(3, 1:3)
    !   write(*,'(a12, 3es16.8)') 'ANIS of A Dz:', anis(4, 1:3)
    !endif

    ! particular residual estimates

    ieta = 6 ! should be alway 6
    allocate(etasA(1:ieta,1:ndim), source = 0.0)
    allocate(etasR(1:ieta,1:ndim), source = 0.0)
    etasR(1, 1:ndim) = elem%eta(dwrEtaKV, 1:ndim)
    etasR(2, 1:ndim) = elem%eta(dwrEtaKB, 1:ndim)
    etasR(3, 1:ndim) = elem%eta(dwrEtaKD, 1:ndim)
    etasR(4, 1:ndim) = elem%eta(dwrEtaKV_dual, 1:ndim)
    etasR(5, 1:ndim) = elem%eta(dwrEtaKB_dual, 1:ndim)
    etasR(6, 1:ndim) = elem%eta(dwrEtaKD_dual, 1:ndim)

    if(present(iprint) .and. iprint) then
       write(*,'(a12, 30es16.8)') 'etasR:', etasR(1:6,1:ndim)
    endif

    !  verification of anisotropic error estimates
    ! TODO: FR_ANI not needed ??? computes etasA
    call Anis_estim_tot(ieta, deg, anis(1:ndim, 1:4, 1:3), lambda_K, sigma_K, phi_K, &
         etasR(1:ieta, 1:ndim), etasA(1:ieta, 1:ndim), eta, optim, diff_term )
    eta = eta / 2  !VD

!    if (elem%i == 1) then
!      print*, "etasR = ", etasR
!      print*, "etasA = ", etasA
!      print*, "tot eta = ", eta
!    end if



    ! TODO: FR_ANI eta is ALGEBRAIC error ??? ???
    ! estimate of type \bar{\eta}^{II}  (or \eta^{III}   !VD
    if(elem%i == 1)  state%estim( dwrA ,1:ndim) = 0.
    state%estim( dwrA ,1 ) = state%estim( dwrA ,1 ) + eta !!!elem%eta( dwr_aver ,1)

    ! ! test OUTPUTS
    ! if( abs(elem%xc(1) -0.6) < 0.04 .and. abs(elem%xc(2) -0.05) < 0.05 &
    !      .and. state%space%adapt%adapt_level >= 4) then
    !    !if(elem%i <= 10) &
    !    !if( eta/elem%eta( dwr_aver ,1) > 10) &
    !    !if( abs(elem%xc(1) - 1.) < 0.02 .and. abs(elem%xc(2) - 1.) < 0.02 ) &
    !    ! if(optim == -10) then
    !    !    !write(500+state%space%adapt%adapt_level ,'(2es12.4, i5, 3es12.4, a2, 10es12.4)') &
    !    !    write(* ,'(a8, i5, 3es12.4, a2, 10es12.4)') &
    !    !         ' ANIS:', &
    !    !                               !elem%xc(:),
    !    !         elem%i, lambda_K, sigma_K, phi_K ,'|', &
    !    !         eta, elem%eta( dwr_aver ,1), eta/elem%eta( dwr_aver ,1), &
    !    !                               !etasA(1) / elem%eta(dwrWeightKV, 1), &
    !    !                               !etasA(2) / elem%eta(dwrWeightKB, 1), &
    !    !                               !etasA(3) , elem%eta(dwrWeightKD, 1), &
    !    !                               !etasA(3) / elem%eta(dwrWeightKD, 1) , &
    !    !                               !etasA(4) / elem%eta(dwrWeightKV_dual, 1), &
    !    !                               !etasA(5) / elem%eta(dwrWeightKB_dual, 1), &
    !    !                               !etasA(6) / elem%eta(dwrWeightKD_dual, 1)
    !    !         etasA(6) , elem%eta(dwrWeightKD_dual, 1), &
    !    !         etasA(6) / elem%eta(dwrWeightKD_dual, 1) !, &

    !    !    if( etasA(6) / elem%eta(dwrWeightKD_dual, 1) > 1E+8) then
    !    !       dof2 = size ( elem%wST_LS, 2)
    !    dof2 = DOFtriang(elem%deg +1)
    !    dof1 = DOFtriang(elem%deg)
    !    print*,'###############################################################'
    !    print*,'dof2 = ', dof2, dof1
    !    !write(*,'(30es12.4)') elem%wST_LS(1, 1:, 1)
    !    write(*,'(30es12.4)') elem%wST_LS(1, 1:dof2, 1)
    !    write(*,'(30es12.4)') elem%zST_LS(1, 1:dof2, 1)
    !    allocate(wR(1:dof2) )
    !    wR(1:dof2) = elem%wST_LS(1, 1:dof2, 1)
    !    call PlotElemFunction3D(99, elem, dof2, wR(1:dof2)  )
    !    !       call PlotElemFunction3D(199, elem, dof2, wR(1:dof2)  )

    !    wR = 0.
    !    wR(1:dof1) =  elem%wST_LS(1, 1:dof1, 1)
    !    call PlotElemFunction3D(98, elem, dof2, wR(1:dof2)  )
    !    !       call PlotElemFunction3D(198, elem, dof2, wR(1:dof2)  )

    !    wR = 0.
    !    wR(1:dof2) = elem%wST_LS(1, 1:dof2, 1)
    !    wR(1:dof1) = 0.
    !    call PlotElemFunction3D(97, elem, dof2, wR(1:dof2)  )
    !    !       call PlotElemFunction3D(197, elem, dof2, wR(1:dof2)  )

    !    Qdof = elem%Qdof
    !    allocate(xi (1:Qdof, 1:2), Fxi(1:Qdof, 1:2), yi(1:2) )
    !    call Set_Coeffs_Gradient(degP, degQ, Du(0:degP), mA(1:2, 1:2), ai(0:degQ) )
    !    call FindAnisotropy(degQ, ai(0:degQ), anis(2, 1:3), .true. )
    !    call ComputeF(elem, Qdof, state%space%V_rule(elem%Qnum)%lambda(1:Qdof,1:nbDim), &
    !         Fxi(1:Qdof,1:nbDim) )
    !    if( elem%i == 2093 ) then
    !       call Draw_deriv(71, degQ, ai(0:degQ) )
    !       call Draw_anis(72, degQ, anis(2, 1:3) )
    !    endif
    !    print*,'degQ = ', degQ, Qdof, ', elem%i = ', elem%i
    !    write(*, '(a12,50es12.4)' ) 'ai(0:degQ) = ', ai(0:degQ)
    !    write(*, '(a12,50es12.4)' ) 'ANIS 1:', anis(1, 1:3)
    !    write(*, '(a12,50es12.4)' ) 'ANIS 2:', anis(2, 1:3)
    !    write(*, '(a12,50es12.4)' ) 'ANIS 3:', anis(3, 1:3)
    !    write(*, '(a12,50es12.4)' ) 'ANIS 4:', anis(4, 1:3)
    !    !       !do i=1, Qdof
    !    !       !xi(i, 1:2) = Fxi(i, 1:2) - elem%xc(1:2)
    !    do i = 1, 7
    !       do l = 0, 100
    !          !xi(i, 1:2) = Fxi(i, 1:2) - elem%xc(1:2)
    !          xi(i, 1) = cos( 2.*pi * l /100) * elem%diam / 5 * (i-1) /2
    !          xi(i, 2) = sin( 2.*pi * l /100) * elem%diam / 5 * (i-1) /2
    !          Fxi(i, 1:2) = xi(i, 1:2) + elem%xc(1:2)
    !          val = 0.
    !          do j=0, degQ
    !             val = val + ai(j) * xi(i, 1)**(degQ-j) * xi(i, 2)**j
    !          enddo
    !          yi(1) =  cos( anis(2, 3) ) * xi(i, 1) + sin(anis(2,3) ) * xi(i,2)
    !          yi(2) = -sin( anis(2, 3) ) * xi(i, 1) + cos(anis(2,3) ) * xi(i,2)
    !          val1 = anis(2, 1) * &
    !               ( yi(1)**2 + yi(2)**2 * anis(2,2)**(-2./degQ) )**(degQ/2.)
    !          val2 =  anis(2, 1) * ( dot_product(yi(1:2), yi(1:2) ))**(degQ/2.)
    !          write(89, '(50es12.4)' ) Fxi(i, 1:2), val, val1, val2, & !xi(i, 1:2), &
    !               anis(2, 1:2),  dot_product(yi(1:2), yi(1:2)) , &
    !               yi(1)**2 + yi(2)**2 * anis(2,2)**(-2./degQ)
    !       enddo
    !    enddo
    !    !       write(81, '(50es12.4)' ) elem%xc(1:2), 1E-5, 1E-5
    !    !       stop "73yd3hd3ijdswijsiw"
    !    !    endif
    ! endif


    ! seeking of the optimal anisotropy, minimum of the functional
    eta_min = 1E+20

    !call SeekMinimumGlobally(ieta, deg, anis(1:4, 1:3), lambda_K, sigma_opt, phi_opt, &
    !     etasR(1:ieta), etasR(1:ieta), eta_min, diff_term)
    !print*
    !print*,'               eta            sigma              phi '
    !print*,'-----------------------------------------------------------'
    !!write(*, '(a10, 60es16.8)') 'globaly:',eta_min, sigma_opt, phi_opt
    !write(99, '(60es16.8)') sigma_opt, phi_opt, eta_min
    !write(99, '(60es16.8)') sigma_opt, phi_opt+pi, eta_min


    ! SEEKING of the global minimum iteratively
    sigma_opt = sigma_K   ! intial approximations
    phi_opt = phi_K       ! intial approximations
    call SeekMinimumIteratively(ieta, deg, anis(1:ndim, 1:4, 1:3), lambda_K, sigma_opt, phi_opt, &
         etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim), eta_min, elem%i, diff_term)

    eta_min = eta_min /2.

    ! shape restrictions
    sigma_opt = min(sigma_opt, sigma_K * sigma_increase_max)

!    if (elem%i == 1) then
!      print*, 'after seek:', elem%i, eta_min/2 , sigma_opt, phi_opt
!    end if

    !print*
    !!write(98, '(60es16.8)') sigma_opt, phi_opt, eta_min
    !!write(98, '(60es16.8)') sigma_opt, phi_opt+pi, eta_min


    !if( abs(elem%xc(1) -0.6) < 0.04 .and. abs(elem%xc(2) -0.05) < 0.05 &
    !     .and. state%space%adapt%adapt_level >= 4) then
    !   write(*,'(a10, 20es12.4)') 'iterative seeks:', &
    !        eta_min, lambda_K, lambda,sigma_opt, phi_opt
    !endif


    ! SEEKING of the global minimum using Newton-like method
    !call SeekMinimumNewton(ieta, deg, anis(1:4, 1:3), lambda_K, sigma_opt, phi_opt, &
    !     etasR(1:ieta), eta_min)

    !write(*,'(a10, 20es12.4)') 'Newton:', eta_min, sigma_opt, phi_opt
    !print*

    !write(*,'(a8, i5, 3(3es12.4, a2))') &
    !     ' ANIS3:', elem%i, lambda_K, sigma_K, phi_K ,'|', lambda, sigma_opt, phi_opt,'|',&
    !     lambda_K / lambda


    deallocate(anis, mA, Du, Dz, ai, etasA, etasR)


  end subroutine Set_DWR_anisotropy



  !> seeking of the minimum by an iterative process
subroutine SeekMinimumIteratively(ieta, deg, anis, lambda, sigma_opt, phi_opt, &
     etasR, etasA, eta_min, ie, diff_term)
  integer, intent(in) :: ieta    ! number of estimator
  integer, intent(in) :: deg     ! current polynomial degree
  real, dimension(1:ndim, 1:4, 1:3), intent(in) :: anis  ! anisotropy of the function
  real, intent(in) ::    lambda   ! size of the triangle
  real, intent(inout) :: sigma_opt, phi_opt    ! anisotropy of the triangle
  real, dimension(1:ieta, 1:ndim), intent(in) :: etasR  ! computed primal and dual residuas
  real, dimension(1:ieta, 1:ndim), intent(inout) :: etasA  ! estimates of the particular "weights"
  real, intent(out) :: eta_min    ! error estimate
  integer, intent(in) :: ie   ! element index
  logical, intent(in) :: diff_term  ! problem with nonzero diffusion -> True
  integer :: i, j, k, optim, max_iter, max_iter_one_dir, iter, i_accur, n_accur
  real :: pi, q_sigma, D_phi, eta_tot, D_phi_init, q_sigma_init
  real :: sigma, phi, sigmaP, sigmaM, phiP, phiM, eta_totP, eta_totM
  logical ::  next_step
  integer ::  conver_sigma , conver_phi, ieta_min
  logical :: sigma_stepP, sigma_stepM, phi_stepP, phi_stepM
!  real, dimension(:,:), allocatable  :: params_init, sigma_init, phi_init


  ! TODO: FR_ANI what is this for? weird ordering
  if (ie == 1)  print*, "Find_Params_init is not used anymore!"
  !allocate( params_init(1:ieta, 1:5) )
  !call Find_Params_init(ieta, deg, anis(1:4, 1:3), params_init(1:ieta, 1:3) )
  ! array for the initial guess and the resulting minima for each local function
  ! call Find_init_params_new
  !allocate( sigma_init(1:ieta), source = 1.0 )
  !allocate( phi_init(1:ieta), source = 0.0 )

!  optim = 0
!  do i =1, ieta
!      ! TODO: FR_ANI What is this for? etasA is rewritten each time
!      call Anis_estim_tot(ieta, deg, anis(1:ndim,1:4, 1:3), lambda, &
!          sigma_init(i), phi_init(i), etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim), eta_min, optim)
!      !write(80,*) params_init(i, 1:2), eta_min
!  enddo
  ! FR_ANI: What is the output of the above??? COMMENTED


  pi =  asin(1.0)*2.0

  optim = 0
  max_iter = 150         ! maximal number of inner loops
  max_iter_one_dir = 20   ! maximal number of steps in one direction

  ! initial "sigma" and "phi steps
  D_phi_init = pi / 10
  q_sigma_init = 1.5

  ! number of "levels of accuracy"  ==> accuracy is the D_phi(_init)/ 2**n_accur
  ! corresponds to the number of outer loops
  n_accur = 5

  !ieta_min = 1
  !ieta_min = 3
  !ieta_min = 5
  ieta_min = 6

  ! TODO: FR_ANI - what is this cycle for - it seems that is not needed
  ! FR_ANI Commented the cycle!
  ! we start at each local minimum
  !do k= ieta_min, ieta  !!  k==1 <==> k==3  &&  k==2 <=> k==4
   ! initial guess for iterative process
   ! TAKEN from sigma_K and phi_K
   !sigma_opt = params_init(k, 1)
   !phi_opt = params_init(k, 2)
  ! we take the original sigma and phi as the initial guess
  ! FR_ANI etasA are already computed - COMMENTED
!  call Anis_estim_tot(ieta, deg, anis(1:ndim,1:4, 1:3), lambda, sigma_opt, phi_opt,  &
!        etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim), eta_min, optim )

   ! initial "sigma" and "phi steps
   D_phi = D_phi_init
   q_sigma = q_sigma_init

   do i_accur = 1, n_accur  !

      D_phi = D_phi / 2
      q_sigma = sqrt(q_sigma)

      ! loop over phi optimization: in each iter, phi is fixed,
      ! for the fixed phi, we seek for the best sigma
      do iter = 1, max_iter
         conver_sigma = 0
         conver_phi = 0

         !varying phi
         sigma = sigma_opt

         phi_stepP = .true.
         phi_stepM = .true.

         do j=1,  max_iter_one_dir

            ! step phi := phi + D_phi
            if( phi_stepP ) then
               phiP = phi_opt + D_phi
               if(phiP > pi) phiP = phiP - pi

               call Anis_estim_tot(ieta, deg, anis(1:ndim, 1:4, 1:3), lambda, sigma, phiP,  &
                    etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim), eta_totP, optim, diff_term )
            endif

            ! step phi := phi - D_phi
            if( phi_stepM ) then
               phiM = phi_opt - D_phi
               if(phiM < 0) phiM = phiM + pi

               call Anis_estim_tot(ieta, deg, anis(1:ndim, 1:4, 1:3), lambda, sigma, phiM,  &
                    etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim), eta_totM, optim, diff_term )
            endif

            if(eta_min <= eta_totP .and. eta_min <= eta_totM) then
               ! minimum was found
               !conver_sigma = conver_sigma
               phi_stepP = .false.
               phi_stepM = .false.

            elseif (eta_totP < eta_min .and. eta_totP <= eta_totM) then
               phi_opt = phiP
               eta_min = eta_totP
               phi_stepM = .false.
               conver_phi = conver_phi + 1

            elseif (eta_totM < eta_min) then
               phi_opt = phiM
               eta_min = eta_totM
               phi_stepP = .false.
               conver_phi = conver_phi + 1
            else
               write(*,'(a8,6es16.8)') 'etas:',eta_min, eta_totP, eta_totM
               stop "strange 33ju3 A"
            endif
            !!write(80+k, '(60es12.4)') sigma_opt, phi_opt, eta_min

            if(.not. phi_stepP .and. .not. phi_stepM) goto 25
         end do  ! do j=1,20
25       continue ! phi_opt was found for the given step size D_phi and sigma

         ! varying of sigma
         phi = phi_opt

         sigma_stepP = .true.
         sigma_stepM = .true.

         do j=1,  max_iter_one_dir

            ! step sigma = sigma * q_sigma
            if( sigma_stepP ) then
               sigmaP = sigma_opt * q_sigma
               !!if(sigmaP > 5) then  ! direct limitation of size of sigma
               !!   eta_totP = 2 * eta_min
               !!else

               call Anis_estim_tot(ieta, deg, anis(1:ndim, 1:4, 1:3), lambda, &
                        sigmaP, phi, etasR(1:ieta,1:ndim), &
                        etasA(1:ieta,1:ndim), eta_totP, optim, diff_term )

               !!endif
            endif

            ! step sigma = sigma / q_sigma
            if( sigma_stepM ) then
               sigmaM = sigma_opt / q_sigma
               if(sigmaM < 1) then   ! sigma must be at least one
                  eta_totM = 2 * eta_min
               else
                  call Anis_estim_tot(ieta, deg, anis(1:ndim,1:4, 1:3), lambda, &
                        sigmaM, phi, etasR(1:ieta,1:ndim), etasA(1:ieta,1:ndim),&
                        eta_totM, optim, diff_term)
               endif

            endif

            if(eta_min <= eta_totP .and. eta_min <= eta_totM) then
               ! minimum was found
               !conver_sigma = conver_sigma
               sigma_stepP = .false.
               sigma_stepM = .false.

            elseif (eta_totP < eta_min .and. eta_totP <= eta_totM) then
               sigma_opt = sigmaP
               eta_min = eta_totP
               sigma_stepM = .false.
               conver_sigma = conver_sigma + 1

            elseif (eta_totM < eta_min) then
               sigma_opt = sigmaM
               eta_min = eta_totM
               sigma_stepP = .false.
               conver_sigma = conver_sigma + 1
            else
               write(*,'(a8,6es16.8)') 'etas:',eta_min, eta_totP,  eta_totM
               stop "strange 33ju3 B"
            endif
            !!write(80+k, '(60es12.4)') sigma_opt, phi_opt, eta_min, q_sigma, D_phi

            if(.not. sigma_stepP .and. .not. sigma_stepM) goto 15

            if(sigma <= 1.) then
               sigma_opt = 1
               phi_opt = phi ! arbitrary
               goto 15
            endif

         enddo   ! do j=1, max_iter_one_dir
15       continue
         ! FR_ANI what is this for?
         if(conver_phi == 0 .or. conver_sigma == 0) goto 20
          !print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~', iter
      enddo  !  iter = 1, max_iter
20    continue

      !!write(70+k, '(3es12.4, i5)') sigma_opt, phi_opt, eta_min, i_accur


   enddo  ! do i_accur = 1, n_accur


!   !endif
!   params_init(k, 3) = sigma_opt
!   params_init(k, 4) = phi_opt
!   params_init(k, 5) = eta_min

  ! FR_ANI commented the cycle
  !enddo  ! k=1, ieta

!  ! seeking of the minimum from several candidates
!  eta_min = 1E+30
!  do k= ieta_min,ieta
!     if(params_init(k, 5) < eta_min) then
!        eta_min   = params_init(k, 5)
!        sigma_opt = params_init(k, 3)
!        phi_opt   = params_init(k, 4)
!     endif
!
!     ! ! not necessary for the computation, only the value at the initial gues
!     ! ! can be commented
!     ! call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, &
!     !      params_init(k, 1), params_init(k, 2),  &
!     !      etasR(1:ieta), eta_tot, optim )
!
!
!     ! write(*,'(a5, i3, 6es12.4, a2, 2es12.4)') &
!     !      'resu:', k, eta_tot, params_init(k, 1:5), '|', q_sigma, D_phi
!     !write(900+state%space%adapt%adapt_level, *)ie, params_init(k, 3:5)
!  enddo
!  !write(900+state%space%adapt%adapt_level, '(x)')

!  deallocate( params_init, sigma_init, phi_init)

  !print*,'~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~', iter, optim


  !!! OLD VERSION
!   ! TODO: FR_ANI - what is this cycle for - it seems that is not needed
!  ! we start at each local minimum
!  do k= ieta_min, ieta  !!  k==1 <==> k==3  &&  k==2 <=> k==4
!
!     ! initial guess for iterative process
!     ! TAKEN from sigma_K and phi_K
!     !sigma_opt = params_init(k, 1)
!     !phi_opt = params_init(k, 2)
!
!     call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigma_opt, phi_opt,  &
!          etasR(1:ieta), etasA(1:ieta), eta_min, optim )
!
!     !!write(80+k, '(60es12.4)') sigma_opt, phi_opt, eta_min
!
!     ! initial "sigma" and "phi steps
!     D_phi = D_phi_init
!     q_sigma = q_sigma_init
!
!     do i_accur = 1, n_accur  !
!
!        D_phi = D_phi / 2
!        q_sigma = sqrt(q_sigma)
!
!
!        ! loop over phi optimization: in each iter, phi is fixed,
!        ! for the fixed phi, we seek for the best sigma
!        do iter = 1, max_iter
!           conver_sigma = 0
!           conver_phi = 0
!
!           !varying phi
!           sigma = sigma_opt
!
!           phi_stepP = .true.
!           phi_stepM = .true.
!
!           do j=1,  max_iter_one_dir
!
!              ! step phi := phi + D_phi
!              if( phi_stepP ) then
!                 phiP = phi_opt + D_phi
!                 if(phiP > pi) phiP = phiP - pi
!
!                 call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigma, phiP,  &
!                      etasR(1:ieta), etasA(1:ieta), eta_totP, optim )
!              endif
!
!              ! step phi := phi - D_phi
!              if( phi_stepM ) then
!                 phiM = phi_opt - D_phi
!                 if(phiM < 0) phiM = phiM + pi
!
!                 call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigma, phiM,  &
!                      etasR(1:ieta), etasA(1:ieta), eta_totM, optim )
!              endif
!
!
!              if(eta_min <= eta_totP .and. eta_min <= eta_totM) then
!                 ! minimum was found
!                 !conver_sigma = conver_sigma
!                 phi_stepP = .false.
!                 phi_stepM = .false.
!
!              elseif (eta_totP < eta_min .and. eta_totP <= eta_totM) then
!                 phi_opt = phiP
!                 eta_min = eta_totP
!                 phi_stepM = .false.
!                 conver_phi = conver_phi + 1
!
!              elseif (eta_totM < eta_min) then
!                 phi_opt = phiM
!                 eta_min = eta_totM
!                 phi_stepP = .false.
!                 conver_phi = conver_phi + 1
!              else
!                 write(*,'(a8,6es16.8)') 'etas:',eta_min, eta_totP, eta_totM
!                 stop "strange 33ju3 A"
!              endif
!              !!write(80+k, '(60es12.4)') sigma_opt, phi_opt, eta_min
!
!              if(.not. phi_stepP .and. .not. phi_stepM) goto 25
!           end do  ! do j=1,20
!25         continue
!
!           ! varying of sigma
!           phi = phi_opt
!
!           sigma_stepP = .true.
!           sigma_stepM = .true.
!
!           do j=1,  max_iter_one_dir
!
!              ! step sigma = sigma * q_sigma
!              if( sigma_stepP ) then
!                 sigmaP = sigma_opt * q_sigma
!                 !!if(sigmaP > 5) then  ! direct limitation of size of sigma
!                 !!   eta_totP = 2 * eta_min
!                 !!else
!
!                 call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigmaP, phi,  &
!                      etasR(1:ieta), etasA(1:ieta), eta_totP, optim )
!
!                 !!endif
!              endif
!
!              ! step sigma = sigma / q_sigma
!              if( sigma_stepM ) then
!                 sigmaM = sigma_opt / q_sigma
!                 if(sigmaM < 1) then   ! sigma must be at least one
!                    eta_totM = 2 * eta_min
!                 else
!                    call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigmaM, phi,  &
!                         etasR(1:ieta), etasA(1:ieta), eta_totM, optim )
!                 endif
!
!              endif
!
!              if(eta_min <= eta_totP .and. eta_min <= eta_totM) then
!                 ! minimum was found
!                 !conver_sigma = conver_sigma
!                 sigma_stepP = .false.
!                 sigma_stepM = .false.
!
!              elseif (eta_totP < eta_min .and. eta_totP <= eta_totM) then
!                 sigma_opt = sigmaP
!                 eta_min = eta_totP
!                 sigma_stepM = .false.
!                 conver_sigma = conver_sigma + 1
!
!              elseif (eta_totM < eta_min) then
!                 sigma_opt = sigmaM
!                 eta_min = eta_totM
!                 sigma_stepP = .false.
!                 conver_sigma = conver_sigma + 1
!              else
!                 write(*,'(a8,6es16.8)') 'etas:',eta_min, eta_totP,  eta_totM
!                 stop "strange 33ju3 B"
!              endif
!              !!write(80+k, '(60es12.4)') sigma_opt, phi_opt, eta_min, q_sigma, D_phi
!
!              if(.not. sigma_stepP .and. .not. sigma_stepM) goto 15
!
!              if(sigma <= 1.) then
!                 sigma_opt = 1
!                 phi_opt = phi ! arbitrary
!                 goto 15
!              endif
!
!           enddo   ! do j=1, max_iter_one_dir
!15         continue
!
!           if(conver_phi == 0 .or. conver_sigma == 0) goto 20
!           !print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~', iter
!        enddo  !  iter = 1, max_iter
!20      continue
!
!        !!write(70+k, '(3es12.4, i5)') sigma_opt, phi_opt, eta_min, i_accur
!
!
!     enddo  ! do i_accur = 1, n_accur
!     !endif
!     params_init(k, 3) = sigma_opt
!     params_init(k, 4) = phi_opt
!     params_init(k, 5) = eta_min
!
!  enddo  ! k=1, ieta
!
!  ! seeking of the minimum from several candidates
!  eta_min = 1E+30
!  do k= ieta_min,ieta
!     if(params_init(k, 5) < eta_min) then
!        eta_min   = params_init(k, 5)
!        sigma_opt = params_init(k, 3)
!        phi_opt   = params_init(k, 4)
!     endif
!
!     ! ! not necessary for the computation, only the value at the initial gues
!     ! ! can be commented
!     ! call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, &
!     !      params_init(k, 1), params_init(k, 2),  &
!     !      etasR(1:ieta), eta_tot, optim )
!
!
!     ! write(*,'(a5, i3, 6es12.4, a2, 2es12.4)') &
!     !      'resu:', k, eta_tot, params_init(k, 1:5), '|', q_sigma, D_phi
!     !write(900+state%space%adapt%adapt_level, *)ie, params_init(k, 3:5)
!  enddo
!  !write(900+state%space%adapt%adapt_level, '(x)')
!
!  deallocate( params_init)
!  print*,'~~~~~~~~~~~ after SeekMinimumIteratively', iter, optim


end subroutine SeekMinimumIteratively


  !> seeking of the minimum by evaluating of all possible combinations
subroutine SeekMinimumGlobally(ieta, deg, anis, lambda, sigma_opt, phi_opt, &
            etasR, etasA, eta_min, diff_term)
  integer, intent(in) :: ieta    ! number of estimator
  integer, intent(in) :: deg     ! current polynomial degree
  real, dimension(1:4, 1:3), intent(in) :: anis  ! anisotropy of the function
  real, intent(inout) :: lambda, sigma_opt, phi_opt    ! anisotropy of the triangle
  real, dimension(1:ieta), intent(in) :: etasR  ! computed primal and dual residuas
  real, dimension(1:ieta), intent(inout) :: etasA  ! estimates of the particular "weights"
  real, intent(out) :: eta_min    ! error estimate
  logical, intent(in) :: diff_term ! True for problems with nonzero diffusion
  integer :: N_phi, N_sigma, i, j, optim
  real :: pi, phi, sigma, q_sigma, D_phi, eta_tot

  optim = 0

  q_sigma = 1.05 !1.005

  N_phi = 50 !200
  N_sigma = 20 !200

  sigma = 1.
  pi =  asin(1.0)*2.0

  D_phi = pi / N_phi

  eta_min = 1E+15

  do j=0,N_phi
     phi = j * D_phi

     do i=0, N_sigma
        sigma = q_sigma**i

        call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda, sigma, phi,  &
             etasR(1:ieta), etasA(1:ieta), eta_tot, optim, diff_term )

        write(100, '(60es12.4)') sigma, phi, eta_tot
        write(101, '(60es12.4)') sigma, phi+pi, eta_tot

        if(eta_tot < eta_min) then
           eta_min = eta_tot
           sigma_opt = sigma
           phi_opt = phi
        endif
     enddo

     !sigma = sigma * q_sigma

     write(100,'(x)')
     write(101,'(x)')
  enddo

end subroutine SeekMinimumGlobally

!> seeking of the anisotropy of the given function by its coeffs
subroutine FindAnisotropy(degP, ai, anis, upper )
  integer, intent(in) :: degP
  real, dimension(0:degP), intent(in) :: ai ! coeffs of interpolation functins
  real, dimension(1:3), intent(out) :: anis  ! components, A, \rho, \vp
  logical, intent(in) :: upper  ! really upper bound using expansion
  integer:: i, n, j
  real :: pi, phi, f, f_max, a_max, f_per, rho, area_min, a, f_der, min_f, area
  real :: f1, f2, f_bound, f_max_opt, rho_opt, scale
  real, dimension(:), allocatable :: xi
  real, dimension(:,:), allocatable :: D, R, RT, M

  n = 100
  pi =  asin(1.0)*2.0


  ! seeking of maximal derivatives
  f_max = -10.


  do i=0, n, 1  ! equidistance 2 degrees is enough
     phi =  pi * i / n   ! angle


     f = Eval_Dir_Derivative(degP, ai(0:degP), cos(phi), sin(phi))

     if( f  > f_max )then
        f_max = f
        a_max = phi

     endif
     !print*,'???', i, f, f_max
  enddo

  ! size in the perpendicular direction
  phi = a_max + pi /2

  ! derivative in the perpendicular direction
  f_per =  Eval_Dir_Derivative(degP, ai(0:degP), cos(phi), sin(phi))

  rho = f_max/max(1E-15, f_per)

  anis(1) = f_max
  anis(2) = rho
  anis(3) = a_max

  !print*,'a_max =', f_max, f_per


  if(upper) then
     ! setting of the upper estimate of type
     ! DirectDeriv M= max_f ( x^T Q D Q^T x )**(degP / 2)
     ! D = diag (1, rho**(-2/degP) )

     allocate(D(1:2, 1:2), R(1:2, 1:2), RT(1:2, 1:2), M(1:2, 1:2), xi(1:2) )

     R(1,1) = cos(a_max)
     R(1,2) = -sin(a_max)
     R(2,1) = sin(a_max)
     R(2,2) = cos(a_max)

     RT(1:2, 1:2) = transpose(R(1:2, 1:2) )

     D(2,1) = 0.
     D(1,2) = 0.

     f1 = f_max
     f2 = rho

     D(1,1) = 1.
     if(rho >0) then
        D(2,2) = rho ** (-2./degP)
     else
        D(2,2) = 0.
     endif

     M(1:2, 1:2) = matmul(R(1:2, 1:2), matmul(D(1:2, 1:2), RT(1:2, 1:2)))

     !area = (rho / f_max**2)**(-1./degP)
     !write(*,'(a6, 3i5,8es12.4)') &
     !     '???',elem%i, 0, 0, f_max, min_f, rho, area

     area_min = 1E+30

     ! seeking of optimal values f_max and rho satisfying the estimate
     ! for several values of f_max, we seek the minimal rho
     do j = 0, 10
        f_max = f1 * (1.01**j)
        !rho = f2
        rho = 1E+6

        if(rho >0) then
           D(2,2) = rho ** (-2./degP)
        else
           D(2,2) = 0.
        endif

        M(1:2, 1:2) = matmul(R(1:2, 1:2), matmul(D(1:2, 1:2), RT(1:2, 1:2)))

        do i=0,n
           !a = 2 * pi * i / n
           a =  pi * i / n
           xi(1) = cos(a)
           xi(2) = sin(a)

           ! directional derivative in direction a
           f_der = Eval_Dir_Derivative(degP, ai(0:degP), cos(a), sin(a) )

           !f  = M(1,1) * cos(a)**2 + ( M(1,2) + M(2,1) )*cos(a) *sin(a) &
           !     + M(2,2)*sin(a)**2
           f = dot_product(xi(1:2), matmul(M(1:2, 1:2), xi(1:2) ) )

           f_bound =  f_max * abs(f)**(degP/2.)

           if(f_der > 1.001 * f_bound ) then !.and. a /=a_max) then
              ! estimate is violated,
              ! we have to modify the anisotropy of the interpolation error function

              !! decrease of the ratio of the interpolation error function
              scale = ((f_der / f_max)**(2./degP) - (cos(a - a_max))**2) &
                   /(sin(a - a_max))**2
              scale = scale**(degP/2.)
              min_f = f_max * scale

              !! increase of the size of the interpolation error function
              !scale = f_der / f_bound
              !f_max = f_max * scale
              !min_f = min_f * scale

              ! refreshing of the matrix on the right-hand-side
              rho = f_max / min_f
              D(2,2) = rho ** (-2./degP)
              M(1:2, 1:2) = matmul(R(1:2, 1:2), matmul(D(1:2, 1:2), RT(1:2, 1:2)))

           endif
        enddo  ! i=1,n

        if(rho > 0. .and. f_max > 0.) then
           !write(*,'(a8,20es14.6)') '>><><><', &
           !     rho, f_max, min_f, f_max/ min_f, (rho / f_max**2)**(-1./degP), &
           !     f_max * min_f, (f_max * min_f)**(1./degP)
           area = (rho / f_max**2)**(-1./degP)
           !!area = (f_max * min_f)**(1./degP)
        else
           area = 1.
        endif

        !write(*,'(a6, 3i5,8es12.4)') &
        !     '???',elem%i, j, i, f_max, min_f, rho, area

        if(area < area_min) then
           area_min = area
           f_max_opt = f_max
           rho_opt = rho
        else
           ! we expect that the dependence of area on f_max is strictly convex
           ! so when are starts to increase, the optimal value was found
           goto 20
        endif

        !if(elem%i == itest)  then
        !   ifig = ifig + 1
        !   call DrawEstimate(ifig, elem%xc(1:2), f_max, degP, M(1:2, 1:2) )
        !endif

     enddo   ! j=1,5

     deallocate(D, R, RT, M, xi)

20   continue
       ! the optimal anisotropy of interpolation error function found
     !f_max = f_max_opt
     !  rho = rho_opt
     !  min_f = f_max / rho
     !  area = (rho / f_max**2)**(-1./degP)

     anis(1) = f_max_opt
     anis(2) = rho_opt
     anis(3) = a_max


     !write(*,'(a6, 3i5,8es12.4)') &
     !     '???',elem%i, -99, -99, f_max, min_f, rho, area
     !write(*,*) '(__________________________________)'

     !print*,'a_max =', anis(1:3)

  endif  ! if(upper)


end subroutine FindAnisotropy




! seeking of the anisotropy of the given function by its coeffs
subroutine  Draw_deriv(ifile, degP, ai)
  integer, intent(in) :: ifile, degP
  real, dimension(0:degP), intent(in) :: ai ! coeffs of interpolation functins
  integer:: i, j, n, m
  real :: pi, phi, f, f_max, a_max, f_per, a1, a2, ratio

  n = 100
  pi =  asin(1.0)*2.0

  m = 0

  ratio = 1.
  do j = 0, m

     do i=0, 2*n, 1  ! equidistance 2 degrees is enough
        phi =  pi * i / n   ! angle

        a1 = cos(phi) !* ratio
        a2 = sin(phi) !* ratio

        f = Eval_Dir_Derivative(degP, ai(0:degP), a1, a2)
        !write(ifile, *)  f * a1  , f* a2, f
        write(ifile, *)  f * cos(phi)  , f* sin(phi), f
     end do
     write(ifile, '(x)')

     ratio = ratio * 0.1

  enddo

end subroutine Draw_deriv

! seeking of the anisotropy of the given function by its coeffs
subroutine  Draw_anis(ifile, degP, anis)
  integer, intent(in) :: ifile, degP
  real, dimension(1:3), intent(in) :: anis ! coeffs of interpolation functins
  integer:: i, n
  real :: pi, phi, f, f_max, a_max, f_per, f1, rho, a1, a2, ratio


  a_max = anis(3)
  rho = anis(2)
  f_max = anis(1)

  n = 100
  pi =  asin(1.0)*2.0

  !ratio = 100.

  !print*,'###', f_max, rho, a_max, degP
  ! estimates for each direction

  do i=0, 2*n, 1  ! equidistance 2 degrees is enough
     phi =  pi * i / n   ! angle

     a1 = cos(phi-a_max) !* ratio
     a2 = sin(phi-a_max) !* ratio
     f = f_max * (a1**2 + rho**(-2./degP) * a2**2 )**(1.*degP/2)

     f1 = f_max * (a1**2 + rho**(-2.) * a2**2 )**(1./2)
     write(ifile, *) f * cos(phi) , f*sin(phi), f, f1
  end do

end subroutine Draw_anis

  !> compute the directional derivative in the direction ( cos(a), sin(a) )
  function Eval_Dir_Derivative(degP, der, a1, a2)
    real :: Eval_Dir_Derivative
    integer, intent(in) :: degP
    real, dimension(0:degP), intent(in) :: der  ! partial derivatives of degree degP
    real, intent(in) :: a1, a2
    real :: f, fac
    integer :: j

    f = 0.

    do j=0, degP
       f = f + a2**j * a1**(degP - j) * der(j)
    enddo

    Eval_Dir_Derivative =  max( abs(f), 1E-50)   ! limitation for minimal derivatives

  end function Eval_Dir_Derivative

  subroutine Eval_products( deg, ai, phi, prod)
    integer, intent(in) :: deg
    real, dimension(1:deg), intent(in) :: ai
    real, intent(in) :: phi
    real, intent(out) :: prod
    integer :: i, j

    prod = 1.
    do j=1,deg

       prod = prod * ( cos(phi) + ai(j) * sin(phi) )
    enddo

  end subroutine Eval_products


  ! function factor(n)
  !   integer :: factor
  !   integer, intent(in) :: n
  !   integer :: i, fac

  !   fac = 1.
  !   do i = 2, n
  !      fac = fac * i
  !   enddo
  !   factor = fac

  ! end function factor



  subroutine Set_Coeffs_Function(degP, Du, ai)
    integer, intent(in) :: degP   ! p+1
    real, dimension(0:degP), intent(in) :: Du  ! degP-th order derivatives
    real, dimension(0:degP), intent(out) :: ai ! coeffs of the degP-function
    integer :: i

    ! evaluation of coeffs of interpolation error
    do i=0, degP
       ai(i) =  1./factorial(i)/factorial(degP -i) * Du(i)
       !print*,'e33',factorial(i), factorial(degP -i), ai(i)
    enddo
  end subroutine Set_Coeffs_Function



  subroutine Set_Coeffs_Gradient(degP, degQ, Du, mA,  ai)
    integer, intent(in) :: degP   ! p+1
    integer, intent(in) :: degQ   ! 2p
    real, dimension(0:degP), intent(in) :: Du  ! degP-th order derivatives
    real, dimension(1:2, 1:2), intent(in) :: mA  ! matrix coeffs
    real, dimension(0:degQ), intent(out) :: ai ! coeffs of the degQ-function
    real, dimension(:), allocatable ::  Dw
    real, dimension(:,  :), allocatable ::  beta
    integer :: i, j, deg

    deg = degP - 1

    if( 2*deg /= degQ) then
       print*,' Inconsistency in Set_Coeffs_Gradient:', degP, degQ, deg
       stop
    endif

    allocate(beta(1:2, 0:deg), Dw(0:degP) )

    Dw = 0.
    ! evaluation of coeffs of interpolation error
    do i=0, degP
       Dw(i) =  1./factorial(i)/factorial(degP -i) * Du(i)
       !print*,'e44',factorial(i), factorial(degP -i), Dw(i)
    enddo

    beta(:,:) = 0.
    ! evaluation of coeffs of gradient \mA \nabla e_I
    do i=0, deg
       beta(1, i) = mA(1,1) * ( i+1) * Dw(i+1) + mA(1,2) *(degP -i) * Dw(i)
       beta(2, i) = mA(2,1) * ( i+1) * Dw(i+1) + mA(2,2) *(degP -i) * Dw(i)
    enddo

    ! \evalation of square of magnitude of  gradient \mA \nabla e_I
    ai(:) = 0
    do i=0, deg

       do j=0, i

          ai(i) = ai(i) + beta(1, j)*beta(1, i-j) + beta(2, j)*beta(2, i-j)

          if(degQ-i > i) &
               ai(degQ-i) = ai(degQ-i) + beta(1, deg-j) *beta(1, deg -i+j)  &
               +  beta(2, deg-j) *beta(2, deg-i+j)

          !print*,'i,j ===', i,  j, i-j, deg-j, deg-i+j

          !write(*,'(a10, 20i5)') '####'
          !write(*,'(a10, 20i5)') '####',deg
          !write(*,'(a10, 20i5)') '####',deg, degP
          !write(*,'(a10, 20i5)') '####',deg, degP,degQ
          !write(*,'(a10, 20i5)') '####',deg, degP,degQ, i
          !write(*,'(a10, 20i5)') '####',deg, degP,degQ, i,j, degQ-i, i-j, deg-j, deg-i+j
       enddo
    enddo

    ! do i=0, -deg

    !    do j=0, i

    !       !ai(i) = ai(i) + beta(1, j)*beta(1, i-j) + beta(2, j)*beta(2, i-j)

    !       if(degQ-i > i) &
    !            ai(degQ-i) = ai(degQ-i) + beta(1, deg-j) *beta(1, deg -i+j)  &
    !            +  beta(2, deg-j) *beta(2, deg-i+j)

    !       print*,'i,j ===', degQ-i, j, i-j, deg-j, deg-i+j

    !       !write(*,'(a10, 20i5)') '####'
    !       !write(*,'(a10, 20i5)') '####',deg
    !       !write(*,'(a10, 20i5)') '####',deg, degP
    !       !write(*,'(a10, 20i5)') '####',deg, degP,degQ
    !       !write(*,'(a10, 20i5)') '####',deg, degP,degQ, i
    !       !write(*,'(a10, 20i5)') '####',deg, degP,degQ, i,j, degQ-i, i-j, deg-j, deg-i+j
    !    enddo
    ! enddo

    !do i=0, degQ
    !   print*,'e77',i,  ai(i)
    !enddo


    deallocate(beta)

    deallocate(dw)

  end subroutine Set_Coeffs_Gradient


  !> seeks the initial values for the optimalization as the local minimas
  !> itype = 1   volume estimate of E_int
  !> itype = 2   boundary estimate of E_int
  !> itype = 3   face  estimate of A \nabla E_int
  subroutine  Find_Params_init(ieta, deg, anis, params_init)
    integer, intent(in) :: ieta    ! number of estimator
    integer, intent(in) :: deg     ! current polynomial degree
    real, dimension(1:4, 1:3), intent(in) :: anis  ! anisotropy of the function
    real, dimension(1:ieta, 1:2), intent(inout) :: params_init
!    real, dimension(1:ieta, 1:ndim), intent(inout) :: sigma_init
!    real, dimension(1:ieta, 1:ndim), intent(inout) :: phi_init
    real :: lambda, sigma, phi    ! anisotropy of the triangle
    real :: etasA    ! error estimate
    integer :: optim  ! <0 => returns the optimal values
    integer :: k

    etasA = 1.
    optim = -1
    lambda = 1

    ! TODO: FR_ANI sigma and phi are not initialized here !!!

    ! TODO: control the ordering, now it should be same as in Anis_estim_tot
!    do k = 1, ndim
!      ! primal solution
!      call Anis_estim(1, deg, anis(k,1, 1:3), lambda, sigma, phi,  etasA, optim )
!      sigma_init(4,k) = sigma
!      phi_init(4,k) = phi
!      ! boundary
!      call Anis_estim(2, deg, anis(1, 1:3), lambda, sigma, phi,  etasA, optim )
!      sigma_init(5,k) = sigma
!      phi_init(5,k) = phi
!    end do

    ! primal solution
    ! volume
    call Anis_estim(1, deg, anis(1, 1:3), lambda, sigma, phi,  etasA, optim )
    params_init(2,1) = sigma
    params_init(2,2) = phi
    ! boundary
    call Anis_estim(2, deg, anis(1, 1:3), lambda, sigma, phi,  etasA, optim )
    params_init(4,1) = sigma
    params_init(4,2) = phi
    ! boundary gradient
    call Anis_estim(3, deg, anis(2, 1:3), lambda, sigma, phi,  etasA, optim  )
    params_init(6,1) = sigma
    params_init(6,2) = phi


    ! dual solution
    call Anis_estim(1, deg, anis(3, 1:3), lambda, sigma, phi,  etasA, optim  )
    params_init(1,1) = sigma
    params_init(1,2) = phi

    call Anis_estim(2, deg, anis(3, 1:3), lambda, sigma, phi,  etasA, optim  )
    params_init(3,1) = sigma
    params_init(3,2) = phi

    call Anis_estim(3, deg, anis(4, 1:3), lambda, sigma, phi,  etasA, optim  )
    params_init(5,1) = sigma
    params_init(5,2) = phi



  end subroutine Find_Params_init

!  !> seeks the initial values for the optimalization as the local minimas
!  !> itype = 1   volume estimate of E_int
!  !> itype = 2   boundary estimate of E_int
!  !> itype = 3   face  estimate of A \nabla E_int
!  subroutine  Find_init_params_new(ieta, deg, anis, sigma_init, phi_init)
!    integer, intent(in) :: ieta    ! number of estimator
!    integer, intent(in) :: deg     ! current polynomial degree
!    real, dimension(1:ndim, 1:4, 1:3), intent(in) :: anis  ! anisotropy of the function
!    real, dimension(1:ieta), intent(out) :: sigma_init, phi_init
!!    real, dimension(1:ieta, 1:ndim), intent(inout) :: sigma_init
!!    real, dimension(1:ieta, 1:ndim), intent(inout) :: phi_init
!!    real :: lambda, sigma, phi    ! anisotropy of the triangle
!!    real :: etasA    ! error estimate
!!    integer :: optim  ! <0 => returns the optimal values
!    integer :: k
!
!
!    phi = anis(3) + pi/2
!    if(phi > pi) phi = phi - pi
!    !sigma =  anis(2)**(1./degP)
!    ! here sigma is the square root of sigma from APNUM 14 !!!!x
!    sigma =  sqrt( anis(2)**(1./degP) )
!
!
!  end subroutine Find_init_params_new

  !> evaluate the total the anisotropic error estimate
  !> itype = 1   volume estimate of E_int
  !> itype = 2   boundary estimate of E_int
  !> itype = 3   face  estimate of A \nabla E_int
  subroutine Anis_estim_tot(ieta, deg, anis, lambda, sigma, phi, etasR, etasA, &
                            eta_tot, optim, diff_term )
    integer, intent(in) :: ieta    ! number of estimator (always 6?)
    integer, intent(in) :: deg     ! current polynomial degree
    real, dimension(1:ndim, 1:4, 1:3), intent(in) :: anis  ! anisotropy of the function
    real, intent(in) ::    lambda   ! size of the triangle
    real, intent(inout) :: sigma, phi    ! anisotropy of the triangle
    real, dimension(1:ieta, 1:ndim), intent(in)    :: etasR  ! computed primal and dual residuas
    real, dimension(1:ieta, 1:ndim), intent(inout) :: etasA  ! estimates of the particular "weights"
    real, intent(inout) :: eta_tot    ! error estimate
    integer, intent(inout) :: optim  ! <0 => returns the optimal values
    logical, intent(in) :: diff_term  ! true for problems with nonzero diffusion
    integer :: k

    etasA(:,:) = 0.0

    ! FR_ANI what is this for?
    if(optim >= 0) optim = optim+1

    ! NEW - anis(1:ndim, 1:4, 1:3)
    do k = 1,ndim
      ! primal solution
      ! volume anizotropy of w is the weight of the dual indicator -> 4
      call Anis_estim(1, deg, anis(k, 1, 1:3), lambda, sigma, phi,  etasA(4,k), optim )
      call Anis_estim(2, deg, anis(k, 1, 1:3), lambda, sigma, phi,  etasA(5,k), optim )

      ! dual solution
      call Anis_estim(1, deg, anis(k,3, 1:3), lambda, sigma, phi,  etasA(1,k), optim  )
      call Anis_estim(2, deg, anis(k,3, 1:3), lambda, sigma, phi,  etasA(2,k), optim  )

      if(diff_term) then
         ! 2 - grad(w) -> 6 part od the estimate
         call Anis_estim(3, deg, anis(k, 2, 1:3), lambda, sigma, phi,  etasA(6,k), optim  )
         ! 4 - grad(z) -> 3rd part of the estimate
         call Anis_estim(3, deg, anis(k, 4, 1:3), lambda, sigma, phi,  etasA(3,k), optim  )
      end if
    end do

    ! eta_tot = \sum_{k=1}^ndim \sum_{i=1}^6 etasR(i,k)*etasA(i,k)
    eta_tot = 0.0
    do k = 1, ndim
      eta_tot = eta_tot + dot_product(etasR(1:ieta,k), etasA(1:ieta,k) )
    end do

    ! TODO: this has to be changed when p+1 is used for both primal and dual
    ! global p+1 problem computed only for the dual solution  MI6
    if (state%space%estim_space == 'DWR' .and. DWR%deg_plus) then
       eta_tot = 0.0
      do k = 1, ndim
        eta_tot = eta_tot + 2*dot_product(etasR(1:3,k), etasA(1:3,k) )
      end do
    endif

  end subroutine Anis_estim_tot


  !> evaluate the current type of the anisotropic error estimate
  !> itype = 1   volume estimate of E_int
  !> itype = 2   boundary estimate of E_int
  !> itype = 3   face  estimate of A \nabla E_int
  ! TODO: FR_ANI why is sigma, phi also OUT parameter, seems that Find_params_init uses it that way???
  subroutine Anis_estim(itype, deg, anis, lambda, sigma, phi,  eta, optim)
    integer, intent(in) :: itype   ! type of error estimate
    integer, intent(in) :: deg     ! current polynomial degree
    real, dimension(1:3), intent(in) :: anis  ! anisotropy of the function
    real, intent(in) ::    lambda   ! size of the triangle
    real, intent(inout) :: sigma, phi    ! anisotropy of the triangle
    real, intent(out) :: eta    ! error estimate
    integer, intent(in) :: optim  ! <0 => returns the optimal values

    real, dimension(:,:), allocatable :: mG
    integer :: degP, degQ
    real :: fac, tau, integ, rho_pow, pi

    pi =  asin(1.0)*2.0

    ! RHS member in (32a-c) without G
    select case (itype)
    case(1)   ! volume estimate of E_int
       degP = deg + 1
       degQ = degP
       fac = anis(1)**2 * lambda**(2*deg+4)/(2*deg+4)

    case(2)   ! boundary estimate of E_int
       degP = deg + 1
       degQ = degP
       fac = anis(1)**2 * lambda**(2*deg+3) * sigma

    case(3)   ! boundary estimate of A \nabla E_int
       degP = deg
       degQ = 2*deg
       fac = anis(1) * lambda**(2*deg+1)* sigma

    case default
       stop 'UNKNOWN itype in Anis_estim'
    end select

    !print*,'#D#E#', deg, degP, degQ, 2*deg+4, fac

    rho_pow = anis(2)**(-2./degQ)  ! degQ
    !! ALT
    !rho_pow = anis(2)**(-2.)

    ! coefficients of the estmate (under (33) in ESCO article)
    allocate(mG(1:2, 1:2) )

    if( optim >= 0 .or. optim == -10) then
       tau = phi - anis(3)
       mG(1,1) = sigma*sigma*(cos(tau)**2 + rho_pow*sin(tau)**2)
       mG(1,2) = -sin(tau)*cos(tau)*(1- rho_pow )
       mG(2,1) = mG(1,2)
       mG(2,2) =  1./sigma/sigma*(sin(tau)**2 + rho_pow*cos(tau)**2)

       !print*,'----------------------------------'
       !write(*,'(a12, 3es12.4)') 'mG(1,:):',   mG(1,:)
       !write(*,'(a12, 3es12.4)') 'mG(2,:):',   mG(2,:)
       !print*,'----------------------------------'
       ! integral (33) in Lemma 4.6
       call IntegrateCircle(degP, mG(1:2, 1:2), integ)

       eta = fac* integ

       !write(*,'(a8,4i5, 20es12.4)') &
       !     '#D#E#', deg, degP, degQ, 2*deg+4, fac, integ, eta, sqrt(eta)


       !if(itype == 3 .and. optim == -10)  then
       !   print*
       !   write(*,'(a8,2i5, 30es12.4)')'87yu3s:',itype, deg, fac, anis(1), lambda**(deg+1), integ, rho_pow, eta
       !endif
    ! FR_ANI what is done here ???
    elseif(optim < 0 ) then
       phi = anis(3) + pi/2
       if(phi > pi) phi = phi - pi

       !sigma =  anis(2)**(1./degP)
       ! here sigma is the square root of sigma from APNUM 14 !!!!x
       sigma =  sqrt( anis(2)**(1./degP) )

       !print*,'###',phi, sigma
       !tau = phi - anis(3)
       !mG(1,1) = sigma*(cos(tau)**2 + rho_pow*sin(tau)**2)
       !mG(1,2) = -sin(tau)*cos(tau)*(1- rho_pow )
       !mG(2,1) = mG(1,2)
       !mG(2,2) =  1./sigma*(sin(tau)**2 + rho_pow*cos(tau)**2)

       !write(*,'(a12, 3es12.4)') 'mG(1,:):',   mG(1,:)
       !write(*,'(a12, 3es12.4)') 'mG(2,:):',   mG(2,:)
       !print*,'M11:', 1./sigma
       !stop "7833333333333333333"

       mG= 0.
       mG(1,1) = 1./sigma/sigma
       mG(2,2) = 1./sigma/sigma


       call IntegrateCircle(degP, mG(1:2, 1:2), integ)

       eta = fac* integ

       !write(200+optim, '(60es12.4)') sigma, phi, eta, 10., 1.


    endif

    ! VD
    eta = sqrt(eta)

    deallocate( mG)



  end subroutine Anis_estim

  !> integration of  (x mG x)^deg over circle
  subroutine IntegrateCircle(degP, mG, integ)
    !use G_integrate
    type(Gauss_rule), pointer  :: G_rule
    integer, intent(in) :: degP
    real, dimension(1:2, 1:2), intent(in) :: mG
    real, intent(inout) :: integ
    integer :: i, j, N
    real :: t1, t2, t, dt, pi, r1, f

    pi =  asin(1.0)*2.0

    G_rule => state%space%G_rule(7)  ! SHOULD BE ENOUGH (?)

    N = 10 !20

    dt = 2*pi / N
    t1 = 0.
    !t2 = dt

    integ = 0.

    do i=1, N
       r1 = 0.
       t1 = (i-1)*dt

       do j=1,G_rule%Qdof
          t = t1 + dt*G_rule%lambda(j)
          f = mG(1,1)*cos(t)**2 + 2*mG(1,2)*cos(t)*sin(t) + mG(2,2)*sin(t)**2

          f = f**degP
          !ALT
          ! f = f

          r1 = r1 + f*G_rule%weights(j)
          !print*,'d4343',j, f  !G_rule%weights(j)
          !stop
       enddo
       integ = integ + r1*dt
       !print*,'d422', integ, r1, dt

       !t1 = t1 + dt
       !t2 = t2+dt
    enddo

    !write(*,'(a8, 20es22.14)') 'integ =', integ

  end subroutine IntegrateCircle



  !> evaluate the derivatives of the "highest order" of the primal and dual solutions"
  !> input : primal solution in elem%wST_LS(1:ndim, 1:dofP , j), j=0,1,2; dofP = DOFtriang(elem%deg+j)
  !> input : dual   solution in elem%zST_LS(1:ndim, 1:dofP , j), j=0,1,2; dofP = DOFtriang(elem%deg+j)
  !> output elem%wSS(1: 2*ndim, j, 0:degP) degP = elem%deg + j,  j=0,1,2
  subroutine Eval_High_Order_Deriv_Elem(elem)
    class(element), target, intent(inout) :: elem
    type(volume_rule), pointer :: V_rule
    real, dimension(:,:,:,:), allocatable :: Dw, Dwi  !  working array
    real, dimension(:,:), pointer :: phi
    real, dimension(:,:,:), allocatable :: Dphi
    real, dimension(:), allocatable :: vec
    real, dimension(:,:), allocatable :: Mass, MassInv
    integer :: i,j, k, l, ideg, ndimD, ix, jx
    integer :: degM, dofM, degP, dofP, Qdof, Qnum, dof_sm

    ndimD = 2*ndim

    degM = elem%deg + 2
    dofM = DOFtriang( degM)

    if( allocated(elem%wSS) ) deallocate (elem%wSS)
    allocate( elem%wSS(1: ndimD, 0:2, 0:degM), source = 0.0 )

    ! quadrature rule
    Qnum = state%space%Qdeg(elem%deg+1, 1)
    if(elem%deg == MaxDegreeImplemented ) Qnum = elem%Qnum

    V_rule => state%space%V_rule(Qnum)
    Qdof = V_rule%Qdof

    allocate(Dphi(1:dofM, 1:nbDim, 1:Qdof ) )
    call  Eval_Dphi_plus(elem, V_rule, dofM,  Dphi(1:dofM, 1:nbDim, 1:Qdof) )

    phi => V_rule%phi(1:dofM, 1:Qdof)

    ! temporary array for evaluation of derivatives
    allocate( Dwi(1:ndimD, 0:degM, 0:degM, 1:Qdof), source = 0.0 )
    allocate( Dw (1:ndimD, 0:degM, 0:degM, 1:dofM), source = 0.0 )


    ! matrices and vector for a projection
    allocate(vec(1:dofM), MassInv(1:dofM, 1:dofM) )
    allocate(Mass(1:dofM, 1:dofM), source = 0.0)

    ! Mass matrix
    call IntegrateBlockBBmass(elem, Qnum, dofM, dofM, Mass(1:dofM,1:dofM) )


    !write(99,*)'?????????????????????????????????????'
    !do l=1, dofM
    !   write(99,'(a8, i5,a2,40es12.4)') 'Mx',l,'|', Mass(l, 1:dofM)
    !enddo
    !write(99,*)'____________________________________________________________'



    do ideg=0, 2  ! derivatives of the reconstruction of degree elem%deg+j
       degP = elem%deg + ideg
       dofP = DOFtriang( degP)
!       print*, 'ideg === ', ideg, dofP, state%p_mod_max
!       print*, 'wST === ', elem%wST_LS(1:ndim, 1:dofP, ideg)

       if(degP > MaxDegreeImplemented) then
          ! this degree is not implemented, hence too high derivative
          elem%wSS(1:ndimD, ideg, 0:degP) = 1.E+60

       else

          ! inversion of the mass matrix
          if(dofP > dofM) stop "dimension false in 39u39j3i  ama-hp_DWR.f90"
          MassInv(1:dofP, 1:dofP) = Mass(1:dofP, 1:dofP)
          call MblockInverse(dofP, MassInv(1:dofP, 1:dofP) )

          Dw = 0.
          Dwi = 0.

          ! setting of the reconstructed solution
          Dw(     1:ndim,  0, 0, 1:dofP) = elem%wST_LS(1:ndim, 1:dofP, ideg)
          Dw(ndim+1:ndimD, 0, 0, 1:dofP) = elem%zST_LS(1:ndim, 1:dofP, ideg)

          ! we nullify the lower degree terms, its degP-th derivative = 0
          !dof_sm = DOFtriang( degP - 1)
          !Dw(     1:ndimD,  0, 0, 1:dof_sm ) = 0.

          !if(elem%i == 727 ) write(*,'(a8, 40es12.4)') 'wsTplus', Dw(1:ndim,  0, 0, 1:dofP)

          do k=1,ndimD

             do i=1,degP   ! i = i-th derivative
                do j=0,i  !   d^i w /dx^j dy^(i-j)

                   ! evaluation of the partial derivative in integ nodes
                   do l=1, Qdof  ! integ node

                      ! jx which (i-1)th derivative has to be integrate
                      ! ix according which variable: ix = 1 according x_1, ix = 2 according x_2
                      if(j==0) then
                         jx = 0
                         ix = 1
                      elseif(j==1 .and. degP >2) then
                         jx = 0
                         ix = 2
                      else
                         jx = j-1
                         ix = 2
                      endif

                      Dwi(k, i, j, l) = dot_product(Dw(k,i-1, jx, 1:dofP), Dphi(1:dofP, ix, l) )

                      ! if(ideg == 1 .and. k==1 .and. l <= 3) then
                      !    write(99,'(a5, 5i5, es12.4,a1,40es12.4)') 'Dw  :',i,j,l, ix, jx,  Dwi(k, i, j, l),'|',Dw(k,i-1, jx, 1:dofP)
                      !    write(99,'(a5, 5i5, es12.4,a1,40es12.4)') 'Dphi:',i,j,l, ix, jx,  Dwi(k, i, j, l),'|',Dphi(1:dofP, ix, l)
                      ! end if


                   enddo ! l

                   ! if(ideg >= 1 .and. k==1) then
                   !    if(i == 1 .and. j== 0 ) write(99,*)'##################################################'
                   !    if(j == 0) write(99,*)'----------------------'
                   !    write(99,'(x)')
                   !    write(99,'(a5, 2i5, 40es12.4)') 'derf:',i,j,  Dwi(k, i, j, :)
                   !    !write(99, '(x)')
                   ! end if

                   ! evaluation of the basis coefficients
                   if(state%space%V_rule(Qnum)%Qdof /= Qdof) stop "trouble s3ye83hi3"
                   call IntegrateVectorBplus(elem, Qnum, dofP, Dwi(k, i, j, 1:Qdof), vec(1:dofP) )

                   ! if(ideg >= 0 .and. k==1) then
                   !    write(99,*)'?????????????????????????????????????'
                   !    do l=1, dofP
                   !       write(99,'(a8, i5, es12.4,a2,40es12.4)') 'Mx',l, vec(l),'|', MassInv(l, 1:dofP)
                   !    enddo
                   !    write(99,*)'?????????????????????????????????????'
                   ! endif

                   Dw(k, i,j, 1:dofP) = matmul(MassInv(1:dofP,1:dofP), vec(1:dofP) )

                   ! if(ideg >= 0 .and. k==1) then
                   !    write(99,'(a5, 2i5, 40es12.4)') 'DW##:',i,j,  Dw(k, i, j, :)
                   !    write(99,'(a5, 2i5, 40es12.4)') 'DW00:',i,j,  Dw(k, 0, 0, :)
                   !    write(99,'(a5, 2i5, 40es12.4)') 'DW10:',i,j,  Dw(k, 1, 0, :)
                   !    write(99,'(a5, 2i5, 40es12.4)') 'DW11:',i,j,  Dw(k, 1, 1, :)
                   !    write(99, '(x)')
                   ! end if

                   !MassS(:,:) = Mass(:,:)
                   !call SolveLocalMatrixProblem(dofP, MassS(1:dofP, 1:dofP), 1, vec(1:dofP) )
                   !Dw(k, i,j, 1:dofP) =  vec(1:dofP)
                enddo ! j
             enddo ! i
          enddo ! k

          if(elem%i == -727) then
             !if(elem%xc(1) > 0.95 .and. elem%xc(2) > 0.95) then
             !if( abs(elem%xc(1) -20 ) < 0.5 ) then
             print*,'####  degree =', degP,  dofP,"   ######################"
             write(*,*) '-----', elem%i, '  ----', elem%xc
             k = 1
             do i=1,degP   ! i = i-th derivative
                do j=0,i  !   d^i w /dx^j dy^(i-j)
                   write(*,'(a6,3i5,500es12.4)') 'Dwi:',i,j,Qdof,Dwi(k, i, j, 1:Qdof)
                enddo
             enddo

             write(*,*) '-----------------------------------'

             do i=0,degP   ! i = i-th derivative
                do j=0,i  !   d^i w /dx^j dy^(i-j)
                   write(*,'(a6,3i5,500es12.4)') 'Dw:',i,j,dofP,Dw(k, i, j, 1:dofP)
                enddo
             enddo

             write(*,*) 'DF ########################################'
             !stop 'u39j3o'
          endif

          do k=1, ndimD
             do i=0, degP
                ! elem%wSS contains only the derivatives of maximal order
                !Dwi(k,degP, i, 1:Qdof) is almost constant
                elem%wSS(k, ideg, i) = sum(Dwi(k,degP, i, 1:Qdof) )/Qdof !/ phi(1,1)

                !if( sqrt( (elem%xc(1) - 0.5)**2 + (elem%xc(2) -0.5)**2 ) < -0.05 ) &
                !write(*,'(2es12.4,4i5, 2es12.4,a2, 120es12.4)') &
                !     elem%xc(:), elem%i, n, degP, i, &
                !     Dwi(k,degP, i, 1:2) /  phi(1,1),'|',elem%wSS(k, n, i)

             enddo
             ! if( sqrt( (elem%xc(1) - 0.5)**2 + (elem%xc(2) -0.5)**2 ) < 0.05 ) print*
             !write(100 + 10*state%space%adapt%adapt_level + ideg,'(120es12.4)') &
             !if( elem%xc(1) > 19 .and. elem%xc(1) < 21) &
             !if(elem%i == -599) &
             ! if( abs(elem%xc(1) -0.6) < 0.04 .and. abs(elem%xc(2) -0.05) < 0.05) then
             ! !if(k == 1) &
             !    if(state%space%adapt%adapt_level >= 4) then
             !       write(299,'(2es12.4,3i5, a8, 120es12.4)') &
             !         elem%xc(:), elem%i, ideg, degP, &
             !         '| w: |', Dw(k, 0,0, 0: dofP) !*phi(1, 1)

             !       write(299,'(2es12.4,3i5, a8, 120es12.4)') &
             !         elem%xc(:), elem%i, ideg, degP, &
             !         '| nd |', elem%wSS(k, ideg, 0: degP) !*phi(1, 1)
             !    endif

             ! endif

          enddo  ! k=1, ndimD

       endif  ! if(degP > MaxDegreeImplemented


    enddo  ! n=0,2

    deallocate(Dw, Dwi,  vec, Dphi)
    deallocate( Mass, MassInv)


  end subroutine Eval_High_Order_Deriv_Elem


  !> smoothing of the anisotropy, i.e., of size (lambda) and aspect ratio (sigma)
  !> stored in anisot(:, 1:2)
  subroutine SmoothAnisotropy(grid, nelem,  anisotropy )
    class( mesh ), intent(inout) :: grid
    integer, intent(in) :: nelem
    real, dimension(1: nelem, 1:4), intent(inout) :: anisotropy
    class(element), pointer :: elem
    real, dimension(:,:), allocatable :: V_anisotropy
    real :: weight, rlen
    integer :: i, j, k, l, ipoc, ival,  npoin
    !real :: sigma_limit_max = 50. ! 50. correspondance with AMA%pos: sigma <= 1/(1.5*pos)
    real :: sigma_limit_max = 500. ! 50. correspondance with AMA%pos: sigma <= 1/(1.5*pos)

    if(nelem /= grid%nelem) stop "troubles in variable nelem in SmoothAnisotropy"
    npoin = grid%npoin

    ! limitation of sigma
    do i=1, nelem
       if(anisotropy(i, 2) > sigma_limit_max) &
            write(*,'(a10,2es12.4)') 'LIMIT:', anisotropy(i, 2) , sigma_limit_max

       anisotropy(i, 2) = min(anisotropy(i, 2) , sigma_limit_max)
    enddo


    !components for limitations: lambda (1), sigma (2), degP (4)
    ival = 4

    ! number of smoothing cycles
    ipoc = 1

    !do i=1, 5
    !   write(*,'(a8, 30es12.4)') 'anisot' , anisotropy(i, 1:3)
    !enddo


    print*, "Anisotropy smoothing: ipoc = ", ipoc

    allocate( V_anisotropy(1:npoin, 0:4) )

    do l=1, ipoc
       V_anisotropy(:,:) = 0.

       ! interpolation from elements to vertices
       do i=1,nelem
          elem => grid%elem(i)

          do j = 1, elem%flen
             k = elem%face(idx, j)

             ! distance between the vertex and the barycenter
             rlen = sqrt(dot_product( elem%xc(:) - grid%x(k,:), elem%xc(:) - grid%x(k,:) ) )

             V_anisotropy(k, 1:ival) = V_anisotropy(k, 1:ival) + anisotropy(i, 1:ival) * elem%area !/rlen
             V_anisotropy(k, 0) = V_anisotropy(k,0) + 1. * elem%area  !/rlen
          enddo
       enddo

       ! weighting
       do j=1, ival
          V_anisotropy(1:npoin, j) = V_anisotropy(1:npoin, j)  / V_anisotropy(1:npoin, 0)
       enddo

       ! backward interpolation into elements
       do i=1,grid%nelem
          elem => grid%elem(i)

          ! only lambda, sigma and degP
          anisotropy(i, 1 ) = sum( V_anisotropy(elem%face(idx,:) ,1 ) ) /  elem%flen
          anisotropy(i, 2 ) = sum( V_anisotropy(elem%face(idx,:) ,2 ) ) /  elem%flen
          anisotropy(i, 4 ) = sum( V_anisotropy(elem%face(idx,:) ,4 ) ) /  elem%flen

       enddo
    enddo  ! do l=1, ipoc

    !print*
    !do i=1, 5
    !   write(*,'(a8, 30es12.4)') 'anisotropy' , anisotropy(i, 1:3)
    !enddo
    !stop "eu33i"

    ! setting of the metric matrices
    do i=1,grid%nelem
       elem => grid%elem(i)
       call Set_DWR_metric(elem, 1, anisotropy(i,1), anisotropy(i,2), anisotropy(i,3), &
            anisotropy(i,4))

       ! for UNIFORM refinement
       !call Set_DWR_metric(elem, 1, 4E-2 , 1., 0. )
    enddo

    deallocate( V_anisotropy )


  end subroutine SmoothAnisotropy

  !> test of anisotropic estimates and their optimality
  subroutine TEST_ANIS_ESTIM( )
    class(element), pointer :: elemK
    integer, dimension(:), allocatable :: iloc
    real, dimension(:), allocatable :: Du, Dz, ai, etasR, etasA
    real, dimension(:,:), allocatable :: x, xx, wi, ww, xi, anis
    real, dimension(:,:), allocatable :: Q, M, D, err, mA, r_ieff
    real :: phi, psi, sigma, area, lambda, rl1, rl2 , pi, val, val1, etas, vp
    integer :: deg, dof, degP, dofP, dofPP, Qdof, degQ
    type(volume_rule), pointer :: V_rule
    integer :: i,j,k, k1, ieta, optim, is, ip, Ns, Np, kk
    real :: Dp, Ds
    logical ::  upper = .true.

    print*
    print*
    print*,'###########################################################'
    print*

    pi =  asin(1.0)*2.0


    deg = 1

    degP = deg + 1
    degQ = 2*deg

    ! unit matrix
    allocate(mA( 1:2, 1:2), source = 0.0 ) ! diffusion matrix TODO real diffusion
    mA(1,1) = 1.;  mA(2,2) = 1.


    allocate(iloc(1:3) )
    !iloc(1)= 1; iloc(2)= 2; iloc(3)= 3

    allocate(x(1:3, 1:2) , xx(1:3, 1:2) )

    allocate(err(1:3, 1:3), source = 0.0 ) ! real error

    allocate(r_ieff(1:3, 1:2) )  ! maximal and minimal effectivity index

    ! reference triangle
    xx(1, 1:2) = (/0, 1 /)
    xx(2, 1:2) = (/-sqrt(3.)/2 , -0.5 /)
    xx(3, 1:2) = (/ sqrt(3.)/2, -0.5 /)

    Np = 20  ! number of phi spliting
    Ns = 30   ! number of sigma spliting

    Dp =  2*pi / Np
    Ds = 1.2

    do is =  0, Ns
       r_ieff(:, 1) = 0.    !maximal effectivity index
       r_ieff(:, 2) = 1000. !minimal

       do ip = 0,  Np

          err = 0.

          area = 1.
          psi = 0.

          phi = ip * Dp
          sigma = Ds**is

          rl1 = area * sigma
          rl2 = area/ sigma
          lambda = sqrt( rl1 * rl2)

          allocate(elemK)
          call elemK%init(0, 3, 3, iloc(1:3) )


          call set_coordinates(psi, phi, rl1, rl2, xx(1:3, 1:2), x(1:3, 1:2) )

          allocate(elemK%F)
          elemK%type = 3; elemK%flen = 3;

          !! geometry
          elemK%F%deg = 1
          elemK%F%dof = 3
          allocate(elemK%F%F(1:elemK%F%dof, 1:2) )
          call SetF(elemK, 3, x(1:3, 1:2) )

          ! computing of outer normals
          allocate( elemK%n(1:elemK%flen,1:nbDim))
          allocate( elemK%dn(1:elemK%flen))
          elemK%face(neigh,:) = 0

          do k=1,elemK%flen
             k1 = mod(k,elemK%flen) + 1
             elemK%n(k,1) = x(k1,2) - x(k,2)

             elemK%n(k,2) = x(k,1)  - x(k1,1)

             elemK%dn(k) = sqrt( dot_product(elemK%n(k,:),elemK%n(k,:)) )
             !print*,'k::-', elemK%n(k, :), elemK%dn(k)
          enddo


          !print*,'###', elemK%F%iFlin
          !stop

          do i= 0, 3
             j = mod(i, 3) + 1
             !write(11, *) xx(j, 1:2), 0.
             write(100+is, *) x(j, 1:2), 0.
          enddo
          write(100+is, '(x)')

          !write(*,'(a8, 12es12.4)') 'Fx:', elemK%F%F(1:elemK%F%dof, 1)
          !write(*,'(a8, 12es12.4)') 'Fy:', elemK%F%F(1:elemK%F%dof, 2)

          !! integ nodes
          elemK%deg = deg
          elemK%dof  =  DOFtriang(elemK%deg)
          elemK%dof_plus  =  DOFtriang(elemK%deg + state%p_mod_max )

          elemK%face(fdeg,:) = elemK%deg
          elemK%face(fdof,:) = elemK%dof

          call SetElementQuadraturesDegrees( elemK )

          call ComputeIntegNode(elemK)

          call PrepareOneElement( elemK )

          call ComputeLocalMassMatrix(elemK )

          V_rule => state%space%V_rule(elemK%Qnum)
          Qdof = elemK%Qdof
          dof = elemK%dof
          dofP = elemK%dof_plus
          dofPP  =  DOFtriang(elemK%deg + 2)

          !print*,'Qnum = ;', elemK%Qnum, elemK%Qdof, elemK%dof,'dof:',dof,dofP,dofPP

          allocate( xi(1:Qdof, 1:2), wi(1:Qdof, 1:3), ww(1:dofP, 1:2) )

          ! function in integ nodes
          call ComputeF(elemK, Qdof, V_rule%lambda(1:Qdof, 1:2),  xi(1:Qdof, 1:2) )

          do i=1,Qdof
             wi(i, 1) = FFF(degP, xi(i, 1:2) )
             write(33,*) , xi(i, 1:2), wi(i, 1)
             !!write(*,'(a8, 6es12.4)') 'wi:', xi(i, 1:2), wi(i, 1)
          enddo

          allocate(elemK%wST_LS(1:1, 1:dofPP, -1:2), source = 0.0 )
          allocate(elemK%zST_LS(1:1, 1:dofPP, -1:2), source = 0.0  )

          ! function in DG basis
          call IntegrateVectorB(elemK, dofP, wi(1:Qdof, 1), ww(1:dofP, 1) )

          call SolveLocalMatrixProblem(dofP, elemK%Mass%Mb(1:dofP, 1:dofP), 1, ww(1:dofP, 1) )

          elemK%wST_LS(1, 1:dofP, 1) = ww(1:dofP, 1)

          call PlotElemFunction3D(34, elemK, dofP, elemK%wST_LS(1, 1:dofP, 1)  )

          ! integral of the function
          call IntegrateDGFunction2(elemK, dofP, elemK%wST_LS(1, 1:dofP, 1), val)
          val = sqrt(val)
          err(1, 1) = val ! error in the L^2-norm

          ! integral of the function over \partial K
          call IntegrateBoundarySquareFunction(elemK, dofP, elemK%wST_LS(1, 1:dofP,1), val)
          val = sqrt(val)
          err(2, 1) = val ! error in the L^2-norm over boundary

          ! integral of the function over \partial K
          call IntegrateBoundarySquareGradFunction(elemK,dofP,elemK%wST_LS(1,1:dofP,1),val)
          val = sqrt(val)
          err(3, 1) = val ! error of the gradient in the L^2-norm over boundary

          ! anisotropic estimate
          call Eval_High_Order_Deriv_Elem(elemK)

          allocate(anis(1:4, 1:5) )
          allocate(Du( 0:degP) )  ! derivatives of order deg
          allocate(ai( 0:2*degP), source = 0.0 )  ! coefficients

          Du(0:degP) =  elemK%wSS(1, 1, 0:degP )
          call Set_Coeffs_Function(degP, Du(0:degP), ai(0:degP) )

          call Draw_deriv(40, degP, ai(0:degP) )

          call FindAnisotropy(degP, ai(0:degP), anis(1, 1:3), upper )
          call Draw_anis(41, degP, anis(1, 1:3) )

          !  gradient

          call Set_Coeffs_Gradient(degP, degQ, Du(0:degP), mA(1:2, 1:2), ai(0:degQ) )
          !call Draw_deriv(45, degQ, ai(0:degQ) )

          call FindAnisotropy(degQ, ai(0:degQ), anis(2, 1:3), upper )
          !call Draw_anis(46, degQ, anis(2, 1:3) )


          !write(*,'(a8, 40es12.4)') 'deriv:', Du(:)
          !write(*,'(a8, 40es12.4)') 'coeffs:', ai(:)
          !write(*,'(a8, 40es12.4)') 'anisot:', anis(1, 1:3)

          allocate( Q(1:2, 1:2), M(1:2, 1:2), D(1:2, 1:2) )
          D = 0.
          vp = anis(1, 3)
          Q(1,1)  = cos(vp);  Q(1,2) = -sin(vp);  Q(2,1)  =  sin(vp);   Q(2,2) = cos(vp)
          D(1,1) = 1.;  D(2,2) = anis(1,2)**(-2./degP)
          M(1:2, 1:2) = matmul(Q(1:2, 1:2), matmul(D(1:2, 1:2), transpose( Q(1:2, 1:2))))

          do i=1,Qdof
             ! anisotropy
             val1 = dot_product( xi(i, 1:2), matmul( M(1:2, 1:2), xi(i, 1:2) ))
             wi(i, 2) = anis(1,1) * val1**(degP/2.)

             ! function
             wi(i, 3) = Eval_Dir_Derivative(degP, ai(0:degP), xi(i,1), xi(i,2))
             write(35,*) , xi(i, 1:2), abs(wi(i, 1:2)), wi(i,3)
          enddo
          ! exact integral of the anisotropic estimate
          val1 = sqrt( dot_product(V_rule%weights(1:Qdof), wi(1:Qdof,2)**2)*elemK%F%JF0/2)
          err(1, 2) = val1

          deallocate(Q, D, M)

          ieta = 6
          allocate(etasA(1:ieta), source = 0.0 )
          allocate(etasR(1:ieta), source = 0.0 )
          etasR(1)  = 1.

          !call Anis_estim_tot(ieta, deg, anis(1:4, 1:3), lambda_K, sigma_K, phi_K, &
          !     etasR(1:ieta), etasA(1:ieta), eta, optim )
          optim = 0
          call Anis_estim(1, deg, anis(1, 1:3), lambda, sigma, phi,  etasA(4), optim )
          err(1, 3) = etasA(4)

          call Anis_estim(2, deg, anis(1, 1:3), lambda, sigma, phi,  etasA(5), optim )
          err(2, 3) = etasA(5)

          call Anis_estim(3, deg, anis(2, 1:3), lambda, sigma, phi,  etasA(6), optim  )
          err(3, 3) = etasA(6)

          !write(*,'(a8, 2es12.4, 3(a10, 3es12.4))') 'sim,phi:',sigma, phi, &
          !     ' L^2(K):', err(1, 1), err(1,3), err(1,3)/err(1,1), &
          !     ' L^2(dK):', err(2, 1), err(2,3), err(2,3)/err(2,1), &
          !     ' H^1(dK):', err(3, 1), err(3,3), err(3,3)/err(3,1)

          do kk = 1, 3
             r_ieff(kk, 1) = max(r_ieff(kk, 1),  err(kk,3)/err(kk,1))
             r_ieff(kk, 2) = min(r_ieff(kk, 2),  err(kk,3)/err(kk,1))
          enddo

          write(80,*) sigma, phi, val, etasA(4), etasA(4) / val

          ! OVER-WRITES sigma
          !if(ip == Np .and. is == Ns) then
          !   call Anis_estim(1, deg, anis(1, 1:3), lambda, sigma, phi,  etasA(4), -1 )
          !   write(*,'(a10, 3es12.4, a10, 2es12.4)') &
          !        'OPTIM:', lambda, sigma, phi, ' estims:', val, etasA(4)
          !write(81,*) sigma, phi, val, etasA(4)
          !endif


          deallocate(xi, ww , wi, etasA, etasR, anis, Du, ai )
       enddo ! do ip = 1, Np

       if(sigma <= 50) then
          kk = 177
          open(kk, file = 'eta_III.dat', status='unknown', position='append')
          write(kk, '(i5, 40es12.4)') is, sigma, r_ieff(1:3, 1), r_ieff(1:3, 2)
          close(kk)
       endif

    end do  ! do is = 1, Ns


    print*
    print*
    print*
    print*
    stop "end subroutine TEST_ANIS_ESTIM"

  end subroutine TEST_ANIS_ESTIM

  function FFF(deg, xx)
    real :: FFF
    integer :: deg
    real, dimension(1:2), intent(in) :: xx
    real :: x, y

    x = xx(1)
    y = xx(2)

    if(deg == 2 ) then
       FFF = x**2 + y**2
    else if(deg == 3) then
       !FFF = (x**3 + 3*x*x*y + 3*x*y*y + y**3)
       FFF = (x**3 + 2*x*x*y + x*y*y - y**3)
       !FFF = (10*x**3 + 2*x*x*y + x*y*y - 0.1*y**3)
       !FFF = (x**2 + 3*x*y - 0.5* y**2)

    else if( deg == 4) then
       FFF = (4*x**4 + 2*x*x*y*y + x*y*y*y - y**4)
    else
       stop 'UNKNOWN FFF'
    endif

  end function FFF

  !> setting of physical coordinates
  subroutine set_coordinates(psi, phi, rl1, rl2, xx, x )
    real, intent(in):: psi, phi, rl1, rl2
    real, dimension(1:3, 1:2), intent(in) :: xx
    real, dimension(1:3, 1:2), intent(inout) :: x
    real, dimension(:,:), allocatable :: QT, Q, M, L
    integer :: i

    allocate(QT(1:2,1:2), Q(1:2, 1:2), M(1:2, 1:2), L(1:2, 1:2)  )
    QT = 0.
    Q = 0.
    L = 0.
    QT(1,1) = cos(psi);  QT(1,2) = sin(psi);  QT(2,1) = -sin(psi);  QT(2,2) = cos(psi)
    Q(1,1)  = cos(phi);  Q(1,2) = -sin(phi);  Q(2,1)  =  sin(phi);   Q(2,2) = cos(phi)
    L(1,1) = rl1;  L(2,2) = rl2

    !do i=1,2
    !   write(*,'(3(a5, 2es12.4))') 'Q:',  Q(i, :),'L:',  L(i, :),'QT:',  QT(i, :)
    !enddo

    M = matmul(Q, matmul(L, QT) )
    do i=1,3
       x(i, 1:2) = matmul( M(1:2, 1:2), xx(i, 1:2) )
    enddo

    deallocate(Q, QT, L, M)

  end subroutine set_coordinates
end module ama_hp_DWR
