c
c  *** 	Calcuration of Response Spectra  ***
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 Response Acceleration (gal)
c   sv: Maximum Response Velocity     (cm/sec)
c   sd: Maximum Response Displacement (cm)
c
c  Dimension Data:
c   ddy: Ground Acceleration (gal)
c   acc: Response Acceleration (gal)
c   vel: Response Velocity     (cm/sec)
c   dis: Response Displacement (cm)
c
      program main
      parameter (mm=66000)
      dimension ddy1(mm),ddy2(mm),ddy3(mm),time(mm)
      dimension acc(mm),vel(mm),dis(mm)
      CHARACTER INFILE*30,otfile*30,head1*20,head2*20,head3*20,head4*20
c
      write(*,*)
      write(*,*)'Enter File Name for Ground Motion Acceleraration:'
      read(*,1000) infile
 1000 format(a20)
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
c
      write(*,*)
      write(*,*)'Enter Damping Factor'
      read(*,*)h

      write(*,*)
      write(*,*)'Enter Output File Name (*.csv)'
      read(*,1000) otfile

      open(11,file=otfile)
      write(11,*)'period(s), acc 1, acc 2, acc 3, vel 1, vel 2, vel 3, d
     *is 1, dis 2, dis 3'
      do 100 j=1,300
        period=0.05*j 
	omega=6.283185/period
	call resp(h,omega,dt,nn,ddy1,acc,vel,dis,mm,sa1,sv1,sd1)
	call resp(h,omega,dt,nn,ddy2,acc,vel,dis,mm,sa2,sv2,sd2)
	call resp(h,omega,dt,nn,ddy3,acc,vel,dis,mm,sa3,sv3,sd3)
	write(11,2000) period,sa1,sa2,sa3,sv1,sv2,sv3,sd1,sd2,sd3
  100 continue
      close(11)
 2000 format(f8.4,9(',',f8.2))

      stop
      end
c
c
c	Subroutine for Response of Single-Degree of System
c
	subroutine resp(h,w,dt,nn,ddy,acc,vel,dis,nd,sa,sv,sd)
c
	dimension ddy(nd),acc(nd),vel(nd),dis(nd)
c
	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
c
	do 110 m=2,nn
	dxf=dx
	xf=x
	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

