
module misc

  implicit none

contains

  !  
  ! 領域の面積を計算.
  ! 
  real(8) function area(dkx, dky) result(retval)
    
    real(8), intent(in) :: dkx, dky
    real(8) dpi
    
    dpi = 2.0d0*acos(-1.0d0)
    
    retval = dpi**2/(dkx*dky)
  end function area

  ! 
  ! 最小値と離散化幅からインデックスを計算
  ! 
  integer function coord_to_index(x, xmin, dx) result(retval)
    
    real(8), intent(in) :: x, xmin, dx
    
    retval = nint(aint((x-xmin)/dx))
    
  end function coord_to_index

  !
  ! スペクトルで計算した任意の点の関数値を返す.
  !
  ! \sum psi(kx, ky) exp(i (kx*x+ky*y)) 
  !   
  real(8) function funcval2d(f, x, y, nx, ny, dkx, dky, ky)
    
    complex(8), intent(in) :: f(0:nx/2, 0:ny-1)
    real(8),    intent(in) :: x, y, dkx, dky
    integer,    intent(in) :: nx, ny, ky(0:ny-1)
    
    integer i,j
    real(8) s, theta
    
    s = 0
    do j=0,ny-1
       theta = dky*dble(ky(j))*y
       s = s + dble(f(0,j))*cos(theta)-dimag(f(0,j))*sin(theta)
       do i=1,nx/2
          theta =  dkx*dble(i)*x+dky*dble(ky(j))*y
          s = s + 2.0d0*(dble(f(i,j))*cos(theta)-dimag(f(i,j))*sin(theta))
       end do
    end do
    funcval2d = s
  end function funcval2d

  !
  ! 2次元場の線形補間. 
  ! 点 (x0, y0) を囲む矩形領域の頂点での関数値から, 
  ! 点 (x0, y0)での関数値を線形補間する.
  ! 
  ! 矩形領域の上端と下端で, x軸に沿って線形補間をそれぞれ実行し, 
  ! その後で 点 x0 上でy軸に沿って線形補間する.
  !
  real(8) function InterpolateLinear2d(f, x0, y0, x, y, nx, ny, dkx, dky) result(retval)
    
    real(8), intent(in) :: f(0:nx-1, 0:ny-1)
    real(8), intent(in) :: x0, y0, dkx, dky, x(0:nx-1), y(0:ny-1)
    integer, intent(in) :: nx, ny
    
    integer ix0,ix1,iy0,iy1
    real(8) fy0,fy1,dx,dy,dpi
    
    dpi = 2.0d0*acos(-1.0d0)
    dx = dpi/dkx/dble(nx)
    dy = dpi/dky/dble(ny)

    ix0 = lower_index_in_axis(x0, nx, dkx)
    ix1 = upper_index_in_axis(x0, nx, dkx)
    iy0 = lower_index_in_axis(y0, ny, dky)
    iy1 = upper_index_in_axis(y0, ny, dky)
!    print *, ix0,ix1,iy0,iy1
    fy0 = InterpolateLinear1d(f(ix0,iy0), f(ix1,iy0), (x0-x(ix0))/dx)
    fy1 = InterpolateLinear1d(f(ix0,iy1), f(ix1,iy1), (x0-x(ix0))/dx)
    retval = InterpolateLinear1d(fy1, fy0, (y0-y(iy0))/dy)
  end function InterpolateLinear2d

  ! 
  ! f0, f1の線形補間. 
  ! r = (x - x0)/dx とする.
  ! 
  real(8) function InterpolateLinear1d(f0, f1, r) result(retval)
    
    real(8), intent(in) :: f0, f1, r
    
    retval = (f1-f0)*r + f0
  end function InterpolateLinear1d
  
  !
  ! 与えられた点 x 以下の最大の座標のインデックスを返す.
  ! 点 x は 0 <= x < 2\pi/dkx を満たすとする.
  !
  integer function lower_index_in_axis(x, nx, dkx) result(retval)
    
    real(8), intent(in) :: x, dkx
    integer, intent(in) :: nx
    
    real(8) dpi
    
    dpi = 2.0d0*acos(-1.0d0)
    
    retval = nint(x/dpi*dkx*nx)
    
  end function lower_index_in_axis
  
  !
  ! 周期境界の場合で与えられた点 x より大きい最小の座標のインデックスを返す.
  ! 点 x は 0 <= x < 2\pi/dkx を満たすとする.
  !
  integer function upper_index_in_axis(x, nx, dkx) result(retval)
    
    real(8), intent(in) :: x, dkx
    integer, intent(in) :: nx
    
    integer i
    i = lower_index_in_axis(x, nx, dkx)
    if(i==nx-1)then
       retval = 0
    else
       retval = i+1
    end if
  end function upper_index_in_axis
  
end module misc
