c
c  *** 	Calcuration of Seismic Response Wave and Energy ***
c
c  Parameters:
c   mm: Number of Delta Time for Main Dimension (> nn) 
c   nn: Actual Number of Delta Time
c   dt: Delta Time (dt*nn = duration)
c    h: Dampling Factor
c   sa: Maximum Absolute Response Acceleration (gal)
c   sv: Maximum Relative Response Velocity     (cm/sec)
c  asv: Maximum Absolute Response Velocity     (cm/sec)
c   sd: Maximum Relative Response Displacement (cm)
c
c  Dimension Data:
c   ddy: Ground Acceleration (gal)
c    dy: Ground Velocity     (cm/s)
c   acc: Absolute Response Acceleration (gal)
c   vel: Relative Response Velocity     (cm/sec)
c  avel: Absolute Response Velocity     (cm/sec)
c   dis: Relative Response Displacement (cm)
c
      program main
      parameter (mm=131100)
      IMPLICIT REAL*8(A-H,O-Z)
      dimension ddy1(mm),ddy2(mm),ddy3(mm),time(mm)
      dimension dy1(mm),dy2(mm),dy3(mm)
      dimension acc1(mm),vel1(mm),dis1(mm)
      dimension acc2(mm),vel2(mm),dis2(mm)
      dimension acc3(mm),vel3(mm),dis3(mm)
      CHARACTER INFILE*40,otfile*40,head1*40,head2*40,head3*40,head4*40
c
      write(*,*)
      write(*,*)'Enter File Name for Ground Motion Acceleraration:'
      read(*,1000) infile
 1000 format(a40)
c
      OPEN(11,FILE=INFILE)
        read(11,*)head1,head2,head3,head4
        i=1
   10   read(11,*,end=11)time(i),ddy1(i),ddy2(i),ddy3(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      nn=i-1
      dt=time(2)-time(1)
      write(*,*)'nn and DT=',nn,dt
      IF(nn.gt.mm)then
        WRITE(*,*)' nn > mm !!!!!'
        write(*,*)'nn and mm=',nn,mm
        read(*,*)dummy
        STOP
      end if
c
c *** Ground Velocity ***
c
      dt2=dt/2.d0
      write(*,*)
      write(*,*)'Ground velocities are seved in "Ground-Vel.csv"'
      write(*,*)'Check the ground velocity waveforms!'
c
      open(11,file='Ground-Vel.csv')
      dy1(1)=0.d0
      dy2(1)=0.d0
      dy3(1)=0.d0
      write(11,*)'time(s),vel1,vel2,vel3'
      write(11,1500)0.d0,dy1(1),dy2(1),dy3(1)
      do i=2,nn
        dy1(i)=dy1(i-1)+(ddy1(i-1)+ddy1(i))*dt2
        dy2(i)=dy2(i-1)+(ddy2(i-1)+ddy2(i))*dt2
        dy3(i)=dy3(i-1)+(ddy3(i-1)+ddy3(i))*dt2
        write(11,1500)time(i),dy1(i),dy2(i),dy3(i)
      end do
      close(11)
 1500 format(f8.4,3(',',f12.4))
c
c *** Calculation of Response Waves ***
c
      write(*,*)
      write(*,*)'Enter Period for the Respose Building(s)'
      read(*,*)Per
      om=6.28318530718d0/Per
      write(*,*)'Enter Damping Factor for the Respose Building'
      read(*,*)h
      omh=dsqrt(1.d0-h*h)*om
      PerH=6.28318530718d0/omh
      write(*,*)
      write(*,*)'T and damping T =',Per,PerH
c
      open(15,file='Energy-Density-wave1.csv')
      write(15,*)'time(s),Ek(Kinematic),Ep(Potential),Ed(damping),Et(Tot
     *al)'
      call resp(h,om,dt,nn,ddy1,dy1,acc1,vel1,dis1,mm,sa1,sv1,asv1,sd1)
      close(15)
