%Process TRIOS data version 4: Read in Dave Bowers radiometer data files
%(.txt files). Calculates KdPAR, spectral Kd (Kd_spectral).

% Need to get details about processing for Reflectance Reflectance PAR (R_PAR) and spectral
% Reflectance (R_spectral) from Dave.
%
%Relevent background/definitions:
%
%Kd=Vertical attenation coefficient for downward irradiance. Kd=the linear
%regression coefficient (slope) of the Ln(Ed) with respect to depth. Most
%popular parameter to quantify the amount of light available for
%photosynthesis.
%
%Reflectance=ratio of upward irradiance(Eu) to downward irradiance (Ed) just below
%the surface. It is generally converted to a percentage and important for
%remote sensing. In natural waters R has been found to be proportional to
%bb/(a+bb), where bb is backscatter and a is absorption.

%Upadated to 'linearleastsquares_kaf.m" to calculate Kd
%and the standard errror, and r-square fit metric.

%KAF July 2017

%Inputs:
%Trioslog.xlsx files created using the log trios log sheets
%Trios.txt data files

%Output-saved as a .mat file
%SurveyTrios.mat: Matlab workspace containing the following variables:
%Survey: Survey name
%Stations: character vector of station labels
%Station_spectral: 2 column vector
%Column 1: Kd from 300 to 900 nm
%Column 2: Standard error of spectral Kd
%wv: Vector of TRIOS wavelengths for plotting
%KdPAR: Vector of calculated attenuation of PAR (400-700 nm)for each station
%SEPAR: Vector standard error of KdPAR for each station

%Also, creates spectral attenuation figure from 400-700 nm.



%% Import the log sheet
pnin='/Users/kelsey/Desktop/Dissertation/Profiler Surveys/York River/Axial Surveys/YR160902/TRIOS/';
Survey='YR160902';

[data,~, ~] = xlsread(fullfile(pnin,[Survey,'TRIOSlog.xlsx']),'Sheet1');

stations=unique(data(:,1));

depth_Corr=0.285;

%%
for ii=1:length(stations)
    
    %ii=6;
    name=['S',num2str(stations(ii))]
    in_Ed=find(data(:,1)==stations(ii)& data(:,4)==1); %use index for depth with downwelling light
    in_Eu=find(data(:,1)==stations(ii) & data(:,4)==2); %use index for depth with upwelling light
    
    depth_Ed=data(in_Ed,5)-depth_Corr;
    depth_Eu=data(in_Eu,5)+depth_Corr;
  
    list_Ed=dir([name, 'd*']);
    w_in_Ed=1:2:length(list_Ed); %water index
    d_in_Ed=2:2:length(list_Ed); %deck index (named with a 0 after the number)
    
    list_Eu=dir([name, 'u*']);
    w_in_Eu=1:2:length(list_Eu); %water index
    d_in_Eu=2:2:length(list_Eu); %deck index (named with a 0 after the number)
    
    
    
    for kk=1:length(w_in_Ed)
        w_tmp2=dlmread(num2str([list_Ed(w_in_Ed(kk)).name]),'',6,0);
        %w_PARindex=find(w_tmp(:,1)>=400 &  w_tmp(:,1)<=700)
        %w_Allindex=find(w_tmp(:,1)>=300 &  w_tmp(:,1)<=900);
        %w_tmpPAR=w_tmp(w_PARindex,:);
        % w_tmp2=w_tmp(w_Allindex,:);
        
        d_tmp2=dlmread(num2str([list_Ed(d_in_Ed(kk)).name]),'',6,0);
        %d_PARindex=find(d_tmp(:,1)>=400 &  d_tmp(:,1)<=700)
        %d_tmp2=d_tmp(w_Allindex,:);
        
        water_tmp(kk,:)=w_tmp2(:,2)';
        deck_tmp(kk,:)=d_tmp2(1:184,2)';
    end
    
    w_wv=w_tmp2(:,1)';
    d_wv=d_tmp2(:,1)';
    
    clear w_tmp2 d_tmp2 d_tmp w_tmp w_PARindex w_index d_index kk
    
    wvbins=[325:5:900];
    [~,~,bin]=histcounts(w_wv,wvbins);
    for i=1:size(water_tmp,1)
        for jj=1:length(midpoint(wvbins))
            inbin=bin==jj;
            water(i,jj)=nanmean(water_tmp(i,inbin));
        end
    end
    clear bin
    
    [~,~,bin]=histcounts(d_wv,wvbins);
    for i=1:size(water_tmp,1)
        for jj=1:length(midpoint(wvbins))
            inbin=bin==jj;
            deck(i,jj)=nanmean(deck_tmp(i,inbin));
        end
    end
    clear bin deck_tmp water_tmp
    
    
    for bb=1:size(water,1)
        % for jj=1:size(water,2)-1
        % Water2(bb,jj)=(water(bb,jj).*(wvbins(jj+1)-wvbins(jj)))./1000;
        % Deck2(bb,jj)=(deck(bb,jj).*(wvbins(jj+1)-wvbins(jj)))./1000;
        Water(bb,:)=(water(bb,:).*(5))./1000;
        Deck(bb,:)=(deck(bb,:).*(5))./1000;
        %end
        % end
    end
    
    clear bb jj
    
    for tt=1:size(Water,1)
        Water(tt,Water(tt,:)<0.0015)=0;
    end
    
    
    w_PARindex=find(wvbins>=400 &  wvbins<=700);
    EdPAR=sum(Water(:,w_PARindex),2); %For KdPar
    DeckPAR=sum(Deck(:,w_PARindex),2);
     
