! a f 'g77 -O3 fouriertest4.f'

      program fouriertest4
      implicit none

      integer numoftiles,m,ipix,jpix
      real*8 ra,dec,pi,thetag,phig,wavelength
      parameter (m=9, ipix=2**m, jpix=2**m,pi=3.141592654)
      character*80 tilefilename,skyfilename
      
      numoftiles=500
      tilefilename='n500r750uniform.dat'
      skyfilename='kangaroo.txt'
	  ra=pi/3.0
	  dec=pi/6.0
	  thetag=1.112
	  phig=0.0
	  wavelength=1.5

      call go(ipix,jpix,numoftiles,tilefilename,skyfilename,
     +ra,dec,thetag,phig,wavelength)
      stop
      end
c----
      subroutine go(ipix,jpix,numoftiles,tilefilename,skyfilename,
     +ra,dec,thetag,phig,wavelength)
      implicit none
      integer ipix,jpix,numoftiles,i,j
      character*80 tilefilename,tempname,skyfilename
      real*8 real1(ipix,jpix), imag1(ipix,jpix), ra, dec, thetag,
     +phig,wavelength,tilearray(500,2)

      call cleararray(ipix,jpix,real1)
      call cleararray(ipix,jpix,imag1)

!---- Test Number 1 -----

c$$$      do i=1,ipix
c$$$         do j=1,jpix
c$$$            real1(i,j)=1.0
c$$$         enddo
c$$$      enddo

!---- Test Number 2 -----

c      real1(2044,2050)=1.0

!---- Test Number 3 -----

      call createcircle(ipix,real1)
      
      tempname='original.pgm'
      call rescalenplot(ipix,jpix,real1,tempname)

!---- Test Number 4 -----

c$$$      call readimage(ipix,jpix,skyfilename,real1)
c$$$      
c$$$      tempname='original.pgm'
c$$$      call rescalenplot(ipix,jpix,real1,tempname)

!---- Test Number 5 -----
c$$$      call readtilepositions(tilefilename,tilearray)
c$$$      call calcFFTW(ipix,real1,ra,dec,tilearray,thetag,phig,wavelength)
c$$$      tempname='points.pgm'
c$$$      call rescalenplot(ipix,jpix,real1,tempname)

!---- FFT section ----
      call FFTpack(ipix,jpix,real1,imag1,+1)
c      call FFTpack(ipix,jpix,real1,imag1,-1)
c      call FFTpackmax(ipix,real1)

      tempname='realprofile.dat'
      call createprofiledata(ipix,jpix,real1,tempname)
      tempname='imagprofile.dat'
      call createprofiledata(ipix,jpix,imag1,tempname)
      tempname='real.pgm'
      call rescalenplot(ipix,jpix,real1,tempname)
      tempname='imag.pgm'
      call rescalenplot(ipix,jpix,imag1,tempname)

      return
      end
c----
      subroutine readtilepositions(fname,tilearray)
      implicit none
      integer u,numoftiles,i
      real*8 tilearray(500,2)
      parameter (u=2)
      character*80 fname
      
      open (u, file=fname, status='old')
      
      read (u,*) numoftiles
      
      do i=1, numoftiles
         read(u,*) tilearray(i,1), tilearray(i,2)
      enddo
      return
      end
c----

      subroutine calcFFTW(ipix,FFTW,ra,dec,tilearray,thetag,phig,wavelength)
      implicit none
      integer i,j,ipix
      real*8 FFTW(ipix,ipix),ra,dec,tilearray(500,2),thetazero,phizero,
     +thetag,phig,wavelength,uvpositions(250*499,2)

      call calcpointdirec(thetag,phig,ra,dec,thetazero,phizero)
      call calcuvpositions(tilearray,uvpositions,wavelength,
     +thetazero,phizero,phig)
      call createdeltgrid(FFTW,ipix,ipix,uvpositions)
      return
      end