c
      open(15,file='Energy-Density-wave2.csv')
      write(15,*)'time(s),Ek(Kinematic),Ep(Potential),Ed(damping),Et(Tot
     *al)'
      call resp(h,om,dt,nn,ddy2,dy2,acc2,vel2,dis2,mm,sa2,sv2,asv2,sd2)
      close(15)
c
      open(15,file='Energy-Density-wave3.csv')
      write(15,*)'time(s),Ek(Kinematic),Ep(Potential),Ed(damping),Et(Tot
     *al)'
      call resp(h,om,dt,nn,ddy3,dy3,acc3,vel3,dis3,mm,sa3,sv3,asv3,sd3)
      close(15)
c
c *** Output of Response Waves ***
c
      write(*,*)
      write(*,*)'Enter Output File Name for Response Wave (*.csv)'
      read(*,1000) otfile
c
      open(12,file=otfile)
      write(12,*)'time(s),InAcc1,Acc1,Vel1,Dis1,InAcc2,Acc2,Vel2,Dis2
     *,InAcc3,Acc3,Vel3,Dis3'
      do 100 i=1,nn
        write(12,2000)time(i),ddy1(i),acc1(i),vel1(i),dis1(i),ddy2(i)
     *      ,acc2(i),vel2(i),dis2(i),ddy3(i),acc3(i),vel3(i),dis3(i)
  100 continue
 2000 format(f9.4,12(',',e15.7))
      close(12)
c
      stop
      end
c
c
c     Subroutine for Response of Single-Degree of System using the Nigam-Jennings method ***
c
      subroutine resp(h,w,dt,nn,ddy,dy,acc,vel,dis,nd,sa,sv,asv,sd)
c
      IMPLICIT REAL*8(A-H,O-Z)
      dimension ddy(nd),dy(nd),acc(nd),vel(nd),dis(nd)
c
	w2=w*w
	hw=h*w
	wd=w*dsqrt(1.d0-h*h)
	wdt=wd*dt
	e=dexp(-hw*dt)
	cwdt=dcos(wdt)
	swdt=dsin(wdt)
c
	a11= e*(cwdt+hW*swdt/wd)
	a12= e*swdt/wd
	a21=-e*w2*swdt/wd
	a22= e*(cwdt-hw*swdt/wd)
c
	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
c
	acc(1)=2.d0*h*w*ddy(1)*dt
	vel(1)=-ddy(1)*dt
	dis(1)=0.d0
	dx=vel(1)
	x=0.d0
	sa=0.d0
	sv=0.d0
        asv=0.d0
	sd=0.d0
c
        Ek=0.d0
        Ep=0.d0
        Ed=0.d0
        Et=0.d0
        write(15,1500)0.0,Ek,Ep,Ed,Et
 1500   format(f9.4,4(',',f9.4))
c
      do 110 m=2,nn
	dxf=dx
	xf=x
	ddym=ddy(m)
	ddyf=ddy(m-1)
c
c     * Relative Response Dis(x) and Vel(dx), and Absolute Acc(ddx) *
	x = a12*dxf+a11*xf+b12*ddym+b11*ddyf
	dx= a22*dxf+a21*xf+b22*ddym+b21*ddyf
	ddx=-2.d0*hw*dx-w2*x
c
c     * Absolute Vel(adx)
	adx= dx+dy(m)
c
	acc(m)=ddx
	vel(m)=dx
	dis(m)=x
c
	sa=dmax1(sa,dabs(ddx))
	sv=dmax1(sv,dabs(dx))
	asv=dmax1(asv,dabs(adx))
	sd=dmax1(sd,dabs(x))
c
        Ek=dx*dx/2.d0
        Ep=w2*x*x/2.d0
        Ed=Ed+hw*(dxf*dxf+dx*dx)*dt
        Et=Et-(ddyf*dxf+ddym*dx)*dt/2.d0
        write(15,1500)dt*dfloat(m-1),Ek,Ep,Ed,Et
c
  110 continue

      return
      end

