電離層遅延（放送暦）

今日はGPS放送暦の電離層遅延について少し調べた。
IS-GPS-200Eを読んでみただけだけど、混乱したのは、角度の単位が式で異なること。間違っているかもしれないけれども、また混乱しそうなので、一応単位をメモしておく。

F=1+16*(0.53-EL)^3 (EL:semi-circle)
phi=0.0137/(EL+0.11)-0.022 (EL:semi-circle)
cos(AZ), sin(AZ) (AZ:radian)

もちろん三角関数の中身がsemi-circleになるはずはないのだけど、Fなんかはどっちだ?となってしまった。

さて、IS-GPS-200Eにのっているモデルはklobucharモデルというやつなのだけど、元論文を読んだことはない。だけど、式を見ると、「何時にこの座標でこれだけの遅延量がある」という記述というより、「丸っこい遅延の塊が移動する様子」が記述されているように見える。なるほど・・・、電離層遅延は伝搬していくということなのかな?

勉強はこれでやったことにして、国土地理院からダウンロードしたRINEXファイルを確認してみると・・・、なんと、電離層遅延パラメータが載っていないではないですか！RINEXファイルフォーマットを読むと、オプション項目なんですね。

ふーむ。IGS点のデータを持ってくるか・・・。GEONETは定点観測なので、相対測位が前提ということなのかなあ。

まあ、とりあえずIS-GPS-200Eの該当部分を読んでプログラムにしてみたので、まだ動くかはわからないけど、超ベータ版ということで自分のために載せておくことにした。

double brdion(const double az, const double el, const double phi_u, const double rambda_u, const gtime_t time, const double *ion_model) {
    double dt;  /* ionospheric delay (sec) **not correction value** */

    double phi_i0, phi_i, rambda_i, phi_m, psi;
    double amp, per, f, x;
    gtime_t tl;
    const double R2SC=R2D/180.0;
    const double D2SC=1/180.0;
    const double SC2R=PI;
    int i;

    amp=per=0.0;

    /* estimate phi */
    psi = 0.0137/(el*R2SC) - 0.022;
    phi_i0 = phi_u*D2SC + psi*cos(az); /* az(rad) */
    
    if (fabs(phi_i0) <= 0.416) {
        phi_i = phi_i0;
    }
    else if (phi_i0 > 0.416) {
        phi_i = 0.416;
    }
    else {
        phi_i = -0.416;
    }

    rambda_i = rambda_u*D2SC + (psi*sin(az)/cos(phi_i*SC2R));
    phi_m = phi_i + 0.064*cos(rambda_i*SC2R - 1.617);

    f = 1.0 + 16.0*(0.53-el*R2SC)*(0.53-el*R2SC)*(0.53-el*R2SC);

    /* ion_model: iono model parameters {a0,a1,a2,a3,b0,b1,b2,b3} */
    for (i=0;i<4;i++) {
        amp += ion_model[i]*pow(phi_m,i);
        per += ion_model[4+i]*pow(phi_m,i);
    }
    if (amp < 0.0) {
        amp = 0.0;
    }
    if (per < 72000.0) {
        per = 72000.0;
    }
    tl = timeadd(time, 12.0*3600.0*rambda_i);
    x = 2.0*PI*(tl.time+tl.sec-14.0*3600.0)/per;

    if (fabs(x) <= 1.57) {
        dt = f*5.0e-9;
    }
    else {
        dt = f*(5.0e-9 + amp*(1-x*x/2.0+x*x*x*x/24.0));
    }

    return dt;
}