c---  Generates a random density field smoothed over a comoving scale
c     using the BBKS power spectra of dark matter density fluctuations
c     g77 -o gaussian gaussian.f

c---  Made for Vietnamese Summer School students
c---  Aug. 22-25, 2018. 
c---  Changbom Park (Korea Institute for Advanced Study)

      parameter (ng=128)
      parameter (ng2=ng+2,ngh=ng/2,ngc=ng**3)
      external powint

      real*4 den(ng2,ng,ng),npow
      real*8 sqr,ave
      character*3 dark
      common/powsp/powh,omehs,dark
      common/tophat/rth

      
c---  Present matter density parameter
      omep = 0.3
c---  SCDM, HDM, or power-law power spectrum
      dark = 'pow'
      dark = 'hdm'
      dark = 'cdm'
c---  Slope of the primordial density power spectrum
      npow = 1.0
c---  Physical size of box in h^-1 Mpc
      size = 128
c---  Hubble parameter in 100 km/s/Mpc
      h = 0.7
c---  RMS matter density fluctuation over 8h^-1Mpc tophat scale
      sigma8 = 0.8
c---  Gaussian smoothing scale in h^-1Mpc
      gauss = 3

c---  Some various precalculations
      rng = real(ng)
      rngc = rng**3
      knyq = ng/2
      rknyq2 = real(knyq)**2
      rknyq = twopi*real(knyq)
      rkinv = 1./rknyq
      powh = 0.5*npow
      omehs = 1./(omep*h*size)
c---  8h^-1Mpc in units of box size 
      rth = 8./size
      pi = 3.141592654
      twopi = 2.*pi
      iseed = 1
      xxx = ran3(iiseed)

c---  Normalization of power spectrum
      call qsimp(powint,0.0001,100.,pamp)
c---  pamp = RMS of top-hat smoothed density fluctuation
      pamp = sigma8/sqrt(pamp)

c---  e-folding Gaussian smoothing length
      egauss = gauss*sqrt(2.)
      sfac   = -(pi*egauss/size)**2

c***  Generation of a Gaussian random density field
      pamp = pamp*rngc*sqrt(0.5)
      do 10 k=1,ng
      kz=k-1
      if(kz.gt.ngh)kz=kz-ng
      do 10 j=1,ng
      ky=j-1
      if(ky.gt.ngh)ky=ky-ng
      do 10 i=0,ngh
      xksq = i**2+ky**2+kz**2
      if(xksq.eq.0..or.xksq.ge.rknyq2)then
        den(2*i+1,j,k)=0.
        den(2*i+2,j,k)=0.
      else
        rk = twopi*sqrt(xksq)
        sfact = pamp*power(rk)*exp(xksq*sfac)
        den(2*i+1,j,k) = sfact*gasdev(0)
        den(2*i+2,j,k) = sfact*gasdev(0)
      end if
   10 continue

c---  Satisfy complex conjugation relation for a real array
      do 12 k=2,ngh
      den(1,1,ng+2-k) = den(1,1,k)
      den(2,1,ng+2-k) =-den(2,1,k)
      do 12 j=2,ngh
      den(1,ng+2-j,ng+2-k) = den(1,j,k)
      den(2,ng+2-j,ng+2-k) =-den(2,j,k)
   12 continue

      do 13 j=2,ngh
      den(1,ng+2-j,1) = den(1,j,1)
      den(2,ng+2-j,1) =-den(2,j,1)
   13 continue

      do 14 k=2,ngh
      do 14 j=ngh+2,ng
      den(1,ng+2-j,ng+2-k) = den(1,j,k)
      den(2,ng+2-j,ng+2-k) =-den(2,j,k)
   14 continue
      call src3ft(den,ng,ng,ng,ng2,ng,-1,ier)

      ave = 0d0
      sqr = 0d0
      do 30 k=1,ng
      do 30 j=1,ng
      do 30 i=1,ng
      ave = ave+dble(den(i,j,k))
      sqr = sqr+dble(den(i,j,k))**2
   30 continue
      ave = ave/rngc
      std = sqrt(sqr/rngc-aver**2)
      print *,iseed,' ',real(ave),real(std)

c---  Writing the density
      open(11,file='rho.smooth',form='unformatted')
      write(11)(((den(i,j,k),i=1,ng),j=1,ng),k=1,ng)
      close(11)

      stop
      end


      FUNCTION ran3(idum)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
