% inv6f.m - Invert 6-frequency TAPS data from PLOT6TAPS
% Version with Fluorometer plot
%
% INPUTS:
%     Scattering model matrix (from MAKETS.M)
%     Data file (from PLOT6.M)
%
% OUTPUTS:
%     Plots of BV vrs depth and BV vrs depth by ESR
%     Matlab mat-file containing the following:
%         a = ESR's, a 1 X NA vector
%       fxv = TAPS frequencies, a 1 X 6 vector
%    adepth = depth bins, a 1 X NDEPTHS vector
%      biov = BV's, a NA X NDEPTHS matrix
%     atemp = binned temperatures, a 1 X NDEPTHS vector
%
% Orig 12/30/92;    dvh
% Modified 2/5/94   cfg  to work with TAPSPLOT
% Modified 3/24/94  dvh  fix problem with TAijt
% Modified 10/12/94 cfg  for 6-frequency inverse
% Modified 10/28/94 cfg  to plot predicted Sv at more freqs
% Modified 03/14/97 cfg  to use myinv function call for NNLS
% Modified 08/12/98 cfg  improve user interface
% Modified 06/04/01 cfg  to use GUI in/out


fprintf('\nINV6 -- Invert TAPS-6 CAST data to estimate size-abundance\n')
fprintf('                   04 June 2001\n\n')

home=pwd;					% save current path structure

%  Check for model on workspace, load if not found

isTAij = exist('TAij');ans='n';
if isTAij
	ans=input('A model exists on workspace, keep it (y/n?): ','s');
end
if ~isTAij | ans ~= 'y'
	[model, pname] = uigetfile('*.*', 'Select MODEL file');
	cd(pname);
	cmd=['load ' model];
	eval(cmd)
	TAij = 10 .^ ((Aij)/10); 
end

% -------- Input data file from PLOT6.M ----------------------------

[filename, pname] = uigetfile('*.mat', 'Select DATA file');
cd(pname);
cmd=['load ',filename];
eval(cmd)
	
%   Opportunity to delete a specific frequency

figure(1);clf
semilogx(fxv/1000,bx,'o',fxv/1000,bx);grid
title('Volume scattering strength spectra')
xlabel('Frequency in kHz')
ylabel('Sv in dB')

dedfreq=input('Delete freq # (0 => none):');
while dedfreq>0
	[n m]=size(bx);
	for k=1:n;
		bx(k,dedfreq)=NaN;
	end
	dedfreq=input('Delete freq # (0 => none):');
end

%  Enter/accept Levenberq-Marquadt factor
	
sfactr = 1e-4;    % nominal value
fprintf('sfactr = %3.1e : ',sfactr)
ans=input('new value = ');
if length(ans)>0
	sfactr=ans;
end

% ------------------ Main loop on indx (depth) ----------------------------

