c
c           power spectra cal. and plot
c                    by fft (triangle filter)
c
c           2006.10.10　変更(oti) t0の設定(subroutine power)
c           

c
      character fname*50,infile*33,outfile*33,dum*21
      dimension data(70000),idata(5)
      dimension data0(20,70000),denn(70000)
      dimension p(70000),frq(70000)
      dimension pp(20,70000)
      common /trns/ h(70000),df1,np1
c--
c
      open(10,file='powers.05')
      
      read(10,*)
      read(10,*) infile
      read(10,*)
      read(10,*) outfile
c
c      read(10,*)
c      read(10,*) Nmax
c
      read(10,*)
      read(10,*) icase
c
      read(10,*)
      read(10,*) ipass
      read(10,*)
      read(10,*) isw,jsw,dep,dt,dr,nd
      write(6,*) isw,jsw,dep,dt,dr,nd
c      read(10,*) 
c      read(10,'(a33)') infile
c      read(10,'(a6)') outfile
  
c      nd=2400                                                           ! データ個数に変更する
c
      nmax=1
      do n=1,Nmax
      open(20,file=infile,status='old')
      open(30,file=outfile)
      
c      read(20,*) dum,dt
      do i=1,ipass
       read(20,*) 
      end do
       
      do i=1,nd
       read(20,*,end=100) (data0(j,i),j=1,icase)
      end do
  100 continue 
      nd=i-1
      
      do iiii=1,icase

       do i=1,nd
        data(i)=data0(iiii,i)
       end do

       call sqd100(nd,data,0.0,dr,nc)
      
       write(6,*) '  dt   : ',dt
       write(6,*) '  nd   : ',nd

       call trend(data,nd)

       write(6,*) '%%% power start %%%'
       call power(nd,dt,data,p,df,np)
       write(6,*) '%%% power end   %%%'
c
       do i=1,np
       pp(iiii,i)=p(i)
       end do

      end do
      
      write(30,'(a13,20i13)') '周波数(Hz)',(j,j=1,icase)
      do i=1,np
      frq(i)=df*(i)
      write(30,601) frq(i),(pp(j,i),j=1,icase)
      end do
  600 format(5i4,i6,f6.2,2i6)
  601 format(20f13.5)
c
c  100 continue
c
      close(20)
      close(30)
c

      end do   ! Nmax loop

      close(10)
c
      stop
      end
c
c-----------------------------------------------------------
      subroutine ptrans(data,nd,dt,dep1)
c
c
c        pressure to water surface elevation transform
c                           by fft method
c                                         88/11/25
c  
c
      dimension     data(nd),data1(70000)
c
      if(dep1.le.0.) dep1 = 0.645
c
      ave = 0.0
c
      do 20 i = 1 , nd
      ave = ave + data(i)
      write(99,*)
   20 continue
c
      ave = ave / nd
      dep = ave / 1.03 + dep1
      write(6,*) 'dep=',dep
c
      do 30 i = 1 , nd
      data(i) = data(i) - ave
   30 continue
c
      data(nd+1) =  ( data(1) - data(nd) ) * 1. / 3. + data(nd)
      data(nd+2) =  ( data(1) - data(nd) ) * 2. / 3. + data(nd)
c
      nd2    = nd + 2
      call trans(nd2,dt,data,data1,dep,dep1)
c
      do 40 i = 1 , nd
      data(i) = data1(i) + dep
   40 continue
c
      return
      end
c
      subroutine trans(nmax,dt,data,data1,dep,dep1)
c
      complex      x   (70000)
      dimension    data(nmax),data1(nmax)
c
      n = 1
    1 if(n.ge.nmax) go to 2
      n = 2 * n
      go to 1
    2 continue

c      n=n/2
c
      do 10 i=1,n
      if(i.le.nmax) x(i) = cmplx(data(i),0.0)
      if(i.gt.nmax) x(i) = ( 0.0 , 0.0 )
   10 continue

