program bessim
c
c   BESSIM, "Statistical Visualization" of turbulence spectra (such as 
c     measured by BES on TFTR).  See bessim.doc for info.
c
c   by Greg Hammett, Princeton Plasma Physics Lab, September 22, 1992.
c
c This is free software and comes with ABSOLUTELY NO WARRANTY.
c This software is used as a research tool, and is frequently undergoing
c changes.  Beware, some of the options may no longer work properly.
c
c If you think you have found a bug in an important part of the program,
c please notify me before making public statements about it.  I can be 
c reached by email at hammett@princeton.edu.
c
c GWH 30-Jan-1993 modified to calculate the 2-point correlation function,
c both by Fourier transforming the spectrum, and by manually averaging 
c together many pairs from a specific realization.
c

      implicit none

c try to set up general procedure for 2D FFT's for later use:

      integer mx                ! # of x points
      integer my                ! # of y points
      parameter (mx=128,my=128)
      real den(mx,my) ! density as a function of x and y
      real x(mx), y(my)

      real trigsx(2*mx), trigsy(2*my)
      integer ifaxx(13), ifaxy(13)

c for later generalization to 2D FFT's in a 3D space, DENK and DEN
c become larger while WORK stays the same size:

      complex denk(mx,my)
      complex work(mx,my)
      complex amp(mx,my)  ! random complex coefficients
      real damp,damp2,dt
      integer mframe,nframe

      real denkx(mx),denky(my)
      complex corxfft(mx),coryfft(my)
      real rcorxfft(mx),rcoryfft(my)
      real kx(mx), ky(my)
      real x0,y0,sigx,sigy, rmaxr, rmaxi, ampr,ampi
      complex ampc
      real dmin,dmax,davg
      real d2,d4,kurtosis,corrx,corry
      integer idelta,jdelta,npts
      integer i,j,ii,jj, iget, jget, nco

      real areverse(my,mx) ! HDF stupidly reverses index order
      real pi

      integer iret
      integer shape(2)

      real phase

      integer mcomodes
      parameter (mcomodes=100) ! Max # of coherent modes
      real ncomodes    ! actual # of coherent modes
      real xposco(mcomodes)  ! x position of coherent modes
      real yposco(mcomodes)  ! y position of coherent modes
      real ampco(mcomodes)   ! amplitude of coherent modes

      character amp_alg*1 ! gaussian or fixed amplitude algorithm
      character ph_alg*1  ! random or fixed phase algorithm
      character answer*1

c functions:
      real*8 g05caf ! NAG uniform random number generator
      real*8 g05ddf ! NAG Gaussian random number generator
      real fpintrp ! my linear interpolator

      integer dssdast,dsadata

c Sr(kr) and Sth(ktheta) spectra from BES measurements:
c k is in units of cm**(-1):
      integer mkr,mkth
      parameter (mkr=10,mkth=11)
      real akr(mkr),pkr(mkr)
      real akth(mkth),pkth(mkth)

c Sr(kr):
      data akr / 0., .247, .323, .579, .766, 1.04, 1.2, 1.4, 2.0, 2.8/
      data pkr/ 5.7, 5.0, 4.5, 3.0, 2.0, 1.0, 0.5, .171, .137, 0.0/
c Sth(kth):
      data akth/ 0., .38, .766, .851, 1.1, 1.3, 1.6, 1.9, 2.15, 2.51,
     >          3.0/
      data pkth/ .73, 1.0, 2.85, 3.2, 3.46, 3.0, 2.0, 1.19, .69, .27,
     >           0./

c This data is from Ray Fonck's 1992 EPS paper.  I tried to be careful 
c when digitizing this.  Later, when I was doing Fourier transforms of 
c these spectra, I discovered that they give 2-point correlation 
c functions which are about 10-15% narrower than they should be.  
c I need to be more accurate, perhaps especially in the high k 
c components or with a finer k-grid for the digitization.  
c But I think this is good enough for now.

