!
! License-Identifier: GPL
!
! Copyright (C) 2016 The Yambo Team
!
! Authors (see AUTHORS file for details): AM DV
!
! headers
!
#include<y_memory.h>
!
subroutine plot_gnuplot(eval_only)
 !
 use pars,        ONLY:SP,schlen
 use units,       ONLY:BO2ANG
 use com,         ONLY:msg
 use YPPm,        ONLY:nr,v2plot,r_fixed,v2plot2D,l_norm_to_one,WF_multiplier,plot_dim,&
&                      plot_title,output_string,l_wavefunction,l_current,use_xcrysden,&
&                      plot_is_1D,plot_is_2D,l_free_hole,l_average_eh
 use LIVE_t,      ONLY:live_timing
 use y_memory_alloc
 !
 implicit none
 logical    :: eval_only
 ! 
 ! Work Space...
 !
 integer           :: i1,i2,i3,ir,nr2plot(2),dir2plot(2)
 character(schlen) :: titles(4)
 real(SP)          :: v_max,r_plot_center(3)
 !
 if (plot_dim==3) return
 !
 if (len_trim(plot_title)==0) then
   call live_timing('3D Merge',nr(1))
 else
   call live_timing('3D Merge of '//trim(plot_title),nr(1))
 endif
 !
 if (.not.eval_only) then
   !
   titles=" "
   !
   titles(:3)=(/'||a1 [A]','||a2 [A]','||a3 [A]'/)
   if (l_wavefunction) titles(4) = '|wf|^2 [1 at max]'
   if (l_current)      titles(4) = 'J [nAmpere]'
   !
 endif
 !
 ! SET CENTER OF THE PLOT
 !
 r_plot_center=0.
 if (.not.l_free_hole.and..not.l_average_eh) r_plot_center=r_fixed
 !
 ! DIMENSIONs
 !
 select case(plot_dim)
   !
   case(1)
     !
     do i1=1,3 
       if (plot_is_1D(i1)) nr2plot=nr(i1)
       if (plot_is_1D(i1)) dir2plot=i1
       if (plot_is_1D(i1)) call msg(output_string,"#",(/titles(1),titles(4)/),INDENT=0,USE_TABS=.true.)
     enddo
     call msg(output_string,'#')
     YAMBO_ALLOC(v2plot2D,(nr2plot(1),1))
     v2plot2D=0._SP
     ir=0
     do i1 = 0,nr(1)-1
       do i2 = 0,nr(2)-1
         do i3 = 0,nr(3)-1
           ir = 1 + i1 + i2*nr(1) + i3*nr(1)*nr(2)
           if (plot_is_1D(1)) v2plot2D(i1+1,1)=v2plot2D(i1+1,1)+v2plot(ir)
           if (plot_is_1D(2)) v2plot2D(i2+1,1)=v2plot2D(i2+1,1)+v2plot(ir)
           if (plot_is_1D(3)) v2plot2D(i3+1,1)=v2plot2D(i3+1,1)+v2plot(ir)
         enddo
       enddo
       call live_timing(steps=1)
     enddo
     v_max=maxval(v2plot2D(:,1))
     ! 
     if (l_norm_to_one) v2plot2D=v2plot2D/v_max  
     v2plot2D=v2plot2D*WF_multiplier
     !
     do i1=0,nr2plot(1)-1
       call msg(output_string,'',(/projection_(i1,dir2plot(1)),&
&                            v2plot2D(i1+1,1)/),INDENT=-2,USE_TABS=.true.)
     enddo
     !
     YAMBO_FREE(v2plot2D)
     !
   case(2)
     !
     if (.not.use_xcrysden) then
       if (plot_is_2D(1)) nr2plot=(/nr(1),nr(2)/)
       if (plot_is_2D(2)) nr2plot=(/nr(1),nr(3)/)
       if (plot_is_2D(3)) nr2plot=(/nr(2),nr(3)/)
     else
       if (plot_is_2D(1)) nr2plot=(/nr(1)+1,nr(2)+1/)
       if (plot_is_2D(2)) nr2plot=(/nr(1)+1,nr(3)+1/)
       if (plot_is_2D(3)) nr2plot=(/nr(2)+1,nr(3)+1/)
     endif
     if (plot_is_2D(1)) dir2plot=(/1,2/)
     if (plot_is_2D(2)) dir2plot=(/1,3/)
     if (plot_is_2D(3)) dir2plot=(/2,3/)
     YAMBO_ALLOC(v2plot2D,(nr2plot(1),nr2plot(2)))
     v2plot2D=0.
     ir=0
     do i1 = 0,nr(1)-1
       do i2 = 0,nr(2)-1
         do i3 = 0,nr(3)-1
           ir = 1 + i1 + i2*nr(1) + i3*nr(1)*nr(2)
           if (plot_is_2D(1)) v2plot2D(i1+1,i2+1)=v2plot2D(i1+1,i2+1)+v2plot(ir)
           if (plot_is_2D(2)) v2plot2D(i1+1,i3+1)=v2plot2D(i1+1,i3+1)+v2plot(ir)
           if (plot_is_2D(3)) v2plot2D(i2+1,i3+1)=v2plot2D(i2+1,i3+1)+v2plot(ir)
         enddo
       enddo
       call live_timing(steps=1)
     enddo
     if(use_xcrysden) then
       v2plot2D(nr2plot(1),:nr2plot(2)-1)=v2plot2D(1,:nr2plot(2)-1)
       v2plot2D(:nr2plot(1)-1,nr2plot(2))=v2plot2D(:nr2plot(1)-1,1)
       v2plot2D(nr2plot(1),nr2plot(2))=v2plot2D(1,1)
     endif
     !
     v_max=maxval(v2plot2D(:,:))
     if (l_norm_to_one) v2plot2D=v2plot2D/v_max  
     v2plot2D=v2plot2D*WF_multiplier
     !
     if (.not.eval_only) then
       if (plot_is_2D(1)) call msg(output_string,"#",&
&                             (/titles(1),titles(2),titles(4)/),INDENT=0,USE_TABS=.true.)
       if (plot_is_2D(2)) call msg(output_string,"#",&
&                             (/titles(1),titles(3),titles(4)/),INDENT=0,USE_TABS=.true.)
       if (plot_is_2D(3)) call msg(output_string,"#",&
&                             (/titles(2),titles(3),titles(4)/),INDENT=0,USE_TABS=.true.)
       call msg(output_string,"#")
       do i1=0,nr2plot(1)-1
         do i2=0,nr2plot(2)-1
           call msg(output_string,'',(/projection_(i1,dir2plot(1)),&
&                   projection_(i2,dir2plot(2)),v2plot2D(i1+1,i2+1)/),&
&                   INDENT=-2,USE_TABS=.true.)
         enddo
       enddo
     endif
     !
     if (.not.use_xcrysden) then
       !
       YAMBO_FREE(v2plot2D)
       !
     endif
     !
 end select
 !
 call live_timing()
 !
 contains
   !
   real(SP) function projection_(I,dir)
     use D_lattice,   ONLY:a  
     use FFT_m,       ONLY:fft_dim
     use vec_operate, ONLY:v_norm
     integer :: I,dir
     real(SP):: rv(3)
     !
     rv(:)=I*a(dir,:)/fft_dim(dir)-r_plot_center(:) 
     projection_=dot_product(rv,a(dir,:))/v_norm(a(dir,:))*BO2ANG
     !
   end function
   !
end subroutine