program pot 

c Written by S.Sanvito Feb 2001.
c New version by V. M. Garcia-Suarez Jan 2006.

      implicit none

c Internal parameters
      integer maxp,maxort
      real ryd
      parameter ( maxp   = 10000000 )
      parameter ( maxort   = 10000 )
      parameter ( ryd   = 13.6058d0 )

c Internal variables
      integer lengfn
      character
     .  name*20,fname*25,task*10,recir*10,foutav*25
      integer
     .  i,i1,i2,i3,ip,is,j,l,pos,k,period,per,ir,izm,
     .  mesh(3),np,nspin,nt,indi
      real
     .  f1(maxp,2),f(maxp,2), f1max, f1min,
     .  ni,nj,nl,aver(maxort,2),norm,
     .  x,xmin,xmax,y,ymin,ymax,rmin2,rmax2,r2,aa(2)
      double precision
     .  cell(3,3),mcell(3,3)
      
c Read plot data
      open( unit=1, file='pot.dat',
     .      status='old', form='formatted' )
      read(1,*) name
      read(1,*) task
      read(1,*) recir 
      read(1,*) xmin,xmax 
      read(1,*) ymin,ymax 
      read(1,*) period
      close(1)

c Read potential

      if (task.eq.'vt' .or. task.eq.'vh' .or. task.eq.'rho') then
        if (task .eq. 'vt') then
	 lengfn = index(name,' ') - 1
         fname = name(1:lengfn)//'.VT'
         foutav = name(1:lengfn)//'.VT.dat'
        elseif (task .eq. 'vh') then
	 lengfn = index(name,' ') - 1
         fname = name(1:lengfn)//'.VH'
         foutav = name(1:lengfn)//'.VH.dat'
        elseif (task .eq. 'rho') then
	 lengfn = index(name,' ') - 1
         fname = name(1:lengfn)//'.RHO'
         foutav = name(1:lengfn)//'.RHO.dat'
	endif 
        open(unit=1,file=fname,status='old',form='unformatted')
        read(1) cell
        read(1) mesh, nspin
	write(6,'(a)') 'Parameters: '
	write(6,'(a,3(/,3f12.6))') ' cell =', cell
        write(6,'(a,3i6)') ' mesh =', mesh
        write(6,'(a,3i6)') ' spin =', nspin
        np = mesh(1) * mesh(2) * mesh(3)
	do is = 1,nspin
	 indi=0
	 do izm = 1,mesh(3)
	  do ir = 1,mesh(2)
             read(1) (f(indi+ip,is),ip=1,mesh(1))
	     indi = indi + mesh(1)
	  enddo   
	 enddo
	enddo
        close(1)
      elseif (task.eq.'vt-vh') then
        lengfn = index(name,' ') - 1
        fname = name(1:lengfn)//'.VH'
        foutav = name(1:lengfn)//'.VH.dat'
        open(unit=1,file=fname,status='old',form='unformatted')
        read(1) cell
        read(1) mesh, nspin
        np = mesh(1) * mesh(2) * mesh(3)
        do is = 1,nspin
          read(1) (f1(ip,is),ip=1,np)
        enddo
        close(1, status='keep')
          fname = name(1:lengfn)//'.VT'
        open(unit=1,file=fname,status='old',form='unformatted')
        read(1) cell
        read(1) mesh, nspin
        np = mesh(1) * mesh(2) * mesh(3)
        do is = 1,nspin
         read(1) (f(ip,is),ip=1,np)
        enddo
        do is = 1,nspin
         do ip=1,np
          f(ip,is)=f(ip,is)-f1(ip,1)
         enddo
        enddo
      endif
      
c Write functions to plot onto a file, and perform integration

      open( unit=7,file=foutav,status='unknown')
      close(unit=7,status='delete')
      open( unit=7,file=foutav,status='new')

c Average over the plane parallel to the junction
      	
      do is = 1,nspin
        do l=1,maxort
          aver(l,is)=0.d0
         enddo  
      enddo  

      do i=1,3 
       mcell(:,i)=cell(:,i)/mesh(i)
      enddo

      norm=0.0
      xmin=xmin/0.5291
      xmax=xmax/0.5291
      if(recir.eq.'a') then
        do is = 1,nspin
         do l=0,mesh(3)-1
          do j=0,mesh(2)-1
           do i=0,mesh(1)-1
            aver(l+1,is)=aver(l+1,is)+f(pos,is)
            norm=norm+1.0
           enddo
          enddo
         enddo
        enddo
      else if(recir.eq.'r') then
        ymin=ymin/0.5291
        ymax=ymax/0.5291
        do is = 1,nspin
         do l=0,mesh(3)-1
          do j=0,mesh(2)-1
           do i=0,mesh(1)-1
            pos= 1 + i + mesh(1)*j+mesh(1)*mesh(2)*l
            x = mcell(1,1)*i + mcell(1,2)*j + mcell(1,3)*l
            y = mcell(2,1)*i + mcell(2,2)*j + mcell(2,3)*l
            if (x.ge.xmin .and. x.le.xmax) then
            if (y.ge.ymin .and. y.le.ymax) then
             aver(l+1,is)=aver(l+1,is)+f(pos,is)
             norm=norm+1.0
            endif
            endif
           enddo
          enddo  
         enddo  
        enddo
      else if(recir.eq.'c') then
        rmin2=xmin**2
        rmax2=xmax**2
        do is = 1,nspin
         do l=0,mesh(3)-1
          do j=0,mesh(2)-1
           do i=0,mesh(1)-1
            pos= 1 + i + mesh(1)*j+mesh(1)*mesh(2)*l
            x = mcell(1,1)*i + mcell(1,2)*j + mcell(1,3)*l
            y = mcell(2,1)*i + mcell(2,2)*j + mcell(2,3)*l
            r2=x**2+y**2
            if (r2.ge.rmin2 .and. r2.le.rmax2) then
             aver(l+1,is)=aver(l+1,is)+f(pos,is)
             norm=norm+1.0
            endif
           enddo
          enddo  
         enddo  
        enddo
      else
        write(6,*)'Invalid rectangular-square option'
        stop
      endif

      do is = 1,nspin
       do l=1,maxort
        aver(l,is)=aver(l,is)/norm*ryd*mesh(3)*nspin
       enddo  
      enddo  

c Total average

      aa=0.
      write(6,'(a)')'Total average:'
      do is = 1,nspin
       do l=1,mesh(3)
        aa(is)=aa(is)+aver(l,is)
       enddo  
       write(6,'(a,i2,f12.6)')' ispin =',is,aa(is)/mesh(3)
      enddo  

c Write potential

      do per=1,period
       do l=2,mesh(3)
        write(7,'(1f8.3,6f20.8)') cell(3,3)/dble(mesh(3))
     .      *dble(l-1)*0.5291+cell(3,3)*0.5291*dble(per-1),
     .      (aver(l,is),is=1,nspin),
     .      aver(l,1)-aver(l,2)
       enddo
      enddo
      	
      close(3)
      close(8)

      end