地盤の増幅スペクトルを計算します。大崎順彦（1994）のFESPを参考にしました。ですので、大地震時にせん断ひずみが大きくなって剛性率が低下するという「非線形性」にまでは対応できません。
一応、オリジナルのFESPおよび吉田望らのDYNEQ320（弾性材料使用）との比較・検証はやってみました。DYNEQは散乱減衰の設定の仕方が少し異なっているので、近似値を使用しました。
しばらくの間、このスクリプトを使って色々試してみようと思っています。
#!/bin/bash temp=`mktemp ./XXXX` scilab -nw << _SCI_ >$temp function [z]=acc_resp(f,b_cond,obj,ref) // 大崎（1994）§9のFESP.forをSCILAB用に // 書き直したもの // 加速度-加速度応答を求める // [引数] // f : 周波数（実数） // b_cond : 境界条件(文字列） 'inner' or 'expose' // obj : 対象層(整数) // ref : 基準層(整数) // [その他] // nlayer(I),thick(ARRAY),gamma_w(ARRAY),G(ARRAY), // h_alpha(ARRAY),h_beta(ARRAY)はグローバル変数 // として引き渡す // [返値] // 周波数応答関数 Zobj/Zref を返す // 2018年2月28日 作成 // omega=2*%pi*f; // 円振動数 h=h_alpha./omega+h_beta; // 減衰定数の算出 gb=G.*(1+2*h*%i); // 複素剛性の算出 P=sqrt(gamma_w./gb/ACC)*%i; // 伝播定数の算出 for layer=1:nlayer-1 R(layer)=gb(layer)/gb(layer+1)*P(layer)/P(layer+1); // インピーダンス比の算出 end A(1)=1.; B(1)=1.; for layer=1:nlayer-1 ec=exp(P(layer)*omega*thick(layer)); ea=A(layer)*ec; eb=B(layer)/ec; A(layer+1)=((1.+R(layer))*ea+(1.-R(layer))*eb)/2.; B(layer+1)=((1.-R(layer))*ea+(1.+R(layer))*eb)/2.; end if b_cond=='inner' then z=(A(obj)+B(obj))/(A(ref)+B(ref)); // 内部境界面 elseif b_cond=='expose' then z=(A(obj)+B(obj))/(2*A(ref)); // 解放基盤面 else error("境界条件が不適当です"); exit; end endfunction nlayer=4; //層の数 thick=[3.8 3.2 3.9]; // 層厚(m) gamma_w=[1.50 1.67 1.85 1.95]; // 単位体積重量(tf/m3) G=[1200. 2900. 5700. 50000.]; // せん断弾性係数(tf/m2) h_alpha=[2. 2. 2. 2.]; // マクスウェル型散乱減衰（地盤の減衰） h_beta=[0.02 0.02 0.02 0.02]; // 履歴減衰型材料減衰（土の減衰） ACC=9.8; //重力加速度(m/sec2) obj=1; // 対象層 ref=4; // 基準層 if size(thick,'c')+1 <> nlayer then error("層厚の数が層数と一致しません"); end if size(gamma_w,'c') <> nlayer then error("単位体積重量の数が層数と一致しません"); end if size(G,'c') <> nlayer then error("せん断弾性係数の数が層数と一致しません"); end if size(h_alpha,'c') <> nlayer then error("散乱減衰の数が層数と一致しません"); end if size(h_beta,'c') <> nlayer then error("材料減衰の数が層数と一致しません"); end [str,ier]=lasterror(); if ier==10000 then exit; end df=50/128; freq0=[0:df:50]; for i=1:size(freq0,'c') if i==1 then z=1; else f=freq0(i); // 周波数 z=acc_resp(f,'inner',obj,ref); end mfprintf(6,'%f %f %f\n',freq0(i),real(z),imag(z)); end _SCI_ cat fesp_out.txt |\ awk '{print $1,sqrt($2*$2+$3*$3)}' |\ gmt psxy -JX15c/10c -R0/21/0/8 -BnWeS+t"増幅スペクトルの比較" \ -Bxa5f1+l"周波数(Hz)" -Bya2f1+l"増幅率" -W2p -K > resp.ps cat dyneq_out.txt |\ sed -e "1,2d" |\ gmt psxy -J -R -B -W1p,red -O -K >> resp.ps cat $temp |\ sed -e "/^ *$/d" |\ awk '{print $1,sqrt($2*$2+$3*$3)}' |\ gmt psxy -J -R -B -S+0.4c -W1p,blue -O -K >> resp.ps rm $temp gmt gmtset FONT_TITLE 16p,GothicBBB-Medium-UniJIS-UTF8-H gmt gmtset FONT_LABEL 12p,GothicBBB-Medium-UniJIS-UTF8-H gmt gmtset FONT_ANNOT_PRIMARY 10p,GothicBBB-Medium-UniJIS-UTF8-H cat << + |\ gmt pslegend -R -J -Dx11c/9c+w3c/1.6c+jTL -F+p+gwhite -O >> resp.ps N 1 S 0.2c - 0.4c 0/0/0 2p,black 0.8c 大崎FESP G 0.1c S 0.2c - 0.4c 255/0/0 1p,red 0.8c DYNEQ320 G 0.1c S 0.2c + 0.4c 0/0/255 1p,blue 0.8c SCILAB +
