%This is the driver program to invert Rrs into IOP + uncertainties in inversion products
%This an implemention of the strategy suggeted 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 seawifs_inversion.mat to combine with modis_inversion_5wl

% saves seawifs_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_seawifs
     % bb_model_95_seawifs
     % bb_model_5_seawifs
     % bbp_model_median_seawifs
     % bbp_model_95_seawifs
     % bbp_model_5_seawifs

   % seawifs_inverted_chl, containing:
     % adg_model_median_seawifs
     % adg_model_95_seawifs
     % adg_model_5_seawifs
     % aph_model_median_seawifs
     % aph_model_95_seawifs
     % aph_model_5_seawifs
     % a_model_median_seawifs
     % a_model_95_seawifs
     % a_model_5_seawifs

%Emmanuel Boss, U Maine, 2006

clear all
close all

%Read input data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\SeaWIFS\SeaWIFS_spectra_7_5km.mat';
Lwn=spectra_mean;
F0=[173.004, 190.154, 196.473, 188.158, 183.010];
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_SeaWIFS=datenum_seawifs(II);
wavelength=[412   443   490   510   555];

%APEX 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);
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_SeaWIFS,'linear');
Salt=interp1(tt,S,time_SeaWIFS,'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. 
conv_criteria=0.1;

%paramters 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
    [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);
    
    %calculate the statistics of all the possible solutions found:
    if N>0
    adg_model_median_seawifs(i,:)=median(adg_model);
    adg_model_95_seawifs(i,:)=prctile(adg_model,95);
    adg_model_5_seawifs(i,:)=prctile(adg_model,5);
    aph_model_median_seawifs(i,:)=median(aph_model);
    aph_model_95_seawifs(i,:)=prctile(aph_model,95);
    aph_model_5_seawifs(i,:)=prctile(aph_model,5);
    a_model_median_seawifs(i,:)=median(a_model);
    a_model_95_seawifs(i,:)=prctile(a_model,95);
    a_model_5_seawifs(i,:)=prctile(a_model,5);
    bb_model_median_seawifs(i,:)=median(b_model);
    bb_model_95_seawifs(i,:)=prctile(b_model,95);
    bb_model_5_seawifs(i,:)=prctile(b_model,5);
    bbp_model_median_seawifs(i,:)=median(bbp_model);
    bbp_model_95_seawifs(i,:)=prctile(bbp_model,95);
    bbp_model_5_seawifs(i,:)=prctile(bbp_model,5);
    nn_seawifs(i)=length(N);
    else     %flag cases where no solution was found
    adg_model_median_seawifs(i,:)=-[10 10 10 10 10];
    adg_model_95_seawifs(i,:)=-[10 10 10 10 10];
    adg_model_5_seawifs(i,:)=-[10 10 10 10 10];
    aph_model_median_seawifs(i,:)=-[10 10 10 10 10];
    aph_model_95_seawifs(i,:)=-[10 10 10 10 10];
    aph_model_5_seawifs(i,:)=-[10 10 10 10 10];
    a_model_median_seawifs(i,:)=-[10 10 10 10 10];
    a_model_95_seawifs(i,:)=-[10 10 10 10 10];
    a_model_5_seawifs(i,:)=-[10 10 10 10 10];
    bb_model_median_seawifs(i,:)=-[10 10 10 10 10];
    bb_model_95_seawifs(i,:)=-[10 10 10 10 10];
    bb_model_5_seawifs(i,:)=-[10 10 10 10 10];
    bbp_model_median_seawifs(i,:)=-[10 10 10 10 10];
    bbp_model_95_seawifs(i,:)=-[10 10 10 10 10];
    bbp_model_5_seawifs(i,:)=-[10 10 10 10 10];
    nn_seawifs(i)=length(N);
    end
end

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

figure(1);
plot(time_float,chl,'*',time_SeaWIFS(I),aph_model_median_seawifs(I,2)./0.04,'r*');
datetick('x','mmmyy');
xlabel('APEX-chlorophyll');
ylabel('Inversion-chlorophyll');

figure(2);
plot(tt,bb,'b',time_SeaWIFS(I),bbp_model_median_seawifs(I,2),'r*',time_SeaWIFS(I),bbp_model_95_seawifs(I,2),'r+',time_SeaWIFS(I),bbp_model_5_seawifs(I,2),'r+');
xlabel('APEX-bb');
datetick('x','mmmyy');
ylabel('Inversion-bb');
title('SeaWIFS and Float 0005 bb');

A=bbp_model_median_seawifs(I,2)';
bb_apex_SeaWIFS=interp1(tt,bb,time_SeaWIFS(I),'linear');
II=find(isfinite(bb_apex_SeaWIFS) & A<0.01);

figure(3);
plot(bb_apex_SeaWIFS(II),A(II),'*');
corrcoef(bb_apex_SeaWIFS(II),A(II));
polyfit(bb_apex_SeaWIFS(II),A(II),1);
title('Correlation of Apex and SeaWIFS bb');

B=aph_model_median_seawifs(I,2)./0.022;
chl_apex_SeaWIFS=interp1(time_float,chl,time_SeaWIFS(I),'linear');
II=find(isfinite(chl_apex_SeaWIFS) & abs(B')<12);

figure(4);
plot(chl_apex_SeaWIFS(II),B(II),'*');
corrcoef(chl_apex_SeaWIFS(II),B(II));
polyfit(chl_apex_SeaWIFS(II),B(II)',1);
title('Correlation of Apex and SeaWIFS chl');

figure(5);
plot(Salt(I),adg_model_median_seawifs(I,2),'o');

save seawifs_inversion aph_model_5_seawifs aph_model_95_seawifs time_SeaWIFS chl_apex_SeaWIFS nn_seawifs bb_apex_SeaWIFS bbp_model_5_seawifs bbp_model_median_seawifs aph_model_median_seawifs bbp_model_95_seawifs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% %for cases where we have one or no solutions we use Lee's algorithm:
% I=find(nn==1);
% rrs=Rrs(I,:);
% time=time_SeaWIFS(I);
% Rrs_=(0.52.*rrs)./(1-1.7.*rrs);
% T=Temp(I);
% S=Salt(I);
% clear a b
% for i=1:length(I)
%     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(2);
% hold on
% plot(time,bbp_440,'g*');
% hold off

% figure(1)
% subplot(2,2,1)
% loglog(a_pg(I,2),a_model_median(I,2),'o');
% hold on
% for i=1:length(I)
%     plot([a_pg(I(i),2),a_pg(I(i),2)],[a_model_5(I(i),2),a_model_95(I(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(I,2),aph_model_median(I,2),'o');
% hold on
% for i=1:length(I)
%     plot([a_phi(I(i),2),a_phi(I(i),2)],[aph_model_5(I(i),2),aph_model_95(I(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(I,2),adg_model_median(I,2),'o');
% hold on
% for i=1:length(I)
%     plot([a_dg(I(i),2),a_dg(I(i),2)],[adg_model_5(I(i),2),adg_model_95(I(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(I,2),a_pg(I,2),10);
% [N, intercept_t1, intercept_t2, slope_t1, slope_t2, rsq, rms1, bias]=reg_stats(aph_model_median(I,2),a_phi(I,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(I,2),a_dg(I,2),10);

seawifs_inverted_bb=[bb_model_median_seawifs bb_model_95_seawifs bb_model_5_seawifs bbp_model_median_seawifs bbp_model_95_seawifs bbp_model_5_seawifs];
seawifs_inverted_chl=[adg_model_median_seawifs adg_model_95_seawifs adg_model_5_seawifs aph_model_median_seawifs aph_model_95_seawifs aph_model_5_seawifs a_model_median_seawifs a_model_95_seawifs a_model_5_seawifs];

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

% % clear unecessary variables
% clear A B F0 I I2 II L3 L4 Lwn N Rrs S Salt T Temp V a_model a_model_5_seawifs a_model_95_seawifs a_model_median_seawifs;
% clear a_sea_water adg_model adg_model_5_seawifs adg_model_95_seawifs adg_model_median_seawifs ans apg_model aph_model;
% clear aph_model_5_seawifs aph_model_95_seawifs aph_model_median_seawifs b_model bb bb_apex_SeaWIFS bb_model_5_seawifs;
% clear bb_model_95_seawifs bb_model_median_seawifs bb_sea_water bbp_model bbp_model_5_seawifs bbp_model_95_seawifs;
% clear bbp_model_median_seawifs chl chl_apex_SeaWIFS chl_surface conv_criteria datenum_seawifs float0005 float0005SurfaceChl;
% clear h i mean_chl median_chl nn_seawifs spectra_mean spectra_median spectra_std std_chl time_SeaWIFS time_float time_surface;
% clear tt wavelength wl;
% 
% save 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\APEX\dataFiles\seawifs_inverted'



    
