function [d,p,de]=loadrdi(f,p)

% function [d,p,de]=loadrdi(f,p)

% LADCP-2 software 

%

% Martin Visbeck and Gerd Krahmann , LDEO April-2000

% read RDI raw data

% thanks to Christian Mertens who contributed most of the functions

%

% Added pg_save variable which saves percent good data, 6/2003

% Added added up-looking beam coordiante data, 6/2004

%

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

%                    L O A D R D I . M 

%                    doc: Fri Jun 18 18:21:56 2004

%                    (c) 2004 ladcp@

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

% CHANGES BY ANT

%  Jun 18, 2004: - added p.mask_dn_bins, p.mask_up_bins (later moved)

%  Jun 21, 2004: - clarified large-velocity warning message

%  Jun 22, 2004: - changed large-velocity warning to check only

%		   central hour of cast, as large velocities were found

%	           to occur commonly in uplooker near beginning and end

%		   of cast, but not in middle. The max number of allowed

%		   large velocities before warning is issued was reduced

%		   from 25 to zero, on the other hand.

%		 - removed logging messages that are not useful in general

%		   (e.g. number of bytes per ensemble)

%		 - removed `removed' messages if no ensembles were removed

%  Jun 24, 2004: - BUG: new large-velocity warning code failed for casts

%		        shorter than one hour

%  Jul  4, 2004: - BUG: new large-velocity warning code failed because

%		        in some cases l.tim contains NaNs

%  Jul 21, 2004: - moved bin masking to [edit_data.m]

%  Nov 17, 2007: - bad error message in b2earth()

%		 - removed a lot of commented-out code from b2earth()

%  Nov 18, 2007: - added code for p.allow_3beam_solutions, p.ignore_beam

%  Jan  7, 2009: - tightened use of exist()

%  Oct 28, 2009: - modified p.ignore_beam for dual-head systems

% p=setdefv(p,'pg_save',[1 2 3 4]);

% Default =3 for loadctd_whoi.

p=setdefv(p,'drot',nan);

% how many db should the last bin be below bin 1

p=setdefv(p,'ts_signal_min',-5);

p=setdefv(p,'ignore_beam',[nan nan]);

if nargin<2, p.name='unknown'; end

if existf(p,'poss')

 if isfinite(sum(p.poss))

  drot=magdev(p.poss(1)+p.poss(2)/60, p.poss(3)+p.poss(4)/60);

  if existf(p,'drot')

   if isfinite(p.drot)

    disp([' found drot:',num2str(p.drot),' should be ',num2str(drot)])

   else

    p.drot=drot;

   end

  else

   p.drot=drot;

  end

 end

end

if existf(f,'ladcpdo')==0

 error([' need file name f.ladcpdo  !!!!! '])

end

if ~exist(f.ladcpdo,'file') | length(f.ladcpdo)==sum(f.ladcpdo==' ')

 error([' can not find ADCP data file : ',f.ladcpdo])

end

f=setdefv(f,'ladcpup',' ');

disp('LOADRDI:')

%

% set some defaults, in case they have not been set yet

%

p=setdefv(p,'pglim',0); 

p=setdefv(p,'elim',0.5); 

p=setdefv(p,'vlim',2.5); 

p=setdefv(p,'wlim',0.2);

p=setdefv(p,'timoff',0);

p=setdefv(p,'drot',NaN);

p=setdefv(p,'orig',0);

p=setdefv(p,'weighbin1',1);

p=setdefv(p,'tiltmax',[22 4]);

p=setdefv(p,'name','unknown');

p=setdefv(p,'maxbinrange',0);

p=setdefv(p,'ts_save',0);

p=setdefv(p,'cm_save',0);

p=setdefv(p,'pg_save',3);

p=setdefv(p,'xmv_min',0);

%

% load RDI data

%

[l,message,le]=updown(f.ladcpdo,f.ladcpup,p.pglim,p.elim,...

                      p.maxbinrange,p.ts_save,p.cm_save,p.pg_save,p);

p.warn=l.warn;

% store original data

%

if p.orig

 d.l = l;

end

% check for data set

%

de=0;

if le==1

  disp(' !!!  NO DATA !!! ')

  disp(message)

  de=1;

  return

end

%

% remember which bin come from which instrument (up-down) configuration

% get instrument configuration

%

d.izd=1:length(l.zd);

d.zd=l.zd;

p.serial_cpu_d=l.serial_cpu_d;

p.nping_total=l.npng_d*l.nens_d;

d.bbadcp=l.bbadcp;

p.instid(1)=prod(p.serial_cpu_d+1)+sum(p.serial_cpu_d);

%if d.bbadcp

 [dummy,d.down]=rditype(f.ladcpdo);

 if d.down.Up

  warn=(' up looking instrument detected in do-file');

  p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;

  disp(warn)

  d.zd=-d.zd;

 end

 if existf(l,'zu')

   d.zu=l.zu;

   p.serial_cpu_u=l.serial_cpu_u;

   p.instid(2)=prod(p.serial_cpu_u+1)+sum(p.serial_cpu_u);

   p.nping_total(2)=l.npng_u*l.nens_u;

   d.izu=fliplr(1:length(l.zu));

   d.izd=d.izd+length(d.izu);

   [dummy,d.up]=rditype(f.ladcpup);

   if d.up.Up==0

    warn=(' down looking instrument detected in up-file');

    p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;

    disp(warn)

    d.zu=-d.zu;

   end

 else

   d.izu=[];

   d.zu=[];

 end

%end

%

% apply w velocity threshold

%

d.wrange=5;

izr=d.izd(1:d.wrange);

w(d.izd,:)=meshgrid(medianan(l.w(izr,:),1),d.izd,2);

if existf(l,'zu') 

  izr=[izr,d.izu(1:d.wrange)]; 

  w(d.izu,:)=meshgrid(medianan(l.w(d.izu(1:d.wrange),:),1),d.izu,2);

end

p=setdefv(p,'wizr',izr);

% to normal velocity data

j = find(abs(l.w-w) > p.wlim);

l.u(j) = NaN;

l.v(j) = NaN;

l.w(j) = NaN;

if p.orig

 d.l.problem(j) = d.l.problem(j)+1;

end

% Estimate single ping velocity error from std(W)

nmax=min([length(d.izd),6]);

sw=stdnan(l.w(d.izd(2:nmax),:));

ii=find(sw>0); sw=medianan(sw(ii));

d.down.Single_Ping_Err=sw/tan(d.down.Beam_angle*pi/180)/...

                        sqrt(d.down.Pings_per_Ensemble);

p.beamangle=d.down.Beam_angle;

if existf(l,'zu') 

 nmax=min([length(d.izu),6]);

 sw=stdnan(l.w(d.izu(2:nmax),:));

 ii=find(sw>0); sw=medianan(sw(ii));

 d.up.Single_Ping_Err=sw/tan(d.up.Beam_angle*pi/180)/...

                        sqrt(d.up.Pings_per_Ensemble);

end

% to bottom track velocity data

j = find(abs(l.wb-w(d.izd(1),:)) > p.wlim);

l.ub(j) = NaN;

l.vb(j) = NaN;

l.wb(j) = NaN;

if p.orig

 d.l.problemb(j) = d.l.problemb(j)+1;

end

% Horizontal Velocity limit

%

% to normal velocity data

vel=sqrt(l.u.^2+l.v.^2);

j = find(vel > p.vlim);

l.u(j) = NaN;

l.v(j) = NaN;

l.w(j) = NaN;

if length(j) > 0

  disp(sprintf(' removed %d values because of horizontal speed > %g m/s',length(j),p.vlim));

end

% only warn if large velocities occur during middle hour of cast; this

% excludes near-surface effects when large velocities can be common.

% However, reduce number of allowed large velocities before warning is

% issued from 25 to 10.

nens = size(l.u,2);

enstime = 86400 * (max(l.tim(1,:))-min(l.tim(1,:))) / nens;

jstart = floor(nens/2-1800/enstime);

jend = ceil(nens/2+1800/enstime);

if (jstart < 1), jstart = 1; end

if (jend > length(vel)), jend = length(vel); end

jj = find(vel(jstart:jend) > p.vlim);

skipnens = 1200 / enstime;

[j1,j2] = find(vel > p.vlim);

jj = find(j2>skipnens & j2<size(l.u,2)-skipnens);

if length(jj)>10

 warn=['**  found  ',int2str(length(jj)),'  horizontal velocities > ' int2str(p.vlim) 'm/s in middle hour of cast'];

 disp(warn);

 if length(jj)>100

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

 end

 disp('** WARNING  check maximum velocity setting on CMD-file   ** ')

end

if p.orig

 d.l.problem(j) = d.l.problem(j)+1;

end

% to bottom track velocity data

vel=sqrt(l.ub.^2+l.vb.^2);

j = find(vel > p.vlim);

l.ub(j) = NaN;

l.vb(j) = NaN;

l.wb(j) = NaN;

if p.orig

 d.l.problemb(j) = d.l.problemb(j)+1;

end

%

% apply a time offset, if given

%

if p.timoff~=0

 disp([' WARNING adjusted ADCP time by ',num2str(p.timoff),' days']),

 l.tim=l.tim+p.timoff;

end

% cut time range and apply a time offset, if given

if existf(p,'time_start')==0

  disp(' using whole profile since no start time was given')

  it=find(isfinite(l.tim(1,:)));

  p.time_start=gregoria(l.tim(1,it(1)));

  p.time_end=gregoria(l.tim(1,it(end)));

  d.time_jul=l.tim(1,it);

else

  % fix time for NB-ADCP data

  if l.bbadcp==0

   dum=gregoria(l.tim(1));

   l.tim=l.tim-julian(dum)+julian([p.time_start(1) dum(2:end)]);

   disp(' adjust year for NB-ADCP using given start time ')

  end

  d.time_jul=l.tim(1,:);

  it=find(julian(p.time_start)<=d.time_jul & julian(p.time_end)>=d.time_jul);

  it2=find(julian(p.time_start)>d.time_jul | julian(p.time_end)<d.time_jul);

  if p.orig

   d.l.problem(:,it2) = d.l.problem(:,it2)+10;

   d.l.problemb(:,it2) = d.l.problemb(:,it2)+10;

  end

  d.time_jul = d.time_jul(it);

  disp([' extracting ',int2str(length(it)),' ensembles as profile'])

end

%

% check whether the given time ranges have lead to a profile

%

if length(it)==0

  disp(' ')

  disp(' given time range resulted in empty extracted array!')

  disp(' ')

  disp(' check times')

  disp([' given start time :'])

  disp(p.time_start)

  disp([' internal LADCP start time :'])

  disp(gregoria(l.tim(1,1)))

  disp([' given end time :'])

  disp(p.time_end)

  disp([' internal LADCP end time :'])

  disp(gregoria(maxnan(l.tim(1,:))))

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

  disp('Will try to use all data instead')

  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

  pause(2)

  it=find(isfinite(l.tim(1,:)));

  d.time_jul=l.tim(1,it);

%  p.time_start=gregoria(l.tim(it(1)));

%  p.time_end=gregoria(l.tim(it(end)));

end

d.soundc=0;

%

% rotate for magnetic deviation

%

if isfinite(p.drot)

 [d.ru,d.rv]=uvrot(l.u(:,it),l.v(:,it),p.drot);

 [ub,vb]=uvrot(l.ub(it),l.vb(it),p.drot);

 disp(' apply magnetic deviation, rotate bottom track and water velocities')

else

 d.ru=l.u(:,it);

 d.rv=l.v(:,it);

 ub=l.ub(it);

 vb=l.vb(it);

end

%

% save all bottom track data together and remove bad data

%

d.bvel=[ub',vb',l.wb(it)',l.eb(it)'];

ii=find(d.bvel<-30); 

d.bvel(ii)=NaN; 

d.hbot=l.hb(it);

d.hbot4=l.hb4(:,it);

p.btrk_used=l.btrk_used;

if existf(l,'hs')==1

 d.hsurf=l.hs(it);

end

%

% extract the profile from all recorded data

%

d.firstlastindx = [it(1),it(end)];

d.rw=l.w(:,it);

d.re=l.e(:,it);

d.ts=l.ts(:,it);

if p.ts_save(1)~=0

 d.ts_all_d=l.ts_all_d(it,:,:);

 if size(l.pit,1)==2

  d.ts_all_u=l.ts_all_u(it,:,:);

 end

end

if p.cm_save(1)~=0

 d.cm_all_d=l.cm_all_d(it,:,:);

 if size(l.pit,1)==2

  d.cm_all_u=l.cm_all_u(it,:,:);

 end

end

if p.pg_save(1)~=0

 d.pg_all_d=l.pg_all_d(it,:,:);

 if size(l.pit,1)==2

  d.pg_all_u=l.pg_all_u(it,:,:);

 end

end

d.hdg=l.hdg(:,it);

d.xmc=l.xmc(:,it);

d.xmv=l.xmv(:,it);

d.tint=l.tint(:,it);

d.sv=l.sv(:,it);

d.temp=l.t(:,it);

d.weight=l.cm(:,it);

d.weight=d.weight./medianan(maxnan(d.weight));

d.pit = l.pit(:,it);

d.rol = l.rol(:,it);

%d.tilt=sqrt(l.pit(1,it).^2 + l.rol(1,it).^2);

% more accurate calculation

d.tilt=real(asin(sqrt(sin(l.pit(1,it)/180*pi).^2 +...

 sin(l.rol(1,it)/180*pi).^2)))/pi*180;

% compute tilt difference

rold=mean(abs(diff([0,d.rol(1,:);d.rol(1,:),0]'))');

pitd=mean(abs(diff([0,d.pit(1,:);d.pit(1,:),0]'))');

d.tiltd=sqrt(rold.^2+pitd.^2);

% reduce weight for strong tilt difference

ii=find(d.tilt>p.tiltmax(1));

if length(ii) > 0

  disp([' removed ',num2str(length(ii)),...

        ' profiles due to tilt > ',num2str(p.tiltmax(1)) ' degrees'])

end

d.weight(:,ii)=NaN;

if length(ii)>length(d.tilt)*0.1;

  warn=([' ',int2str(length(ii)*100/length(d.tilt)),...

         '%  tilt > ',int2str(p.tiltmax(1)),' ']);

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;

end

ii=find(d.tiltd>p.tiltmax(2));

if length(ii) > 0

  disp([' removed ',num2str(length(ii)),...

        ' profiles due to tilt derivative > ',num2str(p.tiltmax(2)) ' degrees'])

end        

d.weight(:,ii)=NaN;

% reduce weight for strong echos possibly from crosstalk or bottom

d.tsw=d.weight*0+1;

for i=1:size(d.tsw,1)

 tsmed=median(d.ts(i,:));

 ts=d.ts(i,:)-tsmed;

 ii=find(ts>0);

 if length(ii)>0

  d.weight(i,ii)=d.weight(i,ii).*(1-(ts(ii)/max(ts)).^1.5);

  d.tsw(i,ii)=(1-(ts(ii)/max(ts)).^1.5);

 end

end

% save transmit current volt and internal temperature

for j=1:size(d.xmv,1)

 p.xmc(j)=medianan(d.xmc(j,:),size(d.xmc,2)/4);

 p.xmv(j)=medianan(d.xmv(j,:),size(d.xmv,2)/4);

 p.tint(j)=medianan(d.tint(j,:),size(d.tint,2)/4);

end

% warn if battery low

if p.xmv(1)<p.xmv_min

  warn=([' median Xmit-volt ',num2str(p.xmv(1)),' < ',...

           num2str(p.xmv_min),' BATTERY weak ? ']);

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

end

% check for dead instrument (no pings just listen)

%

% remove suspect ensembles

drw=medianan(abs(diff(d.rw(d.izd,:))));

dru=medianan(abs(diff(d.rv(d.izd,:))));

drv=medianan(abs(diff(d.ru(d.izd,:))));

nbad=find(abs(drw)<0.005 & abs(dru)<0.005 & abs(dru)<0.005);

if length(nbad) > 0.2*length(it)

  warn=([' down looker ',int2str(length(nbad)),' ensembles  ',...

           ' have no flow gradient. ']);

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  warn=(['DOWN LOOKER NOT PINGING ?']);

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  disp(' WARNING WARNING ')

end

if length(nbad)>0

 d.weight(d.izd,nbad)=nan;

 disp([' removed ',int2str(length(nbad)),...

        ' suspect non pinging? (low velocity gradient) ensembles from down-looker'])

end

if length(d.izu)>1

 % remove suspect ensembles

 drw=medianan(abs(diff(d.rw(d.izu,:))));

 dru=medianan(abs(diff(d.rv(d.izu,:))));

 drv=medianan(abs(diff(d.ru(d.izu,:))));

 nbad=find(abs(drw)<0.005 & abs(dru)<0.005 & abs(dru)<0.005);

 if length(nbad) > 0.2*length(it)

   warn=(['   up looker ',int2str(length(nbad)),' ensembles  ',...

            ' have no flow gradient.']);

   disp(warn)

   p.warn(size(p.warn,1)+1,1:length(warn))=warn;

   warn=(['UP LOOKER NOT PINGING ?']);

   disp(warn)

   p.warn(size(p.warn,1)+1,1:length(warn))=warn;

   disp(' WARNING WARNING ')

 end

 if length(nbad)>0

  d.weight(d.izu,nbad)=nan;

  disp([' removed ',int2str(length(nbad)),...

        ' suspect non pinging? (low velocity gradient) ensembles from up-looker'])

 end

end

%

% reduce certainty in bin 1

%

idb1=d.izd(1);

if length(idb1)==1

  d.weight(idb1,:)=d.weight(idb1,:)*p.weighbin1;

end

if length(d.izu)>1

  iub1=find(d.izu==1);

  if length(iub1)==1

    d.weight(iub1,:)=d.weight(iub1,:)*p.weighbin1;

  end

end

if p.weighbin1~=1

 disp([' multiply weight of bin 1 by ',num2str(p.weighbin1)])

end

%

% prepare array

%

d.izm=d.weight+NaN;

%

% save mean correlation and echo amp profile

%

if length(size(l.tsd_m))==2

  d.tsd_m=reshape(l.tsd_m,length(l.tsd_m)/4,4);