gresp.for

一个组合程序

**********************************
*       Program res_spec.f
*       Generates the RESPONSE SPECTRA of a LINEAR Undamped
*       SDOF System to any arbitrary applied GROUND ACCELERATION
**********************************
*       The program uses the Newmark-B integration scheme
*       to evaluate the response at each discrete time step
**********************************
*       Program Author : Sanjoy Chakraborty : 09CE
**********************************
*       Compile using : fl -AH -MW gresp.for /link /SEG:1000
**********************************

	program eq_resp
	implicit real *8 (a-h,o-z)
	real *8 k,m
	character *22 fgraccl
	dimension x(100000),x1(100000),x2(100000),ft(5000),tt(5000),
     1  p(100000),damprat(20),sd(5000,20),sv(5000,20),sa(5000,20)
	data pi/3.141592654/

	call input(m,rm,tt,ft,nfor,tend,h,tpstart,tpend,th,
     1  ndamprat,damprat,ndigitdat,fgraccl)

	call inter_f(nfor,tt,ft,p,h,tend)

**********************************
*       Loop for computing the value of maximum DISPLACEMENT 
*       response for successive values of system DAMPING Ratio
*       and TIME-Period
*       ndamprat = number of damping levels at which the spectra
*                  are to be generated
*       damprat(i) = values of damping ratio specified (i=1,ndamprat)
*       t = time period (t=tpstart,tpend, in increments of th)
*       Sv,Sa = values of SPECTRAL velocity and acceleration respectively,
*               computed from the maxm. displacement response, Sd
**********************************

	do 701 j=1,ndamprat

	  write(*,703)j,damprat(j)
703       format(' Damping Step :',i3,'  Damping Ratio =',f10.3)

	  i=1
	  do 701 t=tpstart,tpend,th
	    k=m*(2.*pi/t)**2
	    omega=sqrt(k/m)
	    c=2.*m*omega*damprat(j)
	    call resp(x,x1,x2,m,c,k,h,tend,rm,p,ic)
	    call maxresp(h,ic,x,x1,x2,xm,x1m,x2m)
	    sv(i,j)=omega*abs(xm)
	    sa(i,j)=omega**2*abs(xm)
	    sd(i,j)=abs(xm)
	    i=i+1
701       continue

	write(*,704)
704     format(///,' Finished Generating Response Spectra :',
     1  ' Writing Output',//)

	call output(sd,sv,sa,ndamprat,tpstart,tpend,th)

	stop
	end

**********************************
	subroutine input(m,rm,tt,ft,nfor,tend,h,
     1  tpstart,tpend,th,ndamprat,damprat,ndigitdat,fgraccl)
**********************************
*       System parameters are input from Unit:1 
*       Ground acceleration values input from Unit:2 and then
*       converted to time vs. effective force data pairs
**********************************

	implicit real *8 (a-h,o-z)
	real *8 m
	character *22 fin,fgraccl
	dimension tt(5000),ft(5000),damprat(20)
	rm=.1e12

	write(*,710)
710     format(//,' Program for Generating Earthquake Response Spectra'
     1  ,//,' Program Author : Sanjoy Chakraborty, Auburn University',
     2  ///,' Input File Name :')
	read(*,'(a)')fin
	open(1,file=fin,status='old')

	write(*,'(//a)')' Reading Solution Parameters :'
	read(1,*)m,tend,h
	read(1,*)tpstart,tpend,th
	read(1,*)ndamprat
	read(1,*)(damprat(i),i=1,ndamprat)
	read(1,'(a22)')fgraccl

	open(2,file=fgraccl,status='old')
	write(*,'(//a//)')' Reading Ground Acceleration Records :'

	nfor=1
110     continue
	  read(2,*,err=111)tt(nfor),accl
	  ft(nfor)=m*accl
	  nfor=nfor+1
	goto 110
111     nfor=nfor-1
	ndigitdat=nfor

	close(1)
	close(2)

	return
	end

**********************************
	subroutine inter_f(nf,tt,ft,p,h,tend)
**********************************
*       Used to interpolate for the forcing function values
*       at the times at which the solution is sought
**********************************

	implicit real *8 (a-h,o-z)
	dimension tt(5000),ft(5000),p(100000)

	ic=1
	do 1001 t=0.,tend,h
	  do 1002 i=1,nf-1
	    if(t.ge.tt(i).and.t.le.tt(i+1))then
		p(ic)=ft(i)+(ft(i+1)-ft(i))*(t-tt(i))/(tt(i+1)-tt(i))
		goto 1003
	    endif
1002      continue
1003    ic=ic+1
1001    continue

	return
	end

**********************************
	subroutine resp(x,x1,x2,m,c,k,h,tend,rm,p,ic)
**********************************
*       Computes the system response (displacement, velocity
*       and acceleration) at each time step
**********************************
*       The loop (=i) parameter is used to determine the 
*       exact location of the system on the hysteretic loop
*       describing the resistance function Rm
*       loop =1  -->  elastic loading stage
*       loop =2  -->  plastic loading stage
*       loop =3  -->  elastic rebound stage
*       loop =4  -->  plastic rebound stage
*       xmax, xmin, and xlim are parameters used to control
*       the physical route of the resistance function along
*       the hysteretic loop
*       xmax = max. +ve plastic deformation (loop=2)
*       xmin = max. -ve plastic deformation (loop=4)
**********************************

	implicit real *8 (a-h,o-z)
	real *8 k,m,kelas,kplas
	dimension x(100000),x1(100000),x2(100000),p(100000)

	x(1)=0.
	x1(1)=0.
	x2(1)=p(1)/m

	xel=rm/k
	a1=3./h
	a2=6./h
	a3=h/2.
	a4=6./h**2
	kelas=k+a4*m+a1*c
	kplas=a4*m+a1*c
	xlim=xel
	xmin=-xel
	loop=1
	ic=2

	do 501 t=h,tend,h

	  if(loop.eq.1)then
	    call res1(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=-rm-(xmin-x(ic))*k
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x(ic).ge.xlim)then
	      loop=2
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.2)then
	    call res1(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=rm
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x1(ic).le.0.)then
	      loop=3
	      xmax=x(ic)
	      xlim=x(ic)-2.*xel
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.3)then
	    call res1(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=rm-(xmax-x(ic))*k
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x(ic).le.xlim)then
	      loop=4
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.4)then
	    call res1(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=-rm
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x1(ic).ge.0.)then
	      loop=1
	      xlim=x(ic)+2.*xel
	      xmin=x(ic)
	    endif
	    ic=ic+1
	    goto 501

	  endif

