module lib_proj_sinusoidal
  private
!---------------------------------------------------------------
  public :: convLonLatToSinusoidal
  public :: convSinusoidalToLonLat
!---------------------------------------------------------------
  real(8), parameter :: pi = acos(-1.d0)
!---------------------------------------------------------------
contains
!===============================================================
!
!===============================================================
subroutine convLonLatToSinusoidal(lon, lat, x, y, lon_0, a, e)
  implicit none
  real(8), intent(in)  :: lon, lat
  real(8), intent(out) :: x, y
  real(8), intent(in)  :: lon_0
  real(8), intent(in)  :: a
  real(8), intent(in)  :: e
  real(8) :: b

  call calcParams(a, e, b)

  x = a*0.5d0 * (lon - lon_0) * cos(lat)
  y = b*0.5d0 * lat / pi
end subroutine convLonLatToSinusoidal
!===============================================================
!
!===============================================================
subroutine convSinusoidalToLonLat_0d(x, y, lon, lat, lon_0, a, e)
  implicit none
  real(8), intent(in)  :: x, y
  real(8), intent(out) :: lon, lat
  real(8), intent(in)  :: lon_0
  real(8), intent(in)  :: a
  real(8), intent(in)  :: e
  real(8) :: b

  call calcParams(a, e, b)

  lat = y * pi / (b*0.5d0)
  lon = x / (cos(lat) * a*0.5d0) + lon_0
end subroutine convSinusoidalToLonLat_0d
!===============================================================
!
!===============================================================
subroutine calcParams(a, e, b)
  implicit none
  real(8), intent(in)  :: a  ! Semi-major axis
  real(8), intent(in)  :: e  ! Eccentricity
  real(8), intent(out) :: b  ! Semi-minor axis

  b = a * sqrt(1-e**2)
end subroutine calcParams
!===============================================================
!
!===============================================================
end module lib_proj_sinusoidal