c----
      subroutine createdeltgrid(deltgrid,ipix,jpix,uvpositions)
      implicit none
      integer i,j,ipix,jpix
      real*8 uvpositions(250*499,2),ax,ay,pi,deltgrid(ipix,jpix)
      parameter (pi=3.141592654)
      
      do i=1,ipix
         do j=1, jpix
            deltgrid(i,j)=0
         enddo
      enddo
      
      do i=1,250*499
         ax=uvpositions(i,1)
         ay=uvpositions(i,2)
         deltgrid(int((ax)+real(ipix)/2),int(-(ay)+real(jpix)/2))=1
         deltgrid(int((-ax)+real(ipix)/2),int(-(-ay)+real(jpix)/2))=1
      enddo
      
      return
      end
c----
      subroutine calcuvpositions(tilearray,uvpositions,wavelength,
     +thetazero,phizero,phig)
      implicit none
      integer i,j,k
      real*8 tilearray(500,2), uvpositions(250*499,2),wavelength,pi,
     +xtemp,ytemp,thetazero,phizero,phig
      parameter (pi=3.141592654)

      wavelength=wavelength

      k=1
      do i=1,500
         do j=i+1,500
            xtemp=(tilearray(i,1)-tilearray(j,1))
            ytemp=(tilearray(i,2)-tilearray(j,2))
            uvpositions(k,1)=xtemp*cos(phig)-ytemp*sin(phig)
            uvpositions(k,2)=xtemp*sin(phig)+ytemp*cos(phig)
            k=k+1
         enddo
      enddo
      
      do k=1,250*499
         xtemp=uvpositions(k,1)
         ytemp=uvpositions(k,2)
         uvpositions(k,1)=xtemp*cos(phizero)+ytemp*sin(phizero)
         uvpositions(k,2)=-xtemp*sin(phizero)+ytemp*cos(phizero)
      enddo
      
      do k=1,250*499
         uvpositions(k,1)=(cos(thetazero))*uvpositions(k,1)
      enddo
      
      do k=1,250*499
         xtemp=uvpositions(k,1)
         ytemp=uvpositions(k,2)
         uvpositions(k,1)=xtemp*cos(phizero)-ytemp*sin(phizero)
         uvpositions(k,2)=xtemp*sin(phizero)+ytemp*cos(phizero)
      enddo
      
      do k=1,250*499
         xtemp=uvpositions(k,1)
         ytemp=uvpositions(k,2)
         uvpositions(k,1)=(xtemp)/wavelength
         uvpositions(k,2)=(ytemp)/wavelength
      enddo

      wavelength=wavelength      
      return
      end
c----
      subroutine calcpointdirec(thetag,phig,ra,dec,thetazero,phizero)
      implicit none
      real*8 thetag,phig,ra,dec,thetazero,phizero,alt,lat,pi
      parameter (pi=3.141592654)
      
      lat=(pi/2.0)-thetag

      alt=((sin(dec))*(sin(lat)))+((cos(dec))*(cos(lat))*(cos(phig-ra)))
      alt=asin(alt)

      thetazero=(pi/2.0)-alt
      
      phizero=((sin(dec))-(sin(lat))*(sin(alt)))/((cos(lat))*(cos(alt)))
      if (phizero .LT. -1.0) then
         phizero=pi
      elseif (phizero .GT. 1.0) then
         phizero=0.0
      else
         phizero=acos(phizero)
      endif
      
      if (sin(phig-ra) .LT. 0.0) then
         phizero=phizero
      else
         phizero=2*pi-phizero
      endif
      return
      end
c----
      subroutine createcircle(ipix,array)
      implicit none
      integer i,j,m,ipix
      real*8 r, array(0:ipix-1,0:ipix-1)
      call cleararray(ipix,ipix,array)
      do i=0,ipix-1
         do j=0,ipix-1
            r=sqrt(((i-ipix/2.d0)**2)+((j-ipix/2.d0)**2))
             if (r .LT. ipix/4.d0) then
               array(i,j)=1
            else
c               array(i,j)=exp(-(r-ipix/4.0)**2)
            endif
         enddo
      enddo
      return
      end
