%This is the driver program to invert Rrs into IOP + uncertainties in inversion products
%This an implemention of the strategy suggested in the Wang, Boss and Roesler, 2005, Applied Optics paper.
%Here it is used to invert in MODIS data and compare to data from the APEX float

%calls the following m.files:
   %IOP_of_sw.m, IOP_inversion.m, v.m, h_Array.m, and reg_stats.m

% saves modis_inverted.mat to the dataFiles folder containing these two arrays (need to uncomment the clear and save statements to run):
   % seawifs_inverted_bb, containing:
     % bb_model_median
     % bb_model_95
     % bb_model_5
     % bbp_model_median
     % bbp_model_95
     % bbp_model_5

   % modis_inverted_chl, containing:
     % adg_model_median
     % adg_model_95
     % adg_model_5
     % aph_model_median
     % aph_model_95
     % aph_model_5
     % a_model_median
     % a_model_95
     % a_model_5
     
%Emmanuel Boss, U Maine, 2006

clear all
close all

%Read input data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\MODIS\MODIS_spectra_7_5km.mat';
Lwn=spectra_median;
F0=[172.718,187.688,194.973,185.820,186.594];
mm=length(Lwn(:,1));
Rrs=1.38*Lwn(:,1:5)./(ones(mm,1)*F0);
II=find(isfinite(Rrs(:,1).*Rrs(:,2).*Rrs(:,3).*Rrs(:,4).*Rrs(:,5)));
Rrs=Rrs(II,:);
time_MODIS=datenum_modis(II);
wavelength=[412	443	488	531	551];

%ARGOS data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\ARGOS\0005\float0005Surface.mat';
chl_surface(1)=float0005SurfaceChl(1,2);
time_surface(1)=float0005SurfaceChl(1,1);
I=find(diff(float0005SurfaceChl(:,1))>4); % get first reading of the day
for i=1:length(I)
    chl_surface(1+i)=float0005SurfaceChl(1+I(i),2);
    time_surface(1+i)=float0005SurfaceChl(1+I(i),1);
end
time_float=time_surface;
chl=chl_surface;

load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\APEX\0005\float0005.mat';
I2=find(float0005(:,5)>1 & float0005(:,5)<9);
S=float0005(I2,7);
T=float0005(I2,6);
bb=float0005(I2,9);
chl=chl_surface;
tt=tt(I2);

%user inputs:
%Temperature and salinity of water under consideration [degree C] (to compute water IOPs)
Temp=interp1(tt,T,time_MODIS,'linear');
Salt=interp1(tt,S,time_MODIS,'linear');

%compute IOPs of seawater
%[a_sea_water,bb_sea_water]=IOP_of_sw(wavelength,Temp,Salt);

%convergence criteria for selection of solutions. e.g. if 0.1 only solution
%that are within 10% of Rrs at all wavelengths will be considered.
%can try different values here 0.1-0.3
conv_criteria=0.1;

%parameters of model relating Rrs to IOP (here based on GSM)
L3 = 0.0949;
L4 = 0.0794;

