module OW_hadam1

use OW_type
use OW_parms_areal
use OW_fftw_areal

implicit none

contains

subroutine owhadamard(CHD)

!complex		:: CCHD(:,:)
real		:: CHD(:,:)
integer,allocatable :: HD(:,:),temp(:,:)
!complex,allocatable :: tfld(:)
integer		:: i,i1,i2,j1,j2,hij

      allocate(HD(2*sog%xm%n*sou%w%n,2*sog%z%n))

      do i1=1,2*sog%xm%n*sou%w%n
         do j1=1,2*sog%z%n
             j2=j1
             j2=iand(ieor(ishft(j2,1),ieor(i1,sog%xm%n)),j2)
!             j2=iand(ieor(ishft(j2,1),ieor(i1,2*sog%z%n)),j2)
             do i2=1,2*sog%z%n-1
                j2=ieor(ishft(j2,-1),iand(j2,1))
             end do
             hij=(-1)**j2
!             if(mod(i1,2).eq.0) then
                HD(i1,j1)=hij
!             else
!                HD(i1,j1)=-hij
!             end if 
         end do
      end do

!      CHD=float(HD)

      allocate(temp(2*sog%xm%n*sou%w%n,2*sog%z%n));temp=0.

      do i=1,sog%xm%n
         temp(2*(i-1)*sou%w%n+1:(2*i-1)*sou%w%n,:)=HD((i-1)*sou%w%n+1:i*sou%w%n,:)
         temp((2*i-1)*sou%w%n+1:2*i*sou%w%n,:)=HD((i-1+sog%xm%n)*sou%w%n+1:(i+sog%xm%n)*sou%w%n,:)
      end do

      CHD=float(temp)

!      do i=1,2*sog%xm%n*sou%w%n
!         write(3,8) HD(i,:)
!      end do
!8     format(1024I3)
!
!      deallocate(temp)
!      allocate(tfld(2*sog%z%n))
!
!      call owfftw(tfld)
!
!      do i=1,2*sog%xm%n*sou%w%n
!         tfld=0.
!         tfld=cmplx(HD(i,:),0)
!         call sfftw_execute(planF)
!!          CHD(i,:)=atan(imag(tfld)*real(tfld)/(real(tfld)*real(tfld)+tiny(1.)))
!!          CHD(i,:)=CHD(i,:)/maxval(CHD(i,:))
!          CCHD(i,:)=tfld/(sqrt(1.*sog%z%n+tiny(1.)))
!      end do
!
!      deallocate(tfld,HD)
!
      do i=1,2*sog%xm%n*sou%w%n
         write(2,9) CHD(i,1:512)
!         write(5,9) real(CCHD(i,:))
!         write(4,9) imag(CCHD(i,:))
      end do
9 format(512F10.5)

end subroutine

end module
