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)
さて、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;
}
