module opt_slip_plane
    use parameters 
    use elements
    use functions 
    use subroutines

    implicit none

    integer :: sdim
    real(kind=dp) :: spos
    logical :: efill
    
    public 
    contains

    subroutine run_slip_plane(arg_pos)
        !Main calling function for the slip_plane option
        integer, intent(inout) :: arg_pos 

        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis,&
                    starting_anum
                
        integer, allocatable :: slip_eles(:), temp_int(:)
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
                         maxp, minp, vel_interp(3, max_basisnum*max_esize**3)
        real(kind=dp), allocatable :: vel_tmp(:,:)

        integer :: type_interp(max_basisnum*max_esize**3)
        logical :: lat_points(max_esize,max_esize, max_esize)
    

        print *, '---------------------Option Slip_Plane----------------------'

        !Initialize variables
        efill = .false.
        slip_enum = 0
        old_atom_num = atom_num

        !!Parse the argument
        call parse(arg_pos)


        !If we are running the efill code then we have to initiate some variables
        if(efill) then 
            new_ele_num = 0
        end if
        allocate(slip_eles(1024))
        !Now loop over all elements, find which ones intersect 
        do ie = 1, ele_num
            if( (spos < maxval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. &
                (spos > minval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie)))) then 
                slip_enum = slip_enum + 1
                if (slip_enum > size(slip_eles)) then 
                    allocate(temp_int(size(slip_eles)+1024))
                    temp_int(1:size(slip_eles)) = slip_eles
                    temp_int(size(slip_eles)+1:)  = 0
                    call move_alloc(temp_int, slip_eles)
                end if
                slip_eles(slip_enum) = ie
                
                !If we aren't efilling then just refine the element
                if(.not.efill) then 
                    starting_anum=atom_num
                    call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
                    do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
                        call apply_periodic(r_interp(:,ia))
                        call add_atom(0, type_interp(ia), r_interp(:,ia))
                    end do

                    if(allocated(vel_atom)) then  
                        call interpolate_vel(type_ele(ie), size_ele(ie), lat_ele(ie), vel_node(:,:,:,ie), vel_interp)
                        if(size(vel_atom,2) < size(type_atom)) then 
                            allocate(vel_tmp(3, size(type_atom)))
                            vel_tmp=0.0_dp
                            vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
                            call move_alloc(vel_tmp, vel_atom)
                        end if
                        do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
                            vel_atom(:, starting_anum+ia) = vel_interp(:, ia)
                        end do
                    end if
                !If we are efilling then the code is slightly more complex
                else
                    !First populate the lat points array 
                    lat_points(1:size_ele(ie),1:size_ele(ie), 1:size_ele(ie)) = .true.

                    !Now start trying to remesh the region, leaving the slip plane as a discontinuity
                    esize = size_ele(ie) - 2
                    nump_ele = size_ele(ie)**3
                    do while(esize > min_efillsize)
                        if(nump_ele < esize**3) then 
                            esize = esize - 2
                        else
                            ele = cubic_cell*(esize -1) 
                            do o = 1, size_ele(ie) - esize
                                do n = 1, size_ele(ie) - esize
                                    latloop:do m = 1, size_ele(ie) - esize
                                        do inod = 1, ng_node(lat_ele(ie))
                                            vlat = ele(:,inod) + (/ m, n, o /)
                                            if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop
                                            call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
                                        end do

                                        !Check to make sure all lattice points exist for the current element
                                        if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop 
                                        !Check to see if the plane intersects this element if not then add it 
                                        maxp = maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
                                        minp = minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
                                        if(.not.(spos < maxp).and.(spos > minp))then
                                            nump_ele = nump_ele - esize**3 
                                            lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false.
                                            call add_element(0, type_ele(ie), esize, lat_ele(ie),  rfill) 
                                            new_ele_num = new_ele_num + 1
                                        end if
                                    end do latloop
                                end do
                            end do
                        end if
                        esize= esize-2
                    end do
                    ! Now add the leftover lattice points as atoms
                    do o = 1, size_ele(ie)
                        do n = 1, size_ele(ie)
                            do m = 1, size_ele(ie)
                                if(lat_points(m,n,o)) then 
                                    call get_interp_pos(m,n,o, ie, ratom(:,:))
                                    do ibasis = 1, basisnum(lat_ele(ie))
                                        call apply_periodic(r_atom(:,ibasis))
                                        call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis))
                                    end do
                                end if
                            end do
                        end do
                    end do
                end if
            end if
        end do 

        !Once we finish adding elements delete the old ones
        call delete_elements(slip_enum, slip_eles(1:slip_enum))

        !Output data
        if(.not.efill) then 
            print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms"
        else 
            print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms and ", new_ele_num, " elements"
        end if

    end subroutine run_slip_plane 

    subroutine parse(arg_pos)
        !This subroutine parses the input arguments to the mode
        integer, intent(inout) :: arg_pos 

        integer :: arglen
        character(len = 100) :: textholder

        !First read the dimension
        arg_pos = arg_pos +1
        call get_command_argument(arg_pos,textholder, arglen) 
        if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
        
        !Check to make sure that the dimension is correct
        select case(trim(adjustl(textholder))) 
        case('x','X')
            sdim = 1
        case('y','Y')
            sdim = 2
        case('z','Z')
            sdim = 3
        case default
            print *, "Error: dimension ", trim(adjustl(textholder)), " is not accepted. Please select from x, y, or z"
        end select 

        !Now parse the position of the slip plane
        arg_pos = arg_pos + 1
        call get_command_argument(arg_pos, textholder, arglen)
        if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
        call parse_pos(sdim, textholder, spos)

        !Now check to see if efill was passed 
        arg_pos = arg_pos + 1
        if(.not.(arg_pos > command_argument_count())) then 
            call get_command_argument(arg_pos, textholder, arglen)
            if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
            if(trim(adjustl(textholder)) == "efill") then 
                arg_pos = arg_pos +1 
                efill = .true.
            end if
        end if
    end subroutine parse
end module opt_slip_plane
