% ACT data
% select number of days
% to remove hawaii, scroll down and uncomment the "no Hawaii" statements
% select row to start in row.xls spreadsheet
% select parameters to run
% uncomment xlswrite statements

close all
clear all

load actXLSData_good % (bad data removed, generated from genXLSData_good.m)
% load actXLSData_all % (all data, generated from genXLSData_all.m)

cw=0.4; % water leaving radiance

% get number of days
days=500;

I_cbl=find((cbl(:,2)-cbl(1,2))<days);
I_hawaii=find((hawaii(:,2)-hawaii(1,2))<days);
I_maine=find((maine(:,2)-maine(1,2))<days);
I_michigan=find((michigan(:,2)-michigan(1,2))<days);
I_michigan2=find((michigan2(:,2)-michigan2(1,2))<days);
I_mlml=find((mlml(:,2)-mlml(1,2))<days);
I_skio=find((skio(:,2)-skio(1,2))<days);
I_usf=find((usf(:,2)-usf(1,2))<days);

tss=[cbl_tss(I_cbl);hawaii_tss(I_hawaii);maine_tss(I_maine);michigan_tss(I_michigan);michigan2_tss(I_michigan2);mlml_tss(I_mlml);skio_tss(I_skio);usf_tss(I_usf)];
tss_sd=[cbl_tss_sd(I_cbl);hawaii_tss_sd(I_hawaii);maine_tss_sd(I_maine);michigan_tss_sd(I_michigan);michigan2_tss_sd(I_michigan2);mlml_tss_sd(I_mlml);skio_tss_sd(I_skio);usf_tss_sd(I_usf)];
poc=[cbl_poc(I_cbl);hawaii_poc(I_hawaii);maine_poc(I_maine);michigan_poc(I_michigan);michigan2_poc(I_michigan2);mlml_poc(I_mlml);skio_poc(I_skio);usf_poc(I_usf)];
bac=[cbl_bac4h(I_cbl);hawaii_bac4h(I_hawaii);maine_bac4h(I_maine);michigan_bac4h(I_michigan);michigan2_bac4h(I_michigan2);mlml_bac4h(I_mlml);skio_bac4h(I_skio);usf_bac4h(I_usf)]-cw;
aqua=[cbl_aqua(I_cbl);hawaii_aqua(I_hawaii);maine_aqua(I_maine);michigan_aqua(I_michigan);michigan2_aqua(I_michigan2);mlml_aqua(I_mlml);skio_aqua(I_skio);usf_aqua(I_usf)];
mcvan=[cbl_mcvan(I_cbl);hawaii_mcvan(I_hawaii);maine_mcvan(I_maine);michigan_mcvan(I_michigan);michigan2_mcvan(I_michigan2);mlml_mcvan(I_mlml);skio_mcvan(I_skio);usf_mcvan(I_usf)];
%added nominal PRE sensor to remove negaives. When still negative occurs,
%use min value as cal.
wetlabs=[cbl_wetlabs(I_cbl)-0.3;hawaii_wetlabs(I_hawaii)-min(hawaii_wetlabs(I_hawaii));maine_wetlabs(I_maine)-min(maine_wetlabs(I_maine));michigan_wetlabs(I_michigan);michigan2_wetlabs(I_michigan2)-0.1;mlml_wetlabs(I_mlml)-0.3;skio_wetlabs(I_skio)-0.3;usf_wetlabs(I_usf)];
ysi=[cbl_ysi(I_cbl);hawaii_ysi(I_hawaii);maine_ysi(I_maine);michigan_ysi(I_michigan);michigan2_ysi(I_michigan2);mlml_ysi(I_mlml);skio_ysi(I_skio);usf_ysi(I_usf)];
% add nominal PRE sensor value to all insitu results to remove negatives
insitu=[cbl_insitu(I_cbl)+3.6;hawaii_insitu(I_hawaii)+3.0;maine_insitu(I_maine)+2.8;michigan_insitu(I_michigan);michigan2_insitu(I_michigan2)+0.7;mlml_insitu(I_mlml)+3.2;skio_insitu(I_skio)+0.6;usf_insitu(I_usf)+6.6];

% remove zeros from ysi, insitu, and wetlabs arrays
II=find(ysi==0);
ysi(II)=NaN;
II=find(insitu<=0);
insitu(II)=NaN;
II=find(wetlabs<=0);
wetlabs(II)=NaN;

allData=[tss bac aqua insitu mcvan wetlabs ysi poc tss_sd]; % all sensors
% allData=[tss bac aqua wetlabs]; % backscatter sensors only
% allData=[tss bac insitu mcvan ysi]; % sidescatter only

