%======================================================================
%                    L O A D C T D P R O F . M 
%                    doc: Thu Jun 17 23:17:25 2004
%                    dlm: Fri Mar  5 15:48:54 2010
%                    (c) 2004 ladcp@
%                    uE-Info: 77 13 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================

% CHANGES BY ANT:
%   Jul 17, 2008: - removed dependence on p.nav_start, because that's
%		    not set in [loadnav.m] any more
%   Jul 27, 2008: - nanmean() -> meannan()
%   Jan  7, 2009: - tightened use of exist()

function [d,p]=loadctdprof(f,d,p,ipos,isub)
% function [d,p]=loadctdprof(f,d,p,ipos,isub)
% LADCP-2 software version 7.0 
%
% ========== THIS PART IS SPECIFIC TO LDEO CTD PROFILE DATA ============
%
% you might want to change it to accomodate your own data format
%
%
% load and merges LDEO ctdfile
% ipos(1) : pressure
% ipos(2) : temperature (in situ)
% ipos(3) : salinity
%
% computes sound speed 
% computes static stability
%
% Martin Visbeck, 6/10/99
% revised March 2002
% revised December 2002

% CLIVAR P02 (VANC32) CTD profiles do not contain lat/lon;
% threfore, they are taken from p.nav_start

if nargin<5, isub=1; end
if nargin<4, ipos=[1,2,3]; end

% read ctd data file
disp(['LOADCTDPROF: load CTD profile ',f.ctdprof])
if ~exist(f.ctdprof,'file')
 warn=([' LOADCTDPROF can not find ',f.ctdprof]);
 p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
 disp(warn)
 return
end
[dctd,posctd]=readctd(f.ctdprof,isub,ipos);

if length(posctd) == 2
	d.ctdprof_lon=posctd(2);
	d.ctdprof_lat=posctd(1);
elseif p.navdata
	d.ctdprof_lon = meannan(d.slon);
	d.ctdprof_lat = meannan(d.slat);
else
	error('do not know how to determine CTD position');
end
d.ctdprof_z=p2z(dctd(:,1),d.ctdprof_lat);
d.ctdprof_p=dctd(:,1);
d.ctdprof_t=dctd(:,2);
d.ctdprof_s=dctd(:,3);

% get N^2
if exist('sw_bfrq','file')
 d.ctdprof_N2=sw_bfrq(d.ctdprof_s,d.ctdprof_t,d.ctdprof_p,d.ctdprof_lat);
 d.ctdprof_N2(end+1)=d.ctdprof_N2(end);
 d.ctdprof_ss=sw_svel(d.ctdprof_s,d.ctdprof_t,d.ctdprof_p);
else
 disp(' download SW routines to get soundspeed and N^2 ')
end

% get maximum ctdepth
disp([' CTD max depth : ',int2str(max(d.ctdprof_z))])
if ~isfinite(p.zpar(2));
 p.zpar(2)=fix(max(d.ctdprof_z));
end


%====================================================================

function depth=p2z(p,lat)
% !!!!!! USES Z=0 AT P=0  (I.E. NOT 1ATM AT SEA SURFACE)
%	pressure to depth conversion using
%	saunders&fofonoff's method (deep sea res.,
%	1976,23,109-111)
%	formula refitted for alpha(p,t,s) = eos80
%	units:
%		depth         z        meter
%		pressure      p        dbars  (original in bars, but below
%                                              division by 10 is included)
%		latitude      lat      deg
%	checkvalue:
%		depth =       9712.654  m
%	for
%		p     =         1000.     bars
%		lat   =           30.     deg
%	real lat,p
        if nargin < 2, lat=54; end
        p=p/10.;
	x=sin(lat/57.29578);
	x=x*x;
	gr=9.780318*(1.0+(5.2788e-3+2.36e-5*x)*x)+1.092e-5*p;
	depth=(((-1.82e-11*p+2.279e-7).*p-2.2512e-3).*p+97.2659).*p;
	depth=depth./gr;
%
%
%====================================================================

function [d,pos]=readctd(file,isub,index)
% read CTD data
id=fopen(file,'r');
A=fread(id);
fclose(id);
A=setstr2(A);

pos = [];

%%%ii=find(A(:,1)=='@');
%%%A=A((ii+1):end,:);

[lt,lv]=size(A);

d=[];

if nargin<2, isub=1; end

for i=1:isub:lt
 d=[d;sscanf(A(i,:),'%g')'];
end

if nargin>2
 d=d(:,index);
end


% -----------------------------------------------------
function        AR=setstr2(A,cr)
% function      AR=setstr2(A,cr)
% reshape charater string vector A
% in record structure
% cr = 10 (default)

% remove line feeds
ii=find(A==13);
A(ii)=[];


if nargin<2, cr=10;, end

ii=[0;find(A==cr)];

n=length(ii)-1;
m=max([diff(ii)]);

% fill array with blanks
AR=ones(n,m)*32;

for i=1:n
 j=(ii(i)+1):(ii(i+1)-1);
 if length(j)>0
  j1=j-j(1)+1;
  AR(i,j1)=A(j)';
 end
end

AR=setstr(AR);