c Choosing a 19 cm (7.5 inches) size box allows us to print this 
c out on paper at actual size, but I eventually decided on 40 cm
c to give a nice number of modes on the page at once:

      x0=40.0   ! x (radial) Length of box (cm)
      y0=40.0   ! y (poloidal) Length of box (cm)

      write(6,*) 
     > 'Bessim, January 30, 1993 Version, by Greg Hammett, PPPL.'
      write(6,*) 
     > 'This is free software and comes with ABSOLUTELY NO WARRANTY.'
      write(6,*) ' '
      write(6,*) 'Type a "," to get the default (in parentheses):'
      write(6,*) 'Radial (x) box length (cm) (40.0)?'
      read(5,*) x0
      write(6,*) 'Poloidal (y) box length (cm) (40.0)?'
      read(5,*) y0

      ph_alg='r'
      write(6,*) 'random or fixed phases (r or f) (r)?'
      read(5,*) ph_alg

      if(ph_alg .eq. 'r') then
         amp_alg='f'
         write(6,*) 'random or fixed amplitudes (r or f) (f)?'
         read(5,*) amp_alg
      else
         ncomodes=40
         write(6,*) ' number of coherent modes (0=none) (40)?'
         read(5,*) ncomodes
      endif

      if (ph_alg .eq. 'r' .or. ncomodes .gt. 0) then

         answer='r'
         write(6,*) 
     >       'random or fixed initial random number (r or f) (r)?'
         read(5,*) answer

         if(answer .eq. 'r') then
         call  g05ccf   ! initialize random number generator from clock, or
         else
         call g05cbf(0) ! initialize repeatable random sequence.
         endif

      endif

      mframe=1
c      write(6,*) 'Number of frames to generate (10)?'
c      read(5,*) mframe

      dt=10
c      write(6,*) 
c     >   '(time step dt)/(decorrelation time) between frames (.2)?'
c      read(5,*) dt

c****************************************************************

      pi=2*abs(acos(0.0))

      call cftfax(mx,ifaxx,trigsx)   !initialize x FFT's
      call cftfax(my,ifaxy,trigsy)   !initialize y FFT's

c generate the grids:
      do i=1,mx
         x(i)=x0*(i-1)/(1.0*(mx-1))
         ii=i-1
         if(ii .gt. mx/2) ii=ii-mx
         kx(i)=ii*2*pi/x0
      enddo
      do j=1,my
         y(j)=y0*(j-1)/(1.0*(my-1))
         jj=j-1
         if(jj .gt. my/2) jj=jj-my
         ky(j)=jj*2*pi/y0
      enddo

c generate the spectrum denk(kx,ky) = denkx(kx)*denky(ky)
c
      do i=1,mx
         denkx(i)=fpintrp(abs(kx(i)),pkr,akr,mkr)
      enddo
      do j=1,my
         denky(j)=fpintrp(abs(ky(j)),pkth,akth,mkth)
      enddo

c Gaussian tests:
      sigx=3.  ! x correlation length (cm)
      sigy=3.  ! y correlation length (cm)
c      do i=1,mx
c         denkx(i)=exp(-kx(i)**2*sigx**2)
c      enddo
c      do i=1,my
c         denky(i)=exp(-ky(i)**2*sigy**2)*(1+8*ky(i)**2*sigy**2)
c Gaussian test:
c         denky(i)=exp(-ky(i)**2*sigy**2)
c      enddo

      write(6,*) 'kx grid=',kx
      write(6,*) 'spectrum vs. kx =',denkx
      write(6,*) 'ky grid=',ky
      write(6,*) 'spectrum vs. ky =',denky

c calculate exact 2-point correlation function from the FFT of the power
c spectra:
      do i=1,mx
         corxfft(i)=denkx(i)
       enddo
      do j=1,my
         coryfft(j)=denky(j)
      enddo

      call cfft99(corxfft,work,trigsx,ifaxx,1,1,mx,1,1)
      call cfft99(coryfft,work,trigsy,ifaxy,1,1,my,1,1)

      do i=1,mx
         rcorxfft(i)=corxfft(i)/real(corxfft(1))
       enddo
      do j=1,my
         rcoryfft(j)=coryfft(j)/real(coryfft(1))
      enddo

      if(ncomodes .gt. 0) then