% make an index array so we can filter by sites
index_cbl=size(cbl(I_cbl));
index_hawaii=size(hawaii(I_hawaii));
index_maine=size(cbl(I_maine));
index_michigan=size(michigan(I_michigan));
index_michigan2=size(michigan2(I_michigan2));
index_mlml=size(mlml(I_mlml));
index_skio=size(skio(I_skio));
index_usf=size(usf(I_usf));

% add michigans together
index_michigans=index_michigan+index_michigan2;

% add all together
index_all=index_cbl+index_hawaii+index_maine+index_michigan+index_michigan2+index_mlml+index_skio+index_usf;
index=[ones(index_cbl(1,1),1)*1;ones(index_hawaii(1,1),1)*2;ones(index_maine(1,1),1)*3;ones(index_michigan(1,1),1)*4;ones(index_michigan2(1,1),1)*5;ones(index_mlml(1,1),1)*6;ones(index_skio(1,1),1)*7;ones(index_usf(1,1),1)*8];
I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,6).*allData(:,4))); % one of each
I_all=find(isfinite(allData(:,1).*allData(:,2).*allData(:,6).*allData(:,8))); % all sensors where we have the same # points for tss and bac
index_b=NaN*ones(index_all(1,1),1);index_b(I)=1;
index_b_all=NaN*ones(index_all(1,1),1);index_b_all(I_all)=1;

I1=find(isfinite(index_b(1:index_cbl(1,1)))); % cbl
I2=find(isfinite(index_b(index_cbl(1,1)+1:index_cbl(1,1)+index_hawaii(1,1))))+index_cbl(1,1); % hawaii
I3=find(isfinite(index_b(index_cbl(1,1)+index_hawaii(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1))))+index_cbl(1,1)+index_hawaii(1,1); % maine
I4=find(isfinite(index_b(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1); % michigan+michigan2
I5=find(isfinite(index_b(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1); % mlml
I6=find(isfinite(index_b(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1); % skio
I7=find(isfinite(index_b(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1)+index_usf(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1); % usf

I1_all=find(isfinite(index_b_all(1:index_cbl(1,1)))); % cbl
I2_all=find(isfinite(index_b_all(index_cbl(1,1)+1:index_cbl(1,1)+index_hawaii(1,1))))+index_cbl(1,1); % hawaii
I3_all=find(isfinite(index_b_all(index_cbl(1,1)+index_hawaii(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1))))+index_cbl(1,1)+index_hawaii(1,1); % maine
I4_all=find(isfinite(index_b_all(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1); % michigan+michigan2
I5_all=find(isfinite(index_b_all(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1); % mlml
I6_all=find(isfinite(index_b_all(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1); % skio
I7_all=find(isfinite(index_b_all(index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1)+1:index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1)+index_usf(1,1))))+index_cbl(1,1)+index_hawaii(1,1)+index_maine(1,1)+index_michigans(1,1)+index_mlml(1,1)+index_skio(1,1); % usf

%%%%%%%%%%%%% start no Hawaii
% noHawaii=find(index~=2);
% tss=tss(noHawaii);
% tss_sd=tss_sd(noHawaii);
% poc=poc(noHawaii);
% bac=bac(noHawaii);
% wetlabs=wetlabs(noHawaii);
% ysi=ysi(noHawaii);
% insitu=insitu(noHawaii);
% mcvan=mcvan(noHawaii);
% aqua=aqua(noHawaii);
%  allData=[tss bac aqua insitu mcvan wetlabs ysi poc]; % all sensors
%%%%%%%%%%%%% end no Hawaii

% find std(TSS)/mean(TSS) for 85 points (no Hawaii)
I=find(isfinite(tss));
% hist(tss_sd(I)/mean(tss(I)))

% cb=find(index==1);
% hi=find(index==2);
% me=find(index~=3);
% m1=find(index~=4);
% m2=find(index~=5);
% ml=find(index~=6);
% sk=find(index~=7);
% fl=find(index~=8);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tss correlations, follow step numbers
% I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,3).*allData(:,4).*allData(:,5).*allData(:,6).*allData(:,7))); % remove NaNs from corr calculation
% tss is allData(:,1)
% bac is .*allData(:,2)
% aquatec is .*allData(:,3) backscatter
% insitu is .*allData(:,4) sidescatter
% mcvan is .*allData(:,5) sidescatter
% wetlabs is .*allData(:,6) backscatter
% ysi is .*allData(:,7) sidescatter

% #1 pick sensors to use:
% I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,3).*allData(:,4).*allData(:,5).*allData(:,6).*allData(:,7))); % all sensors
% I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,4).*allData(:,5).*allData(:,7))); % sidescatter
% I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,3).*allData(:,6))); % backscatter
% I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,3))); % one sensor
I=find(isfinite(allData(:,1).*allData(:,2).*allData(:,6).*allData(:,4))); % tss, bac, wetlabs, insitu