C     REAL MBIG,MSEED,MZ
      REAL ran3,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
C     REAL mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3=mj*FAC
      return
      END

      FUNCTION gasdev(idum)
      INTEGER idum
      REAL gasdev
CU    USES ran3
      INTEGER iset
      REAL fac,gset,rsq,v1,v2,ran3
      SAVE iset,gset
      DATA iset/0/
      if (iset.eq.0) then
1       v1=2.*ran3(idum)-1.
        v2=2.*ran3(idum)-1.
        rsq=v1**2+v2**2
        if(rsq.ge.1..or.rsq.eq.0.)goto 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software +)-*1a311.

      SUBROUTINE qsimp(func,a,b,s)
      INTEGER JMAX
      REAL a,b,func,s,EPS
      EXTERNAL func
c     PARAMETER (EPS=1.e-5, JMAX=20)
      PARAMETER (EPS=1.e-5, JMAX=40)
CU    USES trapzd
      INTEGER j
      REAL os,ost,st
      ost=-1.e30
      os= -1.e30
      do 11 j=1,JMAX
        call trapzd(func,a,b,st,j)
        s=(4.*st-ost)/3.
        if (abs(s-os).lt.EPS*abs(os)) return
        os=s
        ost=st
11    continue
      pause 'too many steps in qsimp'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software +)-*1a311.

      SUBROUTINE trapzd(func,a,b,s,n)
      INTEGER n
      REAL a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software +)-*1a311.

      function power(rk)
c---  Calculates sqrt(Power spectrum) 
c---  Power spectrum consists of primordial scale-free power law and
c---  transfer function
      character(3) dark
      common/powsp/powh,omehs,dark

c***  Primordial scale-free power spectrum : power = sqrt(P(k))
c---  rk = 2*pi*(wave number) 
      power=rk**powh

c***  Transfer function (BBKS 1986, p60)
c---  omehs = 1/(omega*h*v[100])
      qq = rk*omehs
c---  Cold dark matter, adiabatic fluctuations with 3 flavors of
c---  relativistic neutrinos plus photons and negligible baryon
      if(dark.eq.'cdm')then

      if(qq.lt.0.002)then
          power=power*(1-1.17*qq+(1.351*qq)**2)*(1+3.89*qq
     &     +(16.1*qq)**2+(5.46*qq)**3+(6.71*qq)**4)**(-0.25)
      else
          power=power*alog(1+2.34*qq)/(2.34*qq)*(1+3.89*qq
     &     +(16.1*qq)**2+(5.46*qq)**3+(6.71*qq)**4)**(-0.25)
      end if

c---  Hot dark matter, adiabatic fluctuations with 2 species of
c---  relativiatic neutrinos
      else if(dark.eq.'hdm')then
          damp=2.6*qq
          power=power*exp(-0.16*damp-0.5*damp**2)/(1+1.6*qq
     &     +(4.*qq)**1.5+(0.92*qq)**2)
      end if
      return
      end


c---  Integrand for calculating (dM/M)**2 in top hat window
      function powint(xk)
      common/tophat/rth
      data twopi/6.28318531/,fopi/12.566371/

c---  xk = wave number
      rk = twopi*xk
      powint = fopi*xk*xk*power1(rk)**2
c---  Fourier transform of top-hat window function
      xk1   = rth*rk
      if(xk1.gt.0.02)then
      powint = powint*9*((sin(xk1)-xk1*cos(xk1))/xk1**3)**2
      else
      powint = powint*(1-xk1*xk1*0.1)**2
      end if
      return
      end


      function power1(rk)
c---  Calculates sqrt(Power spectrum) 
c---  Power spectrum consists of primordial scale-free power law and
c---  transfer function
      character(3) dark
      common/powsp/powh,omehs,dark

c***  Primordial scale-free power spectrum : power = sqrt(P(k))
c---  rk = 2*pi*(wave number) 
      power=rk**powh

c***  Transfer function (BBKS 1986, p60)
c---  omehs = 1/(omega*h*v[100])
      qq = rk*omehs
