APMP Round-robin GPS Campaign 2004
Round 1: NMI Sydney (3)

Off[General::spell1];
Off[General::spell1];
Needs["Graphics`Graphics`"];
Needs["Graphics`MultipleListPlot`"];
Needs["Statistics`LinearRegression`"];

SetOptions[ListPlot, Frame -> True, Axes -> False, 
        PlotRange -> All, 
        TextStyle -> {FontFamily -> "Helvetica", 
            FontSize -> 6}];
desktop = "c:\\\\documents and settings\\rw1899\\desktop\\";
inch = 72;
cm = inch/2.54;

host = "NMI3";
dataPath = "G:\\time\\share\\GPS Cal\\APMP GPS \
Cal\\ReprocessingSep04\\NMI 3\\";
{mjd1, mjd2} = {53180, 53208};
store = {};

Definitions

Read data file (single or dual frequency)

HexString[num_, minlen_] := 
        StringJoin[Join[CharacterRange["0", "9"], 
              CharacterRange["A", "F"]][[(# + 
                    1) & /@ 
              IntegerDigits[num, 16, 
                minlen]]]];

\(\(StoreExtra[str_] := 
        Module[{instr, tag, value, pos}, 
instr = 
            StringToStream[str]; 
tag = 
            Read[instr, Word]; 
pos = 
            Position[{"X", "Y", "Z", "INT", "REF", "CAB"}, tag]; 
If[MatchQ[pos, {}], 
            Return[]]; 
If[
            pos[[1, 1]] > 
              3, 
value = 
              Read[instr, 
                Word]; 
If[\(! StringMatchQ[
                  value, "DLY"]\), Return[]]; 
tag = 
              StringJoin[tag, " ", value]]; 
Read[
            instr, Word]; 
value = 
            Read[instr, Number]; 
Close[
            instr]; 
pos = 
            Position[store, tag]; 
If[MatchQ[pos, {}], 
            AppendTo[store, {tag, {}}]; \ pos = Length[store], 
            pos = pos[[1, 
                1]]]; 
store[[pos, \(-1\)]] = 
            Union[store[[
                pos, \(-1\)]], {value}];
];\)\)

\(\(ReadCCTF[name_, opts___] := 
        Module[{filename, infile, str, checksum, count = 0, fail = 0, 
            flag = False, data = {}, format}, 
format = 
            Table[Number, {i, 
                If[\(DualFrequency /. {opts}\) /. Options[ReadCCTF], 21, 
                  18]}]; 
filename = 
            If[NumberQ[name], StringJoin[ToString[name], ".cctf"], 
              name]; 
If[\(UseDataPath /. {opts}\) /. 
              Options[ReadCCTF], 
            filename = 
              StringJoin[dataPath, filename]]; 
infile = 
            OpenRead[filename]; 
          format[[2]] = \(format[[\(-1\)]] = 
              Word\); 
While[
            True, 
If[\((str = 
                    Read[infile, String])\) == EndOfFile, 
              Break[]]; 
\(count++\); 
If[
              StringPosition[str, " = "] ≠ {}, 
              StoreExtra[str]]; 
If[
              StringPosition[str, "hhmmss"] ≠ {}, 
              flag = True; \ 
              Continue[]]; 
If[\((\(Header /. {opts}\) /. 
                    Options[ReadCCTF])\)\  && \ \(! flag\), 
              Continue[]]; 
instr = 
              StringToStream[str]; 
AppendTo[data, 
              Flatten[ReadList[instr, format]]]; 
Close[
              instr]; 
checksum = 
              HexString[Mod[Plus @@ Drop[ToCharacterCode[str], \(-2\)], 256], 
                2]; 
If[
              data[[\(-1\), \(-1\)]] == EndOfFile, 
data = 
                Drop[data, \(-1\)]; 
If[\(ShowErrors /. {opts}\) /. Options[ReadCCTF], 
                Print["! Bad format for ", filename, " at line ", 
                  count]]; 
Continue[];
]; 
If[\(! StringMatchQ[data[[\(-1\), \(-1\)]], 
                  checksum]\), 
 (*\ \(Print["! ", 
                    data[[\(-1\), \(-1\)]], " does not match ", checksum, " at line ", 
                    count];\)\ *) 
data = 
                Drop[data, \(-1\)]; 
\(fail++\);
];
]; 
Close[
            infile]; 
If[\(Verbose /. {opts}\) /. 
              Options[ReadCCTF], 
Print["> Read ", 
              Length[data], " tracks from ", 
              filename]]; 
If[\((fail > 
                  0)\)\  && \ \((\(ShowErrors /. {opts}\) /. 
                  Options[ReadCCTF])\), 
            Print["! ", 
              fail, " tracks did not match checksum"]]; 
Return[
            data];
];\)\)
Options[ReadCCTF] := {Header -> True, Verbose -> True, 
          ShowErrors -> True, UseDataPath -> True, 
          DualFrequency -> False};

