CCTF analysis
Round xx: LabXX

Original author Bruce Warrington
Updated for Mathematica V8 by MJW

ClearAll;

Off[General::spell1];
Off[General::spell1];
Needs["ErrorBarPlots`"];

SetOptions[ListPlot,Frame->True,Axes->False,PlotRange->All,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];
SetOptions[ErrorListPlot,Frame->True,Axes->False,PlotRange->All,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];
desktop="/home/mjw/Desktop/tmp/";(* defines a path for any output *)
inch=72;
cm=inch/2.54;

To use this notebook,  you need to set the path to the data (dataPath) and the MJD range. Subdirectories for data are assumed to be cctf_host and cctf_data and CCTF files are assumed to be named MJD.cctf. Options on reading CCTF files which need to be set are whether or not the data is dual-frequency, and whether to compute the checksum (see the Initialise section) below.

Another fixup that may be required is indexing of the CCTF data fields as given in the Definitions below. It is assumed that the travelling receiver and host receiver are dual - frequency.  Indices need to change otherwise.

host="KRISS";
dataPath="/home/mjw/tmp/APMP GPS Cal/Round 2 processing/KRISS/";
{mjd1,mjd2}={53437,53443};
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[StringReplace[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];
If[data[[-1,-1]]==EndOfFile,
data=Drop[data,-1];
If[ShowErrors/.{opts}/.Options[ReadCCTF],Print["! Bad format for ",filename," at line ",count]];
Continue[];
];
If[CheckSum/.{opts}/.Options[ReadCCTF],
checksum=HexString[Mod[Plus@@Drop[ToCharacterCode[str],-2],256],2];
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,CheckSum->True};

Indices may need to change here for TrackLength2 etc (reduce by 3 if the travelling receiver is single frequency)

Select data columns

SV[list_]:=list[[1]];
TrackMJD[list_]:=list[[3]];
TrackTime[list_]:=list[[4]];
TrackLength[list_]:=list[[5]];
Elevation[list_]:=list[[6]]/10;
Azimuth[list_]:=list[[7]]/10;
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};

Allan Variance and Allan Deviation

Average[y_,k_,n_]:=(1/n)(∑_(i=k)^(k+n-1))y[[i]]

AllanVarianceNormal[y_,n_]:=Module[{M,m},
		M=Length[y];
		m=Floor[M/n]-1;
	Return[{m,(1/(2 m))Sum[(((Average[y,k+n,n]-Average[y,k,n]))^2),{k,1,n m,n}]}];
	]

AllanDeviationNormal[y_,n_]:=Block[{t},
	t=AllanVarianceNormal[y,n];
		Return[{t[[1]],(Sqrt[t[[2]]])}]
	]

AllanVarianceOverlapping[y_,n_]:=Module[{M},
		M=Length[y];
	Return[{M-2 n+1,(1/(2 (M-2 n+ 1)))(∑_(k=1)^(M-2 n+1))(((Average[y,k+n,n]-Average[y,k,n]))^2)}];
	]

AllanDeviationOverlapping[y_,n_]:=Block[{t},
	t=AllanVarianceOverlapping[y,n];
		Return[{t[[1]],(Sqrt[t[[2]]])}]
	]

Miscellaneous