c
c place the coherent eddies at randomly chosen positions, though
c quantize the positions to lie on grid points to try try to reduce
c Gibb's oscilations.
c Choose the amplitude to be +-1, eddies of either polarity.
c
         do i=1,ncomodes
            xposco(i)=x0/mx*float(int(mx*g05caf()))
            yposco(i)=y0/my*float(int(my*g05caf()))
            ampco(i)=(-1)**(i)
         enddo
      endif

      if(amp_alg .eq. 'r') then

c generate an initial set of random, complex Amplitude coefficients:
c we generate the x and y components independently.  This guarantees
c that the phase will be randomly distributed uniformly from 0 to 2*pi.
      do i=1,mx
         do j=1,my
            ampr=g05ddf(0.d0,1.d0)   ! NAG's Gaussian random number generator
            ampi=g05ddf(0.d0,1.d0)
            amp(i,j)=cmplx(ampr,ampi)

         enddo
      enddo

      else

c generate a set of initial complex coefficients all with unit
c amplitude and either random or fixed phases:

      do i=1,mx
         do j=1,my

            if(ph_alg .eq. 'r') then
               phase=2*pi*g05caf() ! NAG's uniform random number generator
               ampr=cos(phase)
               ampi=sin(phase)
               amp(i,j)=cmplx(ampr,ampi)
            else
               ampc=0.0
               do nco=1,ncomodes
                  ampc=ampc+ampco(nco)*exp(-(0.0,1.0)*
     >                 (kx(i)*xposco(nco) + ky(j)*yposco(nco)))
                  amp(i,j)=ampc
               enddo
            endif

         enddo
      enddo

      endif

c Loop over a number of frames:

      do nframe=1,mframe

c generate one realization of this spectrum, using the set
c of random complex amplitudes in the array amp.  

c look at just the perturbed part, so turn off the kx=0, ky=0 component:
         amp(1,1)=0.0 

         do i=1,mx
         do j=1,my
            denk(i,j)=sqrt(denkx(i)*denky(j))*amp(i,j)
         enddo
      enddo

