現在、fortranにて、変位応答スペクトルを作成しようと頑張っているのですが、
その前段階として、Newmarkのβ法により
変位、速度、加速度を求めるプログラムを作成いたしました。
こんな感じで。
real h,w,dt,Sacc(0:3000),beta,pi,T
real acc(0:3000),vel(0:3000),def(0:3000)
integer i
open(1,file='ELCNN.csv',status='old')
dt =0.02
pi =3.14159265358979
write(*,*) 'beta?'
read(5,*) beta
write(*,*) 'T?'
read(5,*) T
write(*,*) 'h?'
read(5,*) h
w=2*pi/T
open(2,file='acc.csv',status='new')
open(3,file='vel.csv',status='new')
open(4,file='def.csv',status='new')
do 110 i=0,2687
read(1,*) Sacc(i)
110 continue
acc(0)=Sacc(0)
vel(0)=0.
def(0)=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.)
write(2,100) acc(i)
write(3,100) vel(i)
write(4,100) def(i)
100 format(f15.6)
10 continue
stop
end
コピペしたので、多少ずれてますが。
実は、かなり前に、計算結果を出すことはできていたのですが、
久しぶりに回してみたら、うまくできなくて。
よく調べてみても、計算式は合ってるし。
結果的におかしかったのは、「ｉ」に０を含ませていたこと。
だから、最初の条件に0:3000の０を記入がされておりませんでした。
基本的に、
整数の後にはピリオドをつける
分数は使わない
などのことを心がけたほうが良いということを教わりました。