module mode_da
    !This mode is used to calculate the dislocation analysis algorithm for both CG and atomistic regions

    use parameters
    use io
    use elements
    use neighbors
    
    implicit none 

    integer :: num_d_points
    integer, allocatable :: is_slipped(:,:,:) 
    real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)

    private :: write_xyz
    public 
    contains

    subroutine da(arg_pos)
        !This is the main calling subroutine for the dislocation analysis code
        integer, intent(out) :: arg_pos
        real(kind=dp), dimension(6) :: temp_box_bd

        integer :: ppos
        character(len=100) :: outfile

        !Now call the dislocation analysis code for the coarse-grained elements
        call parse_command(arg_pos)

        call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd)

        !Now allocate the necessary variables
        if(allocated(is_slipped)) deallocate(is_slipped)
        allocate(is_slipped(max_basisnum, max_ng_node, ele_num))

        if(allocated(disl_line)) deallocate(disl_line)
        num_d_points=0
        allocate(disl_line(3,1024))
        allocate(disl_data(1024))

        call cga_da

        !Now create the output file num and write out to xyz format
        ppos = scan(trim(infiles(1)),".", BACK= .true.)
        if ( ppos > 0 ) then 
            outfile = 'da_'//infiles(1)(1:ppos)//'vtk'
        else 
            outfile = 'da_'//infiles(1)//'.vtk'
        end if
        call write_da_vtk(outfile)        

        ppos = scan(trim(infiles(1)),".", BACK= .true.)
        if ( ppos > 0 ) then 
            outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
        else 
            outfile = 'da_'//infiles(1)//'.xyz'
        end if
        call write_da_xyz(outfile)
        
        return
    end subroutine da


    subroutine parse_command(arg_pos)
        !This subroutine parses the arguments for mode command
        integer, intent(out)  :: arg_pos
        integer :: i, arglen
        character(len = 100) :: textholder
        logical :: file_exists

        !Read the input file 
        call get_command_argument(2,textholder, arglen)
        if (arglen == 0) stop "Missing file for dislocation analysis"
        call get_in_file(textholder)

        arg_pos = 3
        return
    end subroutine parse_command

    subroutine cga_da

        integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3)

        integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6)
        real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm, &
                           max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, &
                            diff_mag(6), diffvec(3)
        real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:)
        logical :: is_edge

        !Initialize variables
        l = 0
        max_node_dist = 0

        !Now save the close and far neighbors of the nodes. This is done to attempt to figure out the character of the dislocation
        !by comparing how the burgers vector is distributed over the face.
        do j = 1, 6
            close_neighbors(1, 1, j) = cubic_faces(4, j)
            close_neighbors(2, 1, j) = cubic_faces(2, j)
            far_neighbor(1,j) = cubic_faces(3,j)

            close_neighbors(1, 2, j) = cubic_faces(1, j)
            close_neighbors(2, 2, j) = cubic_faces(3, j)
            far_neighbor(2,j) = cubic_faces(4,j)

            close_neighbors(1, 3, j) = cubic_faces(2, j)
            close_neighbors(2, 3, j) = cubic_faces(4, j)
            far_neighbor(3,j) = cubic_faces(1,j)

            close_neighbors(1, 4, j) = cubic_faces(3, j)
            close_neighbors(2, 4, j) = cubic_faces(1, j)
            far_neighbor(4,j) = cubic_faces(2,j)
        end do

        !First calculate all of the face centroids
        do i = 1, ele_num
            do j = 1, 6
                r(:) = 0.0_dp
                rmax= -Huge(1.0)
                rmin= Huge(1.0)
                do k = 1, 4
                    do ibasis = 1, basisnum(lat_ele(i))
                        r = r + r_node(:, ibasis, cubic_faces(k,j), i) 
                        rmax= max(rmax, r_node(:, ibasis, cubic_faces(k,j), i))
                        rmin= min(rmin, r_node(:, ibasis, cubic_faces(k,j), i))
                    end do
                end do  
                
                rnorm=norm2(rmax-rmin)
                max_node_dist = max(max_node_dist, rnorm)
                r = r/(basisnum(lat_ele(i))*4)

                !add the face centroids, the type, and map the elements faces to the face arrays
                l = l + 1
                face_centroids(:, l) = r 
                face_types(l) = j
                face_ele(l) = i


            end do
        end do

        !Now set the cut off distance as half the distance from opposite nodes in the largest element
        rc = 0.5*max_node_dist
        call calc_NN(l, face_centroids(:,1:l), rc)

        !Now loop overall the faces and make sure the nearest neighbor is the opposite type. If it isn't than we dscard it 
        is_slipped = 0
        do i = 1, l
            nei = nn(i)             
            
            !Skip if it's 0 or the closest face belongs to the same element
            if ((nei == 0).or.(face_ele(i) == face_ele(nei))) cycle
            

            !Check the face types, the way that the faces are laid out in the cubic_faces array face 1's opposite is 6 and 
            ! face 2's opposite is 5 and etc
            vnode = 0
            if(face_types(i) == (7-face_types(nei))) then 
                vnorm = 0
                do j = 1, 4
                    do ibasis = 1, basisnum(lat_ele(face_ele(i)))
                        !Compute the vectors between all nodes at the face.
                        vnode(:,ibasis,j) = r_node(:,ibasis, cubic_faces(j,face_types(i)), face_ele(i)) - &
                                       r_node(:,ibasis, cubic_faces(j,face_types(nei)), face_ele(nei))
                    end do
                end do

                !Now calculate the maximum difference between nodes
                vnorm = 0
                do j = 1, 4
                    do k=j,4
                        v1=vnode(:,1,j) - vnode(:,1,k)
                        v1norm=norm2(v1)
                        if (vnorm < v1norm) then
                            vnorm = v1norm
                            diffvec= v1
                        end if 
                    end do
                end do
                !Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it
                !as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter
                !I think 0.5 works ok though.
                if (vnorm>0.5_dp) then 
                    print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
                             " with neighbor ", face_ele(nei), " with max displacement of ", vnorm
                    l=0
                    do j = 1, 4
                        is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1

                        !This portion of the code is used to determine what the character of the dislocation most likely is
                        !effectively how this works is we look for the two nodes with the most similar burgers vector.
                        !The vector between these two nodes is close to the line direction and as a result we can probably estimate
                        !if it's close to edge, screw, or if it's pretty fairly mixed.
                        do k = j+1, 4
                            l=l+1
                            diff_mag(l) = norm2(vnode(:, 1, j)-vnode(:,1,k))
                            diff_pair(1,l)=cubic_faces(j, face_types(i))
                            diff_pair(2,l)=cubic_faces(k, face_types(i))
                        end do
                    end do
                    minindex = minloc(diff_mag,1)
                    !Now figure out if the min difference between nodes is associated with a 112 direction    
                    is_edge = .false.
                    do j = 1,6
                        if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then
                            is_edge=.true.
                            exit
                        else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then
                            is_edge=.true.
                            exit
                        end if
                    end do

                    if(is_edge) then 
                        print *, 'Dislocation has primarily edge character'
                    else
                        print *, 'Dislocation has primarily screw character'
                    end if

                    if( i < nei) then 
                        num_d_points=num_d_points+2
                        if(num_d_points > size(disl_line,2)) then 
                            allocate(disl_line_temp(3,size(disl_line,2)+1024))
                            disl_line_temp=0.0_dp
                            disl_line_temp(:, 1:size(disl_line,2))=disl_line
                            call move_alloc(disl_line_temp, disl_line)
                        end if
                        if(num_d_points/2 > size(disl_data)) then 
                            allocate(disl_data_temp(size(disl_data)+1024))
                            disl_data_temp=0
                            disl_data_temp(1:size(disl_data)) = disl_data
                            call move_alloc(disl_data_temp, disl_data)
                        end if
                        print *, num_d_points/2, vnorm
                        disl_data(num_d_points/2)=vnorm
                        disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i))
                        disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i))
                    end if
                end if
            end if
        end do

    end subroutine cga_da

    subroutine write_da_vtk(outfile)
        character (len=*), intent(in) :: outfile

        integer :: i 
        open(unit=11, file=outfile, action='write', status='replace', position='rewind')

        write(11, '(a)') '# vtk DataFile Version 4.0.1'
        write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis'
        write(11, '(a)') 'ASCII'
        write(11, '(a)') 'DATASET UNSTRUCTURED_GRID'
        write(11, *) 'POINTS', num_d_points, 'float'
       
        !Write the points to the file
        do i = 1, num_d_points
            write(11, '(3f23.15)') disl_line(:, i)
        end do

        write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3
        do i=1, num_d_points, 2
            write(11,'(3i16)') 2, i-1, i
        end do
        
        write(11,*) 'CELL_TYPES', num_d_points/2
        do i=1, num_d_points, 2
            write(11,'(i16)') 3
        end do

        write(11,*) 'CELL_DATA', num_d_points/2
        write(11, '(a)') 'SCALARS burgers_vec, float'
        write(11, '(a)') 'LOOKUP_TABLE default'
        do i=1, num_d_points, 2
            write(11, '(f23.15)') disl_data((i+1)/2)
        end do

        close(11)

    end subroutine write_da_vtk

    subroutine write_da_xyz(outfile)
        !This subroutine write the element positions to a .xyz file and marks whether they are slipped or not 
        character(len=*), intent(in) :: outfile
        integer :: i, ibasis, inod, outn
        
        open(unit=11, file=outfile, action='write', status='replace', position='rewind')

        !Write number of node_atoms
        write(11, '(i16)') node_atoms
        !Write comment
        write(11, '(a)') "is_slipped x y z"
        !Write nodal positions
        outn = 0
        do i = 1, ele_num
            do inod = 1, ng_node(lat_ele(i))
                do ibasis = 1, basisnum(lat_ele(i))
                        write(11, '(1i16, 3f23.15)') is_slipped(ibasis,inod,i), r_node(:,ibasis,inod,i)
                        outn = outn + 1
                end do 
            end do 
        end do 

        if(outn /= node_atoms) then 
            print *, "outn", outn, " doesn't equal node_atoms ", node_atoms
        end if
        
        close(11)
    end subroutine write_da_xyz
end module