module EarthClass
   use sim
   implicit none

   type :: TimeZoneList_
      type(String_) :: ID(69)
      real(real64) :: offset(69)
   contains
      procedure, public :: init => initTimeZoneList
   end type

   type :: Earth_
      real(real64) :: AirPressure = 101.325d0 ! kPa
      real(real64) :: AxisTilt = 23.4392811d0 ! deg. (https://en.wikipedia.org/wiki/Earth)
      real(real64) :: StandardGravity = 9.80665d0 ! m/s^2

      ! my position and timezone
      real(real64) :: MyPosition(2) = [36.380d0, 140.470d0] ! Ibaraki Pref. JAPAN
      real(real64) :: MyTimeZonePosition = 135.0d0! JAPAN
      character(:), allocatable :: TimeZone
      character(:), allocatable :: TimeZoneName
      real(real64) :: TimeZoneOffset

      ! my time
      integer(int32) :: Year = 2022
      integer(int32) :: Month = 1
      integer(int32) :: Day = 1
      integer(int32) :: Hour = 0
      integer(int32) :: Minute = 0
      integer(int32) :: Second = 0
      integer(int32) :: Day_Per_Month(0:13) = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31]

   contains
      procedure, public :: init => initEarth
      procedure, public :: setTime => setTimeEarthClass
      procedure, public :: setPosition => setPositionEarthClass
      procedure, public :: setTimeZone => setTimeZoneEarthClass
      procedure, public :: getSunTime => getSunTimeEarthClass
      procedure, public :: getSunPosition => getSunPositionEarthClass
      procedure, public :: getTimeZone => getTimeZoneEarthClass
      procedure, public :: getTimeZoneOffset => getTimeZoneOffsetEarthClass
   end type