c---  Cold dark matter, adiabatic fluctuations with 3 flavors of
c---  relativistic neutrinos plus photons and negligible baryon
      if(dark.eq.'cdm')then
      if(qq.lt.0.002)then
          power=power*(1-1.17*qq+(1.351*qq)**2)*(1+3.89*qq
     &     +(16.1*qq)**2+(5.46*qq)**3+(6.71*qq)**4)**(-0.25)
      else
          power=power*alog(1+2.34*qq)/(2.34*qq)*(1+3.89*qq
     &     +(16.1*qq)**2+(5.46*qq)**3+(6.71*qq)**4)**(-0.25)
      end if

c---  Hot dark matter, adiabatic fluctuations with 2 species of
c---  relativiatic neutrinos
      else if(dark.eq.'hdm')then
          damp=2.6*qq
          power=power*exp(-0.16*damp-0.5*damp**2)/(1+1.6*qq
     &     +(4*qq)**1.5+(0.92*qq)**2)
      end if

      power1 = power
      return
      end

c--- src3ft -- real-to-complex 3-d FFT
c--- src3ft(x,l1,l2,l3,ldx,mdx,iopt,ier)
c---   * x = data array:  real *4 x(ldx,mdx,l3)
c---   * l1 = number of rows of data, power of 2
c---   * l2 = number of columns of data, power of 2
c---   * l3 = number of planes of data, power of 2
c---   * ldx = leading dimension of x, ldx >= l1+2
c---   * mdx = middle dimension of x, mdx >= l2
c---   * iopt = option flag:
c---      >= 0   compute forward transform (real to half complex)
c---      < 0    compute inverse transform (half complex to real)
c---   * ier = error flag, returns -1 if incorrectly called

      subroutine src3ft(x,l1,l2,l3,ldx,mdx,iopt,ier)
      real x(*)
      integer inn(3)

c--- check requirements on dimensions
      ier=0
      iflag1=0
      iflag2=0
      iflag3=0
      ipow=1
      do 10 i=1,20
	ipow=ipow*2
	if (l1.eq.ipow) iflag1=1
	if (l2.eq.ipow) iflag2=1
	if (l3.eq.ipow) iflag3=1
10    continue
      if ( (iflag1+iflag2+iflag3.ne.3) .or. (ldx.lt.l1+2) .or.
     >	   (mdx.lt.l2) )  then
	ier=-1
	return
      end if

      inn(1)=l1
      inn(2)=l2
      inn(3)=l3
      l1p2=l1+2

c--- four2() assumes that the input data densely fill the first part of the
c--- array, and it returns data in this form, so some reordering before
c--- and after the call to four2() is generally needed, the exception
c--- being a half-complex matrix that completely fills the array x()

c--- forward transform, starting with a real matrix
      if (iopt.ge.0) then
        do 20 i3=0,l3-1
	do 20 i2=0,l2-1
	do 20 i1=0,l1-1
 	  x(1+i1+i2*l1+i3*l1*l2) = x(1+i1+i2*ldx+i3*ldx*mdx)
20      continue
        call four2(x,inn,3,-1,0)
 	if ((ldx.ne.l1+2).or.(mdx.ne.l2)) then
 	  do 30 i3=l3-1,0,-1
 	  do 30 i2=l2-1,0,-1
 	  do 30 i1=l1+1,0,-1
 	    x(1+i1+i2*ldx+i3*ldx*mdx) = x(1+i1+i2*l1p2+i3*l1p2*l2)
30        continue
	end if

c--- inverse transform, starting with a half-complex matrix
      else
	if ((ldx.ne.l1+2).or.(mdx.ne.l2)) then
	  do 40 i3=0,l3-1
	  do 40 i2=0,l2-1
	  do 40 i1=0,l1+1
	    x(1+i1+i2*l1p2+i3*l1p2*l2) = x(1+i1+i2*ldx+i3*ldx*mdx)
40        continue
        end if
        call four2(x,inn,3,1,-1)
	xfac = 2./float(l1*l2*l3)
	do 50 i3=l3-1,0,-1
	do 50 i2=l2-1,0,-1
	do 50 i1=l1-1,0,-1
 	  x(1+i1+i2*ldx+i3*ldx*mdx) = x(1+i1+i2*l1+i3*l1*l2)*xfac
50      continue

      end if

      return
      end

c----------------------------------------------------------------------------
      subroutine four2 (data,n,ndim,isign,iform)