c----
      subroutine rswap(a,b)
      implicit none
      real*8 a,b,c
      c = a
      a = b
      b = c
      return
      end
c----
      subroutine map2rdata(nmap,map,rdata)
      implicit none
      integer nmap, i, j, k
      real*8 map(nmap,nmap), rdata(2*nmap*nmap)
      ! Put the data in the FFT format:
      k = 1
      do i=1,nmap
         do j=1,nmap	  	
	    rdata(k) 	= map(i,j) ! Real part = our image
	    k = k+1
	    rdata(k)    = 0.    ! Imaginary part = 0
	    k = k+1
         end do
      end do	
      return
      end
c----
      subroutine rschwopp(n,rdata)
      ! Translates the origin by n/2 in both directions, modulo n
      implicit none
      integer n,n2,i,j,k,l,kk,ll
      real*8 rdata(2*n*n)
      n2 = n/2
      do i=1,n2
         do j=1,n2
            kk          = n*(i-1) + j
            ll          = n*(n2+i-1) + (n2+j)
            k           = 2*kk-1
            l           = 2*ll-1
            call rswap(rdata(k),rdata(l))
            k           = k+1
            l           = l+1
            call rswap(rdata(k),rdata(l))
            kk          = n*(i-1) + (n2+j)
            ll          = n*(n2+i-1) + j
            k           = 2*kk-1
            l           = 2*ll-1
            call rswap(rdata(k),rdata(l))
            k           = k+1
            l           = l+1
            call rswap(rdata(k),rdata(l))
         end do
      end do
      return
      end
c----
      subroutine schwopp(n,map)
      ! Translates the origin by n/2 in both directions, modulo n
      implicit none
      integer n
      real*8 map(n,n), rdata(2*4096*4096)
      call map2rdata(n,map,rdata)
      call rschwopp(n,rdata)
      call rdata2map(n,rdata,map)
      return
      end
c----
      subroutine fftmap(nmap,map,rdata)
      implicit none
      integer nmap
      real*8 map(nmap,nmap), rdata(2*nmap*nmap), norm
      integer isign, ndim, nn(2), i
      ndim 	= 2
      nn(1)	= nmap
      nn(2)	= nmap
      call map2rdata(nmap,map,rdata) ! Put the data in the FFT format
      print *,'Fourier transforming....'
      isign 	= 1
      call fourn(rdata,nn,ndim,isign) ! Fourier transform
      norm = 1./nmap
      do i=1,2*nmap*nmap
         rdata(i) = norm*rdata(i)
      end do
      return
      end
c----
      subroutine rdata2map(nmap,rdata,map)
      implicit none
      integer nmap, i, j, k
      real*8 map(nmap,nmap), rdata(2*nmap*nmap)
      ! Unpack the data from FFT format:
      k = 1
      do i=1,nmap
         do j=1,nmap	  	
	    map(i,j) = rdata(k)	! Want real part only.
	    k = k+2
         end do
      end do	
      return
      end
c----
      subroutine readimage(ipix,jpix,fname,imagearray)
      implicit none
      integer ipix,jpix,i,j
      character*80 fname
      real*8 imagearray(ipix,jpix)
      
      open (1, FILE=fname, status='old')

      do i=1, ipix
         read(1,*) (imagearray(j,i), j=1,jpix)
      enddo
      return
      end
c----
      subroutine createprofiledata(ipix,jpix,array,filename)
      implicit none
      character*80 filename
      integer i,j,ipix,jpix,band
      real*8 array(ipix,jpix),integrated

      open (1,file=filename)

      do j=1,jpix
         write(1,*) array(j,1+ipix/2) 
      enddo
      return
      end
c----
      subroutine rescalenplot(ipix,jpix,image,fname)
      implicit none
      integer ipix,jpix
      integer i,j,imagearray(ipix,jpix)
      real*8 image(ipix,jpix),temp(ipix,jpix)
      character*80 fname

      do i=1, ipix
         do j=1, jpix
            temp(i,j)=image(i,j)
         enddo
      enddo
      call rescaleimagearrayreal(temp,ipix,jpix)
      call integerize(temp,imagearray,ipix,jpix)
      call writepgm(ipix,jpix,imagearray,fname)
      return
      end