501     continue

	return
	end

**********************************
	subroutine res1(k,p,x,x1,x2,m,c,ic,a2,a3,a1)
**********************************
*       Used to compute displacement and velocity
*       at any time step
**********************************

	implicit real *8 (a-h,o-z)
	real *8 k,m
	dimension p(100000),x(100000),x1(100000),x2(100000)

	dps=p(ic)-p(ic-1)+x1(ic-1)*(a2*m+3.*c)+x2(ic-1)*
     1      (3.*m+a3*c)
	dx=dps/k
	dx1=a1*dx-3.*x1(ic-1)-a3*x2(ic-1)
	x(ic)=x(ic-1)+dx
	x1(ic)=x1(ic-1)+dx1

	return
	end

**********************************
	subroutine maxresp(h,ic,x,x1,x2,xm,x1m,x2m)
**********************************
*       Computes the value of Maximum Response at each value
*       of TIME PERIOD of the system, then writes the output
*       to Unit:3 in the format Td/T vs. (DLF)max
**********************************

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

	dimension x(100000),x1(100000),x2(100000)

	xm=x(1)
	txm=0.
	x1m=x1(1)
	tx1m=0.
	x2m=x2(1)
	tx2m=0.

	do 605 i=2,ic-1
	  if(abs(x(i)).gt.abs(xm))then
	    xm=x(i)
	    txm=float(i-1)*h
	  endif
	  if(abs(x1(i)).gt.abs(x1m))then
	    x1m=x1(i)
	    tx1m=float(i-1)*h
	  endif
	  if(abs(x2(i)).gt.abs(x2m))then
	    x2m=x2(i)
	    tx2m=float(i-1)*h
	  endif
605     continue

	return
	end

**********************************
	subroutine output(sd,sv,sa,ndamprat,tpstart,tpend,th)
**********************************
*       Prints the Response Spectra Output to 3 files
*       Unit:3 -> period vs. maxm. displacement @ each damping level
*       Unit:4 -> period vs. maxm. pseudo-velocity
*                 Sv = omega*Sd
*       Unit:5 -> period vs. maxm. pseudo-accleration
*                 Sa = omega*Sd
**********************************

	implicit real *8 (a-h,o-z)
	dimension sd(5000,20),sv(5000,20),sa(5000,20)

	open(3,file='sd.out',status='unknown')
	open(4,file='sv.out',status='unknown')
	open(5,file='sa.out',status='unknown')

	i=1
	do 651 tp=tpstart,tpend,th
	  write(3,652)tp,(sd(i,j),j=1,ndamprat)
	  write(4,652)tp,(sv(i,j),j=1,ndamprat)
	  write(5,652)tp,(sa(i,j),j=1,ndamprat)
652       format(f9.3,20e15.7)
	  i=i+1
651     continue

	close(3)
	close(4)
	close(5)

	return
	end