c
c     cooley-tukey fast fourier transform in usasi basic fortran.
c     multi-dimensional transform, each dimension a power of two,
c     complex or real data.
c     transform(k1,k2,...) = sum(data(j1,j2,...)*exp(isign*2*pi*sqrt(-1)
c     *((j1-1)*(k1-1)/n(1)+(j2-1)*(k2-1)/n(2)+...))), summed for all
c     j1 and k1 from 1 to n(1), j2 and k2 from 1 to n(2), etc. for all
c     ndim subscripts.  ndim must be positive and each n(idim) must be
c     a power of two.  isign is +1 or -1.  let ntot = n(1)*n(2)...
c     ...*n(ndim).  then a -1 transform followed by a +1 one
c     (or vice versa) returns ntot (ntot/2 if iform = 0 or
c     -1) times the original data.  iform = 1, 0 or -1, as data is
c     complex, real or the first half of a complex array.  transform
c     values are returned to array data.  they are complex, real or
c     the first half of a complex array, as iform = 1, -1 or 0.
c     the transform of a real array (iform = 0) dimensioned n(1) by n(2)
c     by ... will be returned in the same array, now considered to
c     be complex of dimensions n(1)/2+1 by n(2) by ....  note that if
c     iform = 0 or -1, n(1) must be even, and enough room must be
c     reserved.  the missing values may be obtained by complex conju-
c     gation.  the reverse transformation, of a half complex array
c     dimensioned n(1)/2+1 by n(2) by ..., is accomplished setting iform
c     to -1.  in the n array, n(1) must be the true n(1), not n(1)/2+1.
c     the transform will be real and returned to the input array.
c     running time is proportional to ntot*log2(ntot), rather than
c     the naive ntot**2.
c     original by norman brenner of mit lincoln laboratory, june 1968.
c     see-- ieee audio transactions (june 1967), special issue on fft.
      dimension data(*), n(1)
      ntot=1
      do 10 idim=1,ndim
 10   ntot=ntot*n(idim)
      if (iform) 70,20,20
 20   nrem=ntot
      do 60 idim=1,ndim
      nrem=nrem/n(idim)
      nprev=ntot/(n(idim)*nrem)
      ncurr=n(idim)
      if (idim-1+iform) 30,30,40
 30   ncurr=ncurr/2
 40   call bitrv (data,nprev,ncurr,nrem)
      call cool2 (data,nprev,ncurr,nrem,isign)
      if (idim-1+iform) 50,50,60
 50   call fixrl (data,n(1),nrem,isign,iform)
      ntot=(ntot/n(1))*(n(1)/2+1)

 60   continue
      return
 70   ntot=(ntot/n(1))*(n(1)/2+1)
      nrem=1
      do 100 jdim=1,ndim
      idim=ndim+1-jdim
      ncurr=n(idim)
      if (idim-1) 80,80,90
 80   ncurr=ncurr/2
      call fixrl (data,n(1),nrem,isign,iform)
      ntot=ntot/(n(1)/2+1)*n(1)
 90   nprev=ntot/(n(idim)*nrem)
      call bitrv (data,nprev,ncurr,nrem)
      call cool2 (data,nprev,ncurr,nrem,isign)
 100  nrem=nrem*n(idim)
      return
      end
c
c-------------------------------------------------------------------------
c
      subroutine bitrv (data,nprev,n,nrem)
c
c     shuffle the data by @bit reversal@.
c     dimension data(nprev,n,nrem)
c     data(i1,i2rev,i3) = data(i1,i2,i3), for all i1 from 1 to nprev,
c     all i2 from 1 to n (which must be a power of two), and all i3
c     from 1 to nrem, where i2rev-1 is the bitwise reversal of i2-1.
c     for example, n = 32, i2-1 = 10011 and i2rev-1 = 11001.
      dimension data(*)
      ip0=2
      ip1=ip0*nprev
      ip4=ip1*n
      ip5=ip4*nrem
      i4rev=1
      do 60 i4=1,ip4,ip1
      if (i4-i4rev) 10,30,30
 10   i1max=i4+ip1-ip0
      do 20 i1=i4,i1max,ip0
      do 20 i5=i1,ip5,ip4
      i5rev=i4rev+i5-i4
      tempr=data(i5)
      tempi=data(i5+1)
      data(i5)=data(i5rev)
      data(i5+1)=data(i5rev+1)
      data(i5rev)=tempr
 20   data(i5rev+1)=tempi
 30   ip2=ip4/2
 40   if (i4rev-ip2) 60,60,50
 50   i4rev=i4rev-ip2
      ip2=ip2/2
      if (ip2-ip1) 60,40,40
 60   i4rev=i4rev+ip2
      return
      end
