function [header,data] = readKMS(cfile)
%   reads KMS data file
%   Usage: [header,data] = readKMS(cfile);

fid = fopen(cfile,'r','b');
%  read first header block
header = rdKMShd(fid);
%    rewind file ....
fseek(fid,0,-1);  % -1 = bof (beginning of file)

%   then read data, skipping header blocks ...
%     more sophisticated reading routine will be needed, but first need to
%     understand more about what's what in the header,  and what info we
%     need for plotting (or other purposes)


% % original from Gary
% blockSize = 92*512;
% blockSizeData = 90*512;
% nData = 30*512;
% nch = sum(header.chActive);
% npts = nData/nch;

% fluxgate settings
blockSize = 272*512;
blockSizeData = 270*512;
nData = 90*512;
nch1 = sum(header.chActive); % active chan in KMS
nch2 = 3;                    % fluxgate
nch  = nch1 + nch2;          
npts = nData/nch;

% end of modifications
nMax = 1000000; %1280000;  % 10000000;
data = zeros(nch,nMax);
done  = 0;
i1 = 1;
while ~done
    [block,count] = fread(fid,[blockSize,1],'uchar'); 
        done = count < blockSize;
    %hd = block(1:1024);
    nRead = count-1024; 
  %  only complete package can be read, uncomplete at the end will be skipped
  %      if nRead > 3*nch  % was from Gary
    if nRead == blockSizeData;
        nData = floor(nRead/3);
        npts = floor(nData/nch);
        temp = reshape(block(1025:end),[3,nData]);
        d = (temp(1,:)*256+temp(2,:))*256+temp(3,:);
  %      d(d > 8388608) = d(d > 8388608)- 16777216;
        d(d > 8388607) = d(d > 8388607)- 16777216;    % from John to get the right value for -2.5V
        
%         d = reshape(d(1:npts*nch),nch,npts);
%         d = d/3355.4432;
        d1 = reshape(d(1:npts*nch1),nch1,npts);
        d1 = d1/3355.4432;     % analogue data
        d2 = reshape(d(npts*nch1+1:npts*nch),nch2,npts);
        d2 = d2.*7.168e-4;     % 2.8*256*1e-6   digital data
        
        i2 = i1+npts-1;
%         data(:,i1:i2) = d;
        data(1:nch1,i1:i2) = d1;
        data(nch1+1:nch,i1:i2) = d2;
        done = count < blockSize;
    else
        done = true;
    end
    i1 = i2+1;
end
data = data(:,1:i2);
fclose(fid);