% Some cases may not have upwelling observations   
if ~isempty(in_Eu)
    for kk=1:length(w_in_Eu) %only need water for this because I am just comparing this to Ed
        w_tmp2=dlmread(num2str([list_Eu(w_in_Eu(kk)).name]),'',6,0);
        %w_Allindex=find(w_tmp(:,1)>=300 &  w_tmp(:,1)<=900);
        %w_tmp2=w_tmp(w_Allindex,:);
        water_Eu_tmp(kk,:)=w_tmp2(:,2)';
     
        d_tmp2=dlmread(num2str([list_Eu(d_in_Eu(kk)).name]),'',6,0);
        %d_PARindex=find(d_tmp(:,1)>=400 &  d_tmp(:,1)<=700)
        %_tmp2=d_tmp(w_Allindex,:);
        deck_Eu_tmp(kk,:)=d_tmp2(1:184,2)';

    end
         w_wv_u=w_tmp2(:,1)';
          d_wv_u=d_tmp2(:,1)';
          
          [~,~,bin]=histcounts(w_wv_u,wvbins);
          for i=1:size(water_Eu_tmp,1)
              for jj=1:length(midpoint(wvbins))
                  inbin=bin==jj;
                  water_Eu(i,jj)=nanmean(water_Eu_tmp(i,inbin));
              end
          end
          clear bin
          
          [~,~,bin]=histcounts(d_wv_u,wvbins);
          for i=1:size(water_Eu_tmp,1)
              for jj=1:length(midpoint(wvbins))
                  inbin=bin==jj;
                  deck_Eu(i,jj)=nanmean(deck_Eu_tmp(i,inbin));
              end
          end
          
          clear bin
          
    
    