c
c-----------------------------------------------------------------------
      subroutine cool2 (data,nprev,n,nrem,isign)
c
c     fourier transform of length n by the cooley-tukey algorithm.
c     bit-reversed to normal order.
c     dimension data(nprev,n,nrem)
c     complex data
c     data(i1,j2,i3) = sum(data(i1,i2,i3)*exp(isign*2*pi*i*((i2-1)*
c     (j2-1)/n))), summed over i2 = 1 to n for all i1 from 1 to nprev,
c     j2 from 1 to n and i3 from 1 to nrem.  n must be a power of two.
c     factoring n by fours saves about twenty five percent over factor-
c     ing by twos.
c     note--it is not necessary to rewrite this subroutine into complex
c     notation so long as the fortran compiler used stores real and
c     imaginary parts in adjacent storage locations.  it must also
c     store arrays with the first subscript increasing fastest.
      dimension data(*)
      twopi=6.283185307179586*float(isign)
      ip0=2
      ip1=ip0*nprev
      ip4=ip1*n
      ip5=ip4*nrem
      ip2=ip1
      npart=n
 10   if (npart-2) 50,30,20
 20   npart=npart/4
      go to 10
c     do a fourier transform of length two
 30   ip3=ip2*2
      do 40 i1=1,ip1,ip0
      do 40 i5=i1,ip5,ip3
      j0=i5
      j1=j0+ip2
      tempr=data(j1)
      tempi=data(j1+1)
      data(j1)=data(j0)-tempr
      data(j1+1)=data(j0+1)-tempi
      data(j0)=data(j0)+tempr
 40   data(j0+1)=data(j0+1)+tempi
      go to 140
c     do a fourier transform of length four (from bit reversed order)
 50   ip3=ip2*4
      theta=twopi/float(ip3/ip1)
      sinth=sin(theta/2.)
      wstpr=-2.*sinth*sinth
c     cos(theta)-1, for accuracy.
      wstpi=sin(theta)
      wr=1.
      wi=0.
      do 130 i2=1,ip2,ip1
      if (i2-1) 70,70,60
 60   w2r=wr*wr-wi*wi
      w2i=2.*wr*wi
      w3r=w2r*wr-w2i*wi
      w3i=w2r*wi+w2i*wr
 70   i1max=i2+ip1-ip0
      do 120 i1=i2,i1max,ip0
      do 120 i5=i1,ip5,ip3
      j0=i5
      j1=j0+ip2
      j2=j1+ip2
      j3=j2+ip2
      if (i2-1) 90,90,80
c     apply the phase shift factors
 80   tempr=data(j1)
      data(j1)=w2r*tempr-w2i*data(j1+1)
      data(j1+1)=w2r*data(j1+1)+w2i*tempr
      tempr=data(j2)
      data(j2)=wr*tempr-wi*data(j2+1)
      data(j2+1)=wr*data(j2+1)+wi*tempr
      tempr=data(j3)
      data(j3)=w3r*tempr-w3i*data(j3+1)
      data(j3+1)=w3r*data(j3+1)+w3i*tempr
 90   t0r=data(j0)+data(j1)
      t0i=data(j0+1)+data(j1+1)
      t1r=data(j0)-data(j1)
      t1i=data(j0+1)-data(j1+1)
      t2r=data(j2)+data(j3)
      t2i=data(j2+1)+data(j3+1)
      t3r=data(j2)-data(j3)
      t3i=data(j2+1)-data(j3+1)
      data(j0)=t0r+t2r
      data(j0+1)=t0i+t2i
      data(j2)=t0r-t2r
      data(j2+1)=t0i-t2i
      if (isign) 100,100,110
 100  t3r=-t3r
      t3i=-t3i
 110  data(j1)=t1r-t3i
      data(j1+1)=t1i+t3r
      data(j3)=t1r+t3i
 120  data(j3+1)=t1i-t3r
      tempr=wr
      wr=wstpr*tempr-wstpi*wi+tempr
 130  wi=wstpr*wi+wstpi*tempr+wi
 140  ip2=ip3
      if (ip3-ip4) 50,150,150
 150  return
      end
