function [DAT,data,scans,header,rawheader]=read_sbedat(filename)
%READ_SBEDAT - 	read raw CTD output from .DAT file 
% 
% CTD Calibration toolbox.
%
% INPUT
%	.DAT file filename
%	
% OUTPUT
%	DAT, data, scans and header structures with the contents of the .DAT
%	file
%
%USAGE -	[DAT,data,scans,header,rawheader]=read_sbedat(filename)
%
%EXPLANATION -	This module is written in response to the need to access
% .DAT file information without access to SEASOFT. read_sbedat.m loads the
% binary information in the .DAT file and converts it to engineering units 
% (eg. Hz, volts).  parse_sbedat.m reads a .CON file, extracts the necessary
% calibration coefficients from a .CON file and converts to scientific units 
% (eg. dbar, oC, S/m).
%
%SEE ALSO -	
%	Information about the SBE data formats gleaned from Section 12 
%	of SBE 11plus V2 Deck Unit User's Manual (11pV2_010.pdf) 
%	retrieved from Seabird Electronics Inc. website 
%	(http://www.seabird.com/) on or about 2008/02/27.
%	parse_sbedat.m
%

%AUTHORS - 	MATLAB code by c.m.duncombe rae
%
%CREATED -	2008/02/27
%
%VERSION:
%	$Revision: 1.10 $
%	$Date: 2008/03/22 17:22:33 $
%	$Id: read_sbedat.m,v 1.10 2008/03/22 17:22:33 duncombe Exp duncombe $
%
%CHANGELOG -	
%  2008/03/12 CMDR
%	* added comments and provenances
%	* lose the orientation display we used during creation and debugging
%  2008/03/18 CMDR
% 	* parse the header to get the number of bytes per scan and the 
%	number of voltage words
%	* fix an error in the counting of the header 
%	* determine whether system time is appended
%  2008/03/19 CMDR
%	* introduce flags for the data words 
%

% DEPENDENTS
%    Internal functions
%	function X=shape_mat(x,n)
%	function X=msb(x)
% 	function X=lsb(x)
% 

% PROGRAMMING NOTES
% we should parse the con file from this function and add the coefficients to 
% this data structure
%
% NEWS FLASH: the format of the .CON file is explained in detail in
% Seasave Users Manual, seasave_7.16a.pdf, page 108.
%

%	word	byte	description		abbreviation
% 	0	0-2	Primary temperature	T1
% 	1	3-5	Primary conductivity	C1
% 	2	6-8	Pressure		PR
% 	3	9-11	Secondary temperature	T2
% 	4	12-14	Secondary conductivity	C2
% 	5	15-17	A/D channels 0 and 1	AD01
% 	6	18-20	A/D channels 2 and 3	AD23
% 	7	21-23	A/D channels 4 and 5	AD45
% 	8	24-26	A/D channels 6 and 7	AD67
% 	9	27-29	unused and Surface PAR	SPAR
% 	10	30-32	Pressure temperature, 	PTMC
% 			CTD status, modulo count

% for the minute we dont care about the cruiseid etc
%%%  
%%%  if nargin==0, cruiseid=getcruiseid; end;
%%%  if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%  
%%%  % Retrieve the necessary preferences for this cruise
%%%  group = group_name(cruiseid);

% set a default for number of time bytes, may be changed later if we have to
tbytes=0;

 fid=fopen(filename);
% read the entire file as characters
 filecontents=fread(fid,inf,'uchar');
 [M,N]=size(filecontents);
 
% find where the header ends - we are looking for the sequence "*END*"
% (5 bytes) which signifies the end of the header / beginning of the data
%  2008/03/19 should actuallly be looking for [ 13 10 '*END*' 13 10 ]
 I=1:(M-5);
 J=find( filecontents(I)  =='*' &	...
	 filecontents(I+1)=='E' &	...
	 filecontents(I+2)=='N' &	...
	 filecontents(I+3)=='D' &	...
	 filecontents(I+4)=='*'		);

% lines are ended by [ 13 10 ]
rawheader=filecontents(1:J+6);
K=find(rawheader==13);
L=[0;find(rawheader==10)];
for i=1:length(K),
	header{i}=char(rawheader((L(i)+1):(K(i)-1)).');
	if findstr('* Number of Bytes Per Scan',header{i}),
		nbytes=str2num(regexprep(header{i},'[^0-9]*',''));
	elseif findstr('* Number of Voltage Words',header{i}),
		nvwrds=str2num(regexprep(header{i},'[^0-9]*',''));
	elseif findstr('* Append System Time',header{i}),
		tbytes=4;
	elseif findstr('* FileName',header{i}),
		is=findstr(header{i},'.');
		isub=is-2:1:is-1;
		station=sscanf(header{i}(isub),'%d',2);
	end;
end; 

% check our scan structure:
% 	PTMC is always there
%	It is likely that T1 C1 PR are always there
%	nbytes=nPTMC+nPR+nT1+nC1+nvwrds*AD+tbytes+nT2+nC2
PTMC=1; PR=0; T1=0; C1=0; T2=0; C2=0; AD01=0; AD23=0; AD45=0; AD67=0; SPAR=0;
if nbytes==21 & tbytes==0,
	PR=1; T1=1; C1=1; AD01=1; AD23=1; AD45=1;
elseif nbytes==28 & tbytes==4,
	PR=1; T1=1; C1=1; AD01=1; AD23=1; AD45=1; AD67=1;
else
	error(['Problem with scan byte count: nbytes==' num2str(nbytes) ...
		' and tbytes==' num2str(tbytes) ]);
	return;
end;

if nbytes~=(PTMC+PR+T1+C1+T2+C2+AD01+AD23+AD45+AD67+SPAR).*3+tbytes;
	disp(['  nbytes  T1    C1    PR    T2    C2   AD01  AD23  AD45  ' ...
		'AD67  SPAR  PTMC tbytes']);
	disp([nbytes,T1,C1,PR,T2,C2,AD01,AD23,AD45,AD67,SPAR,PTMC,tbytes]);
	error(['Scan bytes checksum failed']);
	return;
end;

data=filecontents(J+7:M);

% we can find a nominal value for nbytes from the headers
% we should parse them. According to header, time is added 
% to the scan. By trial and error I determine that this 
% occupies 4 bytes. These bytes fetch up at the BEGINNING 
% of the scan record as we parse it. We lose these columns.
% nbytes=28; % we scanned the header to find out nbytes above.
% timeb=1:4;
timeb=1:tbytes; 

scans=shape_mat(data,nbytes);

% syncbyte is our byte 33 column. (note it may not be actual byte 33, 
% seasoft eliminates unused bytes, so it may be that only say 13 
% bytes are there depending on the sensors hooked up to the CTD. But the 
% sync byte will be varying monotonically from 0 to 255. This column is the 
% sync byte and we need to lose the first lot of bytes up to and incl. this 
% column.)

% we only need 256 rows to make the identification.
% and it looks like we are going to merely confirm that it is the last
% column and the deck unit or seasoft did our work for us.
% This code will need modification if we ever come across a situation
% where the deck unit has been made to average scans: the counter then 
% will be more than 1 that we assume below.
diffm=diff(scans(1:512,:));
% eliminate the jump from 255 to 0
diffm(diffm==-255)=1;
% we find any column that is varying by one and only one (have to adjust this
% code if we ever get a situation where the deck unit has been set to average
% scans for us)
syncbyte=find(~any(diffm-1));
if syncbyte~=size(scans,2),
	scans=shape_mat(data(syncbyte+1:length(data)),nbytes);
end;

% now we step over the time bytes, adjusting for the way we constructed scans
timerec=scans(2:length(scans),timeb);
scans(length(scans),:)=[];
i=length(timeb);				% 4	0
multiplier=[256; 1; 1./256];

% Frequencies: are encoded directly in the output
%% DATA WORD 0
if T1, ptempfrq=scans(:,[1:3]+i)*multiplier; i=i+3; end;
%% DATA WORD 1
if C1, pcondfrq=scans(:,[1:3]+i)*multiplier; i=i+3; end;
%% DATA WORD 2
if PR, pressfrq=scans(:,[1:3]+i)*multiplier; i=i+3; end;

%% secondary sensors do not exist
% secondary sensor channels may still be live!
%% DATA WORD 3
if T2, stempfrq=scans(:,[1:3]+i)*multiplier; i=i+3; end;
%% DATA WORD 4
if C2, scondfrq=scans(:,[1:3]+i)*multiplier; i=i+3; end;

% Voltages: must be calculated from the AD output
%% DATA WORD 5
if AD01,
	N=scans(:,1+i).*16+msb(scans(:,i+2)); volt0=5.*(1-(N./4095)); i=i+1;
	N=lsb(scans(:,1+i)).*256+scans(:,i+2); volt1=5.*(1-(N./4095)); i=i+2;
end;
%% DATA WORD 6
if AD23,
	N=scans(:,1+i).*16+msb(scans(:,i+2)); volt2=5.*(1-(N./4095)); i=i+1;
	N=lsb(scans(:,1+i)).*256+scans(:,i+2); volt3=5.*(1-(N./4095)); i=i+2;
end;

% rest of the sensors do not exist 6 bytes		27
%% DATA WORD 7
% A/D channels 4 and 5
if AD45,
	N=scans(:,1+i).*16+msb(scans(:,i+2)); volt4=5.*(1-(N./4095)); i=i+1;
	N=lsb(scans(:,1+i)).*256+scans(:,i+2); volt5=5.*(1-(N./4095)); i=i+2;
end;

%% DATA WORD 8
% A/D channels 6 and 7
if AD67,
	N=scans(:,1+i).*16+msb(scans(:,i+2)); volt6=5.*(1-(N./4095)); i=i+1;
	N=lsb(scans(:,1+i)).*256+scans(:,i+2); volt7=5.*(1-(N./4095)); i=i+2;
end;

%% DATA WORD 9
% first byte of DW 9 is unused, but the word is still reported
% SPAR and PAR channels are in data word 9
% seems we dont have a spar after all
if SPAR,
	i=i+1;
	spar=msb(scans(:,1+i)).*256+scans(:,i+2); i+1;
end;
% 
%% DATA WORD 10
% DW10 is always present, it has the pressure sensor temperature, the 
% ctd status and the modulo count
if PTMC,
	presst=scans(:,1+i).*16+msb(scans(:,i+2)); i+1;
	i=i+1;
	ctdstat=lsb(scans(:,1+i)); i=i+1;
	sync=scans(:,i+1); % we should have reached the syncbyte 
end;
%% END OF SCAN

% Not sure where the scan number should start, lets make it at the beginning
scanno=[1:length(scans)].';

%% BACKWARDS
%% % parse_scans(data);
%% i=1;
%%  ptempfrq=scans(:,i)./256+scans(:,i+1)+scans(:,i+2).*256;
%% i=i+3;
%%  pcondfrq=scans(:,i)./256+scans(:,i+1)+scans(:,i+2).*256;
%% i=i+3;
%%  pressfrq=scans(:,i)./256+scans(:,i+1)+scans(:,i+2).*256;
%% i=i+3;
%%  stempfrq=scans(:,i)./256+scans(:,i+1)+scans(:,i+2).*256;
%% i=i+3;
%%  scondfrq=scans(:,i)./256+scans(:,i+1)+scans(:,i+2).*256;

% Create a structure with our data in it:
% We are going to need to read a con file!

% We figured out the data, we can drop the orientation display
% figure; subplot(1,3,1); plot(ptempfrq);
% 	subplot(1,3,2); plot(pcondfrq);
% 	subplot(1,3,3); plot(pressfrq);

DAT=struct(	 'srcfilename',filename ...
		,'station',station ...
		,'timerec',timerec ...
		,'ptempfrq',ptempfrq ...
		,'pcondfrq',pcondfrq ...
		,'pressfrq',pressfrq ...
		,'volt0',volt0 ...
		,'volt1',volt1 ...
		,'volt2',volt2 ...
		,'volt3',volt3 ...
		,'presst',presst ...
		,'ctdstat',ctdstat ...
		,'scanno',scanno ...
		,'sync',sync ...
		);


fclose(fid);
 
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Internal functions are below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function X=shape_mat(x,n)
m=size(x(:));
if rem(m(1),n)~=0, fill=zeros(n-rem(m(1),n),1); x=[x(:);fill]; end;
X=reshape(x,n,[]).'; 
return;

function X=lsb(x)
% x represents a byte, 0<=x<=255. X is the least significant bits of x,
% 0<=X<=127.
X=bitand(x,15);
return;

function X=msb(x)
% x represents a byte, 0<=x<=255. X is the most significant bits of x,
% 0<=X<=127.
X=bitand(x,240)./16;
return;

