変位応答スペクトル。
ここでは、1質点における、ある固有周期に対する
変位の最大値を出す作業を行う。
横軸に固有周期
縦軸に変位の最大値を対応させれば、
変位応答スペクトルのグラフが出せる。
以下のプログラムは
エルセントロのNS地震波のデータ「NLCNN」を読み込んでいる。
それにより、Newmarkのβ法を用いて、変位・速度・加速度を出している。
real h,w,dt,Sacc(0:3000),beta,pi,T
real acc(0:3000),vel(0:3000),def(0:3000),def2(0:3000)
real defmax(0:3000)
integer i
open(1,file='ELCNN.csv',status='old')
acc(0)=Sacc(0)
vel(0)=0.
def(0)=0.
dt =0.02
pi =3.14159265358979
write(*,*) 'beta?'
read(5,*) beta
write(*,*) 'h?'
read(5,*) h
do 110 i=0,2687
read(1,*) Sacc(i)
110 continue
do 200 j=1,250
T=j*0.1
w=2*pi/T
defmax(j)=0
do 10 i=0,2687
acc(i+1)=-((Sacc(i)+(w**2.)*def(i)+(2.*h*w+(w**2.)*dt)*vel(i))
& +(h*w*dt+(0.5-beta)*(w**2.)*(dt**2.))*acc(i))
& /(1.+h*w*dt+beta*(w**2)*(dt**2))
vel(i+1)=vel(i)+(acc(i+1)+acc(i))*dt/2.
def(i+1)=def(i)+vel(i)*dt+(0.5-beta)*acc(i)*(dt**2.)
& +beta*acc(i+1)*(dt**2.)
if(def(i+1).LT.0.) then
def2(i)=-def(i)
else
def2(i)=def(i)
endif
if(def2(i).gt.defmax(j)) then
defmax(j)=def2(i)
endif
10 continue
open(2,file='def.csv',status='new')
write(2,100) defmax(j)
100 format(f15.6)
200 continue
stop
end
また、コピペのため、ずれております。