contains
   subroutine initTimeZoneList(this)
      class(TimeZoneList_), intent(inout) :: this

      this%id(1)%all = "Africa/Johannesburg"
      this%id(2)%all = "Africa/Lagos"
      this%id(3)%all = "Africa/Windhoek"
      this%id(4)%all = "America/Adak"
      this%id(5)%all = "America/Anchorage"
      this%id(6)%all = "America/Argentina_Buenos_Aires"
      this%id(7)%all = "America/Bogota"
      this%id(8)%all = "America/Caracas"
      this%id(9)%all = "America/Chicago"
      this%id(10)%all = "America/Denver"
      this%id(11)%all = "America/Godthab"
      this%id(12)%all = "America/Guatemala"
      this%id(13)%all = "America/Halifax"
      this%id(14)%all = "America/Los_Angeles"
      this%id(15)%all = "America/Montevideo"
      this%id(16)%all = "America/New_York"
      this%id(17)%all = "America/Noronha"
      this%id(18)%all = "America/Phoenix"
      this%id(19)%all = "America/Santiago"
      this%id(20)%all = "America/Santo_Domingo"
      this%id(21)%all = "America/St_Johns"
      this%id(22)%all = "Asia/Baghdad"
      this%id(23)%all = "Asia/Baku"
      this%id(24)%all = "Asia/Beirut"
      this%id(25)%all = "Asia/Dhaka"
      this%id(26)%all = "Asia/Dubai"
      this%id(27)%all = "Asia/Irkutsk"
      this%id(28)%all = "Asia/Jakarta"
      this%id(29)%all = "Asia/Kabul"
      this%id(30)%all = "Asia/Kamchatka"
      this%id(31)%all = "Asia/Karachi"
      this%id(32)%all = "Asia/Kathmandu"
      this%id(33)%all = "Asia/Kolkata"
      this%id(34)%all = "Asia/Krasnoyarsk"
      this%id(35)%all = "Asia/Omsk"
      this%id(36)%all = "Asia/Rangoon"
      this%id(37)%all = "Asia/Shanghai"
      this%id(38)%all = "Asia/Tehran"
      this%id(39)%all = "Asia/Tokyo"
      this%id(40)%all = "Asia/Vladivostok"
      this%id(41)%all = "Asia/Yakutsk"
      this%id(42)%all = "Asia/Yekaterinburg"
      this%id(43)%all = "Atlantic/Azores"
      this%id(44)%all = "Atlantic/Cape_Verde"
      this%id(45)%all = "Australia/Adelaide"
      this%id(46)%all = "Australia/Brisbane"
      this%id(47)%all = "Australia/Darwin"
      this%id(48)%all = "Australia/Eucla"
      this%id(49)%all = "Australia/Lord_Howe"
      this%id(50)%all = "Australia/Sydney"
      this%id(51)%all = "Etc/GMT12"
      this%id(52)%all = "Europe/Berlin"
      this%id(53)%all = "Europe/London"
      this%id(54)%all = "Europe/Moscow"
      this%id(55)%all = "Pacific/Apia"
      this%id(56)%all = "Pacific/Auckland"
      this%id(57)%all = "Pacific/Chatham"
      this%id(58)%all = "Pacific/Easter"
      this%id(59)%all = "Pacific/Gambier"
      this%id(60)%all = "Pacific/Honolulu"
      this%id(61)%all = "Pacific/Kiritimati"
      this%id(62)%all = "Pacific/Majuro"
      this%id(63)%all = "Pacific/Marquesas"
      this%id(64)%all = "Pacific/Norfolk"
      this%id(65)%all = "Pacific/Noumea"
      this%id(66)%all = "Pacific/Pago_Pago"
      this%id(67)%all = "Pacific/Pitcairn"
      this%id(68)%all = "Pacific/Tongatapu"
      this%id(69)%all = "UTC"

      this%offset(1) = dble(2)
      this%offset(2) = dble(1)
      this%offset(3) = dble(1)
      this%offset(4) = dble(-10)
      this%offset(5) = dble(-9)
      this%offset(6) = dble(-3)
      this%offset(7) = dble(-5)
      this%offset(8) = dble(-4.5)
      this%offset(9) = dble(-6)
      this%offset(10) = dble(-7)
      this%offset(11) = dble(-3)
      this%offset(12) = dble(-6)
      this%offset(13) = dble(-4)
      this%offset(14) = dble(-8)
      this%offset(15) = dble(-3)
      this%offset(16) = dble(-5)
      this%offset(17) = dble(-2)
      this%offset(18) = dble(-7)
      this%offset(19) = dble(-4)
      this%offset(20) = dble(-4)
      this%offset(21) = dble(-3.5)
      this%offset(22) = dble(3)
      this%offset(23) = dble(4)
      this%offset(24) = dble(2)
      this%offset(25) = dble(6)
      this%offset(26) = dble(4)
      this%offset(27) = dble(9)
      this%offset(28) = dble(7)
      this%offset(29) = dble(4.5)
      this%offset(30) = dble(12)
      this%offset(31) = dble(5)
      this%offset(32) = dble(5.75)
      this%offset(33) = dble(5.5)
      this%offset(34) = dble(8)
      this%offset(35) = dble(7)
      this%offset(36) = dble(6.5)
      this%offset(37) = dble(8)
      this%offset(38) = dble(3.5)
      this%offset(39) = dble(9)
      this%offset(40) = dble(11)
      this%offset(41) = dble(10)
      this%offset(42) = dble(6)
      this%offset(43) = dble(-1)
      this%offset(44) = dble(-1)
      this%offset(45) = dble(9.5)
      this%offset(46) = dble(10)
      this%offset(47) = dble(9.5)
      this%offset(48) = dble(8.75)
      this%offset(49) = dble(10.5)
      this%offset(50) = dble(10)
      this%offset(51) = dble(-12)
      this%offset(52) = dble(1)
      this%offset(53) = dble(0)
      this%offset(54) = dble(4)
      this%offset(55) = dble(13)
      this%offset(56) = dble(12)
      this%offset(57) = dble(12.75)
      this%offset(58) = dble(-6)
      this%offset(59) = dble(-9)
      this%offset(60) = dble(-10)
      this%offset(61) = dble(14)
      this%offset(62) = dble(12)
      this%offset(63) = dble(-9.5)
      this%offset(64) = dble(11.5)
      this%offset(65) = dble(11)
      this%offset(66) = dble(-11)
      this%offset(67) = dble(-8)
      this%offset(68) = dble(13)
      this%offset(69) = dble(0)

   end subroutine

