real function dis2lat( y, phi0 )
!  基準緯度 phi0 から空間距離 y [m] 離れた位置における緯度 [rad].
!  y は北に計算する際は正値, 南に計算する際は負値を与えれば計算可能.
  use Phys_Const
  use Math_Const
  implicit none
  real, intent(in) :: y     ! 基準緯度 phi0 からの空間距離 [m] (北向きが正)
  real, intent(in) :: phi0  ! 基準緯度 [rad]

  dis2lat=y/radius+phi0

  return
end function

real function dis2lon( x, lon0, phi0 )
!  基準経度 lon0 から空間距離 x [m] 離れた位置における経度 [rad].
!  ただし, 基準緯度 phi0 [rad] から同緯度上で計測した距離を用いる.
!  x は東に計算する際は正値, 西に計算する際は負値を与えれば計算可能.
  use Phys_Const
  use Math_Const
  implicit none
  real, intent(in) :: x     ! 基準経度 lon0 からの空間距離 [m] (東向きが正)
  real, intent(in) :: lon0  ! 基準経度 [rad]
  real, intent(in) :: phi0  ! 基準緯度 [rad]

  dis2lon=x/(radius*cos(phi0))+lon0

  return
end function

real function dis2mlat(y,phi0)
!  基準緯度 phi0 から空間距離 y [m] 離れた位置における緯度 [rad].
!  ただし, x は東方向に正, 西方向に負を与えれば計算可能.
  use Phys_Const
  use Math_Const
  implicit none
  real, intent(in) :: y     ! 基準緯度 phi0 からの空間距離 [rad] (北向き正).
  real, intent(in) :: phi0  ! 基準緯度 [rad]

  dis2mlat=asin(tanh(log(tan(0.25*pi+0.5*phi0))+y/radius))

  return
end function

real function dis2mlon(x,lam0)
!  基準経度 lam0 から地図距離 x [m] 離れた位置における経度 [rad].
!  ただし, x は東方向に正, 西方向に負を与えれば計算可能.
  use Phys_Const
  implicit none
  real, intent(in) :: x     ! 基準経度 lam0 からの空間距離 [m] (東向きが正)
  real, intent(in) :: lam0  ! 基準経度 [rad]

  dis2mlon=x/radius+lam0

  return
end function

real function ll2radi( lon1, lat1, lon2, lat2 )
! ある球面上での 2 点での緯度, 経度を元に, それらを終始点とする円弧の距離を
! 計算する.
  use Phys_Const
  implicit none
  real, intent(in) :: lon1    ! 点 1 での経度 [rad]
  real, intent(in) :: lat1    ! 点 1 での緯度 [rad]
  real, intent(in) :: lon2    ! 点 2 での経度 [rad]
  real, intent(in) :: lat2    ! 点 2 での緯度 [rad]
  real :: tmp

  tmp=sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon2-lon1)
  ll2radi=acos(tmp)*radius

  return
end function ll2radi

real function lonlat2lamdis( lon, phi, lon0, phi1, phi2, phi0 )
!  基準緯度 phi1, phi2, 基準経度 lon0 から経度緯度 lon, phi [rad] 離れた
!  位置までの地図上の距離 [m].
  use Phys_Const
  use Math_Const
  implicit none
  real, intent(in) :: lon     ! 求めたい点の経度 [rad].
  real, intent(in) :: phi     ! 求めたい点の緯度 [rad].
  real, intent(in) :: lon0  ! 基準経度 [rad]
  real, intent(in) :: phi1  ! 基準緯度 1 [rad]
  real, intent(in) :: phi2  ! 基準緯度 2 [rad]
  real, intent(in), optional :: phi0  ! 地図座標基準緯度 0 [rad].
        ! phi0 が設定されている場合, y 座標の距離を求める.
  real :: n
  real :: rho, rho0

  n=log(cos(phi1)/cos(phi2))/log(tan(0.25*pi-0.5*phi1)/tan(0.25*pi-0.5*phi2))

  rho=(cos(phi1)*(tan(0.25*pi-0.5*phi))**n)/(n*(tan(0.25*pi-0.5*phi1))**n)

  if(present(phi0))then
     rho0=(cos(phi1)*(tan(0.25*pi-0.5*phi0))**n)/(n*(tan(0.25*pi-0.5*phi1))**n)
     lonlat2lamdis=rho0-rho*cos(n*(lon-lon0))
  else
     lonlat2lamdis=rho*sin(n*(lon-lon0))
  end if

  return
end function

subroutine rt2ll( r, theta, lon0, lat0, lon, lat )
! 地球の球面座標系において, (lon0, lat0) の位置を原点とした (r,\theta) 座標を
! 展開したとき, 距離 r, 同位角 theta の位置における球面上での経度を
! 計算する. 距離 r は球面上の円弧距離とし, theta = 90, 270 deg が球面の
! 子午面上を通るように theta は配置される.
! 球面三角法を用いて計算を行うが, 角度 theta は球面三角形の接角度と近似する.
  use Phys_Const
  use Math_Const
  implicit none
  real, intent(in) :: r          ! lon0 からの円弧距離 [m]
  real, intent(in) :: theta      ! lon0 での同位角 [rad]
  real, intent(in) :: lon0       ! 極座標原点の経度 [rad]
  real, intent(in) :: lat0       ! 極座標原点の緯度 [rad]
  real, intent(inout) :: lon     ! 計算された経度 [rad]
  real, intent(inout) :: lat     ! 計算された緯度 [rad]
  real :: thetad, lond, latd, tmplon, tmplat, rratio

  thetad=180.0*theta/pi
  lond=180.0*lon0/pi
  latd=180.0*lat0/pi
  rratio=r/radius

  if(thetad==-90.0.or.thetad==270.0)then
     lon=lon0
     lat=lat0-rratio
  else if(thetad==90.0)then
     lon=lon0
     lat=lat0+rratio
  else if((-90.0<thetad.and.90.0>thetad).or. (270.0<thetad.and.360.0>=thetad))then
     tmplat=cos(lat0)*sin(rratio)*sin(theta)+sin(lat0)*cos(rratio)
     lat=asin(tmplat)
     tmplon=sin(rratio)*cos(theta)/cos(asin(tmplat))
     lon=lon0+asin(tmplon)
  else if((90.0<thetad.and.270.0>thetad).or. (-180.0<=thetad.and.-90.0>thetad))then
     tmplat=cos(lat0)*sin(rratio)*sin(theta)+sin(lat0)*cos(rratio)
     lat=asin(tmplat)
     tmplon=-sin(rratio)*cos(theta)/cos(asin(tmplat))
     lon=lon0-asin(tmplon)
  else
     write(*,*) "### ERROR : (rt2ll:Map_Function)"
     write(*,*) "argument 'theta' is not valid : ", theta
     write(*,*) "STOP."
     stop
  end if

end subroutine rt2ll