program search
c*
c* Copyright (C) 1999 by Peter Mueller, MPIfR.  peter@mpifr-bonn.mpg.de
c*
c* Permission to use, copy, modify, and distribute this software and its
c* documentation for any purpose and without fee is hereby granted, provided
c* that the above copyright notice appear in all copies and that both that
c* copyright notice and this permission notice appear in supporting
c* documentation.  This software is provided "as is" without express or
c* implied warranty.

c.....HP 9000/700 
c.....f77 +Oall +e search.f -o search
c
c.....SUN
c.....f77 -fast -O search.f -o search
c
c.....LINUX
c.....g77 -O search.f -o search
c
      implicit real*4 (a-h,o-z)

      parameter (NDIM = 2**20)

      real*4 pmax, xmax
      real*4 signal(10000)

      common /array/ time(NDIM)

      nmax = 2**18
      iseed = secnds(0.0)
      dum = ran(iseed)
      p = 239.0 + ran(iseed)
      rms = 3.5 
      ip = int(p)

      call testpuls(nmax,p,rms) 
      call ffa(signal,nmax,ip,pmax,xmax)

      open(unit=13,file='puls',status='unknown')
      write(13,'(''# Periode :'',g20.10)') pmax
      pm = -999999.0
      do lm = 1, ip
         if (signal(lm) .gt. pm) then
            pm = signal(lm)
            lmax = lm
         endif
      enddo
      lmm = ip/2 + lmax
      do lm = 1, ip
         km = mod(lm + lmm, ip)
         if (km .eq. 0) km = ip
         write(13,*) lm, signal(km)
      enddo
      

      end
c
c************************************************************************
c
      subroutine ffa(signal,mmax,ip,pmax,xmax)
c*
c* Copyright (C) 1995 by Peter Mueller, MPIfR.   peter@mpifr-bonn.mpg.de
c*
c* Permission to use, copy, modify, and distribute this software and its
c* documentation for any purpose and without fee is hereby granted, provided
c* that the above copyright notice appear in all copies and that both that
c* copyright notice and this permission notice appear in supporting
c* documentation.  This software is provided "as is" without express or
c* implied warranty.

      implicit real*4 (a-h,o-z)

      parameter (NDIM = 2**20)
      
      real puls(NDIM)
      real signal(*)
      real*8 p, p1, p2, p3, p4, p5, ps, dp 
      real*4 x, xp, xmax, xmax0, xmax1, pmax

      common /array/ time(NDIM)

      xmax = -99999999.0
      nmax = ip * 2**int(log(float(mmax)/float(ip)) / log(2.0) + 0.95)
      if(nmax .gt. NDIM) nmax = nmax / 2
      lmax = nmax / (ip +1)
      kp = nint(log(float(lmax+1))/log(2.0))
      dp = dfloat(ip + 1) / dfloat(nmax)
      p1 = 1.d0 / dfloat(ip)
      p2 = p1 * p1
      p3 = p2 * p1
      p4 = p3 * p1
      p5 = p4 * p1
      ps = p1 + p2 + p3 + p4 + p5

      do n = 0, lmax-1
         x = dp * n
         xp = (((((p5*x + p4)*x + p3)*x + p2)* x - ps)*x + p1)*x
         p = ip + x - xp
         do l = 0, kp-1
            if (mod(n, 2**l) .eq. 0) kp1 = kp - l 
         enddo
         do k = kp1, kp
            np = 2**(kp-k)
            joff = nmax - nmax / 2**(k-1)
            ioff = nmax - nmax / 2**max(k-2,0)
            ish = mod((n + np) / np / 2, ip)
            do i = 0, np-1
               iip = i * ip
               i0 = iip + joff
               i1 = 2 * iip + ioff
               i2 = i1 + ip
               do j = 1, ip
                  j1 = j + ish
                  if(j1 .gt. ip) j1 = j1 - ip
                  ind = j + i0
                  jnd = j + i1
                  knd = j1 + i2
                  if(kp1.eq.1 .and. k.eq.1) then
                     puls(ind) = time(jnd) + time(knd)
                  else
                     puls(ind) = puls(jnd) + puls(knd)
                  endif
               enddo
            enddo
         enddo

         xmax1 = xmax
         xmax0 = -9999999.0
         do lm = 1, ip
            xmax0 = max(puls(lm+joff),xmax0)
         enddo
         xmax = max(xmax0,xmax)

c        print*, p, xmax0 / float(lmax)

         if(xmax .gt. xmax1) then
            pmax = p
            do lm = 1, ip
               signal(lm) = puls(lm+joff) / float(lmax)
            enddo
         endif
      enddo
         
      return
      end
c
c************************************************************************
c
      subroutine testpuls(nmax,p,rms)
c*
c* Copyright (C) 1995 by Peter Mueller, MPIfR.   peter@mpifr-bonn.mpg.de
c*
c* Permission to use, copy, modify, and distribute this software and its
c* documentation for any purpose and without fee is hereby granted, provided
c* that the above copyright notice appear in all copies and that both that
c* copyright notice and this permission notice appear in supporting
c* documentation.  This software is provided "as is" without express or
c* implied warranty.

      implicit real*4 (a-h, o-z)

      parameter (NDIM = 2**20)

      real*8 dx

      common /array/ time(NDIM)

      iseed = 42

c     open(unit=21,file='puls.dat',status='unknown')
      dx = 2.0 * acos(0.0) / p
      do n = 1, nmax
         rms1 = rms*(ran(iseed)-0.5)
         time(n) = sin(dx*n)**150 + rms1
c        if(n.le.3*p) write(21,*) n, time(n)
      enddo
c     close(21)

      return
      end
c
c==============================================================================
      function ran(idum)
c==============================================================================
c
c     Donald Knuth's portable random number generator, nabbed from p199
c     of numerical recipes. Gives same random numbers on both Sun's &
c     Hp's given the same initial seed. Pass idum down as negative to 
c     reshuffle the random numbers.
c
      parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
      dimension ma(55)
      data iff /0/
      save
      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
      ran=mj*fac

      end