! ###############################################################
   subroutine initEarth(this, Now, time, DateTime)
      class(Earth_), intent(inout) :: this
      integer(int32), optional, intent(in) :: DateTime(6) ! Year, Month, Day, Hour, Minute, Second
      logical, optional, intent(in) :: Now
      type(Time_), optional, intent(in) :: time
      character(:), allocatable :: d_time

      if (present(Now)) then
         call this%setTime(Now=Now)
      elseif (present(DateTime)) then
         call this%setTime(DateTime=DateTime)
      elseif (present(time)) then
         call this%setTime(time=time)
      else
         call this%setTime(Now=.true.)
      end if

      this%TimeZoneName = "Asia/Tokyo" ! JAPAN
      this%MyPosition(1:2) = [36.380d0, 140.470d0] ! Ibaraki Pref. JAPAN
      this%MyTimeZonePosition = 135.0d0! JAPAN
      this%TimeZone = "+0900" ! JAPAN

   end subroutine

! ###############################################################
   subroutine setTimeEarthClass(this, Now, time, DateTime)
      class(Earth_), intent(inout) :: this
      integer(int32), optional, intent(in) :: DateTime(6) ! Year, Month, Day, Hour, Minute, Second
      logical, optional, intent(in) :: Now
      type(Time_), optional, intent(in) :: time
      type(Time_) :: internal_time
      character(:), allocatable :: d_time

      if (present(DateTime)) then
         this%Year = DateTime(1)
         this%Month = DateTime(2)
         this%Day = DateTime(3)
         this%Hour = DateTime(4)
         this%Minute = DateTime(5)
         this%Second = DateTime(6)
      end if

      if (present(Now)) then
         if (Now) then
            d_time = internal_time%DateAndTime()
            this%Year = fint(d_time(1:4))
            this%Month = fint(d_time(5:6))
            this%Day = fint(d_time(7:8))
            this%Hour = fint(d_time(9:10))
            this%Minute = fint(d_time(11:12))
            this%Second = fint(d_time(13:14))
         end if
      end if

      if (present(time)) then
         this%Year = fint(time%date(1:4))
         this%Month = fint(time%date(5:6))
         this%Day = fint(time%date(7:8))
         this%Hour = fint(time%time(1:2))
         this%Minute = fint(time%time(3:4))
         this%Second = fint(time%time(5:6))
         this%timezone = time%zone
      end if
   end subroutine
! ###############################################################

! ###############################################################
   subroutine setTimeZoneEarthClass(this, name)
      class(Earth_), intent(inout) :: this
      character(*), intent(in) :: name
      type(TimeZoneList_)::TimeZoneList

      this%MyTimeZonePosition = this%getTimeZoneOffset(name=name)*15.0d0
      this%TimeZoneOffset = this%getTimeZoneOffset(name=name)

   end subroutine

! ###############################################################

! ###############################################################
   subroutine setPositionEarthClass(this, MyPosition, MyTimeZonePosition, MyTimeZone)
      class(Earth_), intent(inout) :: this
      real(real64), optional, intent(in) :: MyPosition(2)
      real(real64), optional, intent(in) :: MyTimeZonePosition
      character(5), optional, intent(in) :: MyTimeZone

      this%MyPosition(1:2) = MyPosition!e.g.[36.380d0, 140.470d0] ! Ibaraki Pref. JAPAN
      this%MyTimeZonePosition = MyTimeZonePosition!e.g.[135.0d0]! JAPAN
      this%TimeZone = MyTimeZone! e.g. "+0900" ! JAPAN

   end subroutine