biov=zeros(length(adepth),length(a));
rnorm=zeros(1,length(adepth));
for indx = 1: length(adepth);
	sv = bx(indx,:);                     % select one row of input data
	finitesv=find(isfinite(sv));
	if (length(finitesv)>2 & abs(max(sv))>1) % only do inverse if have 3+ freqs
		
		[x,stuff] = nnlsinv(sv,TAij,sfactr); % compute inverse solution
		
		bv = (4*pi/3) * (a .^ 3)' .* x;      % compute BV from abundance
		rnorm(indx) = stuff(1);
		xnorm(indx) = stuff(2);
		
		abund(indx,:) = x';
		biov(indx,:) = bv';
		
		%                      PLOT RESULTS
		
		h1=figure(1);clf
		
		xmax = max(a);
		ymin = 0;
		ymax = 1.2*max(biov(indx,:));
		if ymax<=ymin;ymax=ymin+1;end
		depstr=['Depth = ' num2str(adepth(indx))];
		
		%                      LEFT SIDE PLOT
		subplot(121);
		bar(a',biov(indx,:));
		axis([0 xmax ymin ymax]);
		grid;
		title('BV vs size')
		xlabel('ESR in mm')
		ylabel('BV')
		
		%                     RIGHT SIDE PLOT
		subplot(122)
		zz=paij*x;
		semilogx(fxp/1000,10*log10(zz),fxv/1000,sv,'ro');
		grid;
		axis([1e2 1e4 -90 -20])
		xlabel('Frequency in KHz')
		ylabel('Sv in dB')
		title('Sv vs f')
		text(110,-23,depstr);
		pause(0.1);
	end										% of if length(finitesv)..
end                                         % of indx loop

%---------------------  END OF DEPTH LOOP ----------------------
%             Compute total biovolume vrs depth
	
clear tbv
tbv=sum(biov');
	
%--------------------------------------------------------------
%		PLOT BIOVOLUME/RESIDUAL ERRORS VRS DEPTH
%--------------------------------------------------------------
	
% Set depth scale for plots
	
mxd=max(abs(adepth));
dd=1;
if mxd>24;dd=5;end
if mxd>50;dd=10;end
if mxd>100;dd=20;end
if mxd>200;dd=50;end
mxd = -ceil(mxd/dd+.1)*dd;

% Set BV scales for plots

figure(3);clf
semilogx(tbv,adepth);grid
title('Total Biovolume vrs Depth')
xlabel('Total BV in mm^3/m^3')
ylabel('Depth')

mtbv=mean(tbv);				% mean total biovolume
mxtbv=max(tbv);				% max total biovolume

if mxtbv/mtbv>3
	mx10=fix(log10(mxtbv)+1);
else
	mx10=fix(log10(mtbv)+1);
end
mx=10^mx10;
mn10=mx10-3;
mn=10^mn10;
if mn<1;mn=1;end

fprintf('Min BV = %6.0f  : enter new value ',mn);
num=input('','s');if length(num)>0,mn=str2num(num);end
fprintf('Max BV = %6.0f  : enter new value ',mx);
num=input('','s');if length(num)>0,mx=str2num(num);end


%         Plot TBV vrs depth & RNORM vrs depth

h2=figure(2);clf;
set(h2,'PaperOrientation','landscape','PaperUnits','normalized');
set(h2,'PaperPosition',[.15 .15 .75 .75]);
	
subplot(121)
iv=find(tbv>0);
semilogx(tbv,adepth);grid
axis([mn mx mxd 0]);
xlabel('BV - mm^3 / M^3')
ylabel('Depth - M')
title('Total biomass vrs depth')

subplot(122)
semilogx(rnorm,adepth);grid
axis([.01 10 mxd 0]);
title('Norm of residual error')
xlabel('Rnorm')

%          PLOT BIOV (LHS) AND ESR/BIOV CONTOUR (RHS)
	
h3=figure(3);clf;
set(h3,'PaperOrientation','landscape','PaperUnits','normalized');
set(h3,'PaperPosition',[.15 .15 .75 .75]);
	
subplot(141)
mnt = floor(min(t/2))*2;
mxt = ceil(max(t/2))*2;
plot(t,depth,'-r');grid
axis([mnt mxt mxd 0]);
title('Temperature','FontWeight','bold','FontSize',12)
ylabel('Depth - M','FontWeight','bold','FontSize',12);
xlabel('C','FontWeight','bold','FontSize',12)
set(gca,'FontWeight','bold','LineWidth',2);
	
subplot(142)
plot(afluor,adepth,'-g');grid
axis([mnf mxf mxd 0]);
title('Fluorescence','FontWeight','bold','FontSize',12)
xlabel('ug / l','FontWeight','bold','FontSize',12)
set(gca,'FontWeight','bold','LineWidth',2);
	
subplot(143)
semilogx(tbv,adepth,'-k');grid
axis([mx 1e4*mx+.1 mxd 0]);
title('Biomass','FontWeight','bold','FontSize',12)
xlabel('mm^3 / M^3','FontWeight','bold','FontSize',12)
set(gca,'FontWeight','bold','LineWidth',2);
	
subplot(144)
cmap=jet;
%cmap(1,1:3)=[0 0 0];       		% un-comment to make zero biomass plot white
lbv=log10(biov+eps);
mina=min(a);maxa=max(a);
colormap(cmap);
surf(a,adepth,lbv);
caxis([-2 2]);
shading('flat');
view(2)
axis([mina maxa mxd 0]);
title('LOG(BV) by ESR','FontWeight','bold','FontSize',12)
xlabel('ESR - mm','FontWeight','bold','FontSize',12)
set(gca,'FontWeight','bold','LineWidth',2);
colorbar('vert')
	
ans=input('Save processed data (y/n)? ','s');
if length(ans)>0
   	if ans == 'y'
		[outfile, pname] = uiputfile('*.*', 'Select OUTPUT file');
		if pname ~= pwd
			cd(pname);
		end
		cmd=['save ' outfile ' a fxv adepth biov'];
		eval(cmd)
	end
end
fprintf('\n')
chdir(home)

