module Positions
    use Const
    implicit none
    
    type Pos
        real(mpc), allocatable :: mem(:)
        contains 
            procedure :: x  => pos_x
            procedure :: p  => pos_p
            procedure :: q  => pos_q
            
            procedure :: Mx => pos_set_x
            procedure :: Mp => pos_set_p
            procedure :: Mq => pos_set_q
           
            procedure :: assign_pos
            generic :: assignment(=) => assign_pos
!             final :: clear_pos
    end type Pos

    interface Pos
        module procedure :: init_pos_pq, init_pos_x
    end interface Pos

    interface operator(+) 
        module procedure :: add_pos_pos
    end interface
    interface operator(*)
        module procedure :: mult_real_pos, mult_pos_real
    end interface
!     interface assignment(=)
!         module procedure :: pos_assign
!     end interface
    

contains

    function init_pos_x(x) result(c)
        real(mpc), intent(in), dimension(:) :: x

        type(Pos) :: c
        integer :: memsize

        memsize = size(x)
        allocate(c%mem(memsize))
        c%mem = x
    end function init_pos_x
    
    
    function init_pos_pq(p,q) result(c)
        real(mpc), intent(in), dimension(:) :: p, q

        type(Pos) :: c
        integer :: memsize, dimen

        if (size(p) /= size(q)) stop "incompatible dimensions"
        
        memsize = size(p) + size(q)
        dimen = size(p)
        allocate(c%mem(memsize))

        c%mem(1:dimen) = p
        c%mem(dimen+1:memsize) = q
    end function init_pos_pq

!=============================
!         getters 
!=============================
    pure function pos_x(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem)))

        v = self%mem
    end function pos_x
    

    pure function pos_p(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem))/2)
        integer :: half

        half = (size(self%mem))/2
        v = self%mem(1:half)
    end function pos_p
    
    
    pure function pos_q(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem))/2)
        integer :: half

        half = (size(self%mem))/2
        v = self%mem(half+1:2*half)
    end function pos_q


!=================================
!         setters
!=================================
    
    subroutine pos_set_x(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)

        self%mem(1:size(self%mem)) = v
    end subroutine pos_set_x
    

    subroutine pos_set_p(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)
        integer :: half

        half = (size(self%mem))/2
        self%mem(1:half) = v
    end subroutine pos_set_p
    
    
    subroutine pos_set_q(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)
        integer :: half

        half = (size(self%mem))/2
        self%mem(half+1:2*half) = v
    end subroutine pos_set_q
    
    
    subroutine clear_pos(self)
        type(Pos), intent(inout) :: self
        if (allocated(self%mem)) deallocate(self%mem)
    end subroutine clear_pos

! **************************
! overloaded operators
! **************************
    
    function add_pos_pos(p1, p2) result(p)
        type(Pos), intent(in) :: p1, p2
        type(Pos) :: p
        p = Pos(p1%p() + p2%p(), p1%q() + p2%q())
    end function add_pos_pos
    
    
    function mult_real_pos(r, p) result(p1)
        type(Pos), intent(in) :: p
        real(mpc), intent(in) :: r
        type(Pos) :: p1
        p1 = Pos(r * p%p(), r * p%q())
    end function mult_real_pos 
    
    
    function mult_pos_real(p, r) result(p1)
        type(pos), intent(in) :: p
        real(mpc), intent(in) :: r
        type(pos) :: p1
        p1 = mult_real_pos(r, p)
    end function mult_pos_real
    
    
    subroutine assign_pos(p, p1)
        class(Pos), intent(out) :: p
        class(Pos), intent(in) :: p1
        allocate(p%mem(size(p1%mem)), source=p1%mem) 
    end subroutine assign_pos
    
    subroutine test_pos()
        type(Pos) :: p1, p2, p
        p1 = Pos((/1.0_mpc,1.0_mpc/), (/0.1_mpc, 0.0_mpc/))
        p2 = Pos((/2.0_mpc,2.0_mpc/), (/0.2_mpc, 0.0_mpc/))
        write(*,*) "q =", p1%q()
        write(*,*) "x =", p1%x()
        call p1%Mx([1.0_mpc, 1.1_mpc, 0.0_mpc, 0.1_mpc])
        p = p1 + p2
        write(*,*) "p' =", p%p(), "q' =", p%q()
        p = p1
        write(*,*) "p' =", p%p(), "q' =", p%q()
        call p1%Mp([0.0_mpc, 0.0_mpc])
        write(*,*) "p' =", p1%p(), "p' =", p%p()
    end subroutine test_pos

end module
