log_mesh_mod.F90

      Module log_mesh
        use iso_fortran_env, only: output_unit, error_unit
        use runtime_error_mod
        use files_mod

        Type logmesh
           Real (Kind=Kind(0.d0))  :: lambda, center, log_lambda
           Real (Kind=Kind(0.d0))  :: range
           Real (Kind=Kind(0.d0))  :: om_st, om_en, dom
           Real (Kind=Kind(0.d0))  :: precision
           Integer        :: nom,nw
           Real (Kind=Kind(0.d0)), pointer :: xom(:),dxom(:)
           Character(len=10) :: type
        end Type logmesh

        Interface lookup_log_mesh
           module procedure lookup_log_mesh_r, lookup_log_mesh_c
        end Interface
        Interface inter_log_mesh
           module procedure inter_log_mesh_r, inter_log_mesh_c
        end Interface

      Contains

        

        subroutine make_log_mesh ( Mesh,  Lambda, Center, Rng, Type, Nw_1 )

          Implicit None

          type(logmesh)      :: mesh
          Real (Kind=Kind(0.d0))       :: lambda, center, rng
          Integer, Optional   :: nw_1
          Integer             :: n, nw
          Character(len=10)   :: type

          Real (Kind=Kind(0.d0))  :: dom, om_st, om_en

          mesh%Center     = center
          mesh%Range      = rng
          If (str_to_upper(type) == "LOG" ) Then
             om_st = center - rng
             om_en = center + rng
             mesh%Om_st = om_st
             mesh%Om_en = om_en
             mesh%Lambda     = lambda
             mesh%Type       = "Log"
             if (Present(nw_1) ) then
                nw = nw_1
             else
                nw = nint(10.d0*log(10.d0)/log(lambda))
             endif
             mesh%Nw         = nw
             mesh%Nom        = 2*nw + 3
             mesh%Log_Lambda = log(lambda)
             Allocate   ( mesh%Xom(2*nw + 3), mesh%DXom(2*nw+3) )
             Do n = 0,nw
                mesh%xom (n+1          ) =  center  -   rng * (lambda**(-n))
             enddo
             mesh%xom   (nw+2         ) =  center
             do n = nw,0,-1
                mesh%xom(nw+3 +(nw-n) ) =  center  +   rng * (lambda**(-n))
             enddo
             mesh%Precision = mesh%Lambda**(-mesh%Nw)
          elseif (str_to_upper(type) == "LIN" ) then
             mesh%Type       = "Lin"
             If ( Present(nw_1)  ) then
                nw   = nw_1
                mesh%Nw   = nw
                mesh%Nom  = 2*nw + 1
                mesh%Type = "Lin"
                Allocate   ( mesh%Xom(2*nw + 1), mesh%DXom(2*nw+1) )
                om_st = center - rng
                om_en = center + rng
                dom = rng/dble(nw_1)
                mesh%Dom   = dom
                mesh%Om_st = om_st
                mesh%Om_en = om_en
                do n  = 1,mesh%Nom
                   mesh%xom(n)  = om_st + dble(n-1)*dom
                enddo
             else
                Write(error_unit,*) 'Make_log_mesh: You need to include Nw for the Lin Mesh'
                Call terminate_on_error(error_generic,__file__,__line__)
             endif
          else
             Write(error_unit,*) 'Make_log_mesh: Mesh has no type!!'
             Call terminate_on_error(error_generic,__file__,__line__)
          endif
          do n = 1,mesh%Nom-1
             mesh%DXom(n) = mesh%xom (n+1) - mesh%xom (n )
          enddo

        end subroutine make_log_mesh

        subroutine clear_log_mesh ( Mesh )
          Implicit None

          type(logmesh)      :: mesh

          deallocate   ( mesh%Xom,  mesh%DXom  )

        end subroutine clear_log_mesh

        Integer  Function m_find(X,Mesh)

          Implicit None

          type(logmesh) :: mesh
          Real (Kind=Kind(0.d0))  :: x
          Integer        :: m

          if ( str_to_upper(mesh%Type)  == "LOG" ) then
             if ( x >  (mesh%OM_en) .or.  x < (mesh%Om_st) ) then
                m = 0
             else
                if      ( x < mesh%Xom(mesh%Nw+1) ) then
                   m = 2 - int( log( (mesh%Center - x)/mesh%Range )  / mesh%Log_Lambda )
                   !Write(6,*) 'Hi 1', X
                elseif  ( x >  mesh%Xom(mesh%Nw+3) ) then
                   m = 2*mesh%Nw + 3 + int( log( (x- mesh%Center) /mesh%Range )  / mesh%Log_Lambda )
                   !Write(6,*) 'Hi 2', X, Mesh%Center +  Mesh%Range
                elseif  ( x >  mesh%Center )      then
                   m = mesh%Nw+3
                else
                   m = mesh%Nw+2
                endif
             endif
             m_find = m
          else
             m_find = int((x - mesh%Om_st)/mesh%dom) + 2
             if (m_find > mesh%Nom) m_find=mesh%Nom
             if (m_find < 2 ) m_find=2
          endif


          !Write(6,*)
          !Write(6,*) 'Point: ', X
          !if (  m > 0 ) then
          !   Write(6,*) 'Your point lies inbetween ', Mesh%Xom(m-1), ' and ', Mesh%Xom(m)
          !else
          !   Write(6,*) 'Out of range '
          !endif

        end Function m_find