c
c-----------------------------------------------------------------------
      subroutine fixrl (data,n,nrem,isign,iform)
c
c     for iform = 0, convert the transform of a doubled-up real array,
c     considered complex, into its true transform.  supply only the
c     first half of the complex transform, as the second half has
c     conjugate symmetry.  for iform = -1, convert the first half
c     of the true transform into the transform of a doubled-up real
c     array.  n must be even.
c     using complex notation and subscripts starting at zero, the
c     transformation is--
c     dimension data(n,nrem)
c     zstp = exp(isign*2*pi*i/n)
c     do 10 i2=0,nrem-1
c     data(0,i2) = conj(data(0,i2))*(1+i)/(1-iform)
c     do 10 i1=1,n/4
c     z = (1+(2*iform+1)*i*zstp**i1)/2
c     i1cnj = n/2-i1
c     dif = data(i1,i2)-conj(data(i1cnj,i2))
c     temp = z*dif
c     data(i1,i2) = data(i1,i2)-temp
c     10  data(i1cnj,i2) = data(i1cnj,i2)+conj(temp)
c     if i1=i1cnj, the calculation for that value collapses into
c     a simple conjugation of data(i1,i2).
      dimension data(*)
      twopi=6.283185307179586*float(isign)
      ip0=2
      ip1=ip0*(n/2)
      ip2=ip1*nrem
      if (iform) 10,60,60
 10   j1=ip1+1
      data(2)=data(j1)
      if(nrem-1)60,60,15
 15   j1=j1+ip0
      i2min=ip1+1
      do 50 i2=i2min,ip2,ip1
      data(i2)=data(j1)
      j1=j1+ip0
      if (n-2) 40,40,20
 20   i1min=i2+ip0
      i1max=i2+ip1-ip0
      do 30 i1=i1min,i1max,ip0
      data(i1)=data(j1)
      data(i1+1)=data(j1+1)
 30   j1=j1+ip0
 40   data(i2+1)=data(j1)
 50   j1=j1+ip0
 60   do 80 i2=1,ip2,ip1
      tempr=data(i2)
      data(i2)=data(i2)+data(i2+1)
      data(i2+1)=tempr-data(i2+1)
      if (iform) 70,80,80
 70   data(i2)=data(i2)/2.
      data(i2+1)=data(i2+1)/2.
 80   continue
      if (n-2) 170,170,90
 90   theta=twopi/float(n)
      sinth=sin(theta/2.)
      zstpr=-2.*sinth*sinth
      zstpi=sin(theta)
      zr=(1.-zstpi)/2.
      zi=(1.+zstpr)/2.
      if (iform) 100,110,110
 100  zr=1.-zr
      zi=-zi
 110  i1min=ip0+1
      i1max=ip0*(n/4)+1
      do 160 i1=i1min,i1max,ip0
      do 150 i2=i1,ip2,ip1
      i2cnj=n+ip0-2*i1+i2
      if (i2-i2cnj) 140,120,120
 120  if (isign*(2*iform+1)) 130,150,150
 130  data(i2+1)=-data(i2+1)
      go to 150
 140  difr=data(i2)-data(i2cnj)
      difi=data(i2+1)+data(i2cnj+1)
      tempr=difr*zr-difi*zi
      tempi=difr*zi+difi*zr
      data(i2)=data(i2)-tempr
      data(i2+1)=data(i2+1)-tempi
      data(i2cnj)=data(i2cnj)+tempr
      data(i2cnj+1)=data(i2cnj+1)-tempi
 150  continue
      tempr=zr-.5
      zr=zstpr*tempr-zstpi*zi+zr
 160  zi=zstpr*zi+zstpi*tempr+zi
 170  if (iform) 240,180,180
 180  i2=ip2+1
      i1=i2
      j1=ip0*(n/2+1)*nrem+1
      go to 220
 190  data(j1)=data(i1)
      data(j1+1)=data(i1+1)
      i1=i1-ip0
      j1=j1-ip0
 200  if (i2-i1) 190,210,210
 210  data(j1)=data(i1)
      data(j1+1)=0.
 220  i2=i2-ip1
      j1=j1-ip0
      data(j1)=data(i2+1)
      data(j1+1)=0.
      i1=i1-ip0
      j1=j1-ip0
      if (i2-1) 230,230,200
 230  data(2)=0.
 240  return
      end