c----
      subroutine rescaleimagearrayreal(uvgrid,ipix,jpix)
      implicit none
      integer ipix,jpix,i,j
      real*8 uvgrid(ipix,jpix),maxint,minint,maxabs
      maxint=uvgrid(1,1)
      minint=uvgrid(1,1)
      
      do i=1,ipix
         do j=1,jpix
            if (uvgrid(i,j) .GT. maxint) then
               maxint = uvgrid(i,j)
            endif
            if (uvgrid(i,j). LT. minint) then
               minint = uvgrid(i,j)
            endif
         enddo
      enddo
      
      if (abs(maxint) .GT. abs(minint)) then
         maxabs = abs(maxint)
      else
         maxabs = abs(minint)
      endif
      
      do i=1, ipix
         do j=1, jpix
            uvgrid(i,j)=127.0+((uvgrid(i,j)/maxabs)*127.0)
         enddo
      enddo
      return
      end
c----
      subroutine integerize(realarray,intarray,ipix,jpix)
      implicit none
      
      integer ipix,jpix,i,j
      integer intarray(ipix,jpix)
      real*8 realarray(ipix,jpix)
      
      do i=1, ipix
         do j=1, jpix
            intarray(i,j)=realarray(i,j)
         enddo
      enddo
      return
      end
c----
      subroutine writepgm(nx,ny,image,fname)
      implicit none
!       The array image contains nx*ny pixel values between 0 and 255.
      integer nx, ny, i, image(nx*ny)
      character*80 fname
!       call compress_image(nx*ny,image)
      open(2,file=fname)
      write(2,'(a2)') 'P5'
      write(2,'(1i4,1x,1i4)') nx, ny
      write(2,'(a3)') '255'
      write(2,'(9999999a)') (char(image(i)),i=1,(nx*ny))
      return
      end
c----
      subroutine FFTpackmax(ipix,real1)
      implicit none
      integer ipix
      real*8 real1(ipix,ipix),rdata(2*ipix*ipix)
      call schwopp(ipix,real1)
      call fftmap(ipix,real1,rdata)
      call rdata2map(ipix,rdata,real1)
      call schwopp(ipix,real1)
      return
      end
c----
      subroutine FFTpack(ipix,jpix,real1,imag1,sign)
      implicit none
      integer i,j,ipix,jpix,sign,dimvect(2)
      real*8 real1(ipix,jpix),imag1(ipix,jpix),complx(2*ipix,jpix)

      dimvect(1)=ipix
      dimvect(2)=jpix

      call conv2FTformat(ipix,jpix,real1,imag1,complx)
      call reorderFmatrix(ipix,jpix,complx)
      call check_symmetry(ipix,complx)
      call fourn(complx,dimvect,2,sign)
      call reorderFmatrix(ipix,jpix,complx)
      call conv2normform(ipix,jpix,complx,real1,imag1)
      
      return
      end
c----
      SUBROUTINE fourn(data,nn,ndim,isign)
      INTEGER isign,ndim,nn(ndim)
      REAL*8 data(*)
      INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,
     *k2,n,nprev,nrem,ntot
      REAL*8 tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      ntot=1
      do 11 idim=1,ndim
        ntot=ntot*nn(idim)
11    continue
      nprev=1
      do 18 idim=1,ndim
        n=nn(idim)
        nrem=ntot/(n*nprev)
        ip1=2*nprev
        ip2=ip1*n
        ip3=ip2*nrem
        i2rev=1
        do 14 i2=1,ip2,ip1
          if(i2.lt.i2rev)then
            do 13 i1=i2,i2+ip1-2,2
              do 12 i3=i1,ip3,ip2
                i3rev=i2rev+i3-i2
                tempr=data(i3)
                tempi=data(i3+1)
                data(i3)=data(i3rev)
                data(i3+1)=data(i3rev+1)
                data(i3rev)=tempr
                data(i3rev+1)=tempi