! ###############################################################

   function getSunTimeEarthClass(this, Now) result(times)
      !http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunShineAngle/ssa.html
      class(Earth_), intent(inout) :: this
      integer(int32), intent(in) :: Now(6) ! Year, Month, Day, Hour, Minute, Second
      real(real64) :: times(2)
      type(Math_)  :: math
      real(real64) :: delta
      real(real64) :: omega
      real(real64) :: J, e
      real(real64) :: t, TT, h, A, TT1, TT2, t1, t2, sinA, cosA, Ts, buf
      integer(int32) :: total_days

      this%Year = Now(1)
      this%Month = Now(2)
      this%Day = Now(3)
      this%Hour = Now(4)
      this%Minute = Now(5)
      this%Second = Now(6)

      total_days = sum(this%Day_Per_Month(1:this%Month - 1)) + this%Day - 1
      omega = 2.0d0*math%pi/365.0d0
      J = dble(total_days) + 0.50d0

      delta = 0.33281d0 - 22.984d0*cos(omega*J) - 0.34990d0*cos(2.0d0*omega*J) &
              - 0.13980d0*cos(3.0d0*omega*J) + 3.7872d0*sin(omega*J) &
              + 0.03250d0*sin(2.0d0*omega*J) + 0.07187d0*sin(3.0d0*omega*J)

      e = 0.0072d0*cos(omega*J) - 0.0528d0*cos(2.0d0*omega*J) - 0.0012d0*cos(3.0d0*omega*J) &
          - 0.1229d0*sin(omega*J) - 0.1565d0*sin(2.0d0*omega*J) &
          - 0.0041d0*sin(3.0d0*omega*J)

      Ts = dble(this%Hour) + dble(this%Minute)/60.0d0 + dble(this%Second)/3600.0d0

      TT = Ts + (this%MyPosition(2) - this%MyTimeZonePosition)/15.0d0 + e

      t = 15.0d0*TT - 180.0d0

      buf = sin(radian(this%MyPosition(1)))*sin(radian(delta)) &
            + cos(radian(this%MyPosition(1)))*cos(radian(delta))*cos(radian(t))
      h = asin(buf)

      sinA = cos(radian(delta))*sin(radian(t))/cos(h)
      cosA = (sin(h)*sin(radian(this%MyPosition(1))) &
              - sin(radian(delta)))/cos(h)/cos(radian(this%MyPosition(1)))
      A = atan2(sinA, cosA) + math%pi

      t = acos(-tan(radian(delta))*tan(radian(this%MyPosition(1))))
      TT1 = (-degrees(t) + 180.0d0)/15.0d0
      t1 = TT1 - (this%MyPosition(2) - 135.0d0)/15.0d0 - e

      TT2 = (degrees(t) + 180.0d0)/15.0d0
      t2 = TT2 - (this%MyPosition(2) - 135.0d0)/15.0d0 - e

      times(1) = t1 ! 日の出(hr)
      times(2) = t2 ! 日没(hr)

   end function

   function getSunPositionEarthClass(this, Now, DateTime, dt) result(angles)
      !http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunShineAngle/ssa.html
      class(Earth_), intent(inout) :: this
      integer(int32), optional, intent(in) :: DateTime(6) ! Year, Month, Day, Hour, Minute, Second
      logical, optional, intent(in) :: Now
      integer(int32), optional, intent(in) :: dt
      real(real64) :: angles(2)
      type(Math_)  :: math
      real(real64) :: delta
      real(real64) :: omega
      real(real64) :: J, e
      real(real64) :: t, TT, h, A, TT1, TT2, t1, t2, sinA, cosA, Ts, buf
      integer(int32) :: total_days

      type(Time_) :: time

      character(:), allocatable :: d_time

      if (present(DateTime)) then
         this%Year = DateTime(1)
         this%Month = DateTime(2)
         this%Day = DateTime(3)
         this%Hour = DateTime(4)
         this%Minute = DateTime(5)
         this%Second = DateTime(6)
      end if

      if (present(Now)) then
         if (Now) then
            d_time = time%DateAndTime()
            this%Year = fint(d_time(1:4))
            this%Month = fint(d_time(5:6))
            this%Day = fint(d_time(7:8))
            this%Hour = fint(d_time(9:10))
            this%Minute = fint(d_time(11:12))
            this%Second = fint(d_time(13:14))
         end if
      end if

      if (present(dt)) then
         ! dt (sec.)
         this%Second = this%Second + dt
         do
            if (this%Second >= 60) then
               this%Second = this%Second - 60
               this%Minute = this%Minute + 1
            else
               exit
            end if
         end do

         do
            if (this%Minute >= 60) then
               this%Minute = this%Minute - 60
               this%hour = this%hour + 1
            else
               exit
            end if
         end do

         do
            if (this%hour >= 24) then
               this%hour = this%hour - 24
               this%day = this%day + 1
            else
               exit
            end if
         end do

         do
            if (this%day > this%Day_Per_Month(this%month)) then
               this%day = this%day - this%Day_Per_Month(this%month)
               this%Month = this%Month + 1
            else
               exit
            end if
         end do

         do
            if (this%month >= 13) then
               this%month = this%month - 12
               this%Year = this%Year + 1
            else
               exit
            end if
         end do
      end if

      total_days = sum(this%Day_Per_Month(1:this%Month - 1)) + this%Day - 1
      omega = 2.0d0*math%pi/365.0d0
      J = dble(total_days) + 0.50d0

      delta = 0.33281d0 - 22.984d0*cos(omega*J) - 0.34990d0*cos(2.0d0*omega*J) &
              - 0.13980d0*cos(3.0d0*omega*J) + 3.7872d0*sin(omega*J) &
              + 0.03250d0*sin(2.0d0*omega*J) + 0.07187d0*sin(3.0d0*omega*J)

      e = 0.0072d0*cos(omega*J) - 0.0528d0*cos(2.0d0*omega*J) - 0.0012d0*cos(3.0d0*omega*J) &
          - 0.1229d0*sin(omega*J) - 0.1565d0*sin(2.0d0*omega*J) &
          - 0.0041d0*sin(3.0d0*omega*J)

      Ts = dble(this%Hour) + dble(this%Minute)/60.0d0 + dble(this%Second)/3600.0d0

      TT = Ts + (this%MyPosition(2) - this%MyTimeZonePosition)/15.0d0 + e

      t = 15.0d0*TT - 180.0d0

      buf = sin(radian(this%MyPosition(1)))*sin(radian(delta)) &
            + cos(radian(this%MyPosition(1)))*cos(radian(delta))*cos(radian(t))
      h = asin(buf)

      sinA = cos(radian(delta))*sin(radian(t))/cos(h)
      cosA = (sin(h)*sin(radian(this%MyPosition(1))) &
              - sin(radian(delta)))/cos(h)/cos(radian(this%MyPosition(1)))
      A = atan2(sinA, cosA) + math%pi

      t = acos(-tan(radian(delta))*tan(radian(this%MyPosition(1))))
      TT1 = (-degrees(t) + 180.0d0)/15.0d0
      t1 = TT1 - (this%MyPosition(2) - 135.0d0)/15.0d0 - e

      TT2 = (degrees(t) + 180.0d0)/15.0d0
      t2 = TT2 - (this%MyPosition(2) - 135.0d0)/15.0d0 - e

      angles(1) = degrees(A) ! 方位角(deg.)
      angles(2) = degrees(H) ! 太陽高度(deg.)

      if (Ts < t1) then
         angles(:) = 0.0d0
         return
      end if

      if (Ts > t2) then
         angles(:) = 0.0d0
         return
      end if

   end function