c
c
      call fastft( n,x,-1 )
      call filt(dep,dep1,dt,n,x)
      call fastft( n,x, 1 )
c
      do 20 i=1,nmax
      data1(i) = real(x(i))
   20 continue

c
c
      return
      end
c
      subroutine filt(dep,dep1,dt,n,x)
      complex   x(70000)
      common /trns/ h(70000),df,n2
c
      n2 = n / 2 - 1
      df = 1. / ( dt * n )
c
      do 10 i = 2 , n2
      j = n - i + 2
      f = df * ( i - 1 )
      t = 1. / f
      call wave(9.8,t,dep,wk)
      wkdep =wk*dep
      wkdep1=wk*dep1
      if(wkdep.gt. 89.0) wkdep =89.0
      if(wkdep1.gt.89.0) wkdep1=89.0
      h(i) = cosh(wkdep) / cosh(wkdep1) / 1.03
      write(90,*) h(i)
      if(h(i).gt.24.) h(i) = 24.
      x(i) = h(i) * x(i)
      x(j) = h(i) * x(j)
   10 continue
c
      return
      end



      subroutine power(nmax,dt,data,p,df,np)
      complex      x   (70000)
      dimension    data(1),p(1),pp(70000)
      data ii,jj/4,16/
c%%      data ii,jj/2,8/
c
      nt = 1
      do 10 i = 1 , nmax
      x(i) = cmplx(data(i),0.0)
   10 continue
    1 nt = nt * 2
      if(nt-nmax) 1 , 30 , 2
    2 continue
      do 20 i = nmax + 1 , nt
      x(i) = cmplx(0.0,0.0)
   20 continue
   30 continue
      nh = nt / 2 + 1
      t  = dt * nt
      df = ii / t
c
      call fastft( nt,x,-1 )
c

      t0= nmax * dt

      pp(1) = abs(x(1)) **2 * t*(t/t0)
      do 40 i=2,nh-1
      pp(i) = 2.0 * abs(x(i)) **2 * t*(t/t0)
   40 continue
      pp(nh)= abs(x(nh)) **2 * t*(t/t0)
c
      call trfilt(pp,p,nh,np,ii,jj)
c
      return
      end
      subroutine fastft(n,x,ind)
c
      complex x(1),temp,theta
      j=1
      do 140 i = 1, n
      if(i.ge.j) go to 110
      temp=x(j)
      x(j)=x(i)
      x(i)=temp
  110 m=n/2
  120 if(j.le.m) go to 130
      j=j-m
      m=m/2
      if(m.ge.2) go to 120
  130 j=j+m
  140 continue
      kmax=1
  150 if(kmax.ge.n) go to 180
      istep=kmax*2
      do 170 k = 1 , kmax
      theta=cmplx(0.0,3.141593*float(ind*(k-1))/float(kmax))
      do 160 i = k , n , istep
      j=i+kmax
      temp=x(j)*cexp(theta)
      x(j)=x(i)-temp
      x(i)=x(i)+temp
  160 continue
  170 continue
      kmax=istep
      go to 150
  180 continue
      if(ind.ge.0) return
      do 190 i = 1 , n
      x(i) = x(i) / n
  190 continue
      return
      end
      subroutine trfilt(p,pp,mm,nn,ii,jj)
      dimension p(1),pp(1)
      nn = 0
      do 20 i = 1 , mm , ii
      nn = nn + 1
      s = p(i) * jj
      if(jj.eq.1) go to 20
      do 10 j = 1 , jj-1
      l = i - j
      if(l.le.0) l = 2 - l
      m = i + j
      if(m.ge.mm) m = 2 * mm - m
   10 s = s + float(jj-j) * ( p(l) + p(m) )
   20 pp(nn) = s / float(jj*jj)
      return
      end
      subroutine trend(y,nmax)
      dimension y(nmax)
