clear all
Lres=1; %spectral resolution, try 5, 10, 25 nm
waveH=400:Lres:700;,waveH=waveH';
waveS=[412; 440; 490; 510; 555; 670];
waveM=[412; 440; 490; 530; 550; 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
%%%%
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)
RH=interp1(dat(:,1),dat(:,2:nspec),waveH,'linear');
RS=interp1(dat(:,1),dat(:,2:nspec),waveS,'linear');
RM=interp1(dat(:,1),dat(:,2:nspec),waveM,'linear');

%%%%
%%% call GSM model, output is chl, acdm(443), bbp(443) 
%%% Convert chl to a(443) using a* in model
%%% for first 5 MODIS wavelengths only
%%% assumes Rrs beneath water interface and radiance reflectance
%%% R = 0.095 (bb/a+bb) + 0.079 (bb/a+bb)^2
%%% a = aw + chl*a* + acdm(443)*exp(-0.0206(wave-443))
%%% bb = bbw + bbp(443)*(wave/443)^-1.0337

iopGSM = GSM01_invert(FQ*RM(1:5,:)');%only uses first 5 wavelengths
chl = iopGSM(:,1); ag = iopGSM(:,2); bbp = iopGSM(:,3);
iopGSM(:,1)=0.05582*iopGSM(:,1)% <chl> * a*chl(443) = aph(443)

g1 = 0.0949; g2 = 0.0794; S = 0.0206; S_bbp = 1.0337;
wave = [412 443 490 510 555]';
bbw = [0.00333 0.00239 0.00155 0.00130 0.000925]';	%%%%% Smith/Baker1981
aw =  [0.00469 0.00721 0.01500 0.0325  0.05920]';    %%% Pope/Fry1997
aph = [0.00665 0.05582 0.02055 0.01910 0.01015]';   %%% Maritorena et al (2002)
for n=1:nspec-1
atotal(:,n) = aw + aph*chl(n) + ag(n).*exp(-S*(wave-443*ones(size(wave))));
bbtotal(:,n) = bbw + bbp(n) .* ((443*ones(size(wave)) ./ wave).^S_bbp);
x(:,n) = bbtotal(:,n) ./ (atotal(:,n) + bbtotal(:,n));
Rhat_GSM(:,n) = (g1 .*x(:,n) + g2 .* x (:,n).* x(:,n))/FQ;
figure(n),subplot(4,1,1),plot(waveM,RM(:,n),'k',waveM(1:5),Rhat_GSM(:,n),'c')
figure(n),subplot(4,1,2),plot(waveM(1:5),aph*chl(n),'g--')
figure(n),subplot(4,1,3),plot(waveM(1:5),ag(n).*exp(-S*(waveM(1:5)-443*ones(size(wave(1:5))))),'c--')
figure(n),subplot(4,1,4),plot(waveM(1:5),bbp(n) .* ((443*ones(size(wave(1:5))) ./ wave(1:5)).^S_bbp),'b--')

end
%%%%
%%% call QAA4_MODIS model, output is bbp, aphi, acdm (443)
%%% uses all 6 MODIS 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(waveM,15,33);
for n=1:nspec-1
    n
[b_bp, a_phi, a_cdm]=QAA4_MODIS(RM(:,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');
figure(n),subplot(4,1,1),hold on,plot(waveM,Rhat_QAA/FQ,'m')
figure(n),subplot(4,1,2),hold on,plot(waveM(1:5),a_phi(1:5),'g:')
figure(n),subplot(4,1,3),hold on,plot(waveM,a_cdm,'c:')
figure(n),subplot(4,1,4),hold on,plot(waveM,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=waveM;
N=1;
iopRP95 = RP95_invert(RM'*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(waveM);
    
for n=1:nspec-1
    atotal = a_w + aph.*chl(n) + ag(n).*exp(-Scdom.*(waveM-443.*ones(size(waveM))));
    bbtotal = bb_w + bbps(n) .* ((443.*ones(size(waveM)) ./ waveM).^S_bbps) ;
    x=bbtotal./(bbtotal+atotal)
    Rhat_RP95=(0.095*x) + 0.079.*x.*x;

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

%wavelength=waveH;
%N=3;
%iopRP95 = RP95_invert(RH'*FQ,wavelength,N);
%iopRP95(:,1)=iopRP95(:,1)*0.0222%convert to aphi(443);

('a_phyt(443)   a_cdm(443) b_bp(443)')
iopGSM
iopQAA
iopRP95