%Loop where all the spectra are inverted
for i=1:length(II)
    %compute IOPs of seawater
    [a_sea_water,bb_sea_water]=IOP_of_sw(wavelength,Temp(i),Salt(i));
  
    %Setting up the linear solver:
    i
    [V] = v(L3, L4, Rrs(i,:)');
    [h] = Array_h(a_sea_water, bb_sea_water, V);

    %proceed with inversion (calls IOP_inversion.m)
   [N, adg_model, aph_model, apg_model, a_model, bbp_model, b_model]= IOP_inversion(L3, L4, Rrs(i,:)', V, h, a_sea_water, bb_sea_water, wavelength, conv_criteria);
   % [N, adg_model, aph_model, apg_model, a_model, bbp_model, b_model, eigen_adg_model, eigen_bbp_model]= IOP_inversion(L3, L4, Rrs(i,:)', V, h, a_sea_water, bb_sea_water, wavelength, conv_criteria);
    
    %calculate the statistics of all the possible solutions found:
    if N>0
    adg_model_median(i,:)=median(adg_model);
    adg_model_95(i,:)=prctile(adg_model,95);
    adg_model_5(i,:)=prctile(adg_model,5);
    aph_model_median(i,:)=median(aph_model);
    aph_model_95(i,:)=prctile(aph_model,95);
    aph_model_5(i,:)=prctile(aph_model,5);
    a_model_median(i,:)=median(a_model);
    a_model_95(i,:)=prctile(a_model,95);
    a_model_5(i,:)=prctile(a_model,5);
    bb_model_median(i,:)=median(b_model);
    bb_model_95(i,:)=prctile(b_model,95);
    bb_model_5(i,:)=prctile(b_model,5);
    bbp_model_median(i,:)=median(bbp_model);
    bbp_model_95(i,:)=prctile(bbp_model,95);
    bbp_model_5(i,:)=prctile(bbp_model,5);
%     eigen_adg_model_median(i,:)=median(eigen_adg_model);
%     eigen_adg_model_95(i,:)=prctile(eigen_adg_model,95);
%     eigen_adg_model_5(i,:)=prctile(eigen_adg_model,5);
%     eigen_bbp_model_median(i,:)=median(eigen_bbp_model);
%     eigen_bbp_model_95(i,:)=prctile(eigen_bbp_model,95);
%     eigen_bbp_model_5(i,:)=prctile(eigen_bbp_model,5);
    nn(i)=length(N);
    else     %flag cases where no solution was found
    adg_model_median(i,:)=-[10 10 10 10 10];
    adg_model_95(i,:)=-[10 10 10 10 10];
    adg_model_5(i,:)=-[10 10 10 10 10];
    aph_model_median(i,:)=-[10 10 10 10 10];
    aph_model_95(i,:)=-[10 10 10 10 10];
    aph_model_5(i,:)=-[10 10 10 10 10];
    a_model_median(i,:)=-[10 10 10 10 10];
    a_model_95(i,:)=-[10 10 10 10 10];
    a_model_5(i,:)=-[10 10 10 10 10];
    bb_model_median(i,:)=-[10 10 10 10 10];
    bb_model_95(i,:)=-[10 10 10 10 10];
    bb_model_5(i,:)=-[10 10 10 10 10];
    bbp_model_median(i,:)=-[10 10 10 10 10];
    bbp_model_95(i,:)=-[10 10 10 10 10];
    bbp_model_5(i,:)=-[10 10 10 10 10];
%     eigen_adg_model_median(i,:)=-[10];
%     eigen_adg_model_95(i,:)=-[10];
%     eigen_adg_model_5(i,:)=-[10];
%     eigen_bbp_model_median(i,:)=-[10];
%     eigen_bbp_model_95(i,:)=-[10];
%     eigen_bbp_model_5(i,:)=-[10];
    nn(i)=length(N);
    end
end

% load seawifs inversion data
load seawifs_inversion

I=find(nn'>1 & bbp_model_median(:,1)>0);
IS=find(nn_seawifs>1 & bbp_model_median_seawifs(:,1)'>0);

figure(1);% inversion chl
plot(time_float,chl,'b.',time_MODIS(I),aph_model_median(I,2)./0.04,'g*');
hold on
plot(time_float,chl,'b.',time_SeaWIFS(IS),aph_model_median_seawifs(IS,2)./0.04,'r*');
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',16);
set(gca,'XTick',[732099 732464 732829 733194 733560]);
datetick(gca,'x','mmmyy','keepticks');
xlabel('Float Chlorophyll');
ylabel('MODIS/SeaWIFS Chlorophyll');
legend('Float0005','MODIS', 'Float0005', 'SeaWIFS');
% title('MODIS/Seawifs and Float 0005 chl');
hold off

figure(2);
plot(tt,bb,'b.')
hold on
plot(time_MODIS(I),bbp_model_median(I,2),'go')
plot(time_SeaWIFS(IS),bbp_model_median_seawifs(IS,2),'ro')
for i=1:length(I)
plot([time_MODIS(I(i)) time_MODIS(I(i))],[bbp_model_95(I(i),2) bbp_model_5(I(i),2)],'g-')
end
for i=1:length(IS)
plot([time_SeaWIFS(IS(i)) time_SeaWIFS(IS(i))],[bbp_model_95_seawifs(IS(i),2) bbp_model_5_seawifs(IS(i),2)],'r-')
end
plot(tt,bb,'b.')
hold off
set(gca,'xlim',[datenum('01 jun 2004') datenum('01 jun 2007')]);
set(gca,'XTick',[datenum('01 jun 2004') datenum('01 jun 2005') datenum('01 jun 2006') datenum('01 jun 2007')]);
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',16);
datetick('x', 'mmmyy','keeplimits','keepticks');
% xlabel('APEX-bb');
ylabel('b_b(440)');
% legend('Float0005','MODIS','Float0005', 'SeaWIFS');
% title('MODIS/Seawifs and Float 0005 bb');
hold off

A=bbp_model_median(I,2)';
bb_apex_MODIS=interp1(tt,bb,time_MODIS(I),'linear');
II=find(isfinite(bb_apex_MODIS) & A'<0.01);

figure(3);
plot(bb_apex_MODIS(II),A(II),'r*');
hold on
plot([0 0.008],[0 0.008],'k-')
hold off
corrcoef(bb_apex_MODIS(II),A(II));
polyfit(bb_apex_MODIS(II),A(II)',1);
title('Correlation of Apex and MODIS bb');

B=aph_model_median(I,2)./0.04;
chl_apex_MODIS=interp1(time_float,chl,time_MODIS(I),'linear');
II=find(isfinite(chl_apex_MODIS) & abs(B)<12);

figure(4);
plot(chl_apex_MODIS(II),B(II),'r*');
hold on
plot([0 8],[0 8],'k-')
hold off
% corrcoef(chl_apex_MODIS(II),B(II));
polyfit(chl_apex_MODIS(II),B(II),1);
title('Correlation of Apex and MODIS chl');

figure(6);
plot(Salt(I),adg_model_median(I,2),'o');
title('Correlation of Salinity and ADG');

% figure(6);plot(time_MODIS(I),Salt(I));
% datetick('x','mmmyy');
% ylabel('Salinity');
% 
% figure(7);plot(time_MODIS(I),adg_model_median(I,2));
% datetick('x','mmmyy');
% ylabel('ADG');

figure(7);
[AX,H1,H2] = plotyy(time_MODIS(I),adg_model_median(I,2),time_MODIS(I),Salt(I));
set(get(AX(1),'Ylabel'),'String','ADG')
set(get(AX(2),'Ylabel'),'String','Salinity')

set(AX(1),'xlim',[datenum('01 jun 2004') datenum('01 jun 2007')])
set(AX(2),'xlim',[datenum('01 jun 2004') datenum('01 jun 2007')])
set(AX(1),'XTick',[732099 732464 732829 733194])
set(AX(2),'XTick',[])
datetick(AX(1),'x', 'mmmyy','keeplimits','keepticks')
title({'Salinity and ADG Through Time'})

% size(II);
%0.1 has 175 points
%0.2 has 189 points
%0.3 has 195 points

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% %for cases where we have one or no solutions we use Lee's algorithm:
% III=find(nn==1);
% rrs=Rrs(III,:);
% time=time_MODIS(III);
% Rrs_=(0.52.*rrs)./(1-1.7.*rrs);
% T=Temp(III);
% S=Salt(III);
% clear a b
% for i=1:length(III)
%     rho = log10(max([Rrs_(i,2), Rrs_(i,3), Rrs_(i,4)],[],2)./Rrs_(i,5));
%     Kd(i) = 0.0605+10.^(-1.163-1.969*rho+1.239*rho.^2+0.417*rho.^3-0.984*rho.^4);
%     a(i) = (0.9*Kd(i).*(1-6.8*Rrs_(i,5)))./(1+15.3*Rrs_(i,5));
%     Y(i) = 2.2*(1-1.2*exp(-0.9.*(rrs(i,2)./(rrs(i,5)))));
%     u = (-0.0895+sqrt(0.008+0.499.*rrs(i,5)))./0.249;
%     bb_(i) = (a(i)*u)/(1-u);
% end
% [a_sea_water,bb_sea_water]=IOP_of_sw(555,T,S);
% 
% bbp = bb_ - bb_sea_water';
% bbp_440=bbp.*(551/440).^Y;
% figure(7);
% hold on
% plot(time,bbp_440,'g*');
% datetick('x','mmmyy');
% hold off

% figure(1)
% subplot(2,2,1)
% loglog(a_pg(III,2),a_model_median(III,2),'o');
% hold on
% for i=1:length(III)
%     plot([a_pg(III(i),2),a_pg(III(i),2)],[a_model_5(III(i),2),a_model_95(III(i),2)],'b-')
% end
% hold on
% plot([0.01 5],[0.01 5],'k-')
% axis([0.01 5 0.01 5])
% xlabel('Known a_{pg}(440)')
% ylabel('derived a_{pg}(440)')
% hold off
% 
% subplot(2,2,2)
% loglog(a_phi(III,2),aph_model_median(III,2),'o');
% hold on
% for i=1:length(III)
%     plot([a_phi(III(i),2),a_phi(III(i),2)],[aph_model_5(III(i),2),aph_model_95(III(i),2)],'b-')
% end
% hold on
% plot([0.003 1],[0.003 1],'k-')
% axis([0.003 1 0.003 1])
% xlabel('Known a_{ph}(440)')
% ylabel('derived a_{ph}(440)')
% hold off
% 
% subplot(2,2,3)
% loglog(a_dg(III,2),adg_model_median(III,2),'o');
% hold on
% for i=1:length(III)
%     plot([a_dg(III(i),2),a_dg(III(i),2)],[adg_model_5(III(i),2),adg_model_95(III(i),2)],'b-')
% end
% hold on
% plot([0.0005 3],[0.0005 3],'k-')
% axis([0.0005 3 0.0005 3])
% xlabel('Known a_{dg}(440)')
% ylabel('derived a_{dg}(440)')
% hold off
% 
% 
% %calculate the parametric and nonparametric statistics of all the median solution:
% [N, intercept_t1, intercept_t2, slope_t1, slope_t2, rsq, rms1, bias]=reg_stats(a_model_median(III,2),a_pg(III,2),10);
% [N, intercept_t1, intercept_t2, slope_t1, slope_t2, rsq, rms1, bias]=reg_stats(aph_model_median(III,2),a_phi(III,2),10);
% I=find(adg_model_median(:,2)>0 & nn(:)>1);
% [N, intercept_t1, intercept_t2, slope_t1, slope_t2, rsq, rms1, bias]=reg_stats(adg_model_median(III,2),a_dg(III,2),10);

modis_inverted_bb=[bb_model_median bb_model_95 bb_model_5 bbp_model_median bbp_model_95 bbp_model_5];
modis_inverted_chl=[adg_model_median adg_model_95 adg_model_5 aph_model_median aph_model_95 aph_model_5 a_model_median a_model_95 a_model_5];

% save data to the dataFiles folder (uncomment the clear and save statements first)

% % clear unecessary variables
% clear A AX B F0 H1 H2 I I2 II IS L3 L4 Lwn N Rrs S Salt T Temp V a_model a_model_5 a_model_95 a_model_median;
% clear a_sea_water adg_model adg_model_5 adg_model_95 adg_model_median ans apg_model aph_model;
% clear aph_model_5 aph_model_5_seawifs aph_model_95 aph_model_95_seawifs aph_model_median aph_model_median_seawifs b_model bb;
% clear bb_apex_MODIS bb_apex_SeaWIFS bb_model_5 bb_model_95 bb_model_median bb_sea_water bbp_model bbp_model_5 bbp_model_5_seawifs;
% clear bbp_model_95 bbp_model_95_seawifs bbp_model_median bbp_model_median_seawifs chl chl_apex_MODIS chl_apex_SeaWIFS chl_surface;
% clear conv_criteria datenum_modis float0005 float0005SurfaceChl h i mean_chl median_chl nn nn_seawifs spectra_mean spectra_median;
% clear spectra_std std_chl time_MODIS time_SeaWIFS time_float time_surface tt wavelength wl;
% 
% save 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\APEX\dataFiles\modis_inverted'


