      program wavspe
      parameter(npk=32768,np2=npk/2)
      dimension b(npk),c(2,npk),p(np2)
      ibm=4711
      do i=1,1000
      ibm=ibm*16807
      enddo
      rinv=0.5/(2.**31-1.)
c Genererer virrevandringen
      k=1
      amp=1.
1600  continue
      if(2**k.eq.npk) goto 1700
      amp=amp*0.5
      do i=2**k+1,2**(k+1)
      ibm=ibm*16807
      ran=ibm*rinv
      b(i)=ibm*rinv*amp
      enddo
      k=k+1
      goto 1600
1700  continue
      call wt1(b,npk,-1)
      open(unit=1,file='wavsp1.dat',status='unknown')
      do i=1,npk
      write(1,*) i,b(i)
      c(1,i)=b(i)
      c(2,i)=0.
      enddo
      close(1)
c Analyserer virrevandringen
      call four1(c,npk,1)
      do i=1,np2
      p(i)=0.5*(abs(c(1,i)*c(2,i))+abs(c(1,npk+1-i)*c(2,npk+1-i)))
      enddo
      open(unit=2,file='wavsp2.dat',status='unknown')
      do i=1,np2
      write(2,*) i,p(i)
      enddo
      close(2)
      end
      SUBROUTINE wt1(a,n,isign)
      INTEGER isign,n
      REAL a(n)
      EXTERNAL wtstep
CU    USES wtstep
      INTEGER nn
      if (n.lt.4) return
      if (isign.ge.0) then
        nn=n
1       if (nn.ge.4) then
          call daub4(a,nn,isign)
          nn=nn/2
        goto 1
        endif
      else
        nn=4
2       if (nn.le.n) then
          call daub4(a,nn,isign)
          nn=nn*2
        goto 2
        endif
      endif
      return
      END
      SUBROUTINE daub4(a,n,isign)
      INTEGER n,isign,NMAX
      REAL a(n),C3,C2,C1,C0
      PARAMETER (C0=0.4829629131445341,C1=0.8365163037378079,
     *C2=0.2241438680420134,C3=-0.1294095225512604,NMAX=32768)
      REAL wksp(NMAX)
      INTEGER nh,nh1,i,j
      if(n.lt.4)return
      if(n.gt.NMAX) pause 'wksp too small in daub4'
      nh=n/2
      nh1=nh+1
      if (isign.ge.0) then
        i=1
        do 11 j=1,n-3,2
          wksp(i)=C0*a(j)+C1*a(j+1)+C2*a(j+2)+C3*a(j+3)
          wksp(i+nh)=C3*a(j)-C2*a(j+1)+C1*a(j+2)-C0*a(j+3)
          i=i+1
11      continue
        wksp(i)=C0*a(n-1)+C1*a(n)+C2*a(1)+C3*a(2)
        wksp(i+nh)=C3*a(n-1)-C2*a(n)+C1*a(1)-C0*a(2)
      else
        wksp(1)=C2*a(nh)+C1*a(n)+C0*a(1)+C3*a(nh1)
        wksp(2)=C3*a(nh)-C0*a(n)+C1*a(1)-C2*a(nh1)
        j=3
        do 12 i=1,nh-1
          wksp(j)=C2*a(i)+C1*a(i+nh)+C0*a(i+1)+C3*a(i+nh1)
          wksp(j+1)=C3*a(i)-C0*a(i+nh)+C1*a(i+1)-C2*a(i+nh1)
          j=j+2
12      continue
      endif
      do 13 i=1,n
        a(i)=wksp(i)
13    continue
      return
      END
      SUBROUTINE four1(data,nn,isign)
      INTEGER isign,nn
      REAL data(2*nn)
      INTEGER i,istep,j,m,mmax,n
      REAL tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      n=2*nn
      j=1
      do 11 i=1,n,2
        if(j.gt.i)then
          tempr=data(j)
          tempi=data(j+1)
          data(j)=data(i)
          data(j+1)=data(i+1)
          data(i)=tempr
          data(i+1)=tempi
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
          j=j-m
          m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2
2     if (n.gt.mmax) then
        istep=2*mmax
        theta=6.28318530717959d0/(isign*mmax)
        wpr=-2.d0*sin(0.5d0*theta)**2
        wpi=sin(theta)
        wr=1.d0
        wi=0.d0
        do 13 m=1,mmax,2
          do 12 i=m,n,istep
            j=i+mmax
            tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
            tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
            data(j)=data(i)-tempr
            data(j+1)=data(i+1)-tempi
            data(i)=data(i)+tempr
            data(i+1)=data(i+1)+tempi
12        continue
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
      goto 2
      endif
      return
      end
