# Calculation of response spectrum RSPC<-function(nnn,dt,damp,n,ainp,ndata){ aa <- fft(ainp) a1 <- Re(aa)/n a2 <- Im(aa)/n ww <- numeric(n) for(i in 1:n){ if(i<=n/2+1)ww[i] <- 2*pi*(i-1)/n/dt if(n/2+2<=i)ww[i] <- -2*pi*(n-i+1)/n/dt } tk <- numeric(nnn) dism <- numeric(nnn) velm <- numeric(nnn) accm <- numeric(nnn) for(kkk in 1:nnn){ tk[kkk] <- 10^(log10(2.0*dt)+(log10(10.0)-log10(2.0*dt))/nnn*(kkk-1)) omega0 <- 2*pi/tk[kkk] h1 <- ww*ww-omega0*omega0 h2 <- -2*damp*omega0*ww ah1 <- (a1*h1+a2*h2)/(h1*h1+h2*h2)*n ah2 <- (a2*h1-a1*h2)/(h1*h1+h2*h2)*n dis <- Re(fft(complex(re=ah1,im=ah2),inverse=TRUE))/n vel <- Re(fft(complex(re=-ww*ah2,im=ww*ah1),inverse=TRUE))/n acc <- Re(fft(complex(re=-ww*ww*ah1,im=-ww*ww*ah2),inverse=TRUE))/n dism[kkk] <- max(abs(dis[1:ndata])) velm[kkk] <- max(abs(vel[1:ndata])) accm[kkk] <- max(abs(acc[1:ndata]+ainp[1:ndata])) } return(invisible(list(tk,accm,velm,dism))) } CRAC <- function(ndata,dt,acc){ accmax <- max(abs(acc)) vel <- numeric(ndata) dis <- numeric(ndata) for(i in 2:ndata){ vel[i] <- vel[i-1]+(acc[i-1]+acc[i])*dt/2.0 dis[i] <- dis[i-1]+vel[i-1]*dt+(acc[i-1]/3.0+acc[i]/6.0)*dt*dt } tn <- (ndata-1)*dt tt <- seq(0,ndata-1,by=1) tt <- tt*dt dis <- dis*(3.0*tn-2.0*tt)*tt*tt ss <- ((dis[1]+dis[ndata])/2.0+sum(dis[2:ndata-1]))*dt a1 <- 28.0/13.0/tn/tn*(2.0*vel[ndata]-15.0/tn^5.0*ss) a0 <- vel[ndata]/tn-a1/2.0*tn accr <- acc-a0-a1*tt acmax <- max(abs(accr)) accr <- accr*accmax/acmax return(accr) } nnn <- 200 dt <- 0.01 damp <- 0.05 for(iii in 1:3){ if(iii==1)fnameR="dat_acc_TCGH16_EW2.csv" if(iii==2)fnameR="dat_acc_TCGH16_NS2.csv" if(iii==3)fnameR="dat_acc_TCGH16_UD2.csv" if(iii==1)fnameW="dat_spc_1.txt" if(iii==2)fnameW="dat_spc_2.txt" if(iii==3)fnameW="dat_spc_3.txt" d <- scan(fnameR,skip=3) k <- length(d) acc <- CRAC(k,dt,d) n <- 2 while(n