!*******
        Real(Kind=Kind(0.d0)) Function  lookup_log_mesh_r(f, x,Mesh,m_1)

          Implicit None

          type(logmesh) :: mesh
          Real (Kind=Kind(0.d0)), dimension(:) :: f
          Real (Kind=Kind(0.d0))  :: x
          Integer      , Optional     :: m_1

          Integer ::  m
          Real (Kind=Kind(0.d0)) :: x1,x2,y1,y2,a,b

          m = m_find(x,mesh)
          if (m == 0 ) then
             lookup_log_mesh_r = 0.d0
          else
             x1 = mesh%xom(m-1)
             x2 = mesh%xom(m  )
             y1 = f(m-1)
             y2 = f(m)
             a = (y1-y2)/(x1-x2)
             b = (x1*y2 - x2*y1)/(x1-x2)
             lookup_log_mesh_r = a*x + b
          endif

          If  ( Present(m_1) )  m_1 = m

        end Function lookup_log_mesh_r



!*******
!!$        Complex (Kind=Kind(0.d0)) Function  Lookup_log_mesh_C(f, x,Mesh,m_1)
!!$
!!$          Implicit None
!!$
!!$          Type (logmesh) :: Mesh
!!$          Complex (Kind=Kind(0.d0)), dimension(:) :: f
!!$          Real    (Kind=Kind(0.d0))               :: X
!!$          Integer      , Optional     :: m_1
!!$
!!$
!!$          Integer  ::  n, m
!!$          Complex  (Kind=Kind(0.d0)) :: X1,X2,Y1,Y2,a,b
!!$
!!$          m = m_find(X,Mesh)
!!$          if (m == 0 ) then
!!$             Lookup_log_mesh_C = cmplx(0.d0,0.d0)
!!$          else
!!$             x1 = cmplx( Mesh%xom(m-1),0.d0 )
!!$             x2 = cmplx( Mesh%xom(m  ),0.d0 )
!!$             y1 = f(m-1)
!!$             y2 = f(m  )
!!$             a = (y1-y2)/(x1-x2)
!!$             b = (x1*y2 - x2*y1)/(x1-x2)
!!$             Lookup_log_mesh_C = a*cmplx( x , 0.d0 ) + b
!!$          endif
!!$
!!$          If  ( Present(m_1) )  m_1 = m
!!$
!!$        end Function Lookup_log_mesh_C

        Complex (Kind=Kind(0.d0)) Function  lookup_log_mesh_c(f, x,Mesh,m_1)

          Implicit None

          type(logmesh) :: mesh
          Complex (Kind=Kind(0.d0)), dimension(:) :: f
          Real    (Kind=Kind(0.d0))               :: x
          Integer      , Optional     :: m_1


          Integer  ::  m
          Complex  (Kind=Kind(0.d0)) :: z1,z2, z
          Real     (Kind=Kind(0.d0)) :: x1,x2,t

          m = m_find(x,mesh)
          if (m == 0 ) then
             lookup_log_mesh_c = cmplx(0.d0, 0.d0, kind(0.d0))
          else
             x1 =  mesh%xom(m-1)
             x2 =  mesh%xom(m  )
             t  = (x1 - x)/(x2-x1)
             z1 = f(m-1)
             z2 = f(m  )
             z  = z1 + (z1-z2)*t
             lookup_log_mesh_c = z
          endif

          If  ( Present(m_1) )  m_1 = m

        end Function lookup_log_mesh_c


!******
        Real (Kind=Kind(0.d0)) Function  inter_log_mesh_r(f,Mesh)

          Implicit None

          type(logmesh) :: mesh
          Real     (Kind=Kind(0.d0)), dimension(:) :: f
          Real     (Kind=Kind(0.d0))  :: x
          Integer            :: n

          x = 0.d0
          do n = 1,mesh%Nom-1
             x = x + mesh%DXom(n) * (f(n+1) + f(n) )
          enddo
          inter_log_mesh_r = x / 2.d0

        end Function inter_log_mesh_r

!******
        Complex (Kind=Kind(0.d0)) Function  inter_log_mesh_c(f,Mesh)

          Implicit None

          type(logmesh) :: mesh
          Complex     (Kind=Kind(0.d0)), dimension(:) :: f
          Complex     (Kind=Kind(0.d0))  :: z
          Integer               :: n

          z = cmplx(0.d0, 0.d0, kind(0.d0))
          do n = 1,mesh%Nom-1
             z = z + mesh%DXom(n) * ( f(n+1) + f(n) )
          enddo
          inter_log_mesh_c = z /2.d0

        end Function inter_log_mesh_c



        subroutine print_log_mesh(Mesh)

          Implicit None

          type(logmesh) :: mesh

          Integer :: n

          If (str_to_upper(mesh%Type) == "LOG" ) Then
             Open (unit=10,file="Log_Mesh", status="unknown" )
             Write(10,*) '#    Log Mesh     : '
             Write(10,*) '#    Lambda       : ', mesh%Lambda
             Write(10,*) '#    Range        : ', mesh%Range
             Write(10,*) '#    Center       : ', mesh%Center
             Write(10,*) '#    Nom          : ', mesh%Nom
             Write(10,*) '#    Precision    : ', mesh%Lambda**(-mesh%Nw)
             do n = 1,mesh%Nom
                write(10,"(F16.8)") mesh%xom(n)
             enddo
             close(10)
          endif

          If (str_to_upper(mesh%Type) == "LIN" ) Then
             Open (unit=10,file="Lin_Mesh", status="unknown" )
             Write(10,*) '#    Lin Mesh     : '
             Write(10,*) '#    Range        : ', mesh%Range
             Write(10,*) '#    Center       : ', mesh%Center
             Write(10,*) '#    Nom          : ', mesh%Nom
             Write(10,*) '#    Dom          : ', mesh%dom
             do n = 1,mesh%Nom
                write(10,"(F16.8)") mesh%xom(n)
             enddo
             close(10)
          endif



        end subroutine print_log_mesh


      end Module log_mesh