seismic

地球内部地震波線ソフトと画像 21.Sep.2012

地球内部の，地震の波線を書くソフトはネット上にはあまり見当たらないので，これを作ってみました．
画像は下記画像をクリック．教材としてご活用ください．

１次元地震波速度構造として，iasp91（Kennett & Engdahl,1991）を採用しました．

Traveltimes for global earthquake location and phase identification. B. L. N. Kennett' and E. R. Engdah12.

の中に深さ別の地震波速度の計算式があります．

下のコードの

func_v()

の部分です．

改良版のAK135（AKB135ではありません！）の文献は下記サイトに（ただしTableのみで計算式はありません）．
Download (3.9 MB)
Websiteは下記
http://rses.anu.edu.au/seismology/ak135/intro.html

＜地震波線描画ソフト＞

下記文献を参考に独自のコードで作成しました．変数宣言，ファイル書き出し部分などは省略してあります．
下記，桑原論文では，震波線経路の計算にアークサインを２度使っていますが，私のプログラムではSnellの法則の表示に震波線パラメータpを用いましたので，アークサインを開く回数が１回で済みました．

Ubuntu11.04上のegg（eggxライブラリ使用のコンパイラ）でコンパイルします．

#define R0 6371    /* radius of the Earth */
#define DT 0.0001 /* if increased -> line interupted */
#define I0 0.0       /* I0=0.0 caused bug */
#define V0 5.8

for (i=1;i<=IMAX;i++){

    flg = 0;
    ip=(0.2*i+0.0)*DEGRAD;
    th=0.0;
    vp=V0;
    rp=R0;

    p=rp*sin(ip)/vp;        /*　震波線パラメータ    */

    rpp=rp;

   while( rp <= R0 ){

        l=vp*DT;
        a=rp/R0;

        if (flg==1)
            r=rp+l*cos(ip);
        else
            r=rp-l*cos(ip);

        th0=l*sin(ip)/r;
        th=th+th0;

        func_v();

        if ( v*p/r > 0.99999999999 ){

            flg=1;
            rpp=r;
          r=rp+l*cos(ip);

            func_v();
        }

        in=asin(v*p/r);

        x=r*cos(th);
        y=r*sin(th);
        xx=-x*MX+X00;
        yy=y*MY+Y00;
        pset(w,xx,yy);
    }
｝

// read iasp91_velocity table

func_v()

{
        if (r < 1217.1)
            v=11.24094-4.09689*a*a;
        else if (r < 3482.0)
            v=10.03904+3.75665*a-13.67046*a*a;
        else if (r < 3631.0)
            v=14.49470-1.47089*a;
        else if (r < 5611.0)
            v=25.1486-41.1538*a+51.9932*a*a-26.6083*a*a*a;
        else if (r < 5711.0)
            v=25.96984-16.93412*a;
        else if (r < 5961.0)
            v=29.38896-21.40656*a;
        else if (r < 6161.0)
            v=30.78765-23.25415*a;
        else if (r < 6251.0)
            v=25.41389-17.69722*a;
        else if (r < 6336.0)
            v=8.78541-0.74953*a;
        else if (r < 6351.0)
            v=6.50;
        else
            v=5.80;
}

＜参考文献＞
参考にしたのは，下記の文献ですが，PDF化されていません．

桑原浩二：地震波経路の自動描画（昭和６１年度　東レ理科教育賞受賞作品集）
http://www.toray.co.jp/tsf/rika/hig_003.html