clear all
close all
%read and plot TSG
load('/Volumes/Boss_SSD/Perseverance/prod/Perseverance_InLine_TSG_20260120_20260309_Product_v20260310.mat');
figure(1);clf
scatter(tsg.lat,tsg.sss,8,tsg.sss)
title(sprintf(' SSS [PSU]'))
grid on
ylabel('SSS [PSU]');
xlabel('Latitude [°]');

figure(2);clf
scatter(tsg.lat,tsg.temp_cal,8,tsg.temp_cal)
title(sprintf('SST [°C]'))
grid on
ylabel('SST [°C]');
xlabel('Latitude [°]');

%compute density:
%input Fchl and Turbidity
cd '/Volumes/Boss_SSD/Perseverance/raw/FLNTU999'
B=dir('FLNTU999_2026012*.csv');
C=dir('FLNTU999_2026013*.csv');
D=dir('FLNTU999_202602*.csv');
E=dir('FLNTU999_202603*.csv');
B=[B;C;D;E];
tt=[];
Turb=[];
Fchl=[];
Lon=[];
Lat=[];
for i=1:length(B)
    A=rd_FLNTU(B(i).name);
    Fchl=[Fchl;A.FLNTU_Chla];
    Turb=[Turb;A.FLNTU_Turb];
    Lon=[Lon;A.lon];
    Lat=[Lat;A.lat];
    tt=[tt;datenum(char(A.dt))];
end

figure(3);clf
scatter (Lat,Fchl,4,Fchl)  % Chla  color /lat
title(sprintf(' Chl-a []'))
grid on
ylabel('Fchl [V]');
xlabel('Latitude [deg]');

I=find(Turb<2);
figure(4);clf
scatter (Lat(I),Turb(I),4,Turb(I))  % Chla  color /lat
title(sprintf('Turb []'))
grid on
ylabel('Turb []');
xlabel('Latitude [deg]');

cd '/Volumes/Boss_SSD/Perseverance/prod'
%% input high/low map resolution  
%mresol=input('map/bathy resolution :high map & bathy =2, high bathy & low map =1, low all=0')
%if isempty(mresol)
    mresol=0;
%end

%% plot SST map 
figure(5) 
clf
  maxmap=[5 1 0.5]
  bg=1
  
 minlat=min(tsg.lat)-bg;
 maxlat=max(tsg.lat)+bg;
 minlong=min(tsg.lon)-bg;
 maxlong=max(tsg.lon)+bg;
 minSST=min(tsg.temp_cal);
 maxSST=max(tsg.temp_cal);

 m_proj('Stereographic','longitude',160,'latitude',-60,'radius',20);

 sondes=[-5000 -3000 -2000 -1000 -500 -250 -0];
 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes','linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_gshhs_h('save', 'gumby');
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid;
[lontsg,lattsg]=m_ll2xy(tsg.lon,tsg.lat);
 hold on
 scatter (lontsg,lattsg,8,tsg.temp_cal)
 hold on
 caxis([minSST maxSST])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' SST [°C]'));
 colorbar
 grid on

 print('temprature.png','-dpng')

 %% plot SSS map 
 figure(6)
 clf
mresol= 0
 minSSS=min(tsg.sss);
 maxSSS=max(tsg.sss);
sondes=[-5000 -3000 -2000 -1000 -500 -250 -0];
 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(tsg.lon,tsg.lat);
 hold on
 scatter (lontsg,lattsg,8,tsg.sss)
 hold on
 caxis([minSSS maxSSS])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' SSS [PSU]'));
 colorbar
 grid on
  print('salinity.png','-dpng')

load('/Volumes/Boss_SSD/Perseverance/prod/Perseverance_InLine_ACS_20260120_20260309_Products_full_v20260325');
 
