R

応答スペクトル

Rcpp パッケージ

（データ）気象庁
強震波形（平成12年(2000年)鳥取県西部地震）の
境港市東本町
（参考）
WANtaroHP (振動関係)の「加速度応答スペクトル計算プログラム (大崎先生による) 」
（WANtaroさんは「FFT を用いた地震応答スペクトル計算」のRcodeも公開されています。）
Rcpp 入門

OSはzorin linuxです。

sub01関数とsub02関数を定義。Rcppパッケージを使います。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
library(Rcpp)
###### sub01 #####
cppFunction('
  NumericMatrix sub01(List x) 
  {
    NumericVector acc=x["acc"];
    double dt=x["dt"];
    int nn=x["nn"];
    NumericMatrix veldis(nn,2);           
    for(int i=1; i < nn; i++){
    veldis(i,0)=veldis(i-1,0)+(acc[i-1]+acc[i])*dt/2.0;
    veldis(i,1)=veldis(i-1,1)+veldis(i-1,0)*dt+(acc[i-1]/3.0+acc[i]/6.0)*dt*dt;
    }
    return veldis;
  }
')
###### sub02 #####
cppFunction('
  NumericMatrix sub02(List x) 
  {
    NumericVector acc=x["acc"];
    int nn=x["nn"];
    double A11=x["A11"];
    double A12=x["A12"];
    double A21=x["A21"];
    double A22=x["A22"];
    double B11=x["B11"];
    double B12=x["B12"];
    double B21=x["B21"];
    double B22=x["B22"];
    double xf=x["xf"];
    double dxf=x["dxf"];
    double hw=x["hw"];
    double w2=x["w2"];
    NumericMatrix avdmax(nn,3);
    for(int m=1; m < nn; m++){
                avdmax(m,2)=A12*dxf+A11*xf+B12*acc[m]+B11*acc[m-1];
                avdmax(m,1)=A22*dxf+A21*xf+B22*acc[m]+B21*acc[m-1];
                avdmax(m,0)=-2.0*hw*avdmax(m,1)-w2*avdmax(m,2);
                dxf=avdmax(m,1);
                xf=avdmax(m,2);
	}
    return avdmax;
  }
')

データ取り込み、パラメータ等入力

1
2
3
4
5
6
7
8
9
10
#境港市東本町
seibu1<-read.csv("http://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/001006_tottori-seibu/dat/AA06E9E1.csv",skip=6,header=T)
h=0.05 #減衰定数
dt=0.01 #サンプリング周期
tk<-NULL
nnn=100
for (i in 1:nnn){
        tk[i]=10.0^(log10(2.0*dt)+(log10(10.0)-log10(2.0*dt))/nnn*(i-1))
    }
nt<-length(tk)

メイン

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
for (direction in 1:3){
if (direction==1) {acc<-seibu1$NS}
if (direction==2) {acc<-seibu1$EW}
if (direction==3) {acc<-seibu1$UD}
resacc<-NULL;resvel<-NULL;resdis<-NULL
vel<-NULL;dis<-NULL
nn<-length(acc)
accmax=0;velmax=0;dismax=0
vel[1]=0;dis[1]=0
#############################################
#
#Rcppを利用 sub01関数を定義
#for (i in 2:nn){
#            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
#	}
#
#############################################
veldis<-sub01(list(acc=acc,dt=dt,nn=nn))
vel<-veldis[,1]
dis<-veldis[,2]
accmax<-max(abs(acc))
velmax<-max(abs(vel))
dismax<-max(abs(dis))
tt=(nn-1)*dt
t=seq(0,tt,dt)
dis=dis*(3.0*tt-2.0*t)*t*t
sum=(dis[1]+dis[nn])/2.0
sum=sum+sum(dis[2:(nn-1)])
sum=sum*dt
a1=28.0/13.0/tt/tt*(2.0*vel[nn]-15.0/tt^5.0*sum)
a0=vel[nn]/tt-a1/2.0*tt
acmax=0
acc=acc-a0-a1*t
acmax<-max(abs(acc))
coef=accmax/acmax
acc<-acc*coef
for (i in 1:nt){
            w=2.0*pi/tk[i]
            w2=w*w
            hw=h*w
            wd=w*sqrt(1.0-h*h)
            wdt=wd*dt
            e=exp(-hw*dt)
            cwdt=cos(wdt)
            swdt=sin(wdt)
            A11=e*(cwdt+hw*swdt/wd)
            A12=e*swdt/wd
            A21=-e*w2*swdt/wd
            A22=e*(cwdt-hw*swdt/wd)
            ss=-hw*swdt-wd*cwdt
            cc=-hw*cwdt+wd*swdt
            s1=(e*ss+wd)/w2
            c1=(e*cc+hw)/w2
            s2=(e*dt*ss+hw*s1+wd*c1)/w2
            c2=(e*dt*cc+hw*c1-wd*s1)/w2
            s3=dt*s1-s2
            c3=dt*c1-c2
            B11=-s2/wdt
            B12=-s3/wdt
            B21=(hw*s2-wd*c2)/wdt
            B22=(hw*s3-wd*c3)/wdt
            accmax=2.0*hw*abs(acc[1])*dt
            velmax=abs(acc[1])*dt
            dismax=0.0
            dxf=-acc[1]*dt
            xf=0.0
#############################################
#
#Rcppを利用 sub02関数を定義
#for (m in 2:nn){
#                x=A12*dxf+A11*xf+B12*acc[m]+B11*acc[m-1]
#                dx=A22*dxf+A21*xf+B22*acc[m]+B21*acc[m-1]
#                ddx=-2.0*hw*dx-w2*x
#                if(accmax<=abs(ddx))accmax=abs(ddx)
#                if(velmax<=abs(dx))velmax=abs(dx)
#                if(dismax<=abs(x))dismax=abs(x)
#                dxf=dx
#                xf=x
#            			}
#
#############################################
avdmax<-sub02(list(acc=acc,nn=nn,A11=A11,A12=A12,B11=B11,B12=B12,A21=A21,A22=A22,B21=B21,B22=B22,xf=xf,dxf=dxf,hw=hw,w2=w2))
accmax=max(abs(avdmax[,1]))
velmax=max(abs(avdmax[,2]))
dismax=max(abs(avdmax[,3]))
            resacc[i]=accmax
            resvel[i]=velmax
            resdis[i]=dismax
        }
if (direction==1) {resNS<-data.frame(tk,resacc,resvel,resdis)}
if (direction==2) {resEW<-data.frame(tk,resacc,resvel,resdis)}
if (direction==3) {resUD<-data.frame(tk,resacc,resvel,resdis)}
}

グラフ化

  • インストールされているフォントの一覧表示 fc-list|less
  • plot関数のパラメータ panel.firstにablineが使える！
  • panel.first = abline(v=gridx,h=gridy,col=”lightgray”,lty =”dotted”)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#png("higashi.png",width=1000,height=600)
gridx=c(seq(2*dt,0.09,0.01),seq(0.1,0.9,0.1),seq(1,10,1))
gridy=c(seq(0.001,0.009,0.001),seq(0.01,0.09,0.01),seq(0.1,0.9,0.1),seq(1,9,1),seq(10,90,10),seq(100,900,100),seq(1000,10000,1000))
par(mfrow=c(1,3),family="TakaoPMincho")
#EWが最も大きい
plot(resEW$tk,resEW$resacc,type="l",log="xy",las=1,col="red",xlab="Period (sec)",ylab="",lwd=1.5,
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"))
lines(resNS$tk,resNS$resacc,col="blue",lwd=1.5)
lines(resUD$tk,resUD$resacc,col="green",lwd=1.5)
mtext(text=" h=0.05",side=3,line=0,adj=0)
title("加速度応答スペクトル")
legend(
    "bottomleft",
    c("NS","EW","EW"),
    col = c("blue","red","green"),
    lty = c(1,1,1),
    lwd = c(2,2,2),
    bg  = "white"
)
plot(resEW$tk,resEW$resvel,type="l",log="xy",,las=1,col="red",xlab="Period (sec)",ylab="",lwd=1.5,
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty="dotted"))
lines(resNS$tk,resNS$resvel,col="blue",lwd=1.5)
lines(resUD$tk,resUD$resvel,col="green",lwd=1.5)
mtext(text=" h=0.05",side=3,line=0,adj=0)
title("速度応答スペクトル")
legend(
    "bottomright",
    c("NS","EW","EW"),
    col = c("blue","red","green"),
    lty = c(1,1,1),
    lwd = c(2,2,2),
    bg  = "white"
)
plot(resEW$tk,resEW$resdis,type="l",log="xy",,las=1,col="red",xlab="Period (sec)",ylab="",lwd=1.5,
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty = "dotted"))
lines(resNS$tk,resNS$resdis,col="blue",lwd=1.5)
lines(resUD$tk,resUD$resdis,col="green",lwd=1.5)
mtext(text=" h=0.05",side=3,line=0,adj=0)
title("変位応答スペクトル")
legend(
    "bottomright",
    c("NS","EW","EW"),
    col = c("blue","red","green"),
    lty = c(1,1,1),
    lwd = c(2,2,2),
    bg  = "white"
)
par(mfrow=c(1,1))
#dev.off()