      subroutine sfsqd(natoms,fi,xyz,hkl,sfq)
C  Calculate structure factor from struct.dat file

      include 'arraysize.inc'
c...	passed variables
      integer natoms
      real fi(natoms_max)
c01Dec      real fi(nmax)
      double precision xyz(natoms_max,3)
c01Dec      double precision xyz(nmax,3)
      integer hkl(3)
      real sfq
      
c...	local variables
      integer i
C      complex sf
C      real angle,err
      double precision phase,pi,sfre,sfim

C      sf=cmplx(0.0,0.0)
      sfre=0.0d0
      sfim=0.0d0
      sfq=0.0
      pi=4.0d0*datan(1.0d0)

C      write(*,*) 'natoms ',natoms
      do i=1,natoms
C         write(*,'(a,i2)') 'atom #',i
        phase=0.0d0
        do j=1,3
C          write(*,'(f8.2,i4,f6.3,f8.2)') phase,hkl(j),xyz(i,j),
C     &       phase+dble(hkl(j))*xyz(i,j)
          phase=phase+dble(hkl(j))*xyz(i,j)
        enddo
C        write(*,*) 'added'
C would like to avoid any obvious offsets here so do a
C bit of work on the value of phase
C as a guide consider that the xyz will only be given to
C three decimal places          
C        err=.001*sqrt(float(hkl(1)**2+hkl(2)**2+hkl(3)**2))
C        angle=phase-aint(phase)    ! just look at the actual angle
C        if (abs(angle).lt.err) phase=0.0
C        if (abs(abs(angle)-0.5).lt.err) phase=0.5
C        sf=sf+fi(i)*exp(2.0*pi*cmplx(0,1)*phase)
        sfre=sfre+dble(fi(i))*dcos(2.0d0*pi*phase)
        sfim=sfim+dble(fi(i))*dsin(2.0d0*pi*phase)
C        if ((hkl(1).eq.-20).and.(hkl(2).eq.-20).and.
C     &      (hkl(3).eq.-19)) write(*,*) sf
C         write(*,*) 'next atom'
      enddo

C      sfq=abs(sf)**2.0
C      sfq=sqrt(sfre**2+sfim**2)
      sfq=sngl(sfre*sfre+sfim*sfim)
c0111      sfq=float(sfre*sfre+sfim*sfim)
C      write(*,*) (hkl(i),i=1,3),sfq
C      write(*,*) 'returning'
      return
      end