! ################################################################

! ################################################################
   function getTimeZoneEarthClass(this, name) result(timezone)
      class(Earth_), intent(in) :: this
      character(*), intent(in) :: name
      type(TimeZoneList_) :: TimeZoneList
      character(:), allocatable :: timezone
      integer(int32) :: i

      timezone = ""
      call TimeZoneList%init()
      do i = 1, size(TimeZoneList%ID, 1)
         if (index(TimeZoneList%ID(i)%all, name) == 0) then
            cycle
         else
            timezone = TimeZoneList%ID(i)%all
            return
         end if
      end do

      print *, "[ERROR] no such time-zone name as", name
      stop

   end function
! ################################################################

! ################################################################
   function getTimeZoneOffsetEarthClass(this, name) result(timezoneoffset)
      class(Earth_), intent(in) :: this
      character(*), intent(in) :: name
      type(TimeZoneList_) :: TimeZoneList
      real(real64) :: timezoneoffset
      integer(int32) :: i
      timezoneoffset = 0

      call TimeZoneList%init()
      do i = 1, size(TimeZoneList%ID, 1)
         if (index(TimeZoneList%ID(i)%all, name) == 0) then
            cycle
         else
            timezoneoffset = TimeZoneList%offset(i)
            return
         end if
      end do

      print *, "[ERROR] no such time-zone name as", name
      stop

   end function