% #2 list names of sensors used to copy to xls spreadsheet:
% sensorsI='tss-bac-aqua-insitu-mcvan-wetlabs-ysi';
% sensorsI='tss-bac-wetlabs-aqua-insitu-mcvan-ysi'; % all sensors
% sensorsI='tss-bac-insitu-mcvan-ysi'; % sidescatter
% sensorsI='tss-bac-wetlabs-aqua'; % backscatter
% sensorsI='tss-bac-aqua'; % one sensor
sensorsI='tss-bac-wetlabs-insitu'; % one of each

% #3 enter comments on what was run
comments='all days, good data only, no hawaii';

% #4 select row in xls spreadsheet
n=139;

% #5 pick sensors to compare
% a1=bac; a2=0.02; a3=0.02; a4='bac';
%a1=wetlabs; a2=0.1; a3=0.04; a4='wetlabs';
% a1=aqua; a2=0.3; a3=0.05; a4='aquatec';
% a1=insitu; a2=0.1; a3=0.07; a4='insitu';
% a1=mcvan; a2=0.5; a3=0.19; a4='mcvan';
% a1=ysi; a2=0.2; a3=0.40; a4='ysi';

%b1=bac; b2=0.02; b3=0.02; b4='bac';
b1=wetlabs; b2=0.1; b3=0.04; b4='wetlabs';
% b1=aqua; b2=0.3; b3=0.05; b4='aquatec';
% b1=insitu; b2=0.1; b3=0.07; b4='insitu';
% b1=mcvan; b2=0.5; b3=0.19; b4='mcvan';
% b1=ysi; b2=0.2; b3=0.40; b4='ysi';

% % if running tss, use tss for a, above:
a4='tss';
z1=log(tss(I)/median(tss(I)));
delta_z1=tss_sd(I)./tss(I);
z2=log(b1(I)./median(b1(I)));
delta_z2=max(b2,b1(I)*b3)./b1(I);

% if not running tss:
% z1=log(a1(I)./median(a1(I)));
% delta_z1=max(a2,a1(I)*a3)./a1(I);
% z2=log(b1(I)./median(b1(I)));
% delta_z2=max(b2,b1(I)*b3)./b1(I);

% simmulate data and obtain correlation coefficient
for j=1:length(z1)
     zz1(j,:)=random('norm',z1(j),delta_z1(j),1,1000);
     zz2(j,:)=random('norm',z2(j),delta_z2(j),1,1000);
end
for i=1:1000
     R=corrcoef(zz1(:,i),zz2(:,i));
     rho(i)=R(1,2);
end

mean_rho= mean(rho);
confidence95=3*std(rho);

% get size of dataset
s=size(I); 
s=num2str(s(1,1)); % convert to string for excel entry

ab=[mean_rho confidence95]
cd={a4 b4 s sensorsI comments}

% range to print to in excel
range1=['A' int2str(n)];
range2=['C' int2str(n)];
range3=['D' int2str(n)];

% xlswrite('rho.xls', ab, 'rhodata', range1);
% xlswrite('rho.xls', ' ', 'rhodata', range2); % empty for xls calc
% xlswrite('rho.xls', cd, 'rhodata', range3);

% figure;hist(rho);title('');
% 
% to print on screen:
% confidence95Table=[confidence95_bac confidence95_mcvan confidence95_insitu confidence95_ysi];
% rhoTable=[mean_rho_bac mean_rho_mcvan mean_rho_insitu mean_rho_ysi];
% 
% disp('************************************************************');
% disp('confidence95 ():');
% dsply(confidence95Table,'    bac4h     mcvan(s)  insitu(s) ysi(s)',' ');
% disp('mean rho ():');
% dsply(rhoTable,'    bac4h     mcvan(s)  insitu(s) ysi(s)',' ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% calculate uncertainties
TSS_a=[allData(I1,1);allData(I2,1);allData(I3,1);allData(I4,1);allData(I5,1);allData(I6,1);allData(I7,1)];
cp_a=[allData(I1,2);allData(I2,2);allData(I3,2);allData(I4,2);allData(I5,2);allData(I6,2);allData(I7,2)];
bbp_a=[allData(I1,6);allData(I2,6);allData(I3,6);allData(I4,6);allData(I5,6);allData(I6,6);allData(I7,6)];
bsp_a=[allData(I1,4),;allData(I2,4);allData(I3,4);allData(I4,4);allData(I5,4);allData(I6,4);allData(I7,4)];
std_TSS_a=[allData(I1,9);allData(I2,9);allData(I3,9);allData(I4,9);allData(I5,9);allData(I6,9);allData(I7,9)];

%[p,S]=polyfit(cp_a,TSS_a,1);
%y = polyval(p,cp_a);
%rms_a=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))

