MP線形結合の計算（２）

昨日MP線形結合の計算をしてみたけども、出力してみるとやけに大きな値が出る衛星がある。

なんだろうな、と思ったらGLONASSだった。GLONASSのFDMA波長にきちんと衛星毎に合わせる必要があるのだった。

少し計算を分けてみたのでメモしておく。
例によってrtklibに便利なsatwavelenという関数があることに気づいたので、早速使ってみる。

/* compute mp linear combination for GPS and QZS */
static int mp_gps(const double P1, const double P2, const double L1, const double L2, double *mp1, double *mp2) {
    double a=FREQ1*FREQ1/(FREQ2*FREQ2);
    double lam1=CLIGHT/FREQ1;
    double lam2=CLIGHT/FREQ2;

    *mp1 = P1 - (1+2/(a-1))*L1*lam1 + (2/(a-1))*L2*lam2;
    *mp2 = P2 - (2*a/(a-1))*L1*lam1 + (2*a/(a-1)-1)*L2*lam2;

    return 0;
}

/* compute mp linear combination for GLONASS FDMA */
static int mp_glo(int satnum, const nav_t *nav, const double P1, const double P2, const double L1, const double L2, double *mp1, double *mp2) {
    double a, lam1, lam2;

    lam1 = satwavelen(satnum, 0, nav);
    lam2 = satwavelen(satnum, 1, nav);
    a = lam2*lam2/(lam1*lam1); /* f1^2/f2^2 */

    *mp1 = P1 - (1+2/(a-1))*L1*lam1 + (2/(a-1))*L2*lam2;
    *mp2 = P2 - (2*a/(a-1))*L1*lam1 + (2*a/(a-1)-1)*L2*lam2;

    return 0; 
}

aは、f1^2/f2^2だけど、波長だと分母分子が入れ替わります。
satwavelenがすごく便利でした。このコード例では毎回呼び出すのでちょっと無駄な処理が多いんですが。

出力：
----- epoch start -----
sat num: 17, timetag, PRN, MP1, MP2
 0: 2012/07/09 00:00:00   9       71.8339      100.7876
 1: 2012/07/09 00:00:00  41     -242.7161     -297.4944
 2: 2012/07/09 00:00:00  14      137.7746      184.8548
 3: 2012/07/09 00:00:00  40     -242.3190     -297.2934
 4: 2012/07/09 00:00:00  15      119.1687      160.6320
 5: 2012/07/09 00:00:00  42      377.8123      480.3067
 6: 2012/07/09 00:00:00  31       61.5032       91.2867
 7: 2012/07/09 00:00:00  12       35.9511       58.4904
 8: 2012/07/09 00:00:00  33      183.4078      236.6984
 9: 2012/07/09 00:00:00  22       28.5107       45.0275
10: 2012/07/09 00:00:00  25       21.7292       46.4806
11: 2012/07/09 00:00:00  18       77.3021      105.7750
12: 2012/07/09 00:00:00  43     -527.0718     -645.2571
13: 2012/07/09 00:00:00  51       34.2777       53.0246
14: 2012/07/09 00:00:00  52     -118.9257     -130.6246
15: 2012/07/09 00:00:00  27       98.0804      133.2339
16: 2012/07/09 00:00:00  39       11.8226       19.8136

とりあえずまともそうな値が出たのでよしとする。多分合ってるだろう。
あとは単独測位をして仰角を横軸にしてノイズレベルを調べてみようかな。