% ------------------------------------------------------------------------
% This MATLAB script produces plots for the time offset between two 
% receivers, using the GPS common-view technique. 
% Use this script in conjunction with CCTF.m and MergedCCTF.m
% -> CCTF.m allows the user to load .cctf files and assigns the relevant
% parameters in MATLAB. A filter function (FilterTracks) is also available 
% to filter out bad tracks. 
% -> MergedCCTF.m allows the user to merge the data between 2 receivers,
% and calculates the time offset.
% This script was written by Magnus Hsu (NMI-A) 20 Nov 2012.
% ------------------------------------------------------------------------

% clears MATLAB workspace and memory 
clear all; clc;

% define start and stop MJDs for data analysis
startMJD = 53437;
stopMJD = 53443;

% loads data
sys1 = CCTF(startMJD,stopMJD,'cctf_host\','.cctf');   % for system 1 (e.g. host)
sys2 = CCTF(startMJD,stopMJD,'cctf_data\','.cctf');   % for system 2 (e.g. remote)
 
% filters data to remove bad tracks (i.e. DSG < 999 and TrackLength = 780) 
FilteredSys1 = FilterTracks(sys1,999,780);
FilteredSys2 = FilterTracks(sys2,999,780);

% subtracts the two data sets to yield the time difference
% (i.e. Sys1 - Sys2)
combsys = MergedCCTF(FilteredSys1,FilteredSys2);
combsys.Merge();

% extract time difference data without Mod Iono
refsv = combsys.RefSVNoMDIO(0);

% since data is given in 0.1 ns units, CONVERT data to 1 ns by dividing by
% a factor of 10. 
refsv(:,2) = refsv(:,2)/10;

%% plot time difference data
% note that refsv(:,1) is the cumulative MJD and refsv(:,2) is the time
% difference data between sys1 and sys2.
figure(1)
clf
plot(refsv(:,1),refsv(:,2),'.')
ylabel('Sys1 - Sys2 (ns)');
xlabel('MJD');

%% calculates and plots RefGPS, for instances With and WithOut Iono data
% for System 1 
rx1ModIono = combsys.rx1.Tracks(:,CCTF.REFGPS);
rx1MeasIono = combsys.rx1.Tracks(:,CCTF.REFGPS) + ...
                combsys.rx1.Tracks(:,CCTF.MDIO) - ...
                combsys.rx1.Tracks(:,CCTF.MSIO);
% for System 2
rx2ModIono = combsys.rx2.Tracks(:,CCTF.REFGPS);
rx2MeasIono = combsys.rx2.Tracks(:,CCTF.REFGPS) + ...
                combsys.rx2.Tracks(:,CCTF.MDIO) - ...
                combsys.rx2.Tracks(:,CCTF.MSIO);           

% plot modelled and measured Iono data
figure(2)
clf
% Modelled Iono data for System 1
subplot(2,2,1)
plot(refsv(:,1),rx1ModIono,'.');
title('System 1 (mod iono)');
ylabel('Ref - GPS (ns)');
xlabel('MJD');
% Measured Iono data for System 1
subplot(2,2,2)
plot(refsv(:,1),rx1MeasIono,'.');
title('System 1 (meas iono)');
ylabel('Ref - GPS (ns)');
xlabel('MJD');
% Modelled Iono data for System 2
subplot(2,2,3)
plot(refsv(:,1),rx2ModIono,'.');
title('System 2 (mod iono)');
ylabel('Ref - GPS (ns)');
xlabel('MJD');
% Measured Iono data for System 2
subplot(2,2,4)
plot(refsv(:,1),rx2MeasIono,'.');
title('System 2 (meas iono)');
ylabel('Ref - GPS (ns)');
xlabel('MJD');

%% plot DSG and ISG data
figure(3)
clf
% DSG data for System 1
subplot(2,2,1)
plot(refsv(:,1),combsys.rx1.Tracks(:,CCTF.DSG),'.');
title('System 1 (DSG)');
ylabel('DSG');
xlabel('MJD');
% ISG data for System 1
subplot(2,2,2)
plot(refsv(:,1),combsys.rx1.Tracks(:,CCTF.ISG),'.');
title('System 1 (ISG)');
ylabel('ISG');
xlabel('MJD');
% DSG data for System 2
subplot(2,2,3)
plot(refsv(:,1),combsys.rx2.Tracks(:,CCTF.DSG),'.');
title('System 2 (DSG)');
ylabel('DSG');
xlabel('MJD');
% ISG data for System 2
subplot(2,2,4)
plot(refsv(:,1),combsys.rx2.Tracks(:,CCTF.ISG),'.');
title('System 2 (ISG)');
ylabel('ISG');
xlabel('MJD');

%% calculate means on a day-to-day basis

for ii = 1:(stopMJD - startMJD + 1)
    indx = find(refsv(:,1) > (startMJD + ii - 1) & ....
        refsv(:,1) <= (startMJD + ii));
    % MJD(i):
    day(ii) = fix(refsv(indx(end),1));
    % Mean for MJD(i):
    daymean(ii) = mean(refsv(indx(1):indx(end),2));
    % Standard deviation for MJD(i):
    daystd(ii) = std(refsv(indx(1):indx(end),2));
    % Number of tracks for MJD(i):
    daylength(ii) = length(indx);
end

% plot daily means with errorbars
figure(4)
clf
errorbar(day,daymean,daystd,'.')
xlabel('MJD');
ylabel('Mean of Sys1 - Sys2 (ns)');
    
%% Do linear fit to time difference data between sys1 and sys2
% use polynomial fit of 1st order (since linear)
% however, subtract start MJD since want to get offset directly
P = polyfit(refsv(:,1)-startMJD,refsv(:,2),1);

% define x-vector of equally spaced MJD 
Xfit = (refsv(1,1):(refsv(end,1)-refsv(1,1))/1000:refsv(end,1))-startMJD; 
% resulting y-vector from fit
Yfit = polyval(P,Xfit);

% plot time difference and linear fit
figure(5)
clf
plot(refsv(:,1)-startMJD,refsv(:,2),'.',Xfit,Yfit,'r-')
xlabel('MJD')
ylabel('Sys1 - Sys2 (ns)')
title('Time difference plot of Sys1 vs Sys2')

% display fitted function
disp(['----------------------------------------------------------------']);
disp(['Linear fit is given by:']);
disp(['y = ' num2str(P(1)) '( -' num2str(startMJD) ' + x) + (' num2str(P(2)) ')']); 
disp(['----------------------------------------------------------------']);
disp([num2str(length(refsv(:,2))) ' common-view tracks were analysed between MJD ' num2str(day(1)) 'and MJD ' num2str(day(end)) ]);