%% plot Chl map 
 figure(7)
 clf

 I=find(acs_prod.Chl_lineheight<0.01);
 acs_prod.Chl_lineheight(I)=NaN;
 minChl=prctile(log10(acs_prod.Chl_lineheight),1);
 maxChl=prctile(log10(acs_prod.Chl_lineheight),99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(acs_prod.lon,acs_prod.lat);
 hold on
 scatter (lontsg,lattsg,20,log10(acs_prod.Chl_lineheight),'filled')
 hold on
 caxis([minChl maxChl])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' Log_{10}{Chlorophyll} [mg m^{-3}]'));
 colorbar
 grid on
  print('chlorophyll.png','-dpng')

 %% plot gamma map 
 figure(8)
 clf

 mingamma=prctile(acs_prod.cp_gamma,1);
 maxgamma=prctile(acs_prod.cp_gamma,99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(acs_prod.lon,acs_prod.lat);
 hold on
 scatter (lontsg,lattsg,8,acs_prod.cp_gamma,'filled')
 hold on
 caxis([mingamma maxgamma])
 xlabel('lon');
 ylabel('lat');
 title('\gamma, c_p spectral slope');
 colorbar
 grid on
  print('gamma_cp.png','-dpng')
 

 
 
 %% plot F_Chl map 
 figure(9)
 clf

 I=find(acs_prod.Chl_lineheight<0.01);
 acs_prod.Chl_lineheight(I)=NaN;
 minChl=prctile(log10(Fchl),1);
 maxChl=prctile(log10(Fchl),99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(Lon,Lat);
 hold on
 scatter (lontsg,lattsg,8,log10(Fchl))
 hold on
 caxis([minChl maxChl])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' Log_{10}{F-Chlorophyl} [mg m^{-3}]'));
 colorbar
 grid on
  print('F_{chl}.png','-dpng')

%% plot Fchl/chl map
 [C,IA,IC] = unique(tt);
 Fchl_FB=interp1(tt(IA),Fchl(IA),datenum(acs_prod.dt),'linear');
X=Fchl_FB./acs_prod.Chl_lineheight;
I=find(isfinite(X)==1)
 BB=filtfilt(ones(1,20)/20,1,X(I))
figure(10)
 clf

 min_=prctile(log10(BB),1);
 max_=prctile(log10(BB),99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(acs_prod.lon(I),acs_prod.lat(I));
 hold on
 scatter (lontsg,lattsg,8,log10(BB),'filled')
 hold on
 caxis([min_ max_])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' log10(F_{Chlorophyl}/Chlorophyll) []'));
 colorbar
 grid on
  print('F_{chl}_ch;l.png','-dpng')





  