TimeTag[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]+TimeTag[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)&];
tdata=Last/@tdata;
tmean=Mean[tdata];
tlen=Length[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

Remember to set DualFrequency and CheckSum

store={};
dataHost={};
For[mjd=mjd1,mjd≤mjd2,mjd++,dataHost=Join[dataHost,ReadCCTF[StringJoin["cctf_host/",ToString[mjd],".cctf"],DualFrequency->True,CheckSum->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]]

Separation of the two antennae

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

Workspace

Filter the data

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

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

Match tracks

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

Will need to fix the MDIO offset here

mdio = {DateValue[#1],ModIono[#1]-ModIono2[#1]-154.5}&/@dMerge;

Plot REF - SVN and MDIO difference

ListPlot[{diffdataGPS,mdio},FrameLabel->{"MJD",StringJoin[host," - APMP (ns)"]}]

Examine the modelled and measured ionosphere

{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]];
plot1=Show[GraphicsGrid[{
{ListPlot[{DateValue[#]-tmin,RefGPS[#]}&/@dataHost,title[StringJoin[host," (modelled ionosphere)"]],topts],
ListPlot[{DateValue[#]-tmin,RefGPS[#]+ModIono[#]-MeasIono[#]}&/@dataHost,title[StringJoin[host," (measured ionosphere)"]],topts]},
{ListPlot[{DateValue[#]-tmin,RefGPS[#]}&/@dataTrav,title["APMP (modelled ionosphere)"],topts],
ListPlot[{DateValue[#]-tmin,RefGPS[#]+ModIono[#]-MeasIono[#]}&/@dataTrav,title["APMP (measured ionosphere)"],topts]}
}]]

Examine the residuals for the modelled and measured ionosphere

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[GraphicsGrid[{
{ListPlot[{DateValue[#]-tmin,DSG[#]}&/@dataHost,title[host,"DSG"],topts],
ListPlot[{DateValue[#]-tmin,ISG[#]}&/@dataHost,title[host,"ISG"],topts]},
{ListPlot[{DateValue[#]-tmin,DSG[#]}&/@dataTrav,title["APMP","DSG"],topts],
ListPlot[{DateValue[#]-tmin,ISG[#]}&/@dataTrav,title["APMP","ISG"],topts]}
}]]

Daily means for REF - SVN

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

ErrorListPlot[{Take[#,2],ErrorBar[#[[3]]]}&/@daymeans,Frame->True,Axes->False]

Linear fits, unweighted and DSG weighted

fit={};
AppendTo[fit,LinearModelFit[diffdataGPS,x-53437,x,Weights->Automatic]];
AppendTo[fit,LinearModelFit[diffdataGPS,x-53437,x,Weights->(((DSG[#])^(-2))&/@dMerge)]];

tfit=fit[[1]]
MJDFirst=TrackMJD[dMerge[[1]]];
MJDLast=TrackMJD[dMerge[[-1]]];
MJDMiddle=(DateValue[dMerge[[1]]]+DateValue[dMerge[[-1]]])/2 ;
MeanOffset=tfit[MJDMiddle];
rms=Sqrt[tfit["EstimatedVariance"]];
SEslope=Part[tfit["ParameterConfidenceIntervalTableEntries"],2,2];

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

Print["Summary"];
fit2=fit[[2]]
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," [",fit2["BestFit"]/.{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[fit2["BestFitParameters"]],"] ps/day, with a standard error of ",SEslope*1000," [",1000 Part[fit2["ParameterConfidenceIntervalTableEntries"],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[fit2["BestFit"]/.{x->MJDMiddle},1],"ns"},
{"RMS",Places[rms,1],bar," ","ns"},
{"Slope",SlopeString[slope,SEslope],bar,SlopeString[Last[fit2["BestFitParameters"]],Part[fit2["ParameterConfidenceIntervalTableEntries"],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];

Data quality checks

Allan deviation of REF-SVN

y=(10^(-9))Last/@diffdataGPS;

adev=Table[Prepend[AllanDeviationNormal[y,(2^n)],60 13 (2^n)],{n,0,Floor[Log[2,Length[y]]]-2}];
TableForm[N[adev]]

odev=Table[Prepend[AllanDeviationOverlapping[y,(2^n)],60 13 (2^n)],{n,0,Floor[Log[2,Length[y]]]-1}];
TableForm[N[odev]]

red=Sequence[PlotStyle->{Hue[0],AbsoluteDashing[{6,4}]}];
Show[
ListLogLogPlot[Drop[#,{2}]&/@odev,PlotStyle->{Hue[0.3],AbsolutePointSize[8]},GridLines->Automatic],
ListLogLogPlot[Drop[#,{2}]&/@odev,PlotStyle->Hue[0.3],Joined->True],
ListLogLogPlot[Drop[#,{2}]&/@adev,PlotStyle->{Hue[0.6],AbsolutePointSize[8]},GridLines->Automatic],
ListLogLogPlot[Drop[#,{2}]&/@adev,PlotStyle->Hue[0.6],Joined->True],
LogLogPlot[adev[[1,3]]adev[[1,1]]/t,{t,adev[[1,1]],odev[[-1,1]]}],
LogLogPlot[adev[[1,3]]Sqrt[adev[[1,1]]/t],{t,adev[[1,1]],odev[[-1,1]]}],
Frame->True,Axes->False,DisplayFunction->Identity]

Check: elevation and azimuth plots, including residuals

style={
{GrayLevel[0.5],AbsolutePointSize[6]},
{GrayLevel[0],AbsolutePointSize[6]},
{Hue[0],AbsolutePointSize[6]},
{GrayLevel[0],AbsolutePointSize[3]}
};

Schedule and non - schedule tracks for the travelling receiver

Show[
ListPlot[{Azimuth[#],Elevation[#]}&/@Select[dataTrav,MatchQ[#[[2]],"FF"]&],PlotStyle->style[[1]]],
ListPlot[{Azimuth[#],Elevation[#]}&/@Select[dataTrav,!MatchQ[#[[2]],"FF"]&],PlotStyle->style[[3]]],
TextStyle->{FontSize->18},PlotRange->All
]

Tracks for travelling, host and merged data sets

Show[
ListPlot[{Azimuth[#],Elevation[#]}&/@dataTrav,PlotStyle->style[[1]]],
ListPlot[{Azimuth[#],Elevation[#]}&/@dataHost,PlotStyle->style[[3]]],ListPlot[{Azimuth[#],Elevation[#]}&/@dMerge,PlotStyle->style[[4]]],
TextStyle->{FontSize->18},PlotRange->All
]

Polar plot of above

ToPolar[list_]:=Module[{r,theta},
r=90-Elevation[list];
theta=(90-Azimuth[list]) π/180;
Return[{r Cos[theta],r Sin[theta]}]
];

grid=RGBColor[0, 0, 0.7];
Show[
ListPlot[ToPolar/@dataTrav,PlotStyle->style[[1]]],
ListPlot[ToPolar/@dataHost,PlotStyle->style[[3]]],
ListPlot[ToPolar/@dMerge,PlotStyle->style[[4]]],
Sequence@@(PolarPlot[90-#,{t,0,2 π},PlotStyle->grid]&/@{0,30,60}),
ListPlot[{{-90,0},{90,0}},Joined->True,PlotStyle->grid],
ListPlot[{{0,-90},{0,90}},Joined->True,PlotStyle->grid],
TextStyle->{FontSize->18},PlotRange->{{-90,90},{-90,90}},AspectRatio->1
]

REF-SVN residuals as a function of elevation

tfit=fit[[1]];
tplot[fun_,str_]:=ListPlot[
Table[{fun[dMerge[[i]]],tfit[diffdataGPS[[i,1]]]-diffdataGPS[[i,2]]},{i,Length[dMerge]}],
PlotStyle->AbsolutePointSize[3],GridLines->Automatic,FrameLabel->{str,"residual [ns]","",""},PlotRange->All,
DisplayFunction->Identity];
Show[tplot[Elevation,"elevation [deg]"],TextStyle->{FontSize->18}]

Show[tplot[Azimuth,"azimuth [deg]"],TextStyle->{FontSize->18}]