function [DATA]=parse_sbedat(dat,confilename)
%PARSE_SBEDAT - convert raw CTD output to scientific units using con file info
% 
% CTD Calibration toolbox.
%
%INPUT
%	dat_load structure, configuration filename
%
%OUTPUT
%	data structure, and configuration coefficients structure
%
%USAGE -	DATA=parse_sbedat(dat,confilename)
%
%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 - Seabird Application Note 27 (and others?)
%	SBE 11plusV2 Deck Unit Users Manual (and others?)
%	read_sbedat.m (ctdcal toolbox)
%	polyval.m (Mathworks)
%

%AUTHORS - 	MATLAB code by c.m.duncombe rae
%
%CREATED -	2008/03/19
%
% VERSION:
%	$Revision: 1.9 $
%	$Date: 2008/03/25 15:25:12 $
%	$Id: parse_sbedat.m,v 1.9 2008/03/25 15:25:12 duncombe Exp duncombe $
%
%CHANGELOG -	
%  2008/03/21 CMDR
%	* adding oxygen conversions
%  2008/03/25 CMDR
%	* dat structure was not completely transferred to dat_load structure
%

% DEPENDENTS - 
% 	Internal functions
% 		function C=cond_cal(freq,t,p,g,h,i,j,CTcor,CPcor)
% 		function T=temp_calc(freq,g,h,i,j,f0)
% 		function [P,PT]=press_calc(freq,U,C,D,T,M,B)
%

%
% PROGRAMMING NOTES
% we should move parsing the con file to another function and use this 
% function to do the calculation of the scientific units
%
% NEWS FLASH: the format of the .CON file is explained in detail in
% Seasave Users Manual, seasave_7.16a.pdf, page 108.
%

%
% not sure yet how the best way to treat this is. skip for now.
%
% if nargin==0, cruiseid=getcruiseid; end;
% if isempty(cruiseid), cruiseid=getcruiseid; end;
% 
% % Retrieve the necessary preferences for this cruise
% group = group_name(cruiseid);

% read a con file
%%%% programming notes: investigate functions fileparts, filesep
if ~exist('confilename'),
	error('cant find con file');
	return;
end;

fid=fopen(confilename);
nrec=0;
conrec=fgetl(fid);
while ~feof(fid),
	nrec=nrec+1;
	switch nrec
	  case {1},	% primary conductivity sensor s/n
		cond_sens1_sn=conrec;
	  case {4},	% primary temperature sensor s/n
		temp_sens1_sn=conrec;
	  case {6},	% secondary conductivity sensor s/n
		cond_sens2_sn=conrec;
	  case {9},	% secondary temperature sensor s/n
		temp_sens2_sn=conrec;
	  case {11},	% pressure sensor s/n
		pres_sens_sn=conrec;
	  case {12},	% pressure sensor T coeffs
		T(1:4)=str2num(conrec);
	  case {13},	% pressure sensor C & T coeffs
		tmp=str2num(conrec);
		C=tmp(1:3); T(5)=tmp(4);
	  case {14},	% pressure sensor D M B coeffs
		tmp=str2num(conrec);
		% haven't figured out what tmp(3:5) represent yet
		D=tmp(1:2); M=tmp(6); B=tmp(7);
	  case {52},
		cond1_cal_date=conrec;
	  case {53},
		temp1_cal_date=conrec;
	  % guessing lines 54 and 55 are secondary sensors calibration dates
	  case {56},
		press_cal_date=conrec;
	  case {110},
		tmp=str2num(conrec);
		c1_coeffs=struct('g',tmp(1),'h',tmp(2), ...
			'i',tmp(3),'j',tmp(4),'CTcor',tmp(5), ...
			'CPcor',tmp(6));
	  case {111},
		tmp=str2num(conrec);
		t1_coeffs=struct('f0',tmp(1),'g',tmp(2), ...
			'h',tmp(3),'i',tmp(4),'j',tmp(5));
	  case {163},
		oxy1_sbe43_sn=conrec;
	  case {164},
		oxy1_cal_date=conrec;
	  case {165},
		tmp=str2num(conrec);
		oxy1_coeffs=struct('Soc',tmp(1),'TCor',tmp(2), ...
			'Voffset',tmp(3));
	  case {166},
		tmp=str2num(conrec);
		if (tmp(2)~=0.0 | tmp(3)~=0.0), 
		    error(['This configuration file (' confilename ...
			') has a non-zero ' ...
			'value in Line 166. Check this file ' ...
			'against the Seabird SBE43 calibration sheet. ' ...
			'You may need to edit the appropriate case (166)' ...
			' in ', mfilename ]);
		    return;
		end;
		oxy1_coeffs.PCor=tmp(1);
		oxy1_coeffs.Boc=tmp(3);
	end;
	conrec=fgetl(fid);
