module defpi implicit none real(8),parameter::pi=3.14159265358979323846D0 end module defpi program f90_MAKEWAVE !***************************************************************** ! artificial earthquake motion !***************************************************************** use defpi implicit none integer::ndata,i,k,nn real(8)::accmax,velmax,dismax real(8),allocatable::accinp(:) !Acceleration fo original wave real(8),allocatable::acc(:) !Acceleration real(8),allocatable::vel(:) !Velocity real(8),allocatable::dis(:) !Displacement real(8),allocatable::resaccinp(:) !Response of acceleration of original wave real(8),allocatable::resacc(:) !Response of acceleration real(8),allocatable::resvel(:) !Response of velocity real(8),allocatable::resdis(:) !Response of displacement real(8),allocatable::xr(:) !Complex Fourier series real part real(8),allocatable::xi(:) !Complex Fourier series imaginary part real(8),allocatable::frq(:),Tk(:) !frequency,period real(8),allocatable::rk(:) !spectrum ratio (target wave/calculated wave) real(8)::err1,err2,eps real(8)::sum,tt,x1,y1,x2,y2 integer::iii,maxita !Limit value of iteration real(8)::hh !Damping factor real(8)::dt !Time increment integer::nts !Number of data of target spectrum real(8),allocatable::perts(:) !Period of target spectrum real(8),allocatable::accts(:) !Acceleration of target spectrum (function of period) real(8),allocatable::tsacc(:) !Acceleration of target spectrum (function of frequency) character::fnameR1*50,fnameR2*50,fnameW1*50,fnameW2*50 character::dummy*50,linebuff*300 maxita=1000 !Set limit value of iteration hh=0.05 !Set damping factor eps=0.05 !Set tolerance call getarg(1,fnameR1) !Input of original wave (acceleration) call getarg(2,fnameR2) !Input of target spectrum call getarg(3,fnameW1) !output of calculated wave call getarg(4,fnameW2) !Output of response spectrum open(11,file=fnameR1,status='old') read(11,'(a)') linebuff read(11,*) dummy,dt read(11,*) dummy,ndata nn=2 !set number of data to powers of 2 do nn=nn*2 if(ndata*1<=nn)exit end do allocate(accinp(1:nn)) allocate(acc(1:nn)) allocate(vel(1:nn)) allocate(dis(1:nn)) allocate(xr(1:nn)) allocate(xi(1:nn)) allocate(rk(1:nn/2+1)) allocate(frq(1:nn/2+1)) allocate(Tk(1:nn/2+1)) allocate(resacc(1:nn/2+1)) allocate(resvel(1:nn/2+1)) allocate(resdis(1:nn/2+1)) allocate(tsacc(1:nn/2+1)) allocate(resaccinp(1:nn/2+1)) do k=1,nn/2+1 !set frequency frq(k)=1.0D0/dt/dble(nn)*dble(k-1) end do frq(1)=0.5D0*frq(2) do k=1,nn/2+1 Tk(k)=1.0D0/frq(k) end do do i=1,nn accinp(i)=0.0D0 end do do i=1,ndata read(11,*) accinp(i) end do close(11) !Input of target spectrun open(11,file=fnameR2,status='old') read(11,'(a)') linebuff read(11,*) dummy,nts read(11,*) dummy,dummy allocate(perts(1:nts)) allocate(accts(1:nts)) do i=1,nts read(11,*) perts(i),accts(i) end do close(11) !target acceleration at calculated frequency do k=1,nn/2+1 !frequency tt=1.0D0/frq(k) if(tt=k) j=j-k k=k/2 if(k==0)exit end do j=j+k end do end subroutine FFT subroutine NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,dt,ff) !***************************************************************** !Newmark Beta Method !***************************************************************** real(8),intent(out)::acc2,vel2,dis2 real(8),intent(in)::acc1,vel1,dis1 real(8),intent(in)::am,ac,ak,dt,ff ! real(8)::beta=1.0D0/4.0D0 !Average acceleration method real(8)::beta=1.0D0/6.0D0 !Linear acceleration method acc2=(ff-ac*(vel1+0.5D0*dt*acc1)-ak*(dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1))& /(am+ac*0.5D0*dt+ak*beta*dt*dt) vel2=vel1+0.5D0*dt*(acc1+acc2) dis2=dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1+beta*dt*dt*acc2 end subroutine NEWMARK subroutine CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis) !***************************************************************** !Calculation of response !***************************************************************** integer,intent(in)::ndata,nnn real(8),intent(in)::dt,hh real(8),intent(in)::acc(1:ndata),Tk(1:nnn) real(8),intent(out)::resacc(1:nnn),resvel(1:nnn),resdis(1:nnn) integer::i,j,k,md real(8)::acc2,vel2,dis2 real(8)::acc1,vel1,dis1 real(8)::am,ac,ak,f1,f2,fm,ddt real(8)::racc,rvel,rdis real(8)::acmax,vemax,dimax md=5 ddt=dt/dble(md) do k=1,nnn am=1.0D0 !Mass ak=4.0D0*pi*pi*am/Tk(k)/Tk(k) !Stiffness ac=2.0D0*hh*sqrt(ak*am) !Damping constant acc1=0.0D0 vel1=0.0D0 dis1=0.0D0 acmax=acc1 vemax=vel1 dimax=dis1 do i=1,ndata f1=0.0D0 if(1<=i)f1=-am*acc(i-1) f2=-am*acc(i) do j=1,md fm=f1+(f2-f1)/dble(md)*dble(j) call NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,ddt,fm) acc1=acc2 vel1=vel2 dis1=dis2 end do racc=acc2+acc(i) !Absolute acceleration rvel=vel2 !Relative velocity rdis=dis2 !Relative displacement if(acmax