function polarExB( State )
% Entry point function for Kalman filter inversion of the polar cap
%
% Uses ExB drift from weimer model for the forward projection
%
% Optionally allows prior solution to be included in projection. Prior must be
% half resolution in latitude and longitude.
%
% Requires pre-processed GPS Input data structures and ACE satellite data
% polarExB;                % Create new process
% polarExB(load('State')); % Uses prior State structure as starting point
%------------------------------------------------------------------------------------
if nargin==1; Update(State); return; end

% Define starting State
State.Time   = datenum( 2003,10,30,15,0,0 ):...              % Start time
               10/60/24:...                                 % Time step
               datenum( 2003,10,31,06,0,0 );                 % End time
State.t      = 0;                                           % Start index
State.Grid   = initgrid('SPHERICAL',...                     % Grid type
               [0,100:40:1500]*1E3+6365E3,...                 % Grid radial axis
               [-48:4:48]*pi/180,...                        % Grid latitude axis
               [-48:4:48]*pi/180,...                        % Grid longitude axis
               'ROT',geomag.'*[0,0,-1;0,1,0;1,0,0]);        % Grid rotation
State.Itt    = 1;                                           % Filter itterations
State.kA     = 0.95;                                        % Underlying process const
State.PPath  = []; %'./data/leo/4x4';                        % Prior solution
State.SPath  = './data/leo/4x4';                        % Output path and name

% IGS observations, see loadobs.m
State.Obs(1).Path='./data/{yyyy}/{DOY}/{NAME}.mat';         % Data path
State.Obs(1).Name=asciiread('data/sitelists/IGS_NORTH.sites',[inf,4]); % Data file names
State.Obs(1).Type='TRBIAS';                                 % Data type
State.Obs(1).Sample=30;                                     % Data sampling
State.Obs(1).Error =1; %.5                                  % Data error scaling
State.Obs(1).Time=[-10*60-10,20*60+10]/3600/24;             % Data time window

% LEO observations, see loadobs.m
State.Obs(2).Path='./data/{yyyy}/{DOY}/{NAME}.mat';         % Data path
State.Obs(2).Name=['cham';'sacc'];                          % Data file names
State.Obs(2).Type='TRBIAS';                                 % Data type
State.Obs(2).Sample=10;                                     % Data sampling
State.Obs(2).Error =4; %.5                                  % Data error scaling
State.Obs(2).Time=[-10*60-10,20*60+10]/3600/24;             % Data time window


% Temporary variables
Conv   = {1,[1,2,4,2,1],[1,2,4,2,1]};                       % Convolution kernels
Conv   = convmap(size(State.Grid.X),Conv);                  % Mapped convolution

HPath  = './data/solar/magswe_data_64sec_{BART}.hdf';       % HDF By, Bz and solar wind
Prof   = { [1,250, 80,100],...                              % EOF Chapman profiles
           [1,350,100,160],...
           [1,450,140,250],...
           [1,550,200,400] };


%------------------------------------------------------------------------------------
% Define starting state
State.M  = norm({modelEOF(State.Grid,3,Prof),...     % Radial mapping
           diag(sparse(1:length(State.Grid.Yc))),... % Latitude mapping
           diag(sparse(1:length(State.Grid.Zc)))} ); % Longitude mapping
State.x  = State.M.'*iri(State.Grid,State.Time(1));  % Starting estimate unknowns
State.A  = State.M.'*full(Conv)*State.M;             % Forward projection
[i,j,a]  = find(State.A); [m,n]=size(State.A);
State.A  = sparse(i,j,a,m,n); clear('i','j','a');    % Clean matrix
State.P  = State.A*5.^2;                             % Starting covariance
State.Q  = State.A*5.^2;                             % Process noise
State.R  = [];                                       % Observation variances
State.z  = [];                                       % Current observation
State.H  = [];                                       % Observation matrix, H.x=z
State.B  = zeros(size(State.P,1),1);                 % Underlying process
State.S  = 4;                                        % Data rejection factor
State.ID = repmat('VERTICES',size(State.P,1),1);     % Column ID's for unknowns
State.IR = [];                                       % Row ID's for observations

%------------------------------------------------------------------------------------
% Read HDF data for IMF and solar wind, use day either side for interpolation
HDFfile  = unique(filename(HPath,[State.Time(1)-1,State.Time,State.Time(end)+1]),'rows');
HDFfield = {'year','day','hr','min','sec','B_gsm_x','B_gsm_y','B_gsm_z','Vp'};
By = []; Bz = []; Sv = []; t = [];
for i = 1:size(HDFfile,1)
   d  = hdfread(HDFfile(i,:),'MAGSWE_data_64sec','Fields',HDFfield);
   d  = [double(cell2mat(d(1:4)));double(cell2mat(d(5:9)))].';
   t  = [t;datenum(d(:,1),1,1)+d(:,2)-1+d(:,3)/24+d(:,4)/24/60+d(:,5)/24/3600];
   By = [By;d(:,7)]; Bz = [Bz;d(:,8)]; Sv = [Sv;d(:,9)]; d = [];
end
i = find(By>-9999); State.By = interp1(t(i),By(i),State.Time,'linear',mean(By(i)));
i = find(Bz>-9999); State.Bz = interp1(t(i),Bz(i),State.Time,'linear',mean(Bz(i)));
i = find(Sv>-9999); State.Sv = interp1(t(i),Sv(i),State.Time,'linear',mean(Sv(i)));


savestruct('state.mat',State); clear; pack; Update(load('state'));
return;

%************************************************************************************
% Main loop
%************************************************************************************
function Update( State )

while State.t < length(State.Time)
   State.t = State.t + 1;
   stats(State,State.Time(State.t),'HEAD');
   v  = strmatch('VERTICES',State.ID);

% Underlying solution from prior solution or IRI
   if ~isempty(State.PPath)
      P = load(filename([State.PPath,'{yymmmdd-HHMM}'],State.Time(State.t)));
      N = double(P.x(strmatch('VERTICES',P.ID)));
      N = reshape(P.M*N,[length(P.Grid.Xc),length(P.Grid.Yc),length(P.Grid.Zc)]);
      N = interpn(N,[1:size(N,1)].',[1:0.5:size(N,2)],[1:0.5:size(N,3)].');
      N = State.M.'*N(:); clear('P');
   else N = State.M.' * iri(State.Grid,State.Time(State.t)); end

   State.B(v) = N(:)*(1-State.kA); clear('N');

% Forward projection using ExB drift and model
   [Vx,Vz] = ExB( State.Grid, State.Time(State.t),State.By(State.t),...
                  State.Bz(State.t),State.Sv(State.t) );
   Vx = Vx * (State.Time(2)-State.Time(1))*24*3600/State.Itt; % Scale velocities to distances
   Vz = Vz * (State.Time(2)-State.Time(1))*24*3600/State.Itt;
   State.A = spdiag(State.M.'*spdiffuse(State.Grid,[],Vz,Vx,[1,5,5])*State.M*State.kA,...
                    State.A(max(v)+1:end,max(v)+1:end));

% Load data into State
   State = appendobs(State,State.Grid,[]); Out.Obs = State.Obs;
   for i = 1:length(Out.Obs)
      Out.Obs(i).Time = Out.Obs(i).Time+State.Time(State.t);
      State = appendobs(State,State.Grid,loadobs(Out.Obs(i)));
      stats(State,State.Time(State.t),'RAYS');
   end

% Update filter
   for i = 1:State.Itt
      State = kalman( State );
      stats(State,State.Time(State.t),'ERRS');
   end

% Clamp solution to positive
   N = reshape(State.M*State.x(v),size(State.Grid.X));
   N(find(N<0))=0.001; State.x(v) = State.M.'*N(:); clear('N');

% Save output structure, add fields as required
   Out.Time = State.Time(State.t);
   Out.ID   = State.ID;
   Out.IR   = State.IR;
   Out.M    = State.M;
   Out.Grid = struct('ID',State.Grid.ID,'Xc',State.Grid.Xc,'Yc',State.Grid.Yc,...
                     'Zc',State.Grid.Zc,'R',State.Grid.R);
   Out.x    = State.x;
   Out.P    = round(full(diag(State.P))*1E3)/1E3;
   Out.Vx   = round(Vx*1E3)/1E3;
   Out.Vz   = round(Vz*1E3)/1E3;
   OFile    = filename([State.SPath,'{yymmmdd-HHMM}.mat'],State.Time(State.t));
   savestruct(OFile,Out); clear('Out');;

% Plot results
   show(OFile,'FIG1','NORTH','TEC','LOW','EXB');
   
% Clean and re-load state
   savestruct('state.mat',State); clear; pack; State = load('state');
   OFilebmp=filename([State.SPath,'{yymmmdd-HHMM}.bmp'],State.Time(State.t));
   saveas(1,OFilebmp);
end