end;

fclose(fid);

% disp('we are done with the con file');
% t1_coeffs

% call the calculators

[Press,PressT]=press_calc(dat.pressfrq,dat.presst,C,D,T,M,B);
Temp1=temp_calc(	...
		dat.ptempfrq,	...
		t1_coeffs.g,	...
		t1_coeffs.h,	...
		t1_coeffs.i,	...
		t1_coeffs.j,	...
		t1_coeffs.f0);
Cond1=cond_cal(dat.pcondfrq, ...
		Temp1, ...
		Press, ...
		c1_coeffs.g, ...
		c1_coeffs.h, ...
		c1_coeffs.i, ...
		c1_coeffs.j, ...
		c1_coeffs.CTcor, ...
		c1_coeffs.CPcor);
% problem of this going imaginary if the conductivity goes negative
Salt1=real(sw_salt(10.*Cond1./sw_c3515,Temp1,Press));
Oxy1=oxy_calc_none(dat.volt0,Salt1,Temp1,Press, ...
		oxy1_coeffs.Soc, ...
		oxy1_coeffs.Voffset, ...
		oxy1_coeffs.Boc, ...
		oxy1_coeffs.TCor, ...
		oxy1_coeffs.PCor);

% create the structures

DATA.srcfilename=dat.srcfilename;
DATA.confilename=confilename;
DATA.station=dat.station;
DATA.pressure=Press;
DATA.pressT=PressT;
DATA.temp1=Temp1;
DATA.cond1=Cond1;
DATA.salt1=Salt1;
DATA.oxy1=Oxy1;
DATA.c1_coeffs=c1_coeffs;
DATA.t1_coeffs=t1_coeffs;
DATA.pres_coeffs=struct('C',C,'D',D,'T',T,'M',M,'B',B);
DATA.oxy1_coeffs=oxy1_coeffs;
DATA.scanno=dat.scanno;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following functions convert the sensor frequencies and voltages
% to scientific units using the SBE calibration coefficients. Note
% that the coefficients are laid out in the opposite order in the 
% seabird calibration sheets to the way that polyval uses them, hence
% the use of fliplr().
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [P,PT]=press_calc(freq,U,C,D,T,M,B)
% C, D, T, M and B coefficients provided by SBE
% U pressure temperature digitization, must be converted to oC by M&B consts
% freq is in Hz, to comply with cal sheet, convert to period in usec
% pressure is calibrated in psia (pound square inch absolute): convert to dbar
% (see eg. Seabird Application Note 27)
prd=1e6./freq;
PT=U.*M+B;
C0=polyval(fliplr(C),PT);
D0=polyval(fliplr(D),PT);
T0=polyval(fliplr(T),PT);
arg=1-T0.*T0./(prd.*prd);
Ppsia=C0.*arg.*(1-D0.*arg);
P=(Ppsia-14.7)./1.4503774;
return;

function T=temp_calc(freq,g,h,i,j,f0)
%
[m,n]=size(freq);
if n>1, error([ 'temp_calc is for single column vectors only']); end;
%
arg=log(f0./freq);
T=1./polyval([j,i,h,g],arg)-273.15;
return;

function C=cond_cal(freq,t,p,g,h,i,j,CTcor,CPcor)
%
[m,n]=size(freq);
if n>1, error('cond_calc is for single column vectors only'); end;
% In this case frequency should be in killohertz
freq=freq./1000;
C=polyval([j,i,h,0,g],freq)./(10.*(1+t.*CTcor+p.*CPcor));
return;