! ################################################################

   function JP_Cartesian_Origin(ID) result(lat_lon)
      integer(int32), intent(in) :: ID
      real(real64) :: lat_lon(1:2)

      ! https://www.gsi.go.jp/LAW/heimencho.html

      select case (ID)
      case (1)
         lat_lon(1) = to_DecimalDegree([33., 0., 0.])
         lat_lon(2) = to_DecimalDegree([129., 30., 0.])

      case (2)
         lat_lon(1) = to_DecimalDegree([33., 0., 0.])
         lat_lon(2) = to_DecimalDegree([131., 0., 0.])

      case (3)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([132., 10., 0.])

      case (4)
         lat_lon(1) = to_DecimalDegree([33., 0., 0.])
         lat_lon(2) = to_DecimalDegree([133., 30., 0.])

      case (5)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([134., 20., 0.])

      case (6)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([136., 0., 0.])

      case (7)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([137., 10., 0.])

      case (8)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([138., 30., 0.])

      case (9)
         lat_lon(1) = to_DecimalDegree([36., 0., 0.])
         lat_lon(2) = to_DecimalDegree([139., 50., 0.])

      case (10)
         lat_lon(1) = to_DecimalDegree([40., 0., 0.])
         lat_lon(2) = to_DecimalDegree([140., 50., 0.])

      case (11)
         lat_lon(1) = to_DecimalDegree([44., 0., 0.])
         lat_lon(2) = to_DecimalDegree([140., 15., 0.])

      case (12)
         lat_lon(1) = to_DecimalDegree([44., 0., 0.])
         lat_lon(2) = to_DecimalDegree([142., 15., 0.])

      case (13)
         lat_lon(1) = to_DecimalDegree([44., 0., 0.])
         lat_lon(2) = to_DecimalDegree([144., 15., 0.])

      case (14)
         lat_lon(1) = to_DecimalDegree([26., 0., 0.])
         lat_lon(2) = to_DecimalDegree([142., 0., 0.])

      case (15)
         lat_lon(1) = to_DecimalDegree([26., 0., 0.])
         lat_lon(2) = to_DecimalDegree([127., 30., 0.])

      case (16)
         lat_lon(1) = to_DecimalDegree([26., 0., 0.])
         lat_lon(2) = to_DecimalDegree([124., 0., 0.])

      case (17)
         lat_lon(1) = to_DecimalDegree([26., 0., 0.])
         lat_lon(2) = to_DecimalDegree([131., 0., 0.])

      case (18)
         lat_lon(1) = to_DecimalDegree([20., 0., 0.])
         lat_lon(2) = to_DecimalDegree([136., 0., 0.])

      case (19)
         lat_lon(1) = to_DecimalDegree([26., 0., 0.])
         lat_lon(2) = to_DecimalDegree([154., 0., 0.])

      case default
         print *, "ERROR :: ID should be 1 <= id <= 19"
         stop
      end select
   end function

   function JP_to_World(lat, lon) result(ret)
      real(Real64), intent(in) :: lat, lon
      real(real64) :: ret(1:2)

      ret(2) = lon - lat*0.000046038 - lon*0.000083043 + 0.010040
      ret(2) = lat - lat*0.00010695 + lon*0.000017464 + 0.0046017

   end function

   function to_Cartesian(longitude, latitude, origin) result(xy)
      real(real64), intent(in) :: longitude, latitude
      real(real64), intent(in) :: origin(1:2)
      real(real64) :: xy(1:2)
      real(real64) :: x, y, x1, alp(1:5), A_bar, eta, gma, m0, a, F, phi, lambda &
                      , phi_0, lambda_0, t, lambda_c, lambda_s, sigma, S, A_(0:5), n, t_bar, rho, m, tau, xi
      integer(int32) :: j
      type(Math_) :: math

      ! GRS80
      !https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2xy/bl2xy.htm
      m0 = 0.9999d0
      a = 6377.397155d0*1000.0d0 ! m
      F = 299.152813d0
      n = 1.0d0/(2.0d0*F - 1.0d0)

      phi = radian(latitude)
      lambda = radian(longitude)

      phi_0 = radian(origin(1))
      lambda_0 = radian(origin(2))

      rho = 180.0d0/math%pi

      t = sinh(atanh(sin(phi)) - 2.0d0*sqrt(n)/(1.0d0 + n) &
               *atanh(2.0d0*sqrt(n)/(1.0d0 + n)*sin(phi)))
      t_bar = sqrt(1.0d0 + t*t)

      lambda_c = cos(lambda - lambda_0)
      lambda_s = sin(lambda - lambda_0)

      xi = atan(t/lambda_c)
      eta = atanh(lambda_s/t_bar)

      alp(1) = 1.0d0/2.0d0*n - 2.0d0/3.0d0*n*n + 5.0d0/16.0d0*n**3 &
               + 41.0d0/180.0d0*n**4 &
               - 127.0d0/288.00d0*n**5
      alp(2) = 13.0d0/48.0d0*n**2 - 3.0d0/5.0d0*n**3 &
               + 557.0d0/1440.0d0*n**4 + 281.0d0/630.0d0*n**5
      alp(3) = 61.0d0/240.0d0*n**3 - 103.0d0/140.0d0*n**4 + 15061/26880*n**5
      alp(4) = 49561.0d0/161280.0d0*n**4 - 179.0d0/168.0d0*n**5
      alp(5) = 34729.0d0/80640.0d0*n**5

      A_(0) = 1.0d0 + (n**2)/4.0d0 + (n**4)/64.0d0
      A_(1) = -3.0d0/2.0d0*(n - (n**3)/8.0d0 - (n**5)/64.0d0)
      A_(2) = 15.0d0/16.0d0*(n**2 - (n**4)/4.0d0)
      A_(3) = -35.0d0/48.0d0*(n**3 - 5.0d0/16.0d0*(n**5))
      A_(4) = 315.0d0/512.0d0*(n**4)
      A_(5) = -693.0d0/1280.0d0*(n**5)

      A_bar = m0*a/(1.0d0 + n)*A_(0)

      S = m0*a/(1.0d0 + n)*(A_(0)*phi_0)
      do j = 1, 5
         S = S + m0*a/(1.0d0 + n)*A_(j)*sin(2.0d0*j*phi_0)
      end do

      sigma = 1.0d0
      tau = 0.0d0
      do j = 1, 5
         sigma = sigma + 2.0d0*j*alp(j)*cos(2.0d0*j*xi)*cosh(2.0d0*j*eta)
         tau = tau + 2.0d0*j*alp(j)*sin(2.0d0*j*xi)*sinh(2.0d0*j*eta)
      end do

      gma = atan((tau*t_bar*lambda_c + sigma*t*lambda_s) &
                 /(sigma*t_bar*lambda_c - tau*t*lambda_s))

      m = A_bar/a*sqrt((sigma**2 + tau**2)/(t**2 + lambda_c**2)*(1.0d0 + &
                                                                 ((1.0d0 - n)/(1.0d0 + n)*tan(phi))**2))

      x = A_bar*(xi)
      do j = 1, 5
         x = x + A_bar*alp(j)*sin(2.0d0*j*xi)*cosh(2.0d0*j*eta)
      end do
      x = x - S

      y = A_bar*(eta)
      do j = 1, 5
         y = y + A_bar*alp(j)*cos(2.0d0*j*xi)*sinh(2.0d0*j*eta)
      end do

      ! JP
      xy(1) = x
      xy(2) = y

   end function

   pure function to_DecimalDegree(dd_mm_ss) result(ret)
      real(real32), intent(in) :: dd_mm_ss(1:3)
      real(real64) :: ret

      ret = dd_mm_ss(1) + dd_mm_ss(2)/60.0d0 + dd_mm_ss(3)/60.0d0/60.0d0
   end function

   function to_MomentTensor(phi, delta, lambda, M) result(ret)
      real(real64), intent(in):: phi, delta, lambda, M
      real(real64) :: ret(3, 3)

      ret(1, 1) = -M*(sin(delta)*cos(lambda)*sin(2*phi) &
                      + sin(2*delta)*sin(lambda)*sin(phi)*sin(phi))
      ret(1, 2) = M*(sin(delta)*cos(lambda)*cos(2*phi) &
                     + 0.50d0*sin(2*delta)*sin(lambda)*sin(2*phi))
      ret(1, 3) = -M*(cos(delta)*cos(lambda)*cos(phi) &
                      + cos(2*delta)*sin(lambda)*sin(phi))
      ret(2, 1) = ret(1, 2)
      ret(2, 2) = M*(sin(delta)*cos(lambda)*sin(2*phi) &
                     - sin(2*delta)*sin(lambda)*cos(phi)*cos(phi))
      ret(2, 3) = -M*(cos(delta)*cos(lambda)*sin(phi) &
                      - cos(2*delta)*sin(lambda)*cos(phi))
      ret(3, 1) = ret(1, 3)
      ret(3, 2) = ret(2, 3)
      ret(3, 3) = M*sin(2*delta)*sin(lambda)

   end function to_MomentTensor

end module