[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(cp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(cp_a)+b);
loglog(cp_a,TSS_a,'.',cp_a,y,'-');
rms_a=sqrt(sum(((TSS_a-exp(y))./TSS_a/length(TSS_a)).^2));

% [p,S]=polyfit(bbp_a,TSS_a,1);
% y = polyval(p,bbp_a);
% rms_b=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))

[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(bbp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(bbp_a)+b);
loglog(bbp_a,TSS_a,'.',bbp_a,y,'-');
rms_b=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))


% [p,S]=polyfit(bsp_a,TSS_a,1);
% y = polyval(p,bsp_a);
% rms_c=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))
[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(bsp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(bsp_a)+b);
loglog(bsp_a,TSS_a,'.',bsp_a,y,'-');
rms_c=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))


TSS_a=[allData(I1,1);allData(I3,1);allData(I4,1);allData(I5,1);allData(I6,1);allData(I7,1)];
cp_a=[allData(I1,2);allData(I3,2);allData(I4,2);allData(I5,2);allData(I6,2);allData(I7,2)];
bbp_a=[allData(I1,6);allData(I3,6);allData(I4,6);allData(I5,6);allData(I6,6);allData(I7,6)];
bsp_a=[allData(I1,4);allData(I3,4);allData(I4,4);allData(I5,4);allData(I6,4);allData(I7,4)];
std_TSS_a=[allData(I1,9);allData(I3,9);allData(I4,9);allData(I5,9);allData(I6,9);allData(I7,9)];

% attenuation
%A=log(cp_a./median(cp_a));
%B=log(TSS_a./median(TSS_a));
%[p,S]=polyfit(A,B,1);
%y = polyval(p,A);
[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(cp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(cp_a)+b);
%rms_a=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))
% prctile(abs(B-y),5)
% prctile(abs(B-y),50)
% prctile(abs(B-y),95)
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)

% attenuation at 5%, 50%, and 95%
%prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,5) % 0.0144
%prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,50) % 0.1553
%prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,95) % 0.4488

% backscattering
% A=log(bbp_a./median(bbp_a));
% B=log(TSS_a./median(TSS_a));
% [p,S]=polyfit(A,B,1);
% y = polyval(p,A);
% rms_b=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))
% prctile(abs(B-y),5)
% prctile(abs(B-y),50)
% prctile(abs(B-y),95)
[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(bbp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(bbp_a)+b);
%rms_a=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))
% prctile(abs(B-y),5)
% prctile(abs(B-y),50)
% prctile(abs(B-y),95)
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)
% backscattering at 5%, 50%, and 95%
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,5) % 0.0088
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,50) % 0.0892
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,95) % 0.3525

% side-scattering
% A=log(bsp_a./median(bsp_a));
% B=log(TSS_a./median(TSS_a));
% [p,S]=polyfit(A,B,1);
% y = polyval(p,A);rms_c=sqrt(sum(((TSS_a-y)./TSS_a/length(TSS_a)).^2))
% prctile(abs(B-y),5)
% prctile(abs(B-y),50)
% prctile(abs(B-y),95)
% 
% % side-scattering at 5%, 50%, and 95%
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,5) %0.0322 
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,50) % 0.1929
% prctile(abs(TSS_a-exp(y)*median(TSS_a))./TSS_a,95) % 0.5049
[m,b,r,sm,sb,xc,yc,ct] = lsqcubic(log(bsp_a),log(TSS_a),0.02,std_TSS_a./TSS_a,1e-6);
y=exp(m*log(bsp_a)+b);
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)
% end calculate uncertainties

%multivariate analysis:
[b,bint,r,rint,stats] = regress(log(TSS_a),[log(bsp_a),log(cp_a),ones(85,1)])
y=exp(b(1)*log(bsp_a)+b(2)*log(cp_a)+b(3));
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)

[b,bint,r,rint,stats] = regress(log(TSS_a),[log(bbp_a),log(cp_a),ones(85,1)])
y=exp(b(1)*log(bbp_a)+b(2)*log(cp_a)+b(3));
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)