%% plot FDOM map 
 figure(10)
 clf

 ming=prctile(acs_prod.fdom,1);
 maxg=prctile(acs_prod.fdom,99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(acs_prod.lon,acs_prod.lat);
 hold on
 scatter(lontsg,lattsg,8,acs_prod.fdom,'filled')
 hold on
 caxis([ming maxg])
 xlabel('lon');
 ylabel('lat');
 title(sprintf('F_{DOM}'),'filled');
 colorbar
 grid on
 print('FDOM.png','-dpng')
%below is the analysis of Fchl as function of light and CDOM (Fe)



%turbidity map

 Turb_interp=interp1(tt(IA),Turb(IA),datenum(acs_prod.dt),'linear');

figure(11)
 clf

 min_=prctile(log10(Turb_interp),1);
 max_=prctile(log10(Turb_interp),99);


 if mresol== 0  % low map and low bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_elev('contour',sondes);
 elseif mresol== 1  % low map and high bathy resolution
     m_coast('patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 elseif mresol== 2  % high map and high bathy resolution
     m_usercoast('gumby','patch',[.7 .7 .7],'edgecolor','k');
     m_tbase('contour',sondes,'linecolor','k');
 end
 m_grid
 [lontsg,lattsg]=m_ll2xy(acs_prod.lon,acs_prod.lat);
 hold on
 scatter (lontsg,lattsg,8,log10(Turb_interp),'filled')
 hold on
 caxis([min_ max_])
 xlabel('lon');
 ylabel('lat');
 title(sprintf(' log10(Turbidity) []'));
 colorbar
 grid on
  print('Turbidity.png','-dpng')



%compute SZA along-track
%delta=23.45*sind(360/365*(yeardat-81));



 %compare ship fluorescence with AC-S
figure
 [C,IA,IC] = unique(tt);
 Fchl_FB=interp1(tt(IA),Fchl(IA),datenum(acs_prod.dt),'linear');
 loglog(Fchl_FB,acs_prod.Chl_lineheight,'.')
 xlabel('Fchl')
 ylabel('chl_{lh}')



 load('/Volumes/Boss_SSD/Perseverance/raw/Radiance_data/PERSEVERANCE_RADIOMETERS_20260326_0316.mat')
rad_interp=interp1(rad.time,rad.SMP21_SW_Flux,datenum(acs_prod.dt),'linear');
figure
 semilogy(acs_prod.dt,rad_interp/100,'.k')
 hold on
 semilogy(acs_prod.dt,(Fchl_FB-prctile(Fchl_FB,0.1))./acs_prod.Chl_lineheight,'.')
semilogy(acs_prod.dt,acs_prod.fdom*100,'g')
 hold off
 xlabel('date')
ylabel('Fchl/chl_{lh} - red, SWR/100 - black, FDOM - green')
%ylabel('Fchl/chl_{lh}')
 

figure
 scatter(log10(Fchl_FB),log10(acs_prod.Chl_lineheight),8,rad_interp)
 colorbar
 xlabel('log10(Fchl)')
ylabel('log10(chl_{lh})')
  print('Fchl_to_chl.png','-dpng')

figure
 plot(acs_prod.lat,rad_interp/100,'.k')
 hold on
 plot(acs_prod.lat,acs_prod.fdom*100,'g')
 plot(acs_prod.lat,(Fchl_FB-0.0823)./acs_prod.Chl_lineheight,'r.')
 hold off
 xlabel('Latitude')
ylabel('Fchl/chl_{lh} - red, SWR/100 - black')


 figure
  semilogy((Fchl_FB-0.0823)./acs_prod.Chl_lineheight,rad_interp,'.')

 

%compare POC to chl and turbidity
figure
subplot(2,1,1)
plot(acs_prod.Chl_lineheight,acs_prod.POC_cp,'.')
xlabel('Chl')
ylabel('POC')

 Turb_interp=interp1(tt(IA),Turb(IA),datenum(acs_prod.dt),'linear');

subplot(2,1,2)
plot(Turb_interp,acs_prod.POC_cp,'.')
xlabel('Turb')
ylabel('POC')

figure(17)
subplot(2,1,1)
plot(acs_prod.dt,acs_prod.Chl_lineheight,'.')
ylabel('Chl')
ylabel('Date')

 
subplot(2,1,2)
plot(acs_prod.dt,Fchl_FB,'.')
xlabel('Fchl')
ylabel('Date')



figure;
I=find(rad_interp<prctile(rad_interp,35));
plot((Fchl_FB(I)-0.0823)./acs_prod.Chl_lineheight(I),acs_prod.fdom(I),'.')

 plot(acs_prod.dt(I),rad_interp(I)/100,'.k')
 hold on
 plot(acs_prod.dt(I),(Fchl_FB(I)-0.0823)./acs_prod.Chl_lineheight(I),'.')
 plot(acs_prod.dt(I),acs_prod.fdom(I)*100,'*g')
 hold off
 xlabel('date')
ylabel('Fchl/chl_{lh} - red, SWR/100 - black, FDOM - green')

figure;
I=find(rad_interp<prctile(rad_interp,35));
plot((Fchl_FB(I)-0.0823)./acs_prod.Chl_lineheight(I),acs_prod.fdom(I),'.')

 semilogy(acs_prod.dt(I),rad_interp(I)/100,'.k')
 hold on
semilogy(acs_prod.dt(I),(Fchl_FB(I)-0.0823)./acs_prod.Chl_lineheight(I),'.')
 semilogy(acs_prod.dt(I),acs_prod.fdom(I)*100,'*g')
 hold off
 xlabel('date')
ylabel('Fchl/chl_{lh} - red, SWR/100 - black, FDOM - green')


figure
plot(acs_prod.fdom(I),(Fchl_FB(I)-0.0823)./acs_prod.Chl_lineheight(I),'*')
 xlabel('F_{DOM}')
ylabel('Fchl/chl_{lh}')



