function lw = prcprr(prrFile, cflFile);
% PRCPRR    - Process PRR files into data products

% Brian Schlining
% 27 Jan 2000

close all
%=======================================================================
% Make assumptions about the configuration, This may need to be changed
% for different setups.
%=======================================================================
Ed = 1:7; % This is true for most of the files. This is the profiled downwelled irradiance

Marker = 'sdv^<>ph*x+.';
ScreenSize = get(0, 'ScreenSize');
warning off

%===============
% Read the data
%===============
fprintf(1, '\n***Begin Processing***\n');
eu0 = [];
crd = [];
cfl = [];
switch nargin
case 0
   [eu0 crd cfl] = getprrdata;
case 1
   [eu0 crd cfl] = getprrdata(prrFile);
case 2
   [eu0 crd cfl] = getprrdata(prrFile, cflFile);   
end

if isempty(eu0)
   lw = [];
   return
end

%==================
% Find Lu channels
%==================
channels = strvcat(cfl.info.description);
Lu1 = strmatch('Lu', channels);
Lu2 = strmatch('GND', channels(:,3:end));	% Ignore the ground channel
Lu = setdiff(Lu1, Lu2)';
lambda = str2num(channels(Lu,3:end))';

%==================
% Find Es channels
%==================
% Hide the Ed(profile) channels, they have the same names as the Es channels
% and can make it confusing to figure out which is Es and which is Ed.
channels = channels(max(Ed)+1:end, :);		
Es1 = strmatch('Ed', channels);
Es2 = strmatch('GND', channels(:,3:end));	% ignore the ground channel
Es = [setdiff(Es1, Es2) + max(Ed)]';

%==========================================
% Find the column that contains depth info
%==========================================
fields = lower(strvcat(cfl.info.description));
d      = find(strcmp(lower({cfl.info.description}), 'depth'));
fprintf(1, 'Correcting depth offset');
% Correct the depth offset by assuming the minimum pressure read is the surface pressure
% This is going to change with atmospheric conditions so its mostly true.
offset   = min(eu0(:,d));
eu0(:,d) = eu0(:,d)- offset;
fprintf(1, '...completed\n');