c
      fn    = nmax
      sx    = fn * ( fn - 1.) / 2.
      sxx   = sx * ( 2. * fn - 1.) / 3.
      sxxx  = sx * sx
      sxxxx = sxx * ( 3. * fn **2 - 3. * fn - 1.) / 5.
c
      sy    = 0.
      sxy   = 0.
      sxxy  = 0.
c
      do 10 i = 1 , nmax
      x     = i - 1
      sy    = sy   + y(i)
      sxy   = sxy  + y(i) * x
      sxxy  = sxxy + y(i) * x * x
   10 continue
c
      aa    = sxx   - sx  ** 2 / fn
      ab    = sxxy  - sxx * sy / fn
      ac    = sxxx  - sxx * sx / fn
      ad    = sxy   - sy  * sx / fn
      ae    = sxxxx - sxx ** 2 / fn
      a     = ( aa * ab - ac * ad ) / ( aa * ae - ac * ac )
      b     = ( ab - a * ae ) / ac
      c     = ( sy - a * sxx - b * sx ) / fn
c
      do 20 i = 1 , nmax
      t     = i - 1
      tt    = a * t * t + b * t + c
      y(i)  = y(i) - tt
   20 continue
c
      return
      end
      subroutine wave(g,t,d,wk)
c
      data pi2/6.283185/
c
      wl0 = g * t * t / pi2
      dl0 = d / wl0
      x = pi2 * dl0
       if(dl0.lt.0.2) then
        y = sqrt(x) * ( 1. + x / 6. * ( 1. + x / 5. ) )
       else
        yp = x * ( 1. + 1.26 * exp(-1.84*x) )
        st = exp(-2.*yp)
        y  = x * ( 1. + 2. * st * ( 1. + st ) )
       endif
      wk = y / d
c
      return
      end


c
 
      
c********************************************************************
       subroutine sqd100(n,x,xi,dr,ncn)
c********************************************************************
c
       dimension x(1)
c
c --- noise cut
c
c     author  : n.kishida
c
c     input number of data  : n
c           data x          : x
c           initial setting : xi
c           delta x         : dr
c     output data x         : x
c     dimension x(1)
      ncn  = 0
      a    = xi
      if (xi.eq.0.) a = x(1)
      do 200 i = 1 , n
      b0   = a-x(i)
      b1   = abs(b0)
      if (b1.gt.dr) go to 100
      go to 101
  100 continue
      ncn  = ncn + 1
      x(i) = a
  101 continue
      a    = x(i)
  200 continue
      return
      end

      SUBROUTINE WAVE2 (G,T,DEPTH,WL,RN,RKS)
C
C------------- WAVE --------------------------------------------------
C
      PI  = 3.1415926
      PI2 = PI * 2.0
C
      CON = PI2**2 * DEPTH  / (T**2 * G)
      X1  = CON
c      IF ( X1.GT.87.0 )            GO TO 50
      IF ( X1.GT.10.0 )            GO TO 50
      NKO = 1
   30 TH  = TANH(X1)
      F   = X1 * TH - CON
      F1  = TH + X1 * (1.0-TH**2)
C
      X2  = X1 - F / F1
      IF (ABS(X2-X1).LT.0.0001)    GO TO 10
      IF (NKO.GT.20)              GO TO 20
      NKO = NKO + 1
      X1  = X2
      GO TO 30
   20 X2  = (X1+X2) / 2.0
C
   10 WL  = PI2 * DEPTH / X2
      E   = EXP(X2*2.0)
      SIH = (E-1.0/E) / 2.0
      RN  = (1.0 + X2 * 2.0 / SIH)
      RKS = TANH(X2) * RN
      RKS = 1.0 / SQRT(RKS)
      RN  = RN / 2.0
      RETURN
C   - DEPTH WATER -
   50 CONTINUE
      WL  = G * T**2 / PI2
      RN  = 0.5
      RKS = 1.0
      RETURN
      END


