%  TSplotKMS is a simple driver program for time series plotting.
%  It prompts the user for a data file, in one of several formats,
%   reads this file, and initializes the time series plotting interface.
%
%  First a couple of options which are hard-coded here:
decFac = 1; 

%  return full data set (no decimation);scale into 
scale = 1; %if scale == 1 data are returned with units of nT and mV
           %   note: not divided by dipole length!
%  physical units (but no other filter corrections applied)

tUnits = 'Seconds';
%   tUnits = initial time units for labeling time series plots
%     other choices are 'Seconds','Minutes','Days'

nTunits = 6;
%   Number of Hours (more generally, tUnits) to plot on one pag3e

%  this program just reads one file, obtained by browsing
%   load one NIMS data  ... 
[filename, pathname, filterindex] =     ...
        uigetfile('*','Data File to Plot');

[header,y] = readKMS([pathname filename]);

%    get sampling rate, start time etc. from header
dt = 1./header.sampling_frequency;
%   initialize structure "plotStart"
[plotStart] = plotSetup(tUnits,nTunits,dt);
[nch,nData] = size(y);

%get firmware version
fw = dec2bin(header.firmware_version,8);
fw = (fw*1) - 48;
firmware_version = fw(1)*2^3+fw(2)*2^2+fw(3)*2+fw(4);
fw_dec = fw(5)*2^3+fw(6)*2^2+fw(7)*2+fw(8);
if (fw_dec > 9)
    firmware_version = firmware_version + 0.01*fw_dec;
else
    firmware_version = firmware_version + 0.1*fw_dec;
end
firmware_version 
fpga_version = header.fpga_version
%get gain settings
s_char=dec2bin(header.ch_setting,8);
s_bit=s_char(8:-1:1,8:-1:1);
s_rev=reshape(s_bit',1,64);
s_rev;
for i=1:6
    ii=1+(i-1)*2;
    switch (s_rev(ii:ii+1));
        case '00'
            ch_gain1(i)=1;
        case '10'  %reversed bit order
            ch_gain1(i)=1;
        case '01' %reversed bit order
            ch_gain1(i)=8;
        case '11'
            if (fpga_version >= 64)
                ch_gain1(i)=152;
            else
                ch_gain1(i)=40;
            end
        otherwise
            error('This is impossible')
    end
    ii=13+(i-1)*2;
    switch (s_rev(ii:ii+1));
        case '00'
            ch_gain2(i)=1;
        case '10' %reversed bit order
            ch_gain2(i)=4;
        case '01' %reversed bit order
            ch_gain2(i)=16;
        case '11'
            ch_gain2(i)=64;
        otherwise
            error('This is impossible')
    end
    ch_gain(i)=ch_gain1(i)*ch_gain2(i);
    
end
ch_gain1;
ch_gain2;
ch_gain
scale = ch_gain.*1000;  % correct gain & convert to V
if (fpga_version >= 64)
    scale = scale.*0.4
end
scale;


%   define 2 character strings for labeling channels
chType = ['?','B','E','U','B','B','B','E'];
ch = [];
for k = 1:nch
    ch = [ch; chType(header.sensor_ID{k}+1) lower(header.sensor_component{k})];
end

%   sort out time-date
temp =  strsplit(header.GPS_time_info,',');
temp_time=temp;
hr = str2num(temp{2}(1:2));
mm = str2num(temp{2}(3:4));
ss = str2num(temp{2}(5:10));
ss = floor(ss);   % John is using only integer part
day = str2num(temp{3});
mo = str2num(temp{4});
year = str2num(temp{5});
dataStart = datestr([year mo day hr mm ss]);  % prelimiary start of first file
%   by default let's set clock reference to start of year that data are
%   collected in
dataZero = datestr([year 1 1 0 0 0]);
%   t0 gives sample number of first data point,  referenced to clock zero
%     (i.e., sample 1 is at stroke of clock zero
t0_days = datenum(dataStart)-datenum(dataZero);
first_sample_time = header.first_sample_acq_time;
clk_syn = str2num(header.GPS_syn_clk);
first_sample_time = first_sample_time*1.e-6/clk_syn + 1;
if (header.sampling_frequency == 1000) % delay of A/D
    sigmaADdelay = 38/header.sampling_frequency;
elseif  (header.sampling_frequency == 4000)
    sigmaADdelay = 38/header.sampling_frequency;
else
    sigmaADdelay = 39/header.sampling_frequency;
end
delay_ss = first_sample_time - sigmaADdelay;
delay_days = delay_ss/86400;
t0_days = t0_days + delay_ss/86400;
dataStart = datestr(t0_days + datenum(dataZero)) % real start of this file
delay_ss_dec = delay_ss - floor(delay_ss);

t0 = t0_days*86400;   % seconds from dataZero

temp = strsplit(filename,'.');
sta = [];
for k = 1:nch
    sta =  [sta; temp{1}];
end

%   set up dummy unit system response ... this can use info in header
freq = [1,1e-6]./dt;
[SysTF] = setUnitSysTF(nch,freq);

%  reset nPw (# of points to plot initially on one page) now
%   that dt has been read from data file  (NOTE: dt returned
%   by readTSnims accounts for subsampling by a factor of decFac)
nPw = fix(plotStart.tFac*plotStart.nTunits/dt);
nPw = min(nPw,nData);
plotStart.nPw = nPw;

%   Initialize cell arrays for data, time series metadata, 
%     window metadata, etc.
TSc = {};
TSp = {};
TSd = {};
TSw = {};
TSdata = {};
TSfigInd = zeros(100,1);

%   set up TSd
TSd1 = setTSd(y,t0,dt,ch,sta,SysTF,scale,dataZero);
TSd1.sta = char(TSd1.sta);
TSdata0 = (TSd1.range(:,1)+TSd1.range(:,2))/2;
rng = (TSd1.range(:,2)-TSd1.range(:,1))/2;
rng = plotStart.rangeScale*rng;
TSd1.range(:,1:2) = [ TSdata0-rng TSdata0+rng];

%   are these needed???
TSd1.DataDir = pathname;
TSd1.arrayID = '';
TSd1.dataStart = dataStart;

%   copy into cell arrays
iA = 1;
TSdata{iA} = y;
TSd{iA} = TSd1;

[hTSw,TSw1] = plotTSwin(TSd{1},plotStart);
TSfigInd(hTSw) = 1;
i0 = TSw1.i0; i1 = TSw1.i1; di = TSw1.di;
yFac = TSw1.yFac;
ddt = TSd{1}.dt/TSw1.tFac;
t0 = (TSd{1}.t0)*(1-TSw1.t0Eq0)/TSw1.tFac;
I = [i0:di:i1];
[TSw1] = plotPage(ddt*I+t0,...
                yFac*TSdata{1}(:,I),...
                TSw1,TSd{1},1);
TSw{1} = TSw1;
