function [piddat, vect] = readWorkfile(fname)
% function [piddat, vect] = readWorkfile(fname)
%
% read VERTCON workfile resulting from filter output; 
% return list of PIDS
% and corresponding lat/lon and ortho height difference,
% plus metadata on the datums 
%
% INPUT
% fname name of file; may include path
%
% OUTPUT
% piddat  1 name of point (PID)
%         2 datum1
%         3 datum2
% vect   dimensions (#PID x 4)
%         1-2 lat/long (deg)
%         3   ortho height at first datum (m)
%         4   differential ortho height (m)

% check if *mat file exists, otherwise load from ASCII
[f,n,e] = fileparts(fname);
fnameMat = [n e '.mat'];
% if exist(fnameMat)
%     fprintf('loading from saved *.mat file\n');
%     tic;
%     load(fnameMat);
%     toc
% else
    fprintf('loading from ASCII file\n');
    tic;
    [pid, stateplus, lat, lon, ht, dht, datum1, datum2] = textread(char(fname),...
           '%s %3c %f %f %f %f %s %s %*d %*d %*s %*s %*d %*d %*s %*s');
    toc
    save(fnameMat);
    % columns: pid, state(2col)+exclude_code(1col),
    %    lat, lon, ortho_height (first datum?),
    %    ortho_height_difference, datum1, datum2,
    %    (2x) 4-col quality codes
% end

% do not plot points which have been excluded by workedits/filtering
state = stateplus(:,1:2);
excludeCode = stateplus(:,3);
q = find(double(excludeCode)~=32); % not blank (bad data)

% then parse into output arrays
vect = [lat lon ht dht]; vect(q,:) = [];
pid(q) = []; piddat{1} = pid;
datum1(q) = []; piddat{2} = datum1;
datum2(q) = []; piddat{3} = datum2;

return