Select data columns

SV[list_] := 
        list[[1]];
TrackMJD[list_] := 
        list[[3]];
TrackTime[list_] := 
        list[[4]];
TrackLength[list_] := 
        list[[5]];
Elevation[list_] := 
        list[[6]];
Azimuth[list_] := 
        list[[7]];
RefSV[list_] := 
        list[[8]]/
          10;
SlopeSV[list_] := 
        list[[9]];
RefGPS[list_] := 
        list[[10]]/
          10;
SlopeGPS[list_] := 
        list[[11]];
DSG[list_] := 
        list[[12]];
ModIono[list_] := 
        list[[16]]/
          10;
MeasIono[list_] := 
        list[[18]]/
          10;
SMSIO[list_] := 
        list[[19]];
ISG[list_] := 
        list[[20]];

    
TrackLength2[list_] := 
        list[[5 + 
            21]];
RefSV2[list_] := 
        list[[8 + 21]]/
          10;
RefGPS2[list_] := 
        list[[10 + 21]]/
          10;
ModIono2[list_] := 
        list[[16 + 21]]/
          10;

Sort data (into time order then SV order)

OrderCCTF[a_, b_] := 
        If[a[[3]] == 
            b[[3]], 
          If[a[[4]] == 
              b[[4]], 
            First[a] < First[b], 
            a[[4]] < 
              b[[4]]], 
          a[[3]] < 
            b[[3]]];

