clear all
WL=input('enter wavelength resolution (1) SeaWIFS, (2) MODIS, (3) Hyperspectral')
if WL==1,wave=[412; 440; 490; 510; 555; 670];
elseif WL==2,wave=[412; 440; 490; 530; 550; 670];
elseif WL==3,Lres=input('enter spectral resolution (i.e. 1, 2, 5, 10, 25) in nm')
%Lres=1; %spectral resolution, try 5, 10, 25 nm
wave=400:Lres:700;,wave=wave';
end
N=input('enter number of retrieved parameters (1) aph, adm, bbp; (2) aph, acdom, anap, bbp; (3) aph, acdom, anap, bbplarge, bbpsmall')

%%%%%%%%%%%%%%%%%%%%%%%%%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');

%%%%
%%% 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=wave;

iopRP95 = RP95_invert(R'*FQ,wavelength,N);
if N==1, chl = iopRP95(:,1); ag = iopRP95(:,2);  bbps = iopRP95(:,3); Scdom = 0.015; S_bbps = 1.;
elseif N==2, chl = iopRP95(:,1); ag = iopRP95(:,2);  an = iopRP95(:,3); bbps = iopRP95(:,4); Scdom = 0.02; Snap=0.01; S_bbps = 1.;
elseif N==3, chl = iopRP95(:,1); ag = iopRP95(:,2);  an = iopRP95(:,3); bbps = iopRP95(:,4); bbpl = iopRP95(:,5); Scdom = 0.02; Snap=0.01; S_bbps = 1.; S_bbpl = 0;
end
iopRP95(:,1)=iopRP95(:,1)*0.0222;%convert to aphi(443);
    
%water basis vectors require temperature T and salinity S, output water IOPs at wave resolution
    [a_w,bb_w]=water_iops_PF_TScorr(wave,15,33);
% average spectrum
    aph=phyto_avg_abs(wave);
    
for n=1:nspec-1
    atotal = a_w + aph.*chl(n) + ag(n).*exp(-Scdom.*(wave-443.*ones(size(wave))));
    if N>1, atotal = atotal + an(n).*exp(-Snap.*(wave-443.*ones(size(wave))));,end
    bbtotal = bb_w + bbps(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbps) ;
    if N>2, bbtotal = bbtotal + bbpl(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbpl);,end
    x=bbtotal./(bbtotal+atotal)
    Rhat_RP95=(0.095*x) + 0.079.*x.*x;

figure(n),subplot(4,1,1),hold on,plot(wave,R(:,n),'k',wave,Rhat_RP95/FQ,'g'),legend('Meas','RP95'),ylabel('R')
figure(n),subplot(4,1,2),hold on,plot(wave,aph.*chl(n),'g-'),ylabel('a_{phyt}')
figure(n),subplot(4,1,3),hold on,plot(wave,ag(n).*exp(-Scdom.*(wave-443.*ones(size(wave)))),'c-'),ylabel('a_{cdm}'),
          if N>1,hold on,plot(wave,an(n).*exp(-Snap.*(wave-443.*ones(size(wave)))),'b-',wave,atotal-a_w-aph.*chl(n),'r'),legend('CDOM','NAP','CDM'),end
figure(n),subplot(4,1,4),hold on,plot(wave,bbps(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbps),'b-'),ylabel('b_{bp}')
          if N>2,hold on,plot(wave,bbpl(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbpl),'b-',wave,bbtotal-bb_w,'r'),legend('bbp_{small}','bbp_{large}','bbp_{total}'),end

end

('a_phyt(443)   a_cdm(443) b_bp(443)')

iopRP95
