clear all

waveS=[412; 440; 490; 510; 555; 670];

%%%%%%%%%%%%%%%%%%%%%%%%%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
%%%%

%choose your data file loaded into dat 
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)

RS=interp1(dat(:,1),dat(:,2:nspec),waveS,'linear');
figure(100),plot(waveS,RS),title('Measured spectra at SeaWiFS wavelengths'),ylabel('R(\lambda)'),xlabel('Wavelength(nm)')
%%

%%%%
%%% call QAA4_SeaWIFS model, output is bbp, aphi, acdm (443)
%%% uses all 6 SeaWIFS channels, note the aphi spectra that are output are
%%% contaminated with fluorescence
%%% this is not a regression type model but is an empirical model
%%% requires water IOPs

%water basis vectors require temperature T and salinity S, output water IOPs at wave resolution
    [a_w,bb_w]=water_iops_PF_TScorr(waveS,15,33);
for n=1:nspec-1
    n
[b_bp, a_phi, a_cdm]=QAA4_SeaWIFS(RS(:,n)*FQ,a_w,bb_w);
iopQAA(n,:)=[a_phi(2) a_cdm(2) b_bp(2)]

Rhat_QAA=(0.0949*(bb_w+b_bp')./(bb_w+b_bp'+a_w+a_phi'+a_cdm'))/1.7;%1.7 correction factor
figure(n),subplot(4,1,1),plot(waveS,RS(:,n),'k-',waveS,Rhat_QAA/FQ,'m')
figure(n),subplot(4,1,2),plot(waveS(1:5),a_phi(1:5),'g:')
figure(n),subplot(4,1,3),plot(waveS,a_cdm,'c:')
figure(n),subplot(4,1,4),plot(waveS,b_bp,'b:')

end

%% %%
%%% call RP95_invert model, 
%%% for MODIS and SeaWiFS wavelengths, use 3 outputs, aphi, acdm, bbpsmall (443)
%%% for hyperspectral try 4 aphi, acdom anap bbpsmall (443)
%%% for hyperspectral try 5 aphi, acdom anap bbpsmall bbplarge(443)
%%% 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*a* + ... Scdm = 0.015, Scdom = 0.02, Snap = 0.01
%%% bb = bbw + bbp(443)*(wave/443)^-n, n=-1 small, 0 large

wavelength=waveS;
N=1;
iopRP95 = RP95_invert(RS'*FQ,wavelength,N);
chl = iopRP95(:,1); ag = iopRP95(:,2);  bbps = iopRP95(:,3); Scdom = 0.015; S_bbps = 1.;
iopRP95(:,1)=iopRP95(:,1)*0.0222;%convert to aphi(443);
    
aph=phyto_avg_abs(waveS);
    
for n=1:nspec-1
    atotal = a_w + aph.*chl(n) + ag(n).*exp(-Scdom.*(waveS-443.*ones(size(waveS))));
    bbtotal = bb_w + bbps(n) .* ((443.*ones(size(waveS)) ./ waveS).^S_bbps) ;
    x=bbtotal./(bbtotal+atotal)
    Rhat_RP95=(0.095*x) + 0.079.*x.*x;

figure(n),subplot(4,1,1),hold on,plot(waveS,Rhat_RP95/FQ,'g'),legend('Meas','QAA','RP95'),ylabel('R(\lambda)')
figure(n),subplot(4,1,2),hold on,plot(waveS,aph.*chl(n),'g-'),legend('QAA','RP95'),ylabel('a_{phyt}(\lambda)')
figure(n),subplot(4,1,3),hold on,plot(waveS,ag(n).*exp(-Scdom.*(waveS-443.*ones(size(waveS)))),'c-'),ylabel('a_{cdm}(\lambda)'),legend('QAA','RP95')
figure(n),subplot(4,1,4),hold on,plot(waveS,bbps(n) .* ((443.*ones(size(waveS)) ./ waveS).^S_bbps),'b-'),ylabel('b_{bp}(\lambda)'),legend('QAA','RP95')
end

'Retrieved coefficients at 443 nm'
('a_phyt(443)   a_cdm(443) b_bp(443)')
iopQAA
iopRP95