clear w_tmp2 d_tmp2 d_tmp w_tmp  w_index d_index kk

 
   for bb=1:size(water_Eu,1)
      % for jj=1:size(water_Eu,2)-1
      %      Water_Eu(bb,jj)=(water_Eu(bb,jj).*(w_wv(jj+1)-w_wv(jj)))./1000;
      %      Deck_Eu(bb,jj)=(deck_Eu(bb,jj).*(d_wv(jj+1)-d_wv(jj)))./1000;
      % end
      Water_Eu(bb,:)=(water_Eu(bb,:).*(5))./1000;
      Deck_Eu(bb,:)=(deck_Eu(bb,:).*(5))./1000;
   end
  
   clear bb jj
   
    for tt=1:size(Water_Eu,1)
        Water_Eu(tt,Water_Eu(tt,:)<0.0015)=0;
   end
   
       
    EuPAR=sum(Water_Eu(:,w_PARindex),2);
    DeckPAR_Eu=sum(Deck_Eu(:,w_PARindex),2);
    
    
   surf_Eu=interp1(depth_Eu,Water_Eu./Deck_Eu,[0],'linear','extrap');
   surf_Eu(surf_Eu<0)=0;
   surf_Ed=interp1(depth_Ed,Water./Deck,[0],'linear','extrap');
   surf_Ed(surf_Ed<0)=0;
   %poor=find(surf_Ed<=surf_Eu);
   R(ii,:)=((surf_Eu)./(surf_Ed))';
   %R(poor)=0;
   RPAR(ii)=(sum(surf_Eu(w_PARindex))./sum(surf_Ed(w_PARindex)))*100;
   R665(ii)=R(ii,69)*100;

else
   Water_Eu=NaN;
   Deck_Eu=NaN;
    EuPAR=NaN;
    DeckPAR_Eu=NaN;
end

    
 % Determine KdPAR (400-700nm)  
%Make a figure of ln(Ed/Edeck) versus depth and pick depths to be use in
%calculation of Kd (Similar to Licor) 
figure(100);clf
plot(log(EdPAR./DeckPAR),depth_Ed,'bo-')

set(gca,'Ydir','reverse')
set(gca,'Fontsize',16)
ylabel('Depth (m)')
xlabel('ln(E_d/E_d_e_c_k)')
title({'Click on points to include in K_d calculation','Select Points that maintain constant slope'})

[~,d]=ginput;
include=NaN(length(d),1);

for ll=1:length(d)
         [~, include(ll)]=min(abs(depth_Ed-d(ll)));
         %include(ll)=t;
end
        clear ll x d 

close
  
mdl=fitlm(depth_Ed(include),log(EdPAR(include)./DeckPAR(include)),'linear');
KdPAR(ii,1)=abs(mdl.Coefficients.Estimate(2));
seKDPAR(ii)=mdl.Coefficients.SE(2);
clear mdl

%  %[b, Rsqr,se,~,~]=linearleastsquares_kaf(depth_Ed(in_depth),log(EdPAR(in_depth)./DeckPAR(in_depth)));
% [b, Rsqr,se,~,~]=linearleastsquares_kaf(depth_Ed(include),log(EdPAR(include)./DeckPAR(include))); 
%  KdPAR(ii,1)=abs(b(2));
%  SEPAR(ii,1)=abs(se(2));
%  Rsqr_PAR(ii,1)=Rsqr;
%  %R_PAR(ii,1)=((EuPAR(include(1))./DeckPAR_Eu(include(1))))/(EdPAR(include(1))/DeckPAR(include(1)))*100;
% clear b Rsqr se
% 
% Determine Spectral Kd  (300-900 nm, can change later) 
tol=0.0015; %Set a tolerance for irradiance (Per Dave Bowers). 
for n=1:length(Water)
    %X=Water(Water(in_depth,n)>0,n);     deck=Deck(Water(in_depth,n)>0,n);
    %X=Water(Water(include(1:2),n)>0,n);     deck=Deck(Water(include(1:2),n)>0,n);
     X=Water(include,n);     deck1=Deck(include,n);

