clear all
close all
Lres=input('enter spectral resolution (i.e. 1, 2, 5, 10, 25) in nm')
T=15;,S=33;%or enter your own values for temp and salinity
wave=400:Lres:700;,wave=wave';

%%%%%%%%%%%%%%%%%%%%%%%%%input spectral reflectance and wavelength data into arrays wave and Rmeas
%%% note the data files are
%%% Rrs_E_... or Rrs_L_... to indicate whether the observations are
%%% irradiance reflectance (Eu/Ed) or radiance reflectance (Lu/Ed). Also
%%% note that all of the inversions assume the the upward light (Lu or Eu)
%%% is measured just below the air water interface
%%% It is necessary to adjust the FQ parameter adjustment ~ 0.33/0.09 to
%%% transform the Eu/Ed to Lu/Ed. Not perfect but will work for comparing
%%% models since they will be compared with the same observations
%%%%

%dat=load('Rrs_E_WestCoast.dat');,FQ=0.09/0.33;%from Roesler and Perry 1995 paper
%dat=load('Rrs_E_DabobBay.dat');,FQ=0.09/0.33;%from Roesler and Perry 1995 paper
%dat=load('Rrs_E_GulfMaine.dat');,FQ=0.09/0.33;%from Roesler and Perry 1995 paper
%dat=load('Rrs_E_PugetSound.dat');,FQ=0.09/0.33;%from Roesler and Perry 1995 paper

dat=load('Rrs_L_HL_simulation.dat');, FQ=1;% from Curt's simulated data set, see the xls file for HL inputs

[nl,nspec]=size(dat);
R=interp1(dat(:,1),dat(:,2:nspec),wave,'linear');
figure(100),plot(wave,R),title('Measured spectra'),ylabel('R(\lambda)'),xlabel('Wavelength(nm)')


%%%%
%%% call RB03_invert model, 
%%% Convert chl to a(443) using a* in model
%%% assumes Rrs beneath water interface and irradiance reflectance,
%%% converts to irradiance to radiance reflectance
%%% R = 0.095 (bb/a+bb) + 0.079 (bb/a+bb)^2
%%% a = aw + chl*aph* + acdom443*exp(-0.018) + anap443*exp(-0.011)
%%% bb = bbw + bbpratio*(cp443(L/Lo)^gamma - (chl*aph* + anap443*exp(-0.011))

wavelength=wave;    
%water basis vectors require temperature T and salinity S, output water IOPs at wave resolution
    [a_w,bb_w]=water_iops_PF_TScorr(wave,T,S);
% average spectrum
    aph=phyto_avg_abs(wave);

%%% perform inversion
iopRB03 = RB03_invert(R'*FQ,wavelength);%converting the measured reflectance by FQ factor

chl          = iopRB03(:,1);%output is chl 
acdom443     = iopRB03(:,2);
anap443      = iopRB03(:,3);
cp443        = iopRB03(:,4);
bbpratio     = iopRB03(:,5);
gamma        = iopRB03(:,6);
Scdom=0.018;
Snap =0.011;
%astar443=0.0222;%for the averaged phytoplankton basis vector
%iopRPFG(:,i)=iopRPFG(:,i).*astar443;%convert retrieved chl to aphi(443) using basis vector a* values;
    
%%% now plot results
for n=1:nspec-1
aphi(:,n) = aph(:,1).*chl(n,1);%total phyto abs
acdom(:,n)= acdom443(n).*exp(-Scdom*(wave-443*ones(size(wave))));
anap(:,n)= anap443(n).*exp(-Snap*(wave-443*ones(size(wave))));
bbp(:,n) = bbpratio(n) .*(cp443(n).* ((443*ones(size(wave)) ./ wave).^gamma(n)) - (aphi(:,n) + anap(:,n)));

atotal(:,n) = a_w + aphi(:,n) + acdom(:,n) + anap(:,n);
 bbtotal(:,n) = bb_w + bbp(:,n);
 
g1 = 0.0949; g2 = 0.0794; 
x = bbtotal(:,n) ./ (atotal(:,n) + bbtotal(:,n));
    
    Rhat_RB03(:,n)=(0.095*x) + 0.079.*x.*x;

figure(n),subplot(4,1,1),hold on,plot(wave,R(:,n),'k',wave,Rhat_RB03(:,n)/FQ,'g'),legend('Meas','RB03'),ylabel('R'),title('Roesler Boss 2003')
figure(n),subplot(4,1,2),hold on,plot(wave,aphi(:,n),'k-'),ylabel('a_{phyt} (m^{-1})')
figure(n),subplot(4,1,3),hold on,plot(wave,acdom(:,n),'b',wave,anap(:,n),'c'),ylabel('a_{cdm} (m^{-1})'),legend('CDOM','NAP')
figure(n),subplot(4,1,4),hold on,plot(wave,bbp(:,n),'-k'),ylabel('b_{bp} (m^{-1})'),xlabel('Wavelength (nm)')


end%nspec loop
'Retrieved coefficients at 443 nm'
'chl     acdom443     anap443     cp443     bbpratio     gamma'

iopRB03
