function E = ep_cmpl(Flag)
% EP_CMPL - Compile eqpac files. 
% Reads files in ep1 and ep2, compiles them, and creates derived products
%
% EP_CMPL Reads all the files created by Bob Herlein's programs in directories ep1 and eo2
% After reading the nescessary files, It calculates derived products from them and 
% writes the correct headers
%
% Use As: E = ep_cmpl;
%         E = ep_cmpl(EditFlag);
%
% Inputs: EditFlag can equal 0 or 1
%         1: (default) edit for outliers in pressure, salinity, and temperature
%           If bad data is found the entire row is marked bad.
%         0: Do not edit for outliers
%
% Ouputs: E = Structure with the following fields
%             id: ep id# (ex. 'ep1')
%             dat: matrix of data, refer to INIEP for column names

% Brian Schlining
% 29 Jan 1999

if ~nargin
   EditFlag = 1;
end


%=======================================
% Loop through ep1 and ep 2 directories
%=======================================
for ep = 1:2
   
   %=================================
   % Set path based on computer type
   %=================================
   switch computer
   case 'PCWIN'
      PathS = ['\\tsunami.shore.mbari.org\oasis\eqpac\ep' num2str(ep) '\'];
   otherwise
      PathS      = ['/hosts/tsunami/oasis/eqpac/ep' num2str(ep) '/'];
   end
   
   %=================================================================
   % Combine the noon and 12:30 files. Most of the work is done here
   % using the subfunction COALLATE.
   %=================================================================
   epx = coallate(PathS);   
   
   %=============================================================
   % Locate important variables in epx to simplify the equations
   % If modifications are made to COALLATE these may have to be
   % changes
   %=============================================================
   Pressure     = repmat(epx(:,2),1,7);
   Lu_2m        = epx(:,13:19);  % Actually Lu(1.5m) but rounded for naming purposes, can't use decimals in variable names
   Lu_20m       = epx(:,28:34);
   Ed_20m       = epx(:,21:26);
   Lu490_0m     = epx(:,15);
   Lu555_0m     = epx(:,17);
   Lu490_20m    = epx(:,30);
   Lu555_20m    = epx(:,32);
   Lu443_0m     = epx(:,14);
   Lu443_20m    = epx(:,29);
   Lu683_0m     = epx(:,19);
   Lu683_20m    = epx(:,34);
   Ed_0m        = epx(:,6:11);
   
   % Edit for outliers. This assumes pressure is in column 2, Salinity in column
   if EditFlag
      epx = editep(epx);
   end
   
   %=================================================
   % Calculate Derived products from measured values
   %=================================================
   % Calculate ELw
   % ELw = Ed(0m+) = a*Lu(1.5m) + b
   % Where Lu comes from the OCR100 at 1.5m depth
   [r c]        = size(epx);
   Es           = c + 1;
   Ee           = c + 7;
   aLw          = [0.5574 0.5588 0.5632 0.5778 0.6005 0.6713 0.5904];
   bLw          = [0.014 0.012 0.003 0.003 0.002 0.001 0.001];
   for i = 1:7
      epx(:,Es+i-1) = aLw(i).*Lu_2m(:,i) + bLw(i); % Get these coef
   end
   
   % Calculate LLw
   % LLw = Lu(0m+) = 0.544*exp(((ln(Lu(20m))-ln(Lu(1.5m))/(depth-1.5))*-1.5)+ln(Lu(1.5m)))
   % This is calculated from a linear extrapolation of the 20m and 1.5m Lu data.
   [r c]        = size(epx);
   Ls           = c + 1;
   Le           = c + 7;
   epx(:,Ls:Le) = 0.544*exp((((log(Lu_20m)-log(Lu_2m))./(Pressure-1.5))*-1.5)+log(Lu_2m));
   
   % Calculate KLw
   % Klw = 0.544*(0.98*Ed(0m+))/Ed(20m)*Lu(20m)
   % KLw is Lu(0m+) calculated by back-propagating the 20m Lu data, using K values at each wavelength
   [r c]        = size(epx);
   Ks           = c + 1;
   Ke           = c + 6;
   epx(:,Ks:Ke) = 0.544*(0.98.*Ed_0m)./Ed_20m.*Lu_20m(:,1:6);
   
   % Calculate K, attenuation coeeficeints
   % From p10750 in Morelm 1988
   % K = (ln(0.98*Ed(0m+)/Ed(20m))/Depth
   [r c]        = size(epx);
   ks           = c + 1;
   ke           = c + 6;
   epx(:,ks:ke) = (log(0.98*Ed_0m./Ed_20m))./Pressure(:,1:6);
   
   % Calculated Morel Chl
   % Chl = ((K490 - 0.0217)/0.0690)^(1/0.702)
   % This is a rearrangement of eqn(9) in Morel 1988, using the coeff in table 2 
   [r c]     = size(epx);
   Cs        = c + 1;
   epx(:,Cs) = ((epx(:,ks+2) - 0.0217)/0.0690).^(1/0.702);
   BAD = find(epx(:,Cs) < 0);
   if ~isempty(BAD)
      epx(BAD,Cs) = NaN;
   end
   
   
   % Add OC2V2 Chl
   % chl = a_oc2_c(LU_490,LU_555)
   [r c]      = size(epx);
   C0         = c + 1;
   C20        = c + 2;
   epx(:,C0)  = a_oc2_c(Lu490_0m,Lu555_0m);
   epx(:,C20) = a_oc2_c(Lu490_20m,Lu555_20m);
   
   % Calculate Lu443/Lu555 Chl
   % Chl(0m) = 0.7469*(Lu443(0m)/Lu555(0m))^-1.2797;
   % Chl(20m) = 1.0498*(Lu443(20m)/Lu555(20m))^-1.0948
   [r c]       = size(epx);
   Cx0         = c + 1;
   Cx20        = c + 2;
   epx(:,Cx0)  = 0.7469*(Lu443_0m./Lu555_0m).^-1.2797;
   BAD = find(epx(:,Cx0) < 0);
   if ~isempty(BAD)
      epx(BAD,Cx0) = NaN;
   end
   
   epx(:,Cx20) = 1.0498*(Lu443_20m./Lu555_20m).^-1.0948;
   BAD = find(epx(:,Cx20) < 0);
   if ~isempty(BAD)
      epx(BAD,Cx20) = NaN;
   end
   
   
   
   % Add Fouling Ratio = (Lu490/Lu555)
   [r c]      = size(epx);
   F0         = c + 1;
   F20        = c + 2;
   epx(:,F0)  = Lu683_0m./Lu555_0m;
   epx(:,F20) = Lu683_20m./Lu555_20m;
    
   % Add CZCS Chl: log10(C) = 0.053 + 1.71*log10(Lw550/Lw443)
   [r c] = size(epx);
   Cz    = c + 1;
   epx(:,Cz) = 10.^(0.053 + 1.71*log10(epx(:,47)./epx(:,44)));

   E(ep).id  = ['ep' num2str(ep)];
   E(ep).dat = epx;
   
   
end

%======================================================================================
%======================================================================================
% subfunctions
%======================================================================================
%======================================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COALLATE - Reads the files in ep1 and ep2 and combine into one file
% with alternating noon and 1230 lines the columns are as follows
%  1          -  1200 or 1230 (indicating sampiling at noon or 1230am
%  2          - Pressure
%  3          - DecimalYear
%  4:5        - Temperature Salinity
%  6:19  (0m) - Ed412 Ed443 Ed490 Ed510 Ed555 Ed656 PAR Lu412 Lu443 Lu490 Lu510 Lu555 Lu670 Lu683
% 20          - MCP at 10 m  
% 21:34 (20m) - Ed412 Ed443 Ed490 Ed510 Ed555 Ed670 PAR Lu412 Lu443 Lu490 Lu510 Lu555 Lu670 Lu683
% 35          - MCP at 30 m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Dat = coallate(PathS)
%=========================================
% Loop through each directory ep1 and ep2
%=========================================

%===============================================
% Read in the files needed for processing
%===============================================
% Noon Files ===================================
S00a       = readeqpac_([PathS 'spec.noon']);
[t g]      = unique(S00a(:,1));
S00a       = S00a(g,:);
[r l00a]   = size(S00a);

S20a       = readeqpac_([PathS 'spec.20m']);
[t g]      = unique(S20a(:,1));
S20a       = S20a(g,:);
[r l20a]   = size(S20a);

[M10 M30]  = readmcp([PathS 'mcp']);
[t g]      = unique(M10(:,1));
M10        = M10(g,:);
[r lM10]    = size(M10);

[t g]      = unique(M30(:,1));
M30        = M30(g,:);
[r lM30]    = size(M30);

% 12:30 Files ==================================
S00b       = readeqpac_([PathS 'spec.0m.1230']);
[t g]      = unique(S00b(:,1));                 % Get rid of repetitive dates
S00b       = S00b(g,:);
[r l00b]   = size(S00b);

S20b       = readeqpac_([PathS 'spec.20m.1230']); 
[t g]      = unique(S20b(:,1));
S20b       = S20b(g,:);
[r l20b]    = size(S20b);

% 12:30 and Noon Files =========================
Misc       = readeqpac_([PathS 'spec.misc']);
P          = Misc(:,[1 5 9]);       % DecimalYear PressureNoon Pressure1230

CTD        = readeqpac_([PathS 'temp.0m']);
CTD(:,3)   = salin_(CTD(:,3)/4.2914,CTD(:,2),0);  

%=============================================
% correct the times and convert to local time
%+============================================
StartDate  = min([S00a(1,1) S20a(1,1) M10(1,1) M30(1,1)]);                             % Earliest Date
EndDate    = max([max(S00a(end,1)) max(S20a(end,1)) max(M10(end,1)) max(M30(end,1))]); % Latest Date
nextDate   = StartDate;                                                                % Start here
n = 1;                                                                                 % Increment

%===================================================
% Loop through every single date found in the files
%===================================================
while ~isempty(nextDate) 
   
   % Find the maximum occurence of nextDate in each file. Because files occasionally have two
   % lines with the same date but different data, it was easiest to grab only one of the lines 
   % and ignore the other(s).
   r00a   = max(find(S00a(:,1) == nextDate));
   r00b   = max(find(S00b(:,1) == nextDate));
   rM10   = max(find(M10(:,1)  == nextDate));
   rM30   = max(find(M30(:,1)  == nextDate));
   r20a   = max(find(S20a(:,1) == nextDate));
   r20b   = max(find(S20b(:,1) == nextDate));
   rP     = max(find(P(:,1)    == nextDate));
   rCTD   = max(find(CTD(:,1)  == nextDate));
   
   %==========================================================================================
   % Write the data to a buffer variable. If no data is found use -999 as a missing data code
   %==========================================================================================
   if ~isempty(rP)
      P1230 = P(rP, 3);
      X = [1200 P(rP,2)];
   else
      X = [1200 NaN];
      if exist('P1230') | isempty(P1230)
         P1230 = NaN;
      end
   end   
   
   if ~isempty(r00a)
      X = [X S00a(r00a,1)];
   else
      X = [X nextDate];
   end
   
   if ~isempty(rCTD)
      X = [X CTD(rCTD,2:end)];
   else
      X = [X NaN NaN];
   end
   
   % Noon
   if ~isempty(r00a)
      X = [X S00a(r00a,2:end-1)];
   else
      X = [X ones(1,l00a-2)*NaN];
   end
   
   if ~isempty(rM10)
      X = [X M10(rM10,3:end)];
   else
      X = [X ones(1,lM10-2)*NaN];
   end
   
   if ~isempty(r20a)
      X = [X S20a(r20a,2:end)];
   else
      X = [X ones(1,l20a-1)*NaN];
   end
   
   if ~isempty(rM30)
      X = [X M30(rM30,3:end)];
   else
      X = [X ones(1,lM30-2)*NaN];
   end
   
   % 12:30am  
   if ~isempty(rP)
      Y = [1230 P(rP,3)];
   else
      Y = [1230 P1230];
   end
   
   
   %if isempty(rP)
   %   Y = [1230 lP];
   %else
   %   Y = [1230 NaN];
   %end
   
   
   if ~isempty(r00b)
      Y = [Y S00b(r00b,1) NaN NaN S00b(r00b,2:end-1)];
   else
      Y = [Y nextDate NaN NaN ones(1,l00b-2)*NaN];
   end
   
   if ~isempty(r20b)
      Y = [Y ones(1,lM10-2)*NaN S20b(r20b,2:end) ones(1,lM30-2)*NaN];
   else
      Y = [Y ones(1,lM10-2)*NaN ones(1,l20b-1)*NaN ones(1,lM30-2)*NaN];
   end
   
   %====================================
   % Increment to the next largest date
   %====================================
   i = min(find(S00a(:,1)>StartDate));
   if ~isempty(i);
      nextDate = S00a(i,1); 
   else 
      nextDate = []; 
   end
   
   i = min(find(S00b(:,1)>StartDate));
   if ~isempty(i); 
      nextDate = min([S00b(i,1) nextDate]);
   end
   
   i = min(find(S20a(:,1)>StartDate));
   if ~isempty(i); 
      nextDate = min([S20a(i,1) nextDate]);
   end
   
   i = min(find(S20b(:,1)>StartDate));
   if ~isempty(i); 
      nextDate = min([S20b(i,1) nextDate]); 
   end
   
   i = min(find(M10(:,1)>StartDate));
   if ~isempty(i); 
      nextDate = min([M10(i,1) nextDate]); 
   end
   
   i = min(find(M30(:,1)>StartDate));
   if ~isempty(i);
      nextDate = min([M30(i,1) nextDate]);
   end
   
   %==========================================================
   % Write data from the buffers to the output variable, Dat.
   %==========================================================
   Dat(n,:) = X;      
   n        = n + 1;
   Dat(n,:) = Y;
   n        = n + 1;
   
   StartDate = nextDate;
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EDITEP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function epx = editep(epx);
i = find(epx(:,2) > 30 | epx(:,2) < 10);
if ~isempty(i)
   epx(i,6:end) = NaN;
end

i = find(epx(:,5) > 38 | epx(:,5) < 32);
if ~isempty(i)
   epx(i,6:end) = NaN;
end

i = find(epx(:,4) > 35 | epx(:,4) < 15);
if ~isempty(i)
   epx(i,6:end) = NaN;
end

i = find(isinf(epx));
if ~isempty(i)
   epx(i) = NaN;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READMCP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x10, x30] = readmcp(FileName)
x   = readeqpac_(FileName);
x   = sortrows(x,[2 1]);
i   = min(find(x(:,2) == 30));
x10 = x(1:(i - 1),:);
x30 = x(i:end,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READEQPAC_
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = readeqpac_(FileLoc, NumOfDays)
fid = fopen(FileLoc,'r');
OK  = 1;
while OK		 						    	% Trash the header
   S = fgetl(fid);
   if ~strcmp(S(1),'#')
      break
   end
end
i = 1;
while ~feof(fid)					     	% Read the data
   XS     = fgetl(fid);			    	% Get a line from the file
   X(i,:) = sscanf(XS,'%f')';		   % Read the data from the line
   i      = i + 1;
end
fclose(fid);
X = unique(X,'rows');               % Remove duplicate rows