function [df, dt] = maser_adjust(mjd_range, pps_range, session_date, adjustment_date, polyfit_order, channel)
% function [df, dt] = maser_adjust(mjd_range, pps_range, session_date, adjustment_date, polyfit_order)
%
% Calculates the needed maser frequency adjustment to get nominal frequency at
% needed date. Takes data from mjd_range = [min max] and clips bad points 
% outside pps_range = [min max]. Returns the needed frequency offset (which is 
% *added* to the maser frequency) df and time offset (which is also *added* to 
% to the current time offset.
% by kaj.wiik@hut.fi (29 May 1996)
%
% Added new argument 'channel' and changed the HP GPS 
% to default by kwi (23 May 1997)

data = loadall(mjd_range); % Load the needed data

if nargin == 5
	channel = 3;
end

pps = data(:,channel);
date = data(:,1);
clear data

session_date = session_date - 50000;
adjustment_date = adjustment_date - 50000;

% Cut dates and bad points
ind = find((date > min(mjd_range))&(date < max(mjd_range))...
            &(pps > min(pps_range))&(pps < max(pps_range)));

pps = pps(ind);
date = date(ind) - 50000; % Cut and normalize to MJD 50000

p = polyfit(date, pps, polyfit_order)

rms = sqrt(sum((pps - polyval(p,date)).^2))/length(pps);

figure(1)
plot(date, pps, 0:max(mjd_range-50000), polyval(p, 0:max(mjd_range-50000)));
title(['Polynomial fit result, order ' sprintf('%d', polyfit_order) ', RMS error ' sprintf('%.1e', rms) ]);
xlabel('Time since MJD = 50000');
grid

figure(2)
plot(date, pps, date, polyval(p, date));
title('Result of the adjustment and prediction');
xlabel('Time since MJD = 50000');
grid
hold on

df0 = p(length(p) - 1); % The current frequency offset is the second least
                        % order coefficient

pd = polyder(p); % Polynomial pd is frequency=pd(time)

% Solve for needed frequency offset to get zero at session date
pd(length(pd)) = 0;
dfs = polyval(-pd, session_date );

% Frequency offset at the hydrogen line frequency [Hz]
df = (df0 - dfs)/(1e6*24*60*60)*1420405751;

% Offset that must be *added* to the current offset in microseconds 
% to get 10 us at session date
pn = p;
pn(length(pn) - 1) = dfs;

% PPS phase without offset correction
offset_now = polyval(p, adjustment_date);

pn(length(pn)) = 0;
pn(length(pn)) = offset_now - polyval(pn, adjustment_date);

plotdates = adjustment_date:(session_date + 30);

% PPS phase without offset correction
plot(plotdates, polyval(pn, plotdates),'r')

% Offset correction
dt = 10 - polyval(pn, session_date);

% PPS phase with offset correction
plot(plotdates, polyval(pn, plotdates) + dt,'g')


ax = axis;
line([session_date session_date], [ax(3) ax(4)])
text(adjustment_date, offset_now, 'A')
hold off

%keyboard

