clear all
Lres=input('enter spectral resolution (i.e. 1, 2, 5, 10, 25) in nm')
T=15;,S=33;%or enter your own values
%Lres=1; %spectral resolution, try 5, 10, 25 nm
wave=400:Lres:700;,wave=wave';
nL=length(wave);
statement=('we are going to allow four phytoplankton components (diatom- dinoflagellate- chlorophyte- crytophyte-like endmembers')
N=input('enter number of retrieved parameters (1) aph4, adm, bbp; (2) aph4, acdom, anap, bbp; (3) aph4, 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');
figure(100),plot(wave,R),title('Measured spectra '),ylabel('R(\lambda)'),xlabel('Wavelength(nm)')


%%%%
%%% call RPFG_invert model, 
%%% for the hyperspectral or reduced hyperspectral resolution try 
%%% allowing the 4 aphi, and a range of options for a_cdm bbp
%%% 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 + SUMchl(i)*aph_i* + ... Scdm = 0.015, Scdom = 0.02, Snap = 0.01
%%% bb = bbw + bbp(443)*(wave/443)^-n, n=-1 small, 0 large

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_species_abs(wave);

%%% perform inversion
iopRPFG = RPFG_invert(R'*FQ,wavelength,N);%converting the measured reflectance by FQ factor

chl = iopRPFG(:,1:4);%output is chl 
ag = iopRPFG(:,5);
if N==1,   bbps = iopRPFG(:,6); Scdom = 0.015; S_bbps = 1.;
elseif N==2, an = iopRPFG(:,6); bbps = iopRPFG(:,7); Scdom = 0.02; Snap=0.01; S_bbps = 1.;
elseif N==3, an = iopRPFG(:,6); bbps = iopRPFG(:,7); bbpl = iopRPFG(:,8); Scdom = 0.02; Snap=0.01; S_bbps = 1.; S_bbpl = 0;
end

astar443=[0.0183    0.0173    0.0170    0.0262];%for the four phytoplankton basis vectors
aB2R=[1.3280    1.2426    1.1819    1.7935];%443:676 ratio for the four phytoplankton basis vectors
for i=1:4,iopRPFG(:,i)=iopRPFG(:,i).*astar443(i);,end%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) + aph(:,2).*chl(n,2) + aph(:,3).*chl(n,3) + aph(:,4).*chl(n,4);%total phyto abs

if N==1,  Scdom = 0.015; S_bbps = 1.;
  acdom(:,n)= ag(n).*exp(-Scdom*(wave-443*ones(size(wave))));
  Bbps(:,n) = bbps(n) .* ((443*ones(size(wave)) ./ wave).^S_bbps) ;
    acdm(:,n) = acdom(:,n);,anap(1:nL,n)=zeros;
    Bbp(:,n)  = Bbps(:,n);,Bbpl(1:nL,n)=zeros;
elseif N==2,  Scdom = 0.02; Snap = 0.01; S_bbps = 1.;
    acdom(:,n)= ag(n).*exp(-Scdom*(wave-443*ones(size(wave))));
    anap(:,n) = an(n).*exp(-Snap*(wave-443*ones(size(wave))));
    Bbps(:,n) = bbps(n) .* ((443*ones(size(wave)) ./ wave).^S_bbps) ;
    acdm(:,n) = acdom(:,n) + anap(:,n);
    Bbp(:,n)  = Bbps(:,n);,Bbpl(1:nL,n)=zeros;
elseif N==3, Scdom = 0.02; Snap = 0.01; S_bbps = 1.; S_bbpl = 0;
    acdom(:,n)= ag(n).*exp(-Scdom*(wave-443*ones(size(wave))));
    anap(:,n) = an(n).*exp(-Snap*(wave-443*ones(size(wave))));
    Bbps(:,n) = bbps(n) .* ((443*ones(size(wave)) ./ wave).^S_bbps) ;
    Bbpl(:,n) = bbpl(n) .* ((443*ones(size(wave)) ./ wave).^S_bbpl) ;
    acdm(:,n) = acdom(:,n) + anap(:,n);
    Bbp(:,n)  = Bbps(:,n)  + Bbpl(:,n);
end
 atotal(:,n) = a_w + aphi(:,n) + acdm(:,n);
 bbtotal(:,n) = bb_w + Bbp(:,n);
 
g1 = 0.0949; g2 = 0.0794; 
x = bbtotal(:,n) ./ (atotal(:,n) + bbtotal(:,n));
    
    Rhat_RPFG(:,n)=(0.095*x) + 0.079.*x.*x;

figure(n),subplot(4,2,1),hold on,plot(wave,R(:,n),'k',wave,Rhat_RPFG(:,n)/FQ,'g'),legend('Meas','RPFG'),ylabel('R')
figure(n),subplot(4,2,2),hold on,plot(wave,R(:,n)-(Rhat_RPFG(:,n)/FQ),'m'),ylabel('Residual (R_{meas}-R_{obs}')
figure(n),subplot(4,2,3),hold on,plot(wave,aphi(:,n),'k-'),ylabel('a_{phyt} (m^{-1})')
figure(n),subplot(4,2,4),hold on,plot(wave,chl(n,1)*aph(:,1),'c',wave,chl(n,2)*aph(:,2),'b',wave,chl(n,3)*aph(:,3),'g',wave,chl(n,4)*aph(:,4),'r'),legend('diatom','dinoflagellate','chlorophyte','cryptophyte')
figure(n),subplot(4,2,5),hold on,plot(wave,acdm(:,n),'-k'),ylabel('a_{cdm} (m^{-1})')
figure(n),subplot(4,2,6),hold on,plot(wave,acdom(:,n),'b',wave,anap(:,n),'c'),ylabel('a_{cdm} (m^{-1})'),legend('CDOM','NAP')
figure(n),subplot(4,2,7),hold on,plot(wave,Bbp(:,n),'-k'),ylabel('b_{bp} (m^{-1})'),xlabel('Wavelength (nm)')
figure(n),subplot(4,2,8),hold on,plot(wave,Bbps(:,n),'-b',wave,Bbpl(:,n),'r'),xlabel('Wavelength (nm)'),ylabel('b_{bp} (m^{-1})'),legend('small','large')

%figure(n),subplot(4,1,2),hold on,plot(wave,aphi(:,n),'k-'),ylabel('a_{phyt}')
%hold on,plot(wave,chl(n,1)*aph(:,1),'c',wave,chl(n,2)*aph(:,2),'b',wave,chl(n,3)*aph(:,3),'g',wave,chl(n,4)*aph(:,4),'r'),legend('SUM','diatom','dinoflagellate','chlorophyte','cryptophyte')
%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(:,n)-a_w(:)-aphi(:,n),'k'),legend('CDOM','NAP','SUM'),end
%figure(n),subplot(4,1,4),hold on,plot(wave,bbps(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbps),'m-'),ylabel('b_{bp}')
%          if N>2,hold on,plot(wave,bbpl(n) .* ((443.*ones(size(wave)) ./ wave).^S_bbpl),'r-',wave,bbtotal(:,n)-bb_w(:),'k'),legend('bbp_{small}','bbp_{large}','SUM'),end

end%nspec loop
'Retrieved coefficients at 443 nm'
if N==1,('  a_Diatom     a_Dino   a_Chloro  a_Crpto   a_cdm     b_bp')
elseif N==2,('  a_Diatom     a_Dino   a_Chloro  a_Crpto   a_cdom    a_nap     b_bp')
elseif N==3,(' a_Diatom     a_Dino   a_Chloro   a_Crpto   a_cdom    a_nap     b_bps     b_bpl');
end

iopRPFG
