parameter(nnn=60000)
      dimension ddy(nnn),acc(nnn),vel(nnn),dis(nnn),vels(nnn)
      real ddy1(nnn),ddy2(nnn),acc1(nnn),acc2(nnn),vel1(nnn),vel2(nnn)
      real ddy3(nnn),acc3(nnn),vel3(nnn)
      real velm(200),velms(200),iv1,iv2,iv3,iv
      character fmt*8,dm*1,eqn*6
      data h/0.05/
      data a1/1.936/b1/2.011/a2/2.030/b2/1.251/a3/2.171/b3/1.002/
      data bs/5.0/us/5.5/bl/5.5/ul/6.0/
c     "\mhp\eqt.dat"
c     "\mhp\imvelcns.eq"
c     "\mhp\imvelcew.eq"
c     "\mhp\imvelcud.eq"
c     "\mhp\tsi.dat"
      open(1,file='eqp.dat')!■計算する地震動の名前を書いたファイル
      open(7,file='tsi.dat')!■提案震度出力ファイル
      write(7,*)'eqn   0.1-1 0.5-1   1-2   Ip'
 1    read(1,'(a6)')eqn
      if(eqn.eq.'000')goto2
      write(6,*)eqn
      open(2,file='\era\eq\'//eqn//'ns.eq')!■地震動データ（ns成分）
      open(3,file='\era\eq\'//eqn//'ew.eq')!■地震動データ（ewk成分）
      open(4,file='\era\eq\'//eqn//'ud.eq')!■地震動データ（ud成分）
      read(2,'(a1)')dm
      read(3,'(a1)')dm
      read(4,'(a1)')dm
      read(2,'(4f8.0)')dt,stt,edt,fct
      read(3,'(4f8.0)')dt,stt,edt,fct
      read(4,'(4f8.0)')dt,stt,edt,fct
c     dt: 時間刻み, stt: スタート時刻, edt: 終了時刻, ■
c     fct: 倍率の1/100（0.01だとそのまま、ということ）■
      read(2,'(a8)')fmt! ■地震動データフォーマット(ns)
      read(3,'(a8)')fmt! ■地震動データフォーマット(ew)
      read(4,'(a8)')fmt! ■地震動データフォーマット(ud)
      read(2,'(i4)')inum!■地震動データの１行のデータ数(ns)
      read(3,'(i4)')inum!■地震動データの１行のデータ数(ew)
      read(4,'(i4)')inum!■地震動データの１行のデータ数(ud)
      nn=(edt-stt)/dt
      jj=nn/inum+1
      do j=1,jj
      read(2,fmt)(ddy1(i),i=j*inum-(inum-1),j*inum)
      read(3,fmt)(ddy2(i),i=j*inum-(inum-1),j*inum)
      read(4,fmt)(ddy3(i),i=j*inum-(inum-1),j*inum)
        do i=j*inum-(inum-1),j*inum
        ddy1(i)=ddy1(i)*fct*100.
        ddy2(i)=ddy2(i)*fct*100.
        ddy3(i)=ddy3(i)*fct*100.
        end do
      end do
      do k=1,200
      t=real(k)/100.
      w=6.283185/t
      call resp(h,w,dt,nn,ddy1,acc1,vel1,dis,nnn,sa1,sv1,sd1)
      call resp(h,w,dt,nn,ddy2,acc2,vel2,dis,nnn,sa2,sv2,sd2)
      call resp(h,w,dt,nn,ddy3,acc3,vel3,dis,nnn,sa3,sv3,sd3)!3
      velm(k)=0.
      velms(k)=0.!3
      do m=2,nn
      acc(m)=sqrt(acc1(m)**2+acc2(m)**2)
      vel(m)=sqrt(vel1(m)**2+vel2(m)**2)
      vels(m)=sqrt(vel1(m)**2+vel2(m)**2+vel3(m)**2)!3
      if(vel(m).gt.velm(k))velm(k)=vel(m)
      if(vels(m).gt.velms(k))velms(k)=vels(m)!3
      end do
      end do!k
c
c     ■短周期（0.1-1秒）
      v1=0.
      do k=10,100
      v1=v1+velms(k)!3
      end do
      v1=v1/91.
      iv1=a1*log(v1)/log(10.)+b1
c     ■中周期（0.5-1秒）
      v2=0.
      do k=50,100
      v2=v2+velm(k)
      end do
      v2=v2/51.
      iv2=a2*log(v2)/log(10.)+b2
c     ■長周期（1-2秒）
      v3=0.
      do k=100,200
      v3=v3+velm(k)
      end do
      v3=v3/101.
      iv3=a3*log(v3)/log(10.)+b3
c
      if(iv3.ge.ul)then
       iv=iv3
       ic=1
       goto 3
      end if
      if(iv3.ge.bl.and.iv3.lt.ul)then
       iv=(iv3*(iv3-bl)+iv2*(ul-iv3))*2!■比例配分の領域の幅が0.5なので*2
       ic=11
       goto 3
      end if
c
      if(iv1.le.bs)then
       iv=iv1
       ic=2
       goto 3
      end if
      if(iv1.gt.bs.and.iv1.le.us)then
       iv=(iv1*(us-iv1)+iv2*(iv1-bs))*2!■比例配分の領域の幅が0.5なので*2
       ic=21
       goto 3
      end if
      iv=iv2
      ic=3
c
 3    write(7,'(a6,4f6.2,i3)')eqn,iv1,iv2,iv3,iv,ic
      goto1
 2    continue
      stop
      end
c
      subroutine resp(h,w,dt,nn,ddy,acc,vel,dis,nd,sa,sv,sd)
      dimension ddy(nd),acc(nd),vel(nd),dis(nd)
      w2=w*w
      hw=h*w
      wd=w*sqrt(1.-h*h)
      wdt=wd*dt
      e=exp(-hw*dt)
      cwdt=cos(wdt)
      swdt=sin(wdt)
      a11=e*(cwdt+hw*swdt/wd)
      a12=e*swdt/wd
      a21=-e*w2*swdt/wd
      a22=e*(cwdt-hw*swdt/wd)
      ss=-hw*swdt-wd*cwdt
      cc=-hw*cwdt+wd*swdt
      s1=(e*ss+wd)/w2
      c1=(e*cc+hw)/w2
      s2=(e*dt*ss+hw*s1+wd*c1)/w2
      c2=(e*dt*cc+hw*c1-wd*s1)/w2
      s3=dt*s1-s2
      c3=dt*c1-c2
      b11=-s2/wdt
      b12=-s3/wdt
      b21=(hw*s2-wd*c2)/wdt
      b22=(hw*s3-wd*c3)/wdt
      acc(1)=2.*h*w*ddy(1)*dt
      vel(1)=-ddy(1)*dt
      dis(1)=0.
      dx=vel(1)
      x=0.
      sa=0.
      sv=0.
      sd=0.
      do 110 m=2,nn
      dxf=dx
      xf=x
c     write(6,*)m,ddy(m)
      ddym=ddy(m)
      ddyf=ddy(m-1)
      x=a12*dxf+a11*xf+b12*ddym+b11*ddyf
      dx=a22*dxf+a21*xf+b22*ddym+b21*ddyf
      ddx=-2.*hw*dx-w2*x
      acc(m)=ddx
      vel(m)=dx
      dis(m)=x
      sa=amax1(sa,abs(ddx))
      sv=amax1(sv,abs(dx))
      sd=amax1(sd,abs(x))
 110  continue
      return
      end