%=================================================================================
% Loop until acceptable Lw's are calculated. This loop is controlled by user input
%=================================================================================
OK = 0;
iteration = 0;		% Count the number of times the user tries different depth limits
while ~OK
   iteration = iteration + 1;
   %=========================================
   % Find the top and bottom of the downcast
   %=========================================
   downcast(eu0, crd, cfl);	
   fprintf(1, '\n*****************************************************************\n');
   fprintf(1, 'Select depth limits on the figure, then press any key to continue\n');
   fprintf(1, '*****************************************************************\n\n');
   pause;
   UserData = get(gcf, 'UserData');		% The limits (selected by depth but are really the index of data)
   close(gcf);									% are stored in the figures UserData. Get this data before 
   lims = sort(UserData.lims);			% closing the figure
   eu = eu0([lims(1):lims(2)],:);
   
   %=================================
   % Edit for large variations in Es
   %=================================
   [r c] = size(eu);
   flag  = ones(r,1);
   
   %fprintf(1, 'Despikeing Surface Irradiance data with BBOPKQ.M.')
   %for i = Es						 			
   %   ddw  = 6;
   %   m1   = nanmean(eu(:,i))*.05;  	% Maxiumum 1st difference should vary less than 5%
   %   mStd = nanstd(eu(:,i))*1.5;   	% Each difference whould be within 2 std deviations
   %   flag = bbopkq(eu(:,i), eu(:,d), 10, m1, mStd).*flag;
   %end
   %fprintf(1, '...completed\n');
   
   good = find(flag == 1);
   eu   = eu(good,:);
   
   %========================
   % Bin data into 1 m bins
   %========================
   dww = 1;
   fprintf(1, 'Binning data into %i meter wide bins with BBOPBIN.M.', dww);
   eu = bbopbin(eu, cfl, dww);
   fprintf(1, '...completed\n');
   
   %===========================================================
   % Apply an optimal curve-fit
   % where y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t)
   %===========================================================
   % This requires an older version of matlabs optimization toolbox. I've included it with
   % This program so don't distribute this outside of MBARI
   OPTIONS = [0 1e-3 0 0 1];	% Parameters for LEASTSQ
   lam0 = [1]';					% This is a first order log fit. To do a 2nd order change to [1 0]'
   [r c] = size(eu);
   zf = [0:dww:max(eu(:,d))]';
   yf = ones(length(zf), c)*NaN;		% Variable to store the fitted y values
   n = 0;
   
   for i = [Lu]   
      
      n = n + 1;
      warning off
      fprintf(1, 'Curve fitting %s (column %i)', deblank(fields(i,:)),i);
      
      % Use an optimal fit to the entire water column for the lower wavelengths
      if lambda(n) < 600 & min(eu(:,d)) < 4 %ignore long lambda and casts without near surface data
         [lam,OPTIONS, F, J] = leastsq('fitfunexp', lam0, OPTIONS, [], eu(:,d), eu(:,i));
         c   = fitfunexpc(lam, eu(:,d), eu(:,i));
         yf(:,i)  = funexp(lam, c, zf); 
         lw(1,n) = funexp(lam, c, 0)*.544;
         
      else
         % Use a simple log-linear fit on data between the surface and 12m for the long wavelengths.
         % There is a signifigant contribution to the wavelengths by flourescence. This creates
         % a subsurface peak in Lu which really confounds curve fitting. This method may not
         % look like it fits the overall profile but it seems to give reasonable Lw's
         good = find(eu(:,d) >= 0 & eu(:,d) <= 15);
         X = [eu(good,d) ones(size(eu(good,d)))];
         a = X\log(eu(good,i));
         yf(:,i) = exp(a(1).*zf + a(2));
         lw(1,n) = exp(a(1)*0 + a(2))*0.544;
         
      end
      
      fprintf(1,'....completed\n');
      warning on
   end
   
   %============================================================
   % Plot Lu profiles with fitted curve to see if they are good
   %============================================================
   fhp = figure;
   Position = get(fhp, 'Position');
   set(fhp, 'Position', [15 ScreenSize(4)-45-Position(4) Position(3) Position(4)]);
   fr = 0;
   for i = [Lu]
      fr = fr + 1;
      subplot(3, 3, fr)
      plot(eu(:,i), eu(:,d), '.','Color', [0.74 0.74 1])
      hold on
      plot(yf(:,i), zf, 'r-');
      hold off
      title(cfl.info(i).description);
      set(gca,'YDir','reverse', 'YScale' ,'linear','Xscale', 'log');
      grid
   end
   
   %===================================================================
   %Plot each round of Lw to track the changes due to depth selections
   %===================================================================
   fh = figure;
   Position = get(fh, 'Position');
   set(fh, 'Position', [ScreenSize(3)-15-Position(3) 15 Position(3) Position(4)]);

   if iteration > 1
      for it = 1:iteration-1
         h(it+1) = plot(lambda, lwTmp(it,:), Marker(it), 'Color', [0.74 0.74 1]);
         if it == 1
            hold on
         end
         plot(lambda, lwTmp(it,:), '-', 'Color', [0.74 0.74 1])
         legendS{it+1} = ['Round ' num2str(it)];
         
      end
      YLim = get(gca,'YLim');
      if max(YLim) > 1
         set(gca,'YLim', [-0.05 1])
      end
      
   end
   
   %========================================
   % Plot the most recently calculated Lw's
   %========================================
   h(1) = plot(lambda, lw, 'ro');
   hold on
   legendS{1} = 'Current';
   xlabel('Wavelength [nm]')
   ylabel('Water-leaving radiance, Lw (\muW cm^-^2 nm^-^1 str^-^1)')
   plot(lambda, lw, 'r-');
   hold off
   if iteration > 1
      legend(h, legendS,0);
   end
   
   
   fprintf(1, '\n****************************************\n');
   yn = input('Do you want to accept these Lw''s(y/n)?  ', 's');
   fprintf(1, '****************************************\n\n');
   if findstr('y', lower(yn))
      OK = 1;
   else
      lwTmp(iteration,:) = lw;
   end
   if ishandle(fh)
      close(fh);
   end
   if ishandle(fhp)
      close(fhp);
   end
   
   
   
end

%====================================
% Calculate K for Lu and Ed channels
%====================================
eu                = sortrows(eu, d);
K                 = ones(size(eu))*NaN;
[r c]             = size(eu);
K(2:end, [Lu Ed]) = calc_k(eu(1:r-1,[Lu Ed]), eu(2:r,[Lu Ed]), diff(repmat(eu(:,d),1, length([Lu Ed])))); 

%========================================
% Write the data to a seabass style file
%========================================
[p filename ext] = fileparts(crd.src);
outFilename = fullfile(p, [filename '.prc']);
fprintf(1, 'Writing data to file, %s, using WRITEPRC.M.', outFilename);
writeprc(outFilename, eu, crd, cfl, K, lw, Es, Ed, Lu, lambda);
fprintf(1, '...completed\n');

lwnFilename = fullfile(p, [filename '.lwn']);
writelwn(lwnFilename, eu, crd, cfl, lw, Es, Lu, lambda)
fprintf(1, 'Writing data to file, %s, using WRITELWN.M.', lwnFilename);
fprintf(1, '...completed\n');

fprintf(1, '***Processing Completed***\n\n');
warning on