%if isempty(in_Eu)
%    R(ii,n)=NaN;
%else
%   R(ii,n)=((Water_Eu(include(1),n)/Deck_Eu(include(1),n))/(Water(include(1),n)/Deck(include(1),n)))*100; %spectral reflectance
%end
    if length(X)>=2
      if X(1)<tol || X(2)<tol
        Kd(ii,n)=NaN;
        
      else 
        %[b(:,n), ~,se(:,n),~,~]=linearleastsquares_kaf(depth_Ed(include(1:2)),log(X./deck));
        %Kd(ii,n)=abs(b(2,n));
        mdl=fitlm(depth_Ed(include),log(X./deck1),'linear');
        Kd(ii,n)=abs(mdl.Coefficients.Estimate(2));
        SE_Kd(ii,n)=(mdl.Coefficients.SE(2));

        %SE_Kd(ii,n)
        %SE_Kd(ii,n)=abs(se(2,n));
      end
    else
         Kd(ii,n)=NaN;
        SE_Kd(ii,n)=NaN;
    end
  clear b  se X Y 
end
% 
%     figure(200);clf
%     plot(w_wv(w_PARindex),Kd(ii,w_PARindex),'b-','linewidth',2);hold on
%     %plot(w_wv(w_PARindex),R(ii,w_PARindex),'r','linewidth',2);grid
%     %ylabel('Kd (m^-^1) or Reflectance R (%)')
%     ylabel('Kd (m^-^1)')
%     xlabel('Wavelength (nm)')
%     title([name,' Kd_P_A_R= ',num2str(KdPAR(ii))])
%     set(gca,'FontSize',16) %Set axis font size and limits
%     xlim([400 700])
%     ylim([0 ceil(max(Kd(ii,:)))]) % lower limit.
%     %legend(['K_d ',num2str(depth_Ed(include(1))),'-',num2str(depth_Ed(include(2))),' m'],['R ',num2str(depth_Ed(include(1))),' m']);
%     pause
%    saveas(gcf,[name,'spectral.fig'],'fig')
%    close
% 
include_kd{ii}=include;


 clear include name index Water Deck water deck SE n list kk depth d_wv in_Eu in_Ed d_in_Ed d_in_Eu deck_Eu Deck_Eu DeckPAR_Eu depth_Ed depth_Eu in_depth w_Allindex w_in_Ed w_in_Eu water_Eu Water_Eu DeckPAR EdPAR EuPAR 
clear  deck_Eu_tmp water_Eu_tmp water_Eu deck_Eu Water_Eu Deck_Eu
     

end
 %%
R(isnan(R))=0;
R(R<0)=0;
R(isinf(R))=0;

Kd(isnan(Kd))=0;
Kd(Kd<0)=0;
Kd(isinf(Kd))=NaN;

%Save in .mat file.
save([Survey,'_Trios_june20.mat'])
%%
% clear ii data pnin list_Ed list_Eu w_PARindex 
% %Save in .mat file.
% save([Survey,'_Trios.mat'])
% 
figure (200);clf
%subplot(311)
j=1:length(stations);
plot(w_wv(w_PARindex),Kd(j,w_PARindex),'linewidth',2);hold on
legend(num2str(stations))
ylabel('Kd (m^-^1)')
xlabel('Wavelength (nm)')
set(gca,'FontSize',16) %Set axis font size and limits
title ([Survey,': TRIOS Radiometer']) 
% % 
% % subplot(312)  
% % plot(w_wv(w_PARindex),Kd_db(j,w_PARindex),'linewidth',2);hold on
% % %legend(stations)
% % ylabel('Kd surface (m^-^1)')
% % xlabel('Wavelength (nm)')
% % set(gca,'FontSize',16) %Set axis font size and limits
% % 
% %  
% % subplot(313)
% % plot(w_wv(w_PARindex),R(j,w_PARindex),'linewidth',2);hold on
% % %legend(stations)
% % ylabel('% R surface')
% % xlabel('Wavelength (nm)')
% % set(gca,'FontSize',16) %Set axis font size and limits
% % 
% % 
% % 
% % 
% 