c fix up the spectrum to satisfy reality contraints (we could have used
c a complex-to-real fft for the ky-->y fft below, and saved a factor of 2 
c in speed and space, but for simplicity we just used complex-to-complex
c fft's):
      do i=2,mx
         do j=2,my/2
            denk(mx+2-i,my+2-j)=conjg(denk(i,j))
         enddo
      enddo

c handle the axes carefully (since we don't have the periodic points on
c the grid):
      do i=2,mx/2
         denk(mx+2-i,1)=conjg(denk(i,1))
      enddo
      do j=2,my/2
         denk(1,my+2-j)=conjg(denk(1,j))
      enddo

c zero out the phase-less modes (which are at very high k anyway,
c and we presume the spectrum goes to zero up there):

      denk(1,1)=0.0  ! a constant off-set doesn't show up in plots
      do i=1,mx
         denk(i,my/2+1)=0.0
      enddo
      do j=1,my
         denk(mx/2+1,j)=0.0
      enddo

c do inverse FFT:
      call cfft99(denk,work,trigsx,ifaxx,1,mx,mx,my,1)
      call cfft99(denk,work,trigsy,ifaxy,mx,1,my,mx,1)

c it should already be real (or very close to it):
      rmaxr=0.0
      rmaxi=0.0
      do i=1,mx
         do j=1,my
            rmaxr=max(rmaxr,abs(real(denk(i,j))))
            rmaxi=max(rmaxi,abs(aimag(denk(i,j))))
         enddo
      enddo

      write(6,*) 'max real and imaginary parts =', rmaxr, rmaxi

      dmin=0.0
      dmax=0.0
      davg=0.0
      do i=1,mx
         do j=1,my
            den(i,j)=real(denk(i,j))
            dmin=min(dmin,den(i,j))
            dmax=max(dmax,den(i,j))
            davg=davg+den(i,j)
         enddo
      enddo
      davg=davg/mx/my
      write(6,*) 'dmin, dmax, davg=', dmin,dmax,davg

      d2=0.0
      d4=0.0
      do i=1,mx
         do j=1,my
            d2=d2+(den(i,j)-davg)**2
            d4=d4+(den(i,j)-davg)**4
         enddo
      enddo
      d2=d2/mx/my
      kurtosis=d4/mx/my/d2**2
      write(6,*) 'sigma, kurtosis+3 =', sqrt(d2), kurtosis

c  calculate 2-point correlation in x-direction:
      idelta=0
      do while ((x(1+idelta)-x(1)) .lt. 10.0)

c skip last bit to ignore periodic wrap-around:
      d2=0.0
      d4=0.0
      davg=0.0
      do i=1,mx-idelta
         do j=1,my
            davg=davg+den(i,j)
            d2=d2+den(i,j)**2
            d4=d4+den(i,j)*den(i+idelta,j)
         enddo
      enddo
      npts=my*(mx-idelta)
      corrx=(d4/npts - (davg/npts)**2) / (d2/npts - (davg/npts)**2)
      write(6,*) ' idelta,x, Cor(x), Cor_fft(x)=', 
     >     idelta, x(1+idelta)-x(1), corrx, rcorxfft(1+idelta)

      idelta=idelta+1
      enddo

c  calculate 2-point correlation in y-direction:
      jdelta=0
      do while ((y(1+jdelta)-y(1)) .lt. 10.0)

c skip last bit to ignore periodic wrap-around:
      d2=0.0
      d4=0.0
      davg=0.0
      do i=1,mx
         do j=1,my-jdelta
            davg=davg+den(i,j)
            d2=d2+den(i,j)**2
            d4=d4+den(i,j)*den(i,j+jdelta)
         enddo
      enddo
      npts=(my-jdelta)*mx
      corry=(d4/npts - (davg/npts)**2) / (d2/npts - (davg/npts)**2)
      write(6,*) ' jdelta,y, Cor(y), Cor_fft(y)=', 
     >     jdelta, y(1+jdelta)-y(1), corry, rcoryfft(1+jdelta)

      jdelta=jdelta+1
      enddo

c generate some test data:
c      do i=1,mx
c         do j=1,my
c            den(i,j)=cos(2*pi*i/x0)*cos(4*pi*j/y0)
c         enddo
c      enddo

c reverse the indexing of the array for stupid HDF:
c also move origin to put zero phase mode in the 
c center of the box:
      do i=1,mx
         do j=1,my
            iget=mod(i+mx/2-1,mx)+1
            jget=mod(j+my/2-1,my)+1
            areverse(j,i)=den(iget,jget)
         enddo
      enddo

      shape(2)=mx
      shape(1)=my

c write out density(x,y) to the an HDF file:
      iret = dssdast('BES simulation','density','e11.3','cartesian')
      iret = dsadata('bes.hdf',2,shape,areverse)

c write out density(x,y) to a simple file:
      open(unit=21,file='bes.dat',status='unknown')
      write(21,*) mx, ' = # of x grid points '
      write(21,*) my, ' = # of y grid points '
      write(21,*) 'x grid:'
      write(21,*) (x(i),i=1,mx)
      write(21,*) 'y grid:'
      write(21,*) (y(i),i=1,my)
      write(21,*) 'Phi(x,y):'
      write(21,*) ((den(i,j),i=1,mx),j=1,my)
      close(unit=21)

c Solve a Langevin-type equation to generate a new set of complex amplitude
c coefficients whose correlations decay with a certain correlation time.
c
c The equivalent Fokker-Planck equation can be solved exactly for a
c finite size time step, and the result then put in terms of a
c drag term (exp(-dt)) and a diffusive random step size
c proportional to (1.0-exp(-2*dt))**0.5.
c
c For small dt, this reduces to the usual Langevin equation.
c

      if(nframe .ne. mframe) then

         write(6,*) 
     > 'Repeating pattern, unless using random amplitude algorithm'

         if(amp_alg .eq. 'r') then

            damp=exp(-dt)
            damp2=(1.0-exp(-2*dt))**0.5
            do i=1,mx
               do j=1,my
                ampr=g05ddf(0.d0,1.d0) ! NAG's Gaussian random number generator
                ampi=g05ddf(0.d0,1.d0)
                amp(i,j)=amp(i,j)*damp+cmplx(ampr,ampi)*damp2
                enddo
             enddo

          endif

      endif

      enddo     ! end NFRAME loop

      stop
      end