12            continue
13          continue
          endif
          ibit=ip2/2
1         if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
            i2rev=i2rev-ibit
            ibit=ibit/2
          goto 1
          endif
          i2rev=i2rev+ibit
14      continue
        ifp1=ip1
2       if(ifp1.lt.ip2)then
          ifp2=2*ifp1
          theta=isign*6.28318530717959d0/(ifp2/ip1)
          wpr=-2.d0*sin(0.5d0*theta)**2
          wpi=sin(theta)
          wr=1.d0
          wi=0.d0
          do 17 i3=1,ifp1,ip1
            do 16 i1=i3,i3+ip1-2,2
              do 15 i2=i1,ip3,ifp2
                k1=i2
                k2=k1+ifp1
                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                data(k2)=data(k1)-tempr
                data(k2+1)=data(k1+1)-tempi
                data(k1)=data(k1)+tempr
                data(k1+1)=data(k1+1)+tempi
15            continue
16          continue
            wtemp=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
17        continue
          ifp1=ifp2
        goto 2
        endif
        nprev=n*nprev
18    continue
      return
      END
c----
      subroutine reorderFmatrix(ipix,jpix,complex)
      implicit none
      integer i,j,ipix,jpix
      real*8 complex(2*ipix,jpix),tempcomplex(2*ipix,jpix)
      
      do i=1,2*ipix
         do j=1,jpix
            tempcomplex(i,j)=complex(i,j)
         enddo
      enddo

      do i=1,ipix
         do j=1,jpix/2
            complex(ipix+i,j+jpix/2)=tempcomplex(i,j)
         enddo
      enddo

      do i=1,ipix
         do j=1,jpix/2
            complex(i,j)=tempcomplex(i+ipix,j+jpix/2)
         enddo
      enddo

      do i=1,ipix
         do j=1,jpix/2
            complex(i+ipix,j)=tempcomplex(i,j+jpix/2)
         enddo
      enddo

      do i=1,ipix
         do j=1,jpix/2
            complex(i,j+jpix/2)=tempcomplex(i+ipix,j)
         enddo
      enddo
      return
      end
c----

	subroutine check_symmetry(ipix,c)
	implicit none
	integer ipix, i, j, i1, j1
	complex c(0:ipix-1,0:ipix-1), c1
	do i=1,ipix-1
	  i1 = ipix-1
	  if (i1.eq.ipix) i1 = 0
	  do j=1,ipix-1
	    j1 = ipix - 1
	    if (j1.eq.ipix) j1 = 0
	    c1 =  conjg(c(i1,j1))
	    if (c(i,j).ne.c1) print *,'Asymmetry: ',i,j,c(i,j),c1 
	    !print *,'Symmetry: ',i,j,c(i,j),c1 
	  end do
	end do
	print *,'Symmetry test done.'
	return
	end
	
      subroutine cleararray(ipix,jpix,array)
      implicit none
      integer i,j,ipix,jpix
      real*8 array(ipix,jpix)
      
      do i=1,ipix
         do j=1,jpix
            array(i,j)=0
         enddo
      enddo
      return
      end
c----
      subroutine conv2FTformat(ipix,jpix,real1,imag1,complx)
      implicit none
      integer i,j,ipix,jpix
      real*8 real1(ipix,jpix),imag1(ipix,jpix)
      real*8 complx(2*ipix,jpix)

      do i=1,ipix
         do j=1,jpix
            complx(2*i-1,j)=real1(i,j)
            complx(2*i,j)=imag1(i,j)
         enddo
      enddo
      return
      end
c----
      subroutine conv2normform(ipix,jpix,complex,real1,imag1)
      implicit none
      integer i,j,ipix,jpix
      real*8 complex(2*ipix,jpix),real1(ipix,jpix),imag1(ipix,jpix)

      do i=1,ipix
         do j=1,jpix
            real1(i,j)=complex(2*i-1,j)
            imag1(i,j)=complex(2*i,j)
         enddo
      enddo
      return
      end