SortCCTF[data_] := Sort[data, OrderCCTF[#1, #2] &];

Merge two CCTF data sets on matching tracks

\(\(MergeCCTF[data1_, data2_, opts___] := 
        Module[{d1, d2, i, j, data = {}}, 
d1 = 
            SortCCTF[data1]; 
d2 = 
            SortCCTF[data2]; 
j = 1; 
For[
            i = 1, i ≤ 
              Length[d1], \(i++\), 
While[\((j <= 
                    Length[d2])\)\  && \ 
\((\((d2[[j, 3]] < 
                        d1[[i, 
                          3]])\)\  || \ 
\((\((d2[[j, 
                              3]] == 
                            d1[[i, 
                              3]])\)\  && \ 
\((\((d2[[j, 4]] < 
                                d1[[i, 
                                  4]])\)\  || \ 
\((\((d2[[j, 
                                      4]] == 
                                    d1[[i, 
                                      4]])\)\  && \ 
\((d2[[j, 1]] < 
                                    d1[[i, 
                                      1]])\))\))\))\))\), 
\(j++\)]; 
If[\((j ≤ 
                    Length[d2])\)\  && 
\((d2[[j, 3]] == 
                    d1[[i, 
                      3]])\) && \((d2[[
                      j, 4]] == 
                    d1[[i, 
                      4]])\)\  && \ \((d2[[j, 1]] == 
                    d1[[i, 
                      1]])\), 
AppendTo[
                data, Flatten[{d1[[i]], 
                    d2[[
                      j]]}]]];
]; 
If[\(Verbose /. {opts}\) /. 
              Options[MergeCCTF], 
Print["> First ", 
              Length[d1], " tracks, second ", 
              Length[d2], " tracks, matching ", 
              Length[data], " tracks"]]; 
Return[
            data];
];\)\)
Options[MergeCCTF] := {Verbose -> True};

Miscellaneous

TimeValue[list_] := 
        Module[{val, h, m, s}, 
val = 
            TrackTime[list]; 
h = 
            Quotient[val, 10000]; 
m = 
            Quotient[Mod[val, 10000], 100]; 
s = 
            Mod[val, 100]; 
Return[
            3600 h  +  60 m  +  s];
];

DateValue[list_] := 
        TrackMJD[list] + TimeValue[list]/86400;

Filter[list_, fn_] := 
        Module[{result}, 
result = 
            Select[list, fn]; 
Print["> ", 
            Length[result], " of ", 
            Length[list], " tracks selected."]; 
Return[
            result];
];

MeanByDay[list_] := 
        Module[{day, tlist, tmean, trms, tlen, 
            result = {}}, 
For[
            day = Floor[list[[1, 1]]], 
            day ≤ 
              Floor[list[[-1, 
                  1]]], day++, 
tdata = 
              Select[list, 
                day ≤ #[[1]] < (day + 1) &]; 
tlen = 
              Length[tdata]; 
If[tlen == 0, 
              Continue[]]; 
tdata = 
              Last /@ tdata; 
tmean = 
              Mean[tdata]; 
trms = 
              Sqrt[Sum[(((tdata[[i]] - tmean))^2), {i, tlen}]/tlen]; 
AppendTo[
              result, {day, tmean, trms, tlen}]]; 
Return[
            result];
];

Tweak[graph_] := 
        Module[{fix = graph, pos, grid, }, 
pos = 
            Position[fix, 
              GridLines]; 
pos[[1, -1]]++; 
grid = 
            AbsoluteOptions[fix, 
              GridLines]; 
fix[[
              Sequence @@ 
                pos[[1]]]] = 
            Replace[grid[[1, -1]], {RGBColor[__] -> RGBColor[0.705882,  1,  1], 
                AbsoluteThickness[__] -> AbsoluteThickness[0.5]}, 
              Infinity];  
Return[
            fix];
];

Initialise

store = {};
dataHost = {};
For[mjd = mjd1, mjd ≤ mjd2, mjd++, 
        dataHost = 
          Join[dataHost, 
            ReadCCTF[
              StringJoin["cctf_host\\", ToString[mjd], ".cctf"], 
              DualFrequency -> True, 
              Verbose -> True]]];
store1 = store

TableForm[{#[[1]], 
            PaddedForm[#[[2, 1]], {10, 
                3}, ExponentFunction -> (Null  &)]} & /@ 
        Take[store1, 3]]

store = {};
dataTrav = {};
For[mjd = mjd1, mjd ≤ mjd2, mjd++, 
        dataTrav = 
          Join[dataTrav, 
            ReadCCTF[
              StringJoin["cctf_data\\", ToString[mjd], ".cctf"], 
              DualFrequency -> True, 
              Verbose -> True]]];
store2 = store

TableForm[{#[[1]], 
            PaddedForm[#[[2, 1]], {10, 
                3}, ExponentFunction -> (Null  &)]} & /@ 
        Take[store2, 3]]

Sqrt[
      Plus @@ (#[[2, 1]]^2 & /@ 
            Take[store1 - store2, 3])]

Workspace

dataHost = 
        Filter[dataHost, (TrackLength[#] == 780) && (DSG[#] < 
                  999) && (MeasIono[#] < 999) && (Abs[SMSIO[#]] < 
                  999) && (ISG[#] < 999) &];

dataTrav = 
        Filter[dataTrav, (TrackLength[#] == 780) && (DSG[#] < 
                  999) && (MeasIono[#] < 999) && (Abs[SMSIO[#]] < 
                  999) && (ISG[#] < 999) &];

dMerge = MergeCCTF[dataHost, dataTrav];
diffdataGPS = {DateValue[#1], 
              RefSV[#1] + ModIono[#1] - RefSV2[#1] - ModIono2[#1]} & /@ 
          dMerge;

{xmin, 
          xmax} = DateValue[
              dMerge[[#]]] & /@ {1, -1};
tmin = Floor[xmin];
title[str_] := 
        Sequence[
          FrameLabel -> {StringJoin["MJD - ", 
                ToString[tmin]], "Ref-GPS (ns)", 
              str, ""}];
topts = 
        Sequence[PlotRange -> {{xmin, xmax} - tmin, All}, 
          PlotStyle -> AbsolutePointSize[1], 
          GridLines -> Automatic, 
          FrameStyle -> AbsoluteThickness[0.5], 
          DisplayFunction -> Identity];
plot1 = 
        Show[GraphicsArray[{
{Tweak[
                  ListPlot[{DateValue[#] - tmin, RefGPS[#]} & /@ 
                      dataHost, 
                    title[StringJoin[host, " (modelled ionosphere)"]], 
                    topts]], 
Tweak[
                  ListPlot[{DateValue[#] - tmin, 
                          RefGPS[#] + ModIono[#] - MeasIono[#]} & /@ 
                      dataHost, 
                    title[StringJoin[host, " (measured ionosphere)"]], 
                    topts]]}, 
{Tweak[
                  ListPlot[{DateValue[#] - tmin, RefGPS[#]} & /@ 
                      dataTrav, title["APMP (modelled ionosphere)"], 
                    topts]], 
Tweak[
                  ListPlot[{DateValue[#] - tmin, 
                          RefGPS[#] + ModIono[#] - MeasIono[#]} & /@ 
                      dataTrav, title["APMP (measured ionosphere)"], 
                    topts]]}
}]];

title[str1_, str2_] := 
        Sequence[
          FrameLabel -> {StringJoin["MJD - ", ToString[tmin]], str2, 
              StringJoin[str1, " ", 
                str2], ""}];
topts = 
        Sequence[PlotRange -> {{xmin, xmax} - tmin, {0, 160}}, 
          PlotStyle -> AbsolutePointSize[1], 
          GridLines -> Automatic, 
          FrameStyle -> AbsoluteThickness[0.5], 
          DisplayFunction -> Identity];
plot2 = 
        Show[GraphicsArray[{
{Tweak[
                  ListPlot[{DateValue[#] - tmin, DSG[#]} & /@ dataHost, 
                    title[host, "DSG"], topts]], 
Tweak[
                  ListPlot[{DateValue[#] - tmin, ISG[#]} & /@ dataHost, 
                    title[host, "ISG"], 
                    topts]]}, 
{Tweak[
                  ListPlot[{DateValue[#] - tmin, DSG[#]} & /@ dataTrav, 
                    title["APMP", "DSG"], 
                    topts]], 
Tweak[
                  ListPlot[{DateValue[#] - tmin, ISG[#]} & /@ dataTrav, 
                    title["APMP", "ISG"], 
                    topts]]}
}]];

daymeans = MeanByDay[diffdataGPS];
TableForm[daymeans // N]

{Plus @@ (Last /@ daymeans), Length[diffdataGPS]}

MultipleListPlot[{Take[#, 2], 
              ErrorBar[#[[3]]]} & /@ 
          daymeans, SymbolStyle -> GrayLevel[0], 
        SymbolShape -> PlotSymbol[Box, 4], 
        TextStyle -> {FontSize -> 12}];

fit = {};
AppendTo[fit, 
        Regress[diffdataGPS, {1, x}, x, Weights -> Automatic, 
          RegressionReport -> {BestFit, BestFitParameters, 
              ParameterCITable, EstimatedVariance}]];
AppendTo[fit, 
        Regress[diffdataGPS, {1, x}, x, 
          Weights -> (DSG[#]^(-2) & /@ dMerge), 
          RegressionReport -> {BestFit, BestFitParameters, 
              ParameterCITable, EstimatedVariance}]];

tfit = 
        fit[[1]];
MJDFirst = 
        TrackMJD[
          dMerge[[1]]];
MJDLast = 
        TrackMJD[
          dMerge[[-1]]];
MJDMiddle = (DateValue[
                dMerge[[1]]] + 
              DateValue[
                dMerge[[-1]]])/
          2;
{intercept, slope} = 
        BestFitParameters /. tfit;
MeanOffset = 
        slope MJDMiddle + intercept;
rms = Sqrt[EstimatedVariance /. tfit];
SEslope = Part[ParameterCITable /. tfit, 1, 2, 2];

popts[hue_] := 
        Sequence[PlotJoined -> True, 
          PlotStyle -> {AbsoluteThickness[0.5], 
              Hue[hue]}];
\(\(plot3 = 
        Show[Tweak[
            DisplayTogether[
ListPlot[\(# - {tmin, 0} &\) /@ 
                  diffdataGPS, 
                PlotStyle -> 
                  AbsolutePointSize[
                    1]], 
ListPlot[\({# - 
                        tmin, \(BestFit /. 
                          tfit\) /. {x -> #}} &\) /@ {xmin, xmax}, 
                popts[0]], 
ListPlot[\({MJDMiddle - tmin, 
                      MeanOffset + #} &\) /@ {\(-rms\), rms}, 
                popts[0]], 
 (*\ \(ListPlot[\({# - 
                          tmin, \(BestFit /. 
                            fit[[2]]\) /. {x -> #}} &\) /@ {tmin, tmax}, 
                  popts[0.6]]\)\(,\)\ *) 
PlotRange -> {{xmin, xmax} - tmin, All}, GridLines\  -> Automatic, 
              Frame -> True, FrameStyle -> AbsoluteThickness[0.5], 
              FrameLabel -> {StringJoin["MJD - ", ToString[tmin]], 
                  StringJoin[host, " - APMP (ns)"]}, 
              DisplayFunction -> Identity
]], 
          DisplayFunction -> $DisplayFunction];\)\)

Print["Summary"];
Print[
      Length[dMerge], " common-view tracks were analysed between MJD ",  
      MJDFirst, " and MJD ", MJDLast]; 
    Print["The mean offset (Host Rx - Travelling Rx) between the two \
receivers was ", 
      MeanOffset, " [", BestFit /. 
          fit[[2]] /. {x -> 
            MJDMiddle}, " weighted] ns, with an RMS deviation of ", 
      rms, " ns, at ", N[MJDMiddle], "."];
Print["The slope of the line of best fit was ", 
        slope*1000, " [", 
        1000 Last[
            BestFitParameters /. 
              fit[[2]]], "] ps/day, \
with a standard error of ", SEslope*1000, " [", 
        1000 Part[
            ParameterCITable /. 
              fit[[2]], 1, 2, 
            2], "] ps/day."];

Output

Export[StringJoin[desktop, StringJoin[host, "1.gif"]], plot1, 
        ImageSize -> 15 cm, ImageResolution -> 300];
Export[StringJoin[desktop, StringJoin[host, "2.gif"]], plot2, 
        ImageSize -> 15 cm, 
        ImageResolution -> 300];
Export[StringJoin[desktop, StringJoin[host, "3.gif"]], plot3, 
        ImageSize -> 7 cm, ImageResolution -> 300];

\(\(head = {
{"Summary", 
            FontWeight -> "Bold"}, 
StringJoin[
            ToString[Length[dMerge]], " common-view tracks, MJD ", \ 
            ToString[MJDFirst], " - ", 
            ToString[
              MJDLast]], 
StringJoin["Mean offset at MJD ", ToString[
              N[MJDMiddle]]], 
" "
};\)\)

Places[number_, places_] := 
        NumberForm[number, {100, places}];
SlopeString[slope_, error_] := 
        StringJoin[ToString[Round[1000 slope]], "±", 
          ToString[Round[1000 error]]];

bar = "/";
\(\(info = {
{" ", "(unweighted", 
            bar, "weighted)", " "}, 
{"Mean offset", Places[MeanOffset, 1], bar, 
            Places[\(BestFit /. 
                  fit[[2]]\) /. {x -> MJDMiddle}, 1], "ns"}, 
{"RMS", 
            Places[rms, 1], 
            bar, " ", "ns"}, 
{"Slope", 
            SlopeString[slope, SEslope], bar, 
            SlopeString[
              Last[BestFitParameters /. 
                  fit[[2]]], 
              Part[ParameterCITable /. 
                  fit[[2]], 1, 2, 
                2]], "ps/day"}
};\)\)

Scribe[arg_, opts___] := 
        If[MatchQ[Head[arg], List], 
Scribe[#, opts] & /@ 
            arg, 
StyleForm[
            arg, 
FontFamily -> (FontFamily /. {opts} /. Options[
                    Scribe]), 
FontSize -> (FontSize /. {opts} /. 
                  Options[
                    Scribe]), 
FontWeight -> (FontWeight /. {opts} /. 
                  Options[
                    Scribe])]
];
Options[Scribe] = {FontFamily -> "Helvetica", 
          FontSize -> 14, FontWeight -> "Plain"};

top = 
        GridBox[{If[MatchQ[Head[#], List], 
                  Scribe[#[[1]], 
                    Sequence @@ Rest[#]], Scribe[#]]} & /@ head, 
          ColumnAlignments -> Left];
body = 
        GridBox[Scribe[info], 
          ColumnAlignments -> {Left, Right, Center, Left, 
              Left}];
table = 
      DisplayForm[
        GridBox[{{top}, {body}}, ColumnAlignments -> Left]]

Export[StringJoin[desktop, StringJoin[host, "4.gif"]], 
        ToBoxes[table], ImageSize -> 10 cm, 
        ImageResolution -> 300];

Export[StringJoin[desktop, StringJoin[host, "mean.dat"]], 
        daymeans // N];