[b,bint,r,rint,stats] = regress(log(TSS_a),[log(bsp_a),log(bbp_a),ones(85,1)])
y=exp(b(1)*log(bsp_a)+b(2)*log(bbp_a)+b(3));
prctile(abs(TSS_a-y)./TSS_a,5)
prctile(abs(TSS_a-y)./TSS_a,50)
prctile(abs(TSS_a-y)./TSS_a,95)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tss plots, all sites
% figure('Units','inches','Resize','off','position',[0.5 0.5 3.5 8.0],'paperposition',[1.5 3.5 3.5 9.0]);
% 
% subplot('Position',[0.19 0.72 0.70 0.24]); %top panel: tss vs. bac
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% loglog(allData(I1,1),allData(I1,2),'*',allData(I2,1),allData(I2,2),'s',allData(I3,1),allData(I3,2),'o',allData(I4,1),allData(I4,2),'.',allData(I5,1),allData(I5,2),'+',allData(I6,1),allData(I6,2),'d',allData(I7,1),allData(I7,2),'p');
% % plot         I1: cbl-bac, cbl-tss            I2: hi-bac, hi-tss              I3: me-bac, me-tss              I4: mi-bac, mi-tss              I5: mlml-bac, mlml-tss         I6: skio-bac, skio-tss           I7: usf-bac, usf-tss 
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% axis([ 1 100 0.1 100]);
% % create ellipse around HI data
% % annotation('ellipse',[0.1336 0.5855 0.09914 0.109]);
% %xlabel('PM (g m^{-3})');
% ylabel('c(660) (m^{-1})');
% 
% subplot('Position',[0.19 0.40 0.70 0.24]); %middle panel: tss vs. wetlabs
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% %loglog(allData(I1,1),allData(I1,6)*0.0163,'*',allData(I2,1),allData(I2,6)*0.0163,'s',allData(I3,1),allData(I3,6)*0.0163,'o',allData(I4,1),allData(I4,6)*0.0163,'.',allData(I5,1),allData(I5,6)*0.0163,'+',allData(I6,1),allData(I6,6)*0.0163,'d',allData(I7,1),allData(I7,6)*0.0163,'p')
% loglog(allData(I1,1),allData(I1,6)*1.0,'*',allData(I2,1),allData(I2,6)*1.0,'s',allData(I3,1),allData(I3,6)*1.0,'o',allData(I4,1),allData(I4,6)*1.0,'.',allData(I5,1),allData(I5,6)*1.0,'+',allData(I6,1),allData(I6,6)*1.0,'d',allData(I7,1),allData(I7,6)*1.0,'p')
% 
% %draw legend
% hold on
% plot(45,.035/0.01,'Marker','*','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% plot(45,.016/0.01,'Marker','s','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
% plot(45,.0085/0.01,'Marker','.','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% plot(45,.0042/0.01,'Marker','+','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% plot(45,.0018/0.01,'Marker','d','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
% axis([1 5 1 5]);
% hold off
% 
% annotation('textbox','Position',[0.78 0.490 0.05 0.05],'LineStyle','none',...
%   'FontName','Times New Roman','FontSize',10,'String',{'CB'},'FitHeightToText','on');
% annotation('textbox','Position',[0.78 0.465 0.05 0.05],'LineStyle','none',...
%   'FontName','Times New Roman','FontSize',10,'String',{'HI'},'FitHeightToText','on');
% annotation('textbox','Position',[0.78 0.441 0.05 0.05],'LineStyle','none',...
%   'FontName','Times New Roman','FontSize',10,'String',{'MI'},'FitHeightToText','on');
% annotation('textbox','Position',[0.78 0.415 0.05 0.05],'LineStyle','none',...
%   'FontName','Times New Roman','FontSize',10,'String',{'ML'},'FitHeightToText','on');
% annotation('textbox','Position',[0.78 0.385 0.05 0.05],'LineStyle','none',...
%   'FontName','Times New Roman','FontSize',10,'String',{'SK'},'FitHeightToText','on');
% 
% axis([ 1 100 0.1 100]);
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% %xlabel('PM (g m^{-3})');
% ylabel('b_{b}(700) (NTU)');
% 
% subplot('Position',[0.19 0.08 0.70 0.24]); %bottom panel: tss vs. insitu
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% loglog(allData(I1,1),allData(I1,4),'*',allData(I2,1),allData(I2,4),'s',allData(I3,1),allData(I3,4),'o',allData(I4,1),allData(I4,4),'.',allData(I5,1),allData(I5,4),'+',allData(I6,1),allData(I6,4),'d',allData(I7,1),allData(I7,4),'p')
% axis([1 100 1 100]);
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% xlabel('PM (g m^{-3})');
% ylabel('b_{s}(860) (NTU)');
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tss plots, all sites
% figure('Units','inches','Resize','off','position',[0.5 0.5 3.5 8.0],'paperposition',[1.5 3.5 3.5 9.0]);
% 
% subplot('Position',[0.19 0.72 0.70 0.24]); %top panel: tss vs. bac
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% plot(allData(I1,1),allData(I1,2),'*',allData(I2,1),allData(I2,2),'s',allData(I3,1),allData(I3,2),'o',allData(I4,1),allData(I4,2),'.',allData(I5,1),allData(I5,2),'+',allData(I6,1),allData(I6,2),'d',allData(I7,1),allData(I7,2),'p');
% % plot         I1: cbl-bac, cbl-tss            I2: hi-bac, hi-tss              I3: me-bac, me-tss              I4: mi-bac, mi-tss              I5: mlml-bac, mlml-tss         I6: skio-bac, skio-tss           I7: usf-bac, usf-tss 
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% axis([ 1 100 0.0 40]);
% % create ellipse around HI data
% % annotation('ellipse',[0.1336 0.5855 0.09914 0.109]);
% %xlabel('PM (g m^{-3})');
% %ylabel('c(660) (m^{-1})');
% 
% subplot('Position',[0.19 0.40 0.70 0.24]); %middle panel: tss vs. wetlabs
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% plot(allData(I1,1),allData(I1,6)*1.0,'*',allData(I2,1),allData(I2,6)*1.0,'s',allData(I3,1),allData(I3,6)*1.0,'o',allData(I4,1),allData(I4,6)*1.0,'.',allData(I5,1),allData(I5,6)*1.0,'+',allData(I6,1),allData(I6,6)*1.0,'d',allData(I7,1),allData(I7,6)*1.0,'p')
% 
% 
% axis([ 1 100 0.0 50]);
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% %xlabel('PM (g m^{-3})');
% %ylabel('b_{b}(700) (m^{-1})');
% 
% subplot('Position',[0.19 0.08 0.70 0.24]); %bottom panel: tss vs. insitu
% set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
% plot(allData(I1,1),allData(I1,4),'*',allData(I2,1),allData(I2,4),'s',allData(I3,1),allData(I3,4),'o',allData(I4,1),allData(I4,4),'.',allData(I5,1),allData(I5,4),'+',allData(I6,1),allData(I6,4),'d',allData(I7,1),allData(I7,4),'p')
% axis([1 100 0 80]);
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
% set(gca,'yMinorTick','off');
% set(gca,'xMinorTick','off');
% xlabel('PM (g m^{-3})');
% %ylabel('b_{s}(860) (NTU)');% subplot('Position',[0.55 0.12 0.29 0.35]); %bottom right panel: legend
% % plot(1.5,4.2,'Marker','*','MarkerSize',[10],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% % hold on
% % plot(1.5,3.6,'Marker','s','MarkerSize',[10],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
% % plot(1.5,3.05,'Marker','.','MarkerSize',[10],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% % plot(1.5,2.45,'Marker','+','MarkerSize',[10],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
% % plot(1.5,1.96,'Marker','d','MarkerSize',[10],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
% % axis([1 5 1 5]);
% % hold off
% % % legend text:
% % annotation('textbox','Position',[0.62 0.380 0.05 0.05],'LineStyle','none',...
% %   'FontName','Times New Roman','FontSize',10,'String',{'Chesapeake Bay'},'FitHeightToText','on');
% % annotation('textbox','Position',[0.62 0.340 0.05 0.05],'LineStyle','none',...
% %   'FontName','Times New Roman','FontSize',10,'String',{'Hawaii'},'FitHeightToText','on');
% % annotation('textbox','Position',[0.62 0.300 0.05 0.05],'LineStyle','none',...
% %   'FontName','Times New Roman','FontSize',10,'String',{'Michigan'},'FitHeightToText','on');
% % annotation('textbox','Position',[0.62 0.260 0.05 0.05],'LineStyle','none',...
% %   'FontName','Times New Roman','FontSize',10,'String',{'Moss Landing'},'FitHeightToText','on');
% % annotation('textbox','Position',[0.62 0.220 0.05 0.05],'LineStyle','none',...
% %   'FontName','Times New Roman','FontSize',10,'String',{'Skidaway'},'FitHeightToText','on');
% % box('off');
% % axis('off');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% above two figures on same plot
% figure('Units','inches','Resize','off','position',[0 0 7.25 8.0],'paperposition',[0 0 7.25 8.0]);
figure('Units','inches','Resize','off','outerposition',[0 0 7.25 7.25],'paperposition',[0 0 7.25 7.25]); %density: TOP LEFT PANEL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tss plots, all sites - left panels
subplot('Position',[0.1 0.72 0.275 0.235]); %top left panel: tss vs. bac
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
loglog(allData(I1,1),allData(I1,2),'*',allData(I2,1),allData(I2,2),'s',allData(I3,1),allData(I3,2),'o',allData(I4,1),allData(I4,2),'.',allData(I5,1),allData(I5,2),'+',allData(I6,1),allData(I6,2),'d',allData(I7,1),allData(I7,2),'p');
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
axis([ 1 100 0.1 100]);
ylabel('c(660) (m^{-1})');

subplot('Position',[0.1 0.41 0.275 0.235]); %middle left panel: tss vs. wetlabs
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
loglog(allData(I1,1),allData(I1,6)*1.0,'*',allData(I2,1),allData(I2,6)*1.0,'s',allData(I3,1),allData(I3,6)*1.0,'o',allData(I4,1),allData(I4,6)*1.0,'.',allData(I5,1),allData(I5,6)*1.0,'+',allData(I6,1),allData(I6,6)*1.0,'d',allData(I7,1),allData(I7,6)*1.0,'p')

% draw legend
hold on
plot(45,.035/0.01,'Marker','*','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(45,.016/0.01,'Marker','s','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(45,.0085/0.01,'Marker','.','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(45,.0042/0.01,'Marker','+','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(45,.0018/0.01,'Marker','d','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
axis([1 5 1 5]);
hold off

annotation('textbox','Position',[0.334 0.505 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'CB'},'FitHeightToText','on');
annotation('textbox','Position',[0.334 0.476 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'HI'},'FitHeightToText','on');
annotation('textbox','Position',[0.334 0.450 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'MI'},'FitHeightToText','on');
annotation('textbox','Position',[0.334 0.425 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'ML'},'FitHeightToText','on');
annotation('textbox','Position',[0.334 0.4 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'SK'},'FitHeightToText','on');

axis([ 1 100 0.1 100]);
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
ylabel('b_{b}(700) (NTU)');

subplot('Position',[0.1 0.1 0.275 0.235]); %bottom left panel: tss vs. insitu
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
loglog(allData(I1,1),allData(I1,4),'*',allData(I2,1),allData(I2,4),'s',allData(I3,1),allData(I3,4),'o',allData(I4,1),allData(I4,4),'.',allData(I5,1),allData(I5,4),'+',allData(I6,1),allData(I6,4),'d',allData(I7,1),allData(I7,4),'p')
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
xlabel('PM (g m^{-3})');
ylabel('b_{s}(860) (NTU)');
axis([1 100 1 100]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tss plots, all sites - right panels 
subplot('Position',[0.46 0.72 0.275 0.235]); %top right panel: tss vs. bac
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
plot(allData(I1,1),allData(I1,2),'*',allData(I2,1),allData(I2,2),'s',allData(I3,1),allData(I3,2),'o',allData(I4,1),allData(I4,2),'.',allData(I5,1),allData(I5,2),'+',allData(I6,1),allData(I6,2),'d',allData(I7,1),allData(I7,2),'p');
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
axis([ 1 100 0.0 40]);

subplot('Position',[0.46 0.41 0.275 0.235]); %middle right panel: tss vs. wetlabs
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
plot(allData(I1,1),allData(I1,6)*1.0,'*',allData(I2,1),allData(I2,6)*1.0,'s',allData(I3,1),allData(I3,6)*1.0,'o',allData(I4,1),allData(I4,6)*1.0,'.',allData(I5,1),allData(I5,6)*1.0,'+',allData(I6,1),allData(I6,6)*1.0,'d',allData(I7,1),allData(I7,6)*1.0,'p')
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
axis([ 1 100 0.0 50]);

subplot('Position',[0.46 0.1 0.275 0.235]); %bottom right panel: tss vs. insitu
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
plot(allData(I1,1),allData(I1,4),'*',allData(I2,1),allData(I2,4),'s',allData(I3,1),allData(I3,4),'o',allData(I4,1),allData(I4,4),'.',allData(I5,1),allData(I5,4),'+',allData(I6,1),allData(I6,4),'d',allData(I7,1),allData(I7,4),'p')
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.75);
set(gca,'yMinorTick','off');
set(gca,'xMinorTick','off');
axis([1 100 0 80]);
xlabel('PM (g m^{-3})');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% histograms
% figure('Units','inches','Resize','off','position',[0.5 0.5 7.5 4.5],'paperposition',[1.5 3.5 7.5 4.5]);
h=figure('Units','inches','Resize','off','position',[0 0 7.25 5.25],'paperposition',[0.5 3 7.25 5.25]);

subplot(2,4,1)
hist(allData(I,1))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('PM','FontName','Times New Roman','FontSize',10)
subplot(2,4,5)
hist(log(allData(I,1)/nanmedian(allData(I,1))))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('log\{PM/median(PM)\}','FontName','Times New Roman','FontSize',10)
subplot(2,4,2)
hist(allData(I,2))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('c','FontName','Times New Roman','FontSize',10)
subplot(2,4,6)
hist(log(allData(I,2)/nanmedian(allData(I,2))))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('log\{c/median(c)\}','FontName','Times New Roman','FontSize',10)
subplot(2,4,3)
hist(allData(I,6))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('b_{b}','FontName','Times New Roman','FontSize',10)
subplot(2,4,7)
hist(log(allData(I,6)/nanmedian(allData(I,6))))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('log\{b_{b}/median(b_{b})\}','FontName','Times New Roman','FontSize',10)
subplot(2,4,4)
hist(allData(I,4))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('b_{s}','FontName','Times New Roman','FontSize',10)
subplot(2,4,8)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
hist(log(allData(I,4)/nanmedian(allData(I,4))))
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[.7 .7 .7],'EdgeColor','k')
set(gca,'LineWidth',0.75);
xlabel('log\{b_{s}/median(b_{s})\}','FontName','Times New Roman','FontSize',10)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% poc plots
% POC/tss vs. bb/c
figure('Units','inches','Resize','on','position',[0.5 0.5 3.5 2.5],'paperposition',[1.5 3.5 3.5 2.5]);
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0; 0 0 0; 0 0 0]);
PT=allData(:,8)./allData(:,1);
bbc=allData(:,6)*0.0163./allData(:,2); %wetlab conversion factor
% plot(PT(I1),bbc(I1),'*',PT(I2),bbc(I2),'s',PT(I3),bbc(I3),'o',PT(I4),bbc(I4),'.',PT(I5),bbc(I5),'+',PT(I6),bbc(I6),'d',PT(I7),bbc(I7),'o');
plot(PT(I1_all),bbc(I1_all),'*',PT(I2_all),bbc(I2_all),'s',PT(I3_all),bbc(I3_all),'o',PT(I4_all),bbc(I4_all),'.',PT(I5_all),bbc(I5_all),'+',PT(I6_all),bbc(I6_all),'d',PT(I7_all),bbc(I7_all),'o');

%draw legend
hold on
plot(0.375,.037,'Marker','*','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(0.375,.034,'Marker','s','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
plot(0.375,.0307,'Marker','.','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(0.375,.028,'Marker','+','MarkerSize',[8],'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0]);
plot(0.375,.025,'Marker','d','MarkerSize',[8],'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]);
% axis([1 5 1 5]);
hold off

% legend('CB','Hi','Mi','ML','Sk');
% legend('boxoff');
set(gca,'FontName','Times New Roman','FontSize',10);
set(gca,'LineWidth',0.75);
axis([0 0.45 0.005 0.04])
% axis([0.005 0.04 0 0.45]); %reverse axes
xlabel('POC [g m^{-3}]/ PM [g m^{-3}]')
ylabel('b_{b} [m^{-1}] / c [m^{-1}]');
% xlabel('b_{b} [m^{-1}] / c [m^{-1}]');
% ylabel('POC [g m^{-3}]/ PM [g m^{-3}]')

%draw legend
annotation('textbox','Position',[0.80 0.86 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'CB'},'FitHeightToText','on');
annotation('textbox','Position',[0.80 0.80 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'HI'},'FitHeightToText','on');
annotation('textbox','Position',[0.80 0.73 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'MI'},'FitHeightToText','on');
annotation('textbox','Position',[0.80 0.67 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'ML'},'FitHeightToText','on');
annotation('textbox','Position',[0.80 0.61 0.05 0.05],'LineStyle','none',...
  'FontName','Times New Roman','FontSize',10,'String',{'SK'},'FitHeightToText','on');

%% Create doublearrow
% annotation('doublearrow',[0.3196 0.7732],[0.7966 0.3608]);
% % Create textbox
% annotation(...
%   'textbox',...
%   'EdgeColor',[1 1 1],...
%   'FontName','Times New Roman',...
%   'FontSize',12,...
%   'Position',[0.3643 0.4334 0.4804 0.1019],...
%   'FitHeightToText','off',...
%   'HorizontalAlignment','center',...
%   'String',{'Change consistent with change in','organic concentration (index of refraction)'});
 
