! #############################
module ReactorClass
   use MathClass
   implicit none

   integer(int32), parameter :: PF_MAX_REACTION_NUM = 20
   integer(int32), parameter :: PF_MAX_INPUT_NUM = 5
   integer(int32), parameter :: PF_MAX_OUTPUT_NUM = 5

   type :: Reactor_
      ! substances
      type(Real64Ptr_), allocatable :: substances(:)
      type(String_), allocatable :: substance_names(:)
      integer :: n_substances = 0

      ! reactions
      integer :: n_reaction = 0
      integer(int32):: input_list(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM) = 0
      integer(int32):: output_list(PF_MAX_REACTION_NUM, PF_MAX_OUTPUT_NUM) = 0
      real(real64)  :: reaction_order(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM + PF_MAX_OUTPUT_NUM) = 0
      real(real64)  :: mol_rate(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM + PF_MAX_OUTPUT_NUM) = 0
      real(real64)  :: constants(PF_MAX_REACTION_NUM) = 0

   contains
      procedure :: Init => InitReactor
      procedure :: put => putReactor
      procedure :: define => defineReactor
      procedure :: SearchSubstanceID => SearchSubstanceIDReactor
      procedure :: run => runReactor
   end type
contains

   ! ---------------------------
   subroutine InitReactor(this, max_num_substances)
      class(Reactor_), intent(inout) :: this
      integer(int32), optional, intent(in) :: max_num_substances
      integer(int32) :: n

      if (present(max_num_substances)) then
         n = max_num_substances
      else
         n = 100
      end if
      this%n_substances = 0

      if (allocated(this%substances)) then
         deallocate (this%substances)
      end if
      if (allocated(this%substance_names)) then
         deallocate (this%substance_names)
      end if

      allocate (this%substances(n))
      allocate (this%substance_names(n))

   end subroutine InitReactor
   ! ---------------------------

   ! ---------------------------
   subroutine putReactor(this, substance_name, substance)
      class(Reactor_), intent(inout) :: this
      character(*), intent(in) :: substance_name
      real(real64), target, intent(in) :: substance
      integer(int32) :: i, num
      num = this%n_substances + 1

      if (num > size(this%substances)) then
         print *, "ERROR :: putReactor num>size(this%substances ) "
      else
         this%substances(num)%ptr => substance
         this%substance_names(num) = substance_name
      end if
      this%n_substances = this%n_substances + 1

   end subroutine putReactor
   ! ---------------------------

   pure function SearchSubstanceIDReactor(this, name) result(ret)
      class(Reactor_), intent(in) :: this
      character(*), intent(in) :: name
      integer(int32) :: i, ret

      ret = 0
      do i = 1, this%n_substances
         if (this%substance_names(i)%all == name) then
            ret = i
            return
         end if
      end do

   end function

   ! ---------------------------
   subroutine defineReactor(this, input_list, output_list, constant, mol_rate, reaction_order)
      class(Reactor_), intent(inout) :: this
      type(String_), intent(in) :: input_list(:)
      type(String_), intent(in) :: output_list(:)
      real(real64), intent(in) :: constant ! reacton constant k
      real(real64), optional, intent(in) :: reaction_order(:), mol_rate(:)
      real(real64), allocatable :: mrate(:), rrate(:)
      integer(int32) :: i, j, n

      this%n_reaction = this%n_reaction + 1

      if (present(mol_rate)) then
         mrate = mol_rate
      else
         n = size(input_list) + size(output_list)
         allocate (mrate(n))
         mrate(:) = 1.0d0
      end if

      if (present(reaction_order)) then
         rrate = reaction_order
      else
         n = size(input_list) + size(output_list)
         allocate (rrate(n))
         rrate(:) = 1.0d0
      end if

      if (size(input_list) > PF_MAX_INPUT_NUM .or. size(output_list) > PF_MAX_OUTPUT_NUM) then
         print *, "ERROR :: defineReactor size(mrate) > PF_MAX_INPUT_NUM+PF_MAX_OUTPUT_NUM"
         stop
      end if
      ! register reaction
      ! INPUT
      do i = 1, size(input_list)
         ! search substance id
         this%input_list(this%n_reaction, i) = this%SearchSubstanceID(input_list(i)%all)
      end do

      do i = 1, size(output_list)
         ! search substance id
         this%output_list(this%n_reaction, i) = this%SearchSubstanceID(output_list(i)%all)
      end do

      this%reaction_order(this%n_reaction, 1:size(mrate)) = mrate(1:size(mrate))
      this%mol_rate(this%n_reaction, 1:size(mrate)) = mrate(1:size(mrate))

      this%constants(this%n_reaction) = constant

   end subroutine
   ! ---------------------------

   subroutine runReactor(this, dt)
      class(Reactor_), intent(inout) :: this
      real(real64), intent(in) :: dt
      integer(int32) :: i, j, num_input, num_output
      real(real64) :: A_B_ratio, xi, xi_n, coef, A0, B0, nA, nB
      real(real64) :: k1, k2, k3, k4

      ! run simulation
      if (this%n_reaction == 0) then
         print *, "ERROR :: runReactor this%n_reaction==0"
         stop
      elseif (this%n_reaction == 1) then
         ! single reaction equation
         num_input = 0
         do i = 1, size(this%input_list, 2)
            if (this%input_list(1, i) /= 0) then
               num_input = num_input + 1
            end if
         end do
         num_output = 0
         do i = 1, size(this%output_list, 2)
            if (this%output_list(1, i) /= 0) then
               num_output = num_output + 1
            end if
         end do
         if (num_input == 1 .and. num_output == 1) then
            ! x*A -> y*B reaction equation
            print *, "[ok] Reactor Started >> "
            A_B_ratio = this%reaction_order(1, 2)/this%reaction_order(1, 1)
            this%substances(this%input_list(1, 1))%ptr = this%substances(this%input_list(1, 1))%ptr*exp(this%constants(1)*dt)
       this%substances(this%output_list(1, 1))%ptr = this%substances(this%output_list(1, 1))%ptr*exp(A_B_ratio*this%constants(1)*dt)
         elseif (num_input == 2 .and. num_output == 1) then
            print *, "[ok] Reactor Started >> "
            ! x*A + y*B -> z*C reaction equation
            ! A + y/x*B -> z/x*C reaction equation
            ! by Runge-Kutta mehod
            ! Reacted A = A0 - xi
            ! d xi/dt = k*(A0-xi)**(nA) * (B0 - xi)**(nB)
            ! 反応比率も必要！！
            ! 反応速度式の次数は反応比率によらないが，
            ! 反応量は反応比率によるため．
            coef = this%constants(1)

            A0 = this%substances(this%input_list(1, 1))%ptr
            B0 = this%substances(this%input_list(1, 2))%ptr
            nA = this%reaction_order(1, 1)
            nB = this%reaction_order(1, 2)
            xi_n = 0.0d0
            k1 = coef*((A0 - xi_n)**(nA))*((B0 - xi_n)**(nB))
            k2 = coef*((A0 - (xi_n + dt*0.50d0*k1))**(nA)) &
                 *((B0 - (xi_n + dt*0.50d0*k1))**(nB))
            k3 = coef*((A0 - (xi_n + dt*0.50d0*k2))**(nA)) &
                 *((B0 - (xi_n + dt*0.50d0*k2))**(nB))
            k4 = coef*((A0 - (xi_n + dt*k3))**(nA)) &
                 *((B0 - (xi_n + dt*k3))**(nB))
            xi = xi_n + 1.0d0/6.0d0*dt*(k1 + k2 + k3 + k4)

            this%substances(this%input_list(1, 1))%ptr = A0 - xi
            this%substances(this%input_list(1, 2))%ptr = &
               B0 - dble(this%mol_rate(1, 2)/this%mol_rate(1, 1))*xi

            this%substances(this%output_list(1, 1))%ptr = &
               this%substances(this%output_list(1, 1))%ptr &
               + dble(this%mol_rate(1, 3)/this%mol_rate(1, 1))*xi
         else
            print *, "ERROR :: runReactor num_input/=1 .or. num_output/=1, now being implemented."
            stop
         end if

      else
         print *, "ERROR :: runReactor this%n_reaction>=2, now being implemented."
         stop
      end if
   end subroutine

end module ReactorClass
! #############################

