www.gusucode.com > Ogive_optimization_1.0.6 > Ogive_densitymap_sub.m
function [runok,densitymap_filename,txtlog]=Ogive_densitymap_sub(savepath0,fname,txtlog,errordir,SecOutput,DataOutput,DAT,DAT15,runix,id)
if nargin==0
load Ogive_densitymap_TESTDELETE
% DAT{1}(8)=250;
% DAT{1}(7)=100;
end
warning('off','MATLAB:MKDIR:DirectoryExists')
try
if DAT{12}(10,2)>0
figurepath=[DAT{9},filesep,'3_FIGURES',filesep,'03 FREQUENCY INTERVALS',filesep];
end
figurepath2=[DAT{9},filesep,'3_FIGURES',filesep,'03B EC_T vs BEA_TX',filesep];
% CONSTANTS
L=40.68*1000; %Latent heat of vaporization [J/mol] (Temperature dependent!) (from wikipedia)
Cp=1004.67; %Cp is the specific heat of dry air 1004.67 J kg-1 K-1
% SETUP
runok=1;
SCRIPT_LEVEL=3;
cm_densitymap_gray=load('cm_densitymap_gray');
cm_densitymap_gray=cm_densitymap_gray.cm_densitymap_gray;
h=DAT{2}(4);
hz=DAT{2}(1);
int=250; %downscaling of Ogives AND output matrix resolution
% if DAT{1}(18)==0 %STANDARD METHOD
N1=DAT{1}(7); %indices of timerange and runmean investigated
% elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
% N1=5000; %indices of timerange and runmean investigated
% end
N2=DAT{1}(8);
sig=DAT{1}(9);
ixes=DAT{14};
% LOOP
lvar=size(SecOutput{6,1},2);
ixmultLOG=nan(lvar,2);
Ogive30=cell(lvar,1);
densitymap_filename=cell(1,lvar);
yax=cell(1,lvar);
for ii=runix
if SecOutput{9,1}(ii)==0
if DAT{1}(18)==0 %STANDARD METHOD
mintime=5;
elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
mintime=1;
end
tmax=length(DataOutput{ii,1});
mp=round(tmax/2);
HFbound=DAT{1}(11);
lsp=linspace(mintime,tmax/(hz*60),N1);
mins=floor(DAT{1}(5)/60);
if any(lsp>=30)
lsp30=abs(lsp-30);
ixm30=find(lsp30==min(lsp30),1,'first');
lsp(ixm30)=30;
else
if any(lsp>=mins)
lsp30=abs(lsp-mins);
ixm30=find(lsp30==min(lsp30),1,'first');
lsp(ixm30)=mins;
else
ixm30=length(lsp);
end
end
twin=round(lsp*60*hz/2)'; %indexes
unitscale=[0 3 6];
secday=[1 86400];
Ylabels={'molm^{-2}s^{-1}';'molm^{-2}d^{-1}';'Wm^{-2}'};
Ylabels2={'';'m';'\mu'};
%ls with lowres-runmean emphasis
%
% ls=round(linspace(0,tmax/2,N2-1));ls=ls(2:end);
% ls0=logspace(-2.2,0,N2-1);
% ls0=1-ls0;
% ls0=ls0*(tmax/2);
% ls=round(fliplr(ls0+(30*hz)));
% SPECTRAL DECOMPOSITION AND MASS OGIVE CALCULATION
% covlog=nan(N1,1);
% fluxlog=nan(N1,1);
Ogive_MAT0=cell(N1,1);
Ogive_end_lindet=nan(N1,1);
rawdat1=DataOutput{ii,1};
rawdat2=DataOutput{ii,2};
ix1=max(1,mp-twin);
ix2=min(tmax,mp+twin);
datalength=diff([ix1,ix2]');
ixl2=round(mean(unique(diff(datalength)))*1.5);
[freq_raw_nat0,~,~,~] = Ogive_calcCospectra(hz,h,SecOutput{3,1}(1),rawdat1(ix1(end):ix2(end)),rawdat2(ix1(end):ix2(end)),0,'none','none',0);
rlog_f=log10([min(freq_raw_nat0) max(freq_raw_nat0)]);
lsf=linspace(rlog_f(1),rlog_f(2),int)';
fr_nat=(10.^lsf);
Flux_lindet=nan(N1,1);
BEA_STATS_LOG=cell(N1,1);
% Flux_lindet0=nan(N1,3);
if DAT{1}(18)==0 %STANDARD METHOD
l_Ogive_end=nan(N1,1);
lsrange=[-1.8 0]; %-2.2
ls0=fliplr((10^lsrange(2))-(10.^linspace(lsrange(1),lsrange(2),N2+1)));
ls=round(ls0*(tmax/2));
ls=ls(2:end);
rls=myrange(ls);
ls=round(interp1(rls,[rls(1) floor(length(DataOutput{ii,1})/2)],ls));
Ogive_end=cell(N1,1);
for ia=1:N1
wdat=rawdat1(ix1(ia):ix2(ia));
sdat=rawdat2(ix1(ia):ix2(ia));
WDAT=nan(length(ix1(ia):ix2(ia)),N2);
SDAT=nan(length(ix1(ia):ix2(ia)),N2);
mU2=mean(DataOutput{ii,6}(ix1(ia):ix2(ia)));
N2ix=find(ls<floor(length(wdat)/2),1,'last');
WDAT(:,1)=detrend(wdat);
SDAT(:,1)=detrend(sdat);
if ~isempty(N2ix)
for ic=1:N2ix
WDAT(:,ic+1)=wdat-my_runmean(wdat,ls(ic),'zero');
SDAT(:,ic+1)=sdat-my_runmean(sdat,ls(ic),[],[mean(sdat(1:ixl2)) mean(sdat(end-ixl2:end))]);
end
else
N2ix=0;
end
% ia=1;
% wdat1=rawdat1(ix1(ia):ix2(ia));
% sdat1=rawdat2(ix1(ia):ix2(ia));
% ia=2;
% wdat2=rawdat1(ix1(ia):ix2(ia));
% sdat2=rawdat2(ix1(ia):ix2(ia));
% x1=1:length(wdat1);
% x2=1:length(wdat2);
% ixl=length(wdat1)-length(wdat2);
% ixl2=round(ixl*1.5);
% x2=x2+round(ixl/2);
% ic=2;
% close all
% figure('visible','off')
% subplot(2,3,1)
% hold on
% grid on
% plot(wdat1)
% plot(wdat2,'-r')
% axis tight
% subplot(2,3,4)
% hold on
% grid on
% plot(sdat1)
% plot(sdat2,'-r')
% axis tight
% subplot(2,3,2)
% hold on
% grid on
% plot(x1,runmean(wdat1,ls(ic)))
% plot(x2,runmean(wdat2,ls(ic)),'-r')
% % wdat1_sm=ksr(1:length(wdat1),wdat1,median(unique(diff(1:length(wdat1))))*ls(ic),length(wdat1));
% % wdat2_sm=ksr(1:length(wdat2),wdat2,median(unique(diff(1:length(wdat2))))*ls(ic),length(wdat2));
% % plot(wdat1_sm.f,'-b','linewidth',2)
% % plot(wdat2_sm.f,'-r','linewidth',2)
% wdat1_sm0=my_runmean(wdat1,ls(ic),'zero');
% plot(x1,wdat1_sm0,':g')
% plot(x2,my_runmean(wdat2,ls(ic),'zero'),'--g')
% axis tight
% subplot(2,3,5)
% hold on
% grid on
% plot(x1,runmean(sdat1,ls(ic)))
% plot(x2,runmean(sdat2,ls(ic)),'-r')
% % sdat1_sm=ksr(1:length(sdat1),sdat1,median(unique(diff(1:length(sdat1))))*ls(ic),length(sdat1));
% % sdat2_sm=ksr(1:length(sdat2),sdat2,median(unique(diff(1:length(sdat2))))*ls(ic),length(sdat2));
% % plot(sdat1_sm.f,'-b','linewidth',2)
% % plot(sdat2_sm.f,'-r','linewidth',2)
% sdat1_sm0=my_runmean(sdat1,ls(ic),[],[mean(sdat1(1:ixl2)) mean(sdat1(end-ixl2:end))]);
% plot(x1,sdat1_sm0,':g')
% plot(x2,my_runmean(sdat2,ls(ic),[],[mean(sdat2(1:ixl2)) mean(sdat2(end-ixl2:end))]),'--g')
% axis tight
% set(gcf,'visible','on')
% subplot(2,3,3)
% hold on
% grid on
% plot(wdat1)
% % plot(wdat1_sm.f,'-r')
% plot(wdat1_sm0,'--r')
% axis tight
% subplot(2,3,6)
% hold on
% grid on
% plot(sdat1)
% % plot(sdat1_sm.f,'-r')
% plot(sdat1_sm0,'--r')
% axis tight
% %%%%% DELETE
% covlog_temp=cov(WDAT(:,N2ix+1),SDAT(:,N2ix+1));covlog(ia)=covlog_temp(2);
% if DAT{5}(DAT{13}(1,ii))==3 %QSENS [default unit: Wm-2]
% rho_d=median(DataOutput{ii,5}(ix1(ia):ix2(ia)));
% fluxlog(ia)=covlog(ia).*Cp.*rho_d;
% elseif DAT{5}(DAT{13}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4) [user-defined unit]
% c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
% scale_mult=unitscale(DAT{6}(5,ixes(2,ii))+1);
% time_mult=secday(DAT{6}(4,ixes(2,ii))-5);
% fluxlog(ia)=(covlog(ia).*c_ds).*(10.^scale_mult).*time_mult;
% elseif DAT{5}(DAT{13}(1,ii))==6 %QLAT [default unit: Wm-2]
% c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
% fluxlog(ia)=covlog(ia).*L.*c_ds;
% end
% %%%%%%
WDAT(:,N2ix+2:end)=[];
SDAT(:,N2ix+2:end)=[];
[freq_raw_nat,spectra_raw_nat,~,~] = Ogive_calcCospectra(hz,h,mU2,WDAT,SDAT,0,'none','none',0); %a tiny amount of spectral energy is invariably lost in this operation
spc1=flipud(log(freq_raw_nat));
spc2=flipud(bsxfun(@times,spectra_raw_nat,freq_raw_nat));
Ogive=-flipud(cumtrapz(spc1,spc2,1));
cc=cov(WDAT(:,1),SDAT(:,1));
Flux_lindet(ia)=cc(2);
% Flux_lindet0(ia,1)=cc(2);
% Flux_lindet0(ia,2)=Ogive(ixval,1);
% Flux_lindet0(ia,3)=trapz(freq_raw_nat,spectra_raw_nat(:,1));
if DAT{5}(DAT{14}(1,ii))==3 %QSENS [default unit: Wm-2]
rho_d=median(DataOutput{ii,5}(ix1(ia):ix2(ia)));
Ogive=Ogive.*Cp.*rho_d;
Flux_lindet(ia)=Flux_lindet(ia).*Cp.*rho_d;
elseif DAT{5}(DAT{14}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4) [user-defined unit]
c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
scale_mult=unitscale(DAT{6}(5,ixes(2,ii))+1);
time_mult=secday(DAT{6}(4,ixes(2,ii))-5);
Ogive=(Ogive.*c_ds).*(10.^scale_mult).*time_mult;
Flux_lindet(ia)=(Flux_lindet(ia).*c_ds).*(10.^scale_mult).*time_mult;
elseif DAT{5}(DAT{14}(1,ii))==6 %QLAT [default unit: Wm-2]
c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
Ogive=Ogive.*L.*c_ds;
Flux_lindet(ia)=Flux_lindet(ia).*L.*c_ds;
end
Ogive_MAT0{ia}=interp1q(freq_raw_nat,Ogive,fr_nat)'; %downscale Ogive
ixnonan=find(~isnan(Ogive_MAT0{ia}(1,:)),1,'first');
Ogive_end{ia}=Ogive_MAT0{ia}(:,ixnonan);
l_Ogive_end(ia)=length(Ogive_end{ia});
% BLOCK ENSEMBLE AVERAGE (BEA)
if SecOutput{10,2}{ii}{2}(ia)==1
BEA_STATS=SecOutput{10,2}{ii}{3}{ia};
ixok=find(BEA_STATS(5,:)==1);
Flux_BEA=nan(3,size(BEA_STATS,2));
mwdat=mean(wdat);
msdat=mean(sdat);
for il=ixok
ix_int=round(linspace(1,length(wdat),BEA_STATS(1,il)+1));
Flux_int=nan(BEA_STATS(1,il),1);
Flux2_int=nan(BEA_STATS(1,il),1);
mwdat_int=nan(BEA_STATS(1,il),1);
msdat_int=nan(BEA_STATS(1,il),1);
for il2=1:BEA_STATS(1,il)
wdat2=wdat(ix_int(il2):ix_int(il2+1));
sdat2=sdat(ix_int(il2):ix_int(il2+1));
Flux_temp=cov(detrend(wdat2,'constant'),detrend(sdat2,'constant'));
mwdat_int(il2)=mean(wdat2)-mwdat;
msdat_int(il2)=mean(sdat2)-msdat;
if DAT{5}(DAT{14}(1,ii))==3 %QSENS [default unit: Wm-2]
rho_d=DataOutput{ii,5}(ix1(ia):ix2(ia));
rho_d=median(rho_d(ix_int(il2):ix_int(il2+1)));
Flux_int(il2)=Flux_temp(2).*Cp.*rho_d;
Flux2_int(il2)=mwdat_int(il2)*msdat_int(il2).*Cp.*rho_d;
elseif DAT{5}(DAT{14}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4) [user-defined unit]
c_ds=DataOutput{ii,4}(ix1(ia):ix2(ia));
c_ds=median(c_ds(ix_int(il2):ix_int(il2+1)));
scale_mult=unitscale(DAT{6}(5,ixes(2,ii))+1);
time_mult=secday(DAT{6}(4,ixes(2,ii))-5);
Flux_int(il2)=(Flux_temp(2).*c_ds).*(10.^scale_mult).*time_mult;
Flux2_int(il2)=(mwdat_int(il2)*msdat_int(il2).*c_ds).*(10.^scale_mult).*time_mult;
elseif DAT{5}(DAT{14}(1,ii))==6 %QLAT [default unit: Wm-2]
c_ds=DataOutput{ii,4}(ix1(ia):ix2(ia));
c_ds=median(c_ds(ix_int(il2):ix_int(il2+1)));
Flux_int(il2)=Flux_temp(2).*L.*c_ds;
Flux2_int(il2)=mwdat_int(il2)*msdat_int(il2).*L.*c_ds;
end
end
Flux_BEA(1,il)=mean(Flux2_int);
Flux_BEA(2,il)=mean(Flux_int);
Flux_BEA(3,il)=mean(Flux2_int)+mean(Flux_int);
end
BEA_STATS_LOG{ia}=[BEA_STATS;Flux_BEA];
end
end
% FIG5=figure;
% set(gcf,'position',get(0,'screensize'),'color','w','Visible','off')
% set(gca,'box','on')
% subplot(3,1,1:2)
% hold on
% grid on
% for iac=1:N1
% dl=datalength(iac)/(60*hz);
% if SecOutput{10,2}{ii}{1}==0
% plot(dl,Flux_lindet(iac),'*','color',[.5 .5 .5])
% else
% plot(dl,Flux_lindet(iac),'*r')
% end
% if ~isempty(BEA_STATS_LOG{iac})
% stdflux=nanstd(BEA_STATS_LOG{iac}(8,:));
% medflux=nanmedian(BEA_STATS_LOG{iac}(8,:));
% plot(dl,medflux,'.b')
% plot([dl dl],[medflux medflux-stdflux],'-b')
% plot([dl dl],[medflux medflux+stdflux],'-b')
% plot(dl,medflux-stdflux,'.b')
% plot(dl,medflux+stdflux,'.b')
% end
% end
% % xlabel('Data length [Min]')
% ylabel([DAT{8}{DAT{13}(3,ii),6},' [Wm^{-2}]'])
% title('Red/gray = stationary/instationary EC_T | Blue = BEA_{TX}')
% ax=axis;
% plot([ax(1) ax(2)],[0 0],'-k')
% set(gca,'box','on')
% subplot(3,1,3)
% hold on
% grid on
% for iac=1:N1
% dl=datalength(iac)/(60*hz);
% if ~isempty(BEA_STATS_LOG{iac})
% stdflux=nanstd(BEA_STATS_LOG{iac}(6,:));
% medflux=nanmedian(BEA_STATS_LOG{iac}(6,:));
% plot(dl,medflux,'.k')
% plot([dl dl],[medflux medflux-stdflux],'-k')
% plot([dl dl],[medflux medflux+stdflux],'-k')
% plot(dl,medflux-stdflux,'.k')
% plot(dl,medflux+stdflux,'.k')
%
% stdflux=nanstd(BEA_STATS_LOG{iac}(7,:));
% medflux=nanmedian(BEA_STATS_LOG{iac}(7,:));
% plot(dl,medflux,'.b')
% plot([dl dl],[medflux medflux-stdflux],'-b')
% plot([dl dl],[medflux medflux+stdflux],'-b')
% plot(dl,medflux-stdflux,'.b')
% plot(dl,medflux+stdflux,'.b')
% end
% end
% set(gca,'xlim',ax(1:2))
% xlabel('Data length [Min]')
% ylabel([DAT{8}{DAT{13}(3,ii),6},' [Wm^{-2}]'])
% title('Black=BEA_{TX} term#1 | Blue=BEA_{TX} term#2 (standard blockwise covariance)')
% plot([ax(1) ax(2)],[0 0],'-k')
% set(gca,'box','on')
% path=[datestr(DAT15(3),'yyyy-mm-dd HH-MM-SS'),' -- ',num2str(DAT15(1)),'-',num2str(DAT15(2)),'.fig'];
% if ~(exist([figurepath2,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep],'dir')==7)
% mkdir([figurepath2,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep])
% end
% hgsave([figurepath2,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep,path])
% close(FIG5)
Flux_lindet(SecOutput{10,2}{ii}{1}==0)=nan;
if ~isnan(ixm30)
Ogive30{ii}=Ogive_MAT0{ixm30}(1,:);
end
Ogive_end_prc=prctile(cell2mat(Ogive_end),[2 13.6 50 86.4 98]);
ml=max(l_Ogive_end);
for ia=1:N1
Ogive_end_lindet(ia)=Ogive_end{ia}(end);
if length(Ogive_end{ia})~=ml
Ogive_end{ia}=[Ogive_end{ia};nan(ml-length(Ogive_end{ia}),1)];
end
end
Ogive_MAT=cell2mat(Ogive_MAT0(:));
elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
l_Ogive_end=nan(int,1);
lsrange=[-1.4 0]; %-2.2
ls0=fliplr((10^lsrange(2))-(10.^linspace(lsrange(1),lsrange(2),N2+1)));
ls=round(ls0*(tmax/2));
ls=ls(2:end);
rls=myrange(ls);
ls=round(interp1(rls,[rls(1) floor(length(DataOutput{ii,1})/2)],ls));
Ogive_MAT0F=cell(1,int);
P=polyfit([mintime tmax/(hz*60)],[.4 1],1);
factor=polyval(P,lsp);
maxtime=lsp.*factor;
maxfreq=1./(maxtime*60);
%determine which ia indeces to use (the ones which actually make a difference)
% freqlog=false(N1,1);
% for ia=1:N1
% if ia==1
% freqlog(ia)=true;
% ixnonan_old=find(fr_nat>=maxfreq(ia),1,'first');
% else
% ixnonan_new=find(fr_nat>=maxfreq(ia),1,'first');
% if ixnonan_new<ixnonan_old
% freqlog(ia)=true;
% end
% ixnonan_old=ixnonan_new;
% end
% end
% iaok=unique([find(freqlog==1)' ixm30]);
%loop
ix_tail=find(lsf>lsf(1)+2,1,'first');
ctr=0;
for ia=1:N1
% for ia=iaok
ctr=ctr+1;
wdat=rawdat1(ix1(ia):ix2(ia));
sdat=rawdat2(ix1(ia):ix2(ia));
WDAT=nan(length(ix1(ia):ix2(ia)),N2);
SDAT=nan(length(ix1(ia):ix2(ia)),N2);
mU2=mean(DataOutput{ii,6}(ix1(ia):ix2(ia)));
N2ix=find(ls<floor(length(wdat)/2),1,'last');
WDAT(:,1)=detrend(wdat);
SDAT(:,1)=detrend(sdat);
if ~isempty(N2ix)
for ic=1:N2ix
WDAT(:,ic+1)=wdat-my_runmean(wdat,ls(ic),'zero');
SDAT(:,ic+1)=sdat-my_runmean(sdat,ls(ic),[],[mean(sdat(1:ixl2)) mean(sdat(end-ixl2:end))]);
end
else
N2ix=0;
end
% ia=1;
% wdat1=rawdat1(ix1(ia):ix2(ia));
% sdat1=rawdat2(ix1(ia):ix2(ia));
% ia=2;
% wdat2=rawdat1(ix1(ia):ix2(ia));
% sdat2=rawdat2(ix1(ia):ix2(ia));
% x1=1:length(wdat1);
% x2=1:length(wdat2);
% ixl=length(wdat1)-length(wdat2);
% ixl2=round(ixl*1.5);
% x2=x2+round(ixl/2);
% ic=2;
% close all
% figure('visible','off')
% subplot(2,3,1)
% hold on
% grid on
% plot(wdat1)
% plot(wdat2,'-r')
% axis tight
% subplot(2,3,4)
% hold on
% grid on
% plot(sdat1)
% plot(sdat2,'-r')
% axis tight
% subplot(2,3,2)
% hold on
% grid on
% plot(x1,runmean(wdat1,ls(ic)))
% plot(x2,runmean(wdat2,ls(ic)),'-r')
% % wdat1_sm=ksr(1:length(wdat1),wdat1,median(unique(diff(1:length(wdat1))))*ls(ic),length(wdat1));
% % wdat2_sm=ksr(1:length(wdat2),wdat2,median(unique(diff(1:length(wdat2))))*ls(ic),length(wdat2));
% % plot(wdat1_sm.f,'-b','linewidth',2)
% % plot(wdat2_sm.f,'-r','linewidth',2)
% wdat1_sm0=my_runmean(wdat1,ls(ic),'zero');
% plot(x1,wdat1_sm0,':g')
% plot(x2,my_runmean(wdat2,ls(ic),'zero'),'--g')
% axis tight
% subplot(2,3,5)
% hold on
% grid on
% plot(x1,runmean(sdat1,ls(ic)))
% plot(x2,runmean(sdat2,ls(ic)),'-r')
% % sdat1_sm=ksr(1:length(sdat1),sdat1,median(unique(diff(1:length(sdat1))))*ls(ic),length(sdat1));
% % sdat2_sm=ksr(1:length(sdat2),sdat2,median(unique(diff(1:length(sdat2))))*ls(ic),length(sdat2));
% % plot(sdat1_sm.f,'-b','linewidth',2)
% % plot(sdat2_sm.f,'-r','linewidth',2)
% sdat1_sm0=my_runmean(sdat1,ls(ic),[],[mean(sdat1(1:ixl2)) mean(sdat1(end-ixl2:end))]);
% plot(x1,sdat1_sm0,':g')
% plot(x2,my_runmean(sdat2,ls(ic),[],[mean(sdat2(1:ixl2)) mean(sdat2(end-ixl2:end))]),'--g')
% axis tight
% set(gcf,'visible','on')
% subplot(2,3,3)
% hold on
% grid on
% plot(wdat1)
% % plot(wdat1_sm.f,'-r')
% plot(wdat1_sm0,'--r')
% axis tight
% subplot(2,3,6)
% hold on
% grid on
% plot(sdat1)
% % plot(sdat1_sm.f,'-r')
% plot(sdat1_sm0,'--r')
% axis tight
% %%%%% DELETE
% covlog_temp=cov(WDAT(:,N2ix+1),SDAT(:,N2ix+1));covlog(ia)=covlog_temp(2);
% if DAT{5}(DAT{13}(1,ii))==3 %QSENS [default unit: Wm-2]
% rho_d=median(DataOutput{ii,5}(ix1(ia):ix2(ia)));
% fluxlog(ia)=covlog(ia).*Cp.*rho_d;
% elseif DAT{5}(DAT{13}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4) [user-defined unit]
% c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
% scale_mult=unitscale(DAT{6}(5,ixes(2,ii))+1);
% time_mult=secday(DAT{6}(4,ixes(2,ii))-5);
% fluxlog(ia)=(covlog(ia).*c_ds).*(10.^scale_mult).*time_mult;
% elseif DAT{5}(DAT{13}(1,ii))==6 %QLAT [default unit: Wm-2]
% c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
% fluxlog(ia)=covlog(ia).*L.*c_ds;
% end
% %%%%%%
WDAT(:,N2ix+2:end)=[];
SDAT(:,N2ix+2:end)=[];
[freq_raw_nat,spectra_raw_nat,~,~] = Ogive_calcCospectra(hz,h,mU2,WDAT,SDAT,0,'none','none',0);
spc1=flipud(log(freq_raw_nat));
spc2=flipud(bsxfun(@times,spectra_raw_nat,freq_raw_nat));
Ogive=-flipud(cumtrapz(spc1,spc2,1));
cc=cov(WDAT(:,1),SDAT(:,1));
Flux_lindet(ia)=cc(2);
% Flux_lindet0(ia,1)=cc(2);
% Flux_lindet0(ia,2)=Ogive(ixval,1);
% Flux_lindet0(ia,3)=trapz(freq_raw_nat,spectra_raw_nat(:,1));
if DAT{5}(DAT{14}(1,ii))==3 %QSENS [default unit: Wm-2]
rho_d=median(DataOutput{ii,5}(ix1(ia):ix2(ia)));
Ogive=Ogive.*Cp.*rho_d;
Flux_lindet(ia)=Flux_lindet(ia).*Cp.*rho_d;
elseif DAT{5}(DAT{14}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4) [user-defined unit]
c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
scale_mult=unitscale(DAT{6}(5,ixes(2,ii))+1);
time_mult=secday(DAT{6}(4,ixes(2,ii))-5);
Ogive=(Ogive.*c_ds).*(10.^scale_mult).*time_mult;
Flux_lindet(ia)=(Flux_lindet(ia).*c_ds).*(10.^scale_mult).*time_mult;
elseif DAT{5}(DAT{14}(1,ii))==6 %QLAT [default unit: Wm-2]
c_ds=median(DataOutput{ii,4}(ix1(ia):ix2(ia)));
Ogive=Ogive.*L.*c_ds;
Flux_lindet(ia)=Flux_lindet(ia).*L.*c_ds;
end
Ogive_downscaled=interp1q(freq_raw_nat,Ogive,fr_nat)'; %downscale Ogive
Ogive_downscaled(:,fr_nat<maxfreq(ia))=nan;
if ia==1
% freqlog2(ia)=true;
ixnonan_old=find(~isnan(Ogive_downscaled(1,:)),1,'first');
for ic=ixnonan_old:int
Ogive_MAT0F{ic}=Ogive_downscaled(:,ic);
l_Ogive_end(ic)=length(Ogive_MAT0F{ic});
end
else
ixnonan_new=find(~isnan(Ogive_downscaled(1,:)),1,'first');
% if ixnonan_new<ixnonan_old
% freqlog2(ia)=true;
% for ic=ixnonan_new:ixnonan_old-1
% Ogive_MAT0F{ic}=Ogive_downscaled(:,ic);
% l_Ogive_end(ic)=length(Ogive_MAT0F{ic});
% end
xx=linspace(1,4,ix_tail+1);
pf=polyfit([1 4],[0 50],1);
yy2=polyval(pf,xx);
prcs=[yy2;100-yy2];
ctr=0;
for ic=ixnonan_new:ixnonan_new+ix_tail;%(int-ixnonan_0)
ctr=ctr+1;
prcs_ic=prctile(Ogive_downscaled(:,ic),prcs(:,ctr)');
if any(Ogive_downscaled(:,ic)>=prcs_ic(1) & Ogive_downscaled(:,ic)<=prcs_ic(2))
ixok=Ogive_downscaled(:,ic)>=prcs_ic(1) & Ogive_downscaled(:,ic)<=prcs_ic(2);
Ogive_MAT0F{ic}=[Ogive_MAT0F{ic};Ogive_downscaled(ixok,ic)];
end
l_Ogive_end(ic)=length(Ogive_MAT0F{ic});
end
% for ic=ixnonan_new:ixnonan_new+ix_tail;%(int-ixnonan_0)
% Ogive_MAT0F{ic}=[Ogive_MAT0F{ic};Ogive_downscaled(:,ic)];
% l_Ogive_end(ic)=length(Ogive_MAT0F{ic});
% end
% end
% ixnonan_old=ixnonan_new;
end
Ogive_MAT0{ia}=Ogive_downscaled; %downscale Ogive
end
Flux_lindet(SecOutput{10,2}{ii}{1}==0)=nan;
if ~isnan(ixm30)
Ogive30{ii}=Ogive_MAT0{ixm30}(1,:);
end
Ogive_end_prc=prctile(Ogive_MAT0F{1},[2 13.6 50 86.4 98]);
ml=max(l_Ogive_end);
Ogive_MAT=nan(ml,int);
for ia=1:int
if length(Ogive_MAT0F{ia})~=ml
Ogive_MAT(:,ia)=[Ogive_MAT0F{ia};nan(ml-length(Ogive_MAT0F{ia}),1)];
else
Ogive_MAT(:,ia)=Ogive_MAT0F{ia};
end
end
end
% FORMATION OF MASS OGIVE DENSITY FIELD
medO=nanmedian(Ogive_MAT);
prcO=cell(3,1);
prcO{1} = prctile(Ogive_MAT,[13.6 100-13.6]);
prcO{2} = prctile(Ogive_MAT,[25 75]);
prcO{3} = prctile(Ogive_MAT,[35 65]);
yax{ii} = prctile(Ogive_MAT,[1 99]);
%determine upper and lower limit of the frequency ranges (FR)
FR0=nan(1,3);
FR0(1)=lsf(1);
FR0(2)=HFbound;
if DAT{2}(5)>1 %signal from closed path IRGA
medO2=myrunmean(medO,3,1,nan);
accrange=max(abs(medO2))/((1/3)*100);
ixok=find(medO2>accrange | medO2<-accrange,1,'last');
lsfix=lsf(ixok);
FR0(3)=max(0,lsfix);
else
FR0(3)=.7; %fixed
end
NUM_INTVL=DAT{1}(10);%5;
FR=nan(NUM_INTVL,2);
FR(1,1)=FR0(1);
FR(end,1)=FR0(2);
FR(:,2)=FR0(3);
intvl0=diff(FR0(1:2))/8;
[~,ixl1]=min(abs(lsf-(FR0(1)+intvl0/2)));
[~,ixl2]=min(abs(lsf-(FR0(2)-intvl0/2)));
% ixc=5:8:70;
% lvsm=length(ixc);
% figure
% ctr=0;
% ymult=1.2;
% for ikn=ixc
% ctr=ctr+1;
%
% subaxis(3,lvsm,ctr,'Spacinghoriz', 0.01,'Spacingvert',0.025, 'Padding',0, 'Margin',.03);
% hold on
% grid on
% medO_sm=Smooth(medO,ikn);
% yi1_sm=Smooth(prcO{2}(1,:),ikn);
% yi2_sm=Smooth(prcO{2}(2,:),ikn);
% plot(lsf,medO_sm,'color',[.7 .7 .7])
% plot(lsf,yi1_sm,'color',[.7 .7 .7])
% plot(lsf,yi2_sm,'color',[.7 .7 .7])
% plot(lsf(ixl1:ixl2),medO_sm(ixl1:ixl2),'-r','linewidth',2)
% plot(lsf(ixl1:ixl2),yi1_sm(ixl1:ixl2),'-m')
% plot(lsf(ixl1:ixl2),yi2_sm(ixl1:ixl2),'-m')
% title(num2str(ikn))
% axis tight
% axis square
% ax=axis;
% plot([lsf(ixl1) lsf(ixl1)],[ax(3) ax(4)],'-g')
% plot([lsf(ixl2) lsf(ixl2)],[ax(3) ax(4)],'-g')
% plot([FR0(1) FR0(1)],[ax(3) ax(4)],'-k')
% plot([FR0(2) FR0(2)],[ax(3) ax(4)],'-k')
% set(gca,'box','on','xticklabel',[],'yticklabel',[])
%
% subaxis(3,lvsm,ctr+lvsm,'Spacinghoriz', 0.01,'Spacingvert',0.025, 'Padding',0, 'Margin',.03);
% hold on
% grid on
% dx=gradient(lsf);
% dy_med=gradient(medO_sm);
% dy_y1=gradient(yi1_sm);
% dy_y2=gradient(yi2_sm);
% curv_med = gradient(atan2(dy_med,dx)) ./ hypot(dx,dy_med);
% curv_y1 = gradient(atan2(dy_y1,dx)) ./ hypot(dx,dy_y1);
% curv_y2 = gradient(atan2(dy_y2,dx)) ./ hypot(dx,dy_y2);
% medc=median(curv_med(3:end-3));
% stdc=std(curv_med(3:end-3));
% plot(lsf,curv_med,'color',[.7 .7 .7])
% plot(lsf,curv_y1,'color',[.7 .7 .7])
% plot(lsf,curv_y2,'color',[.7 .7 .7])
% plot(lsf(ixl1:ixl2),curv_med(ixl1:ixl2),'-r','linewidth',2)
% plot(lsf(ixl1:ixl2),curv_y1(ixl1:ixl2),'-m')
% plot(lsf(ixl1:ixl2),curv_y2(ixl1:ixl2),'-m')
% plot([lsf(1) lsf(end)],[medc medc],'-r')
% plot([lsf(1) lsf(end)],[medc-stdc medc-stdc],'--b')
% plot([lsf(1) lsf(end)],[medc+stdc medc+stdc],'--b')
% plot([lsf(1) lsf(end)],[0 0],'-k')
% axis tight
% axis square
% ylm=[0 max([max(abs(curv_med(ixl1:ixl2)));max(abs(curv_y1(ixl1:ixl2)));max(abs(curv_y2(ixl1:ixl2)));max(abs(medc+stdc));max(abs(medc-stdc))])*ymult];
% set(gca,'ylim',[-max(abs(ylm)) max(abs(ylm))])
% ax=axis;
% plot([lsf(ixl1) lsf(ixl1)],[ax(3) ax(4)],'-g')
% plot([lsf(ixl2) lsf(ixl2)],[ax(3) ax(4)],'-g')
% plot([FR0(1) FR0(1)],[ax(3) ax(4)],'-k')
% plot([FR0(2) FR0(2)],[ax(3) ax(4)],'-k')
% set(gca,'box','on','xticklabel',[],'yticklabel',[])
% title('curvature')
%
% subaxis(3,lvsm,ctr+(2*lvsm),'Spacinghoriz', 0.01,'Spacingvert',0.025, 'Padding',0, 'Margin',.03);
% hold on
% grid on
% d1med=Smooth(diff(curv_med),2);
% d1y1=Smooth(diff(curv_y1),2);
% d1y2=Smooth(diff(curv_y2),2);
% plot(lsf(1:end-1),d1med,'color',[.7 .7 .7])
% plot(lsf(1:end-1),d1y1,'color',[.7 .7 .7])
% plot(lsf(1:end-1),d1y2,'color',[.7 .7 .7])
% plot(lsf(ixl1:ixl2),d1med(ixl1:ixl2),'-r','linewidth',2)
% plot(lsf(ixl1:ixl2),d1y1(ixl1:ixl2),'-m')
% plot(lsf(ixl1:ixl2),d1y2(ixl1:ixl2),'-m')
% plot([lsf(1) lsf(end)],[0 0],'-k')
% axis tight
% axis square
% ylm=[0 max([max(abs(d1med(ixl1:ixl2)));max(abs(d1y1(ixl1:ixl2)));max(abs(d1y2(ixl1:ixl2)))])*ymult];
% set(gca,'ylim',[-max(abs(ylm)) max(abs(ylm))])
% ax=axis;
% plot([lsf(ixl1) lsf(ixl1)],[ax(3) ax(4)],'-g')
% plot([lsf(ixl2) lsf(ixl2)],[ax(3) ax(4)],'-g')
% plot([FR0(1) FR0(1)],[ax(3) ax(4)],'-k')
% plot([FR0(2) FR0(2)],[ax(3) ax(4)],'-k')
% set(gca,'box','on','xticklabel',[],'yticklabel',[])
% title('diff^1')
%
% % subaxis(4,lvsm,ctr+(3*lvsm),'Spacinghoriz', 0.01,'Spacingvert',0.025, 'Padding',0, 'Margin',.03);
% % hold on
% % grid on
% % d2med=Smooth(diff(d1med),8);
% % d2y1=Smooth(diff(d1y1),8);
% % d2y2=Smooth(diff(d1y2),8);
% % plot(lsf(2:end-1),d2med,'color',[.7 .7 .7])
% % plot(lsf(2:end-1),d2y1,'color',[.7 .7 .7])
% % plot(lsf(2:end-1),d2y2,'color',[.7 .7 .7])
% % plot(lsf(ixl1+1:ixl2+1),d2med(ixl1:ixl2))
% % plot(lsf(ixl1+1:ixl2+1),d2y1(ixl1:ixl2),'-m')
% % plot(lsf(ixl1+1:ixl2+1),d2y2(ixl1:ixl2),'-m')
% % plot([lsf(1) lsf(end)],[0 0],'-k')
% % axis tight
% % axis square
% % ylm=[0 max([max(abs(d2med(ixl1:ixl2)));max(abs(d2y1(ixl1:ixl2)));max(abs(d2y2(ixl1:ixl2)))])*ymult];
% % set(gca,'ylim',[-max(abs(ylm)) max(abs(ylm))])
% % ax=axis;
% % plot([lsf(ixl1) lsf(ixl1)],[ax(3) ax(4)],'-g')
% % plot([lsf(ixl2) lsf(ixl2)],[ax(3) ax(4)],'-g')
% % plot([FR0(1) FR0(1)],[ax(3) ax(4)],'-k')
% % plot([FR0(2) FR0(2)],[ax(3) ax(4)],'-k')
% % set(gca,'box','on','xticklabel',[],'yticklabel',[])
% % title('diff^2')
% end
diff1zero_log=zeros(1,size(lsf,1));
ixc=1:30;
for ikn=ixc
dx=gradient(lsf);
medO_sm=ksr(1:length(medO),medO,median(unique(diff(1:length(medO))))*1,length(medO));
% medO_sm=Smooth(medO,ikn);
NUMdat=1+(length(prcO)*2);
curv=cell(NUMdat,1);
d1=cell(NUMdat,1);
dy=cell(NUMdat,1);
medyi=cell(NUMdat,1);
medyi{1}=medO_sm.f;
ctr=0;
for ikn2=1:length(prcO)
for ikn3=1:2
ctr=ctr+1;
medyi{ctr+1}=mySmooth([],prcO{ikn2}(ikn3,:),ikn,[]);
% yy=prcO{ikn2}(ikn3,:);
% medO_sm_ksr=ksr(1:length(yy),yy,median(unique(diff(1:length(yy))))*ikn,length(yy));
% medyi{ctr+1}=medO_sm_ksr.f;
% medyi{ctr+1}=Smooth(prcO{ikn2}(ikn3,:),ikn);
end
end
for ikn2=1:NUMdat
dy{ikn2}=gradient(medyi{ikn2});
curv{ikn2} = gradient(atan2(dy{ikn2}(:),dx)) ./ hypot(dx,dy{ikn2}(:));
%Diff^1
d1{ikn2}=mySmooth([],diff(curv{ikn2}),2,[]);
% yy=diff(curv{ikn2});
% medO_sm_ksr=ksr(1:length(yy),yy,median(unique(diff(1:length(yy))))*2,length(yy));
% d1{ikn2}=medO_sm_ksr.f;
% d1{ikn2}=Smooth(diff(curv{ikn2}),2);
d1_sgn=sign(d1{ikn2});
d1_sgnfilt=filter([.5 .5],1,d1_sgn);
if any(d1_sgnfilt==0)
ixzero=find(d1_sgnfilt==0);
ixzero_ok=find(ixzero>ixl1 & ixzero<ixl2);
if ~isempty(ixzero_ok)
diff1zero_log(ixzero(ixzero_ok))=diff1zero_log(ixzero(ixzero_ok))+sign(curv{ikn2}(ixzero(ixzero_ok)))';
end
end
end
end
smfac=3;
diff1zero_log_sm=mySmooth([],diff1zero_log,smfac,[]);
% yy=diff1zero_log;
% diff1zero_log_sm=ksr(1:length(yy),yy,median(unique(diff(1:length(yy))))*smfac,length(yy));
% diff1zero_log_sm=diff1zero_log_sm.f;
% diff1zero_log_sm=Smooth(diff1zero_log,smfac);
diff2zero_log_sm=diff(diff1zero_log_sm,1);
diff3zero_log_sm=mySmooth([],diff(diff2zero_log_sm),smfac,[]);
% yy=diff(diff2zero_log_sm);
% diff3zero_log_sm=ksr(1:length(yy),yy,median(unique(diff(1:length(yy))))*smfac,length(yy));
% diff3zero_log_sm=diff3zero_log_sm.f;
% diff3zero_log_sm=Smooth(diff(diff2zero_log_sm),smfac);
% figure
% hold on
% plot(lsf,diff1zero_log_sm,'-b')
% plot(lsf(1:end-1),diff2zero_log_sm,'-r')
% plot(lsf(2:end-1),diff3zero_log_sm,'-m')
% plot(lsf,diff1zero_log./max(abs(diff1zero_log))*25,'-g')
% ax=axis;
% plot([ax(1) ax(2)],[0 0],'-k')
% grid on
logok=true(1,length(lsf));
logok(lsf<lsf(ixl1+1) | lsf>lsf(ixl2-1))=false;
diff1zero_log_sm(~logok)=0;
diff3zero_log_sm(~logok(2:end-1))=0;
diff1zero_log_sm_med=nanmedian(diff1zero_log_sm(logok));
diff1zero_log_sm(~logok)=diff1zero_log_sm_med;
diff1zero_log_sm_nonan=diff1zero_log_sm(logok);diff1zero_log_sm_nonan(isnan(diff1zero_log_sm_nonan))=[]; %avoid use of nanstd below which requires the stats toolbox
diff1zero_log_sm_std=std(diff1zero_log_sm_nonan)*1.5;
ix_above_std=find(diff1zero_log_sm>=(diff1zero_log_sm_med+diff1zero_log_sm_std));
ix_below_std=find(diff1zero_log_sm<=(diff1zero_log_sm_med-diff1zero_log_sm_std));
if DAT{12}(10,2)>0
HF=figure('Visible','off');
subplot(3,3,[1 2])
set(gcf,'color','w','position',get(0,'screensize'))
hold on
grid on
plot(lsf,diff1zero_log)
plot(lsf,diff1zero_log_sm,'-r')
axis tight
ax=axis;
ax(1)=ax(1)-.15;
set(gca,'xlim',ax(1:2))
plot([ax(1) ax(2)],[0 0],'-k')
plot([ax(1) ax(2)],[diff1zero_log_sm_med-diff1zero_log_sm_std diff1zero_log_sm_med-diff1zero_log_sm_std],'--r')
plot([ax(1) ax(2)],[diff1zero_log_sm_med diff1zero_log_sm_med],'--r')
plot([ax(1) ax(2)],[diff1zero_log_sm_med+diff1zero_log_sm_std diff1zero_log_sm_med+diff1zero_log_sm_std],'--r')
xtx=get(gca,'xtick');
xtxlab=cell(length(xtx),1);
for ih=1:length(xtx)
xtxlab{ih}=['10^{',num2str(xtx(ih)),'}'];
end
set(gca,'fontsize',12,'fontweight','bold','box','on','xticklabel',xtxlab)
xlabel('Natural frequency','fontsize',14,'fontweight','bold')
ylabel('Ogive curvature analysis','fontsize',14,'fontweight','bold')
end
for ikn=1:NUM_INTVL-2
if ~isempty(ix_above_std) || ~isempty(ix_below_std)
[~,mnix]=max(abs(diff1zero_log_sm-diff1zero_log_sm_med));
dsm_sgnfilt3=filter([.5 .5],1,sign(diff3zero_log_sm));
dxzero3=unique([ixl1;find(dsm_sgnfilt3==0 & logok(2:end-1)');ixl2]);
mnix2=dxzero3(find(dxzero3<mnix,1,'last'));
if isempty(mnix2)
FR(ikn+1,1)=lsf(1);
else
FR(ikn+1,1)=lsf(mnix2);
end
if ikn<(NUM_INTVL-2)
dsm=diff(diff1zero_log_sm);
dsm_sgnfilt=filter([.5 .5],1,sign(dsm));
dxzero=[ixl1;find(dsm_sgnfilt==0 & logok(1:end-1)');ixl2];
[~,mnix_zero]=min(abs(dxzero-mnix));
% figure
% hold on
% plot(lsf,diff1zero_log_sm,'-r')
% axis tight
% ax=axis;
% grid on
% for iku=1:length(dxzero)
% plot([lsf(dxzero(iku)) lsf(dxzero(iku))],[ax(3) ax(4)],'--k')
% end
diff1zero_log_sm(dxzero(max(1,mnix_zero-1)):dxzero(min(length(dxzero),mnix_zero+1)))=diff1zero_log_sm_med;
diff3zero_log_sm(dxzero(max(1,mnix_zero-1)):dxzero(min(length(dxzero),mnix_zero+1)))=0;
ix_above_std=find(diff1zero_log_sm>=(diff1zero_log_sm_med+diff1zero_log_sm_std));
ix_below_std=find(diff1zero_log_sm<=(diff1zero_log_sm_med-diff1zero_log_sm_std));
end
end
end
if any(isnan(FR(:,1)))
ixnan=isnan(FR(:,1));
FR(ixnan,:)=[];
end
sizFR1=size(FR,1);
if DAT{12}(10,2)>0
for ih=1:sizFR1
plot([FR(ih,1) FR(ih,1)],[ax(3) ax(4)],'--k','linewidth',2)
end
end
FR=sortrows(FR);
Ogt=Ogive_MAT(~isnan(Ogive_MAT));
rOgt=[nanmin(Ogt) nanmax(Ogt)];
if DAT{1}(18)==0 %STANDARD METHOD
p=.5;
rlog_Og_r=prctile(Ogt,[p 100-p]);
if sign(rOgt(1))~=sign(rOgt(2))
while sign(rlog_Og_r(1))==sign(rlog_Og_r(2))
p=p*.2;
rlog_Og_r=prctile(Ogt,[p 100-p]);
end
else
rlog_Og_r(abs(rlog_Og_r)==min(abs(rlog_Og_r)))=0;
end
lsOg=linspace(rlog_Og_r(1),rlog_Og_r(2),int);
elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
lsOg=linspace(rOgt(1),rOgt(2),int);
end
hist_scl_x_all=linspace(0,10,int);
hist_scl_all=exp(-(hist_scl_x_all - 6).^2./(2*(2.5)^2));
[mx,mxix]=max(hist_scl_all);
hist_scl_all(1:mxix)=mx;
if DAT{1}(18)==0 %STANDARD METHOD
% hist_scl = gaussmf(linspace(0,10,int),[sig 5]);
% d = [3 2 3 1 1 1];
% w=[.5 .5 1 1 1 1];
% edges=linspace(min(d),max(d),10);
% [count,bin] = histc(d,edges);
% totfreq = accumarray(bin',w',[length(edges) 1]); %4 comes from the number of
% figure %columns in the 'edges' vector
% bar(edges,totfreq, 'histc');
% figure
% ctr=0;
% for sig=.4:.2:10
% ctr=ctr+1;
% ctr=0;
% sigv=2:.4:12;
% nums=nan(length(sigv),1);
% scale=nan(length(sigv),1);
% for sig2=sigv
% ctr=ctr+1;
% hist_scl_x=linspace(0,10,N1);
% hist_scl=exp(-(hist_scl_x - 0).^2./(2*(sig2)^2));
% scale(ctr)=sum(l_Ogive_end)/round(sum(hist_scl.*l_Ogive_end'));
% nums(ctr)=round(sum((scale(ctr)*hist_scl).*l_Ogive_end'));
%
%
%
% % hist_scl=hist_scl*(max(hist_scl_x)/trapz(hist_scl_x,hist_scl));
%
%
%
%
% end
% SCALING BY NORMAL DISTRIBUTION
hist_scl_x=linspace(0,10,N1);
hist_scl=exp(-(hist_scl_x - 0).^2./(2*(sig)^2));
scale=sum(l_Ogive_end)/round(sum(hist_scl.*l_Ogive_end'));
hist_scl=scale*hist_scl;
% hist_scl=hist_scl*(max(hist_scl_x)/trapz(hist_scl_x,hist_scl));
hist_scl_vec=nan(size(Ogive_MAT,1),1);
l_Ogive_end_cs=[1;cumsum(l_Ogive_end)];
for ij=1:N1
hist_scl_vec(l_Ogive_end_cs(ij):l_Ogive_end_cs(ij+1))=hist_scl(ij);
end
MAT=nan(int);
for ik=1:int
MAT(:,ik)=hist_scl_all(ik).*histwv2(Ogive_MAT(:,ik),hist_scl_vec,lsOg);
% MAT(:,ik)=hist_scl(ik).*hist(Ogive_MAT(:,ik),lsOg)';
end
%POSSIBLE FAST SOLUTION. CHECK IF THIS CAN BE SCALED AS WELL:
% MAT=hist(Ogive_MAT,lsOg);
MAT(1,:)=0;
MAT(end,:)=0;
elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
MAT=nan(int);
for ik=1:int
% MAT(:,ik)=hist_scl_all(ik).*histwv2(Ogive_MAT(:,ik),hist_scl_vec,lsOg);
% MAT(:,ik)=hist_scl_all(ik).*hist(Ogive_MAT(~isnan(Ogive_MAT(:,ik)),ik),lsOg)';
MAT(:,ik)=hist(Ogive_MAT(~isnan(Ogive_MAT(:,ik)),ik),lsOg)';
end
%POSSIBLE FAST SOLUTION. CHECK IF THIS CAN BE SCALED AS WELL:
% MAT=hist(Ogive_MAT,lsOg);
MATsum=sum(MAT);
MAT=bsxfun(@times,MAT,(ones(size(MATsum))*mean(MATsum))./MATsum);
MAT=bsxfun(@times,MAT,hist_scl_all);
MAT(isnan(MAT))=0;
end
% OGIVE OPTIMIZATION -->
lsf0=lsf;
ixmed=median(unique(diff(lsf)));
ixmedOg=median(unique(diff(lsOg)));
extendOgperc=20; %percent extra Ogive range
ixmultOg=round((diff(myrange(lsOg))*(extendOgperc/100))/ixmedOg);
if DAT{1}(17)==1
lsf1n=lsf(1)*2;
else
extendfreqLF=1000; %minutes
lsf1n=log10(1/((((1/10^lsf(1))/60)+extendfreqLF)*60));
end
ixmult1=round(diff([lsf1n;lsf(1)])/ixmed);
extendfreqHF=1e5; %hz
extendfreqHF10=log10(extendfreqHF/2);
ixmult2=round(diff([lsf(end);extendfreqHF10])/ixmed);
ixmultLOG(ii,:)=[ixmult1 ixmult2];
medO2=[nan(1,ixmult1),medO,nan(1,ixmult2)];
lsf2=linspace(lsf(1)-(ixmult1*ixmed),lsf(end)+(ixmult2*ixmed),ixmult1+length(lsf)+ixmult2);
fr_nat=(10.^lsf2)';
MAT2=[zeros(int,ixmult1) MAT zeros(int,ixmult2)];
if ~isnan(ixm30)
Ogive302{ii}=[nan(1,ixmult1) Ogive30{ii} nan(1,ixmult2)];
end
int2=ixmult1+int+ixmult2;
MAT3=[zeros(ixmultOg,int2);MAT2;zeros(ixmultOg,int2)];
lsOg2=linspace(lsOg(1)-(ixmedOg*ixmultOg),lsOg(end)+(ixmedOg*ixmultOg),int+2*ixmultOg);
if DAT{1}(18)==0 %STANDARD METHOD
gstd=4;
myfilter = fspecial('gaussian',[10 10], gstd);
elseif DAT{1}(18)==1 %EXPERIMENTAL METHOD (MAINTAIN INSTANTANEOUSNESS)
gstd=3.5;
myfilter = fspecial('gaussian',[20 20], gstd);
end
MAT_filtered = imfilter(MAT3, myfilter, 'replicate');
MAT_filtered_neg=-MAT_filtered;
% gstd=6;
% myfilter = fspecial('gaussian',[60 60], gstd);
% MAT_filtered = imfilter(MAT3, myfilter, 'replicate');
% MAT_filtered_neg2=-MAT_filtered;
% subplot(5,10,ctr)
% imagesc(MAT_filtered_neg)
% title(num2str(sig),'fontsize',7,'fontweight','bold')
% set(gca,'xtick',[],'ytick',[])
% end
lsOg=lsOg2;
medO=medO2;
lsf=lsf2;
Ogive30{ii}=Ogive302{ii};
% determine frequencywise cumulative sum of MAT_filtered_neg
% for each frequency interval to be optimized:
sumFR=nan(sizFR1,1);
mnixFR=nan(sizFR1,1);
MAT_nansum_cumsum=fliplr(cumsum(fliplr(nansum(MAT_filtered_neg))));
for ih=1:sizFR1
[~,mnixFR(ih)]=min(abs(lsf-FR(ih,1)));
sumFR(ih)=MAT_nansum_cumsum(mnixFR(ih));
end
FR=[FR sumFR];
if DAT{12}(10,2)>0
subplot(3,3,[4 8])
imagesc(lsf,lsOg,MAT_filtered_neg)
if DAT{5}(DAT{14}(1,ii))==3 %QSENS
ylabel('Q_{SENS} = c_p\rho_dOg_{w\theta_v} [Wm^{-2}]','fontsize',14,'fontweight','bold')
elseif DAT{5}(DAT{14}(1,ii))==5 %ARBITRARY SCALAR (E.G.: FCO2 OR FCH4)
ylabel([DAT{8}{ixes(3,ii),6},' = c_{ds}Og_{wr_',DAT{8}{ixes(3,ii),4},'} [',[Ylabels2{DAT{6}(5,ixes(2,ii))+1} Ylabels{DAT{6}(4,ixes(2,ii))-5}],']',],'fontsize',14,'fontweight','bold')
elseif DAT{5}(DAT{14}(1,ii))==6 %QLAT
ylabel('Q_{LAT} = L_hc_{ds}Og_{wr_q} [Wm^{-2}]','fontsize',14,'fontweight','bold')
end
xlabel('Natural frequency','fontsize',14,'fontweight','bold')
hold on
colormap(cm_densitymap_gray)
plot(lsf,medO,'-r')
for ih=1:length(prcO)
plot(lsf0,prcO{ih}(1,:),'-c')
plot(lsf0,prcO{ih}(2,:),'-c')
end
axis tight
set(gca,'box','on','xlim',ax(1:2))
ax=axis;
xtx=get(gca,'xtick');%xtx2=unique(round(xtx));xtx2=xtx2(xtx2>=xtx(1) & xtx2<=xtx(end));
ytx=get(gca,'ytick');
for ih=1:length(xtx)
plot([xtx(ih) xtx(ih)],[ax(3) ax(4)],':','color',[.7 .7 .7])
end
for ih=1:length(ytx)
plot([ax(1) ax(2)],[ytx(ih) ytx(ih)],':','color',[.7 .7 .7])
end
for ih=1:sizFR1
plot([FR(ih,1) FR(ih,1)],[ax(3) ax(4)],'--k','linewidth',2)
end
plot([ax(1) ax(2)],[0 0],'-k')
xtxlab=cell(length(xtx),1);
for ih=1:length(xtx)
xtxlab{ih}=['10^{',num2str(xtx(ih)),'}'];
end
set(gca,'ydir','normal','fontsize',12,'fontweight','bold','box','on','xticklabel',xtxlab)
set(gca,'ylim',ax(3:4))
path=[datestr(DAT15(3),'yyyy-mm-dd HH-MM-SS'),' -- ',num2str(DAT15(1)),'-',num2str(DAT15(2)),'.fig'];
if ~(exist([figurepath,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep],'dir')==7)
mkdir([figurepath,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep])
end
hgsave([figurepath,num2str(ii),'_',DAT{8}{ixes(3,ii),3},filesep,path])
close(HF)
end
end
filename_com=[datestr(DAT15(3),'yyyy-mm-dd HH-MM-SS'),' -- ',num2str(DAT15(1)),'-',num2str(DAT15(2)),'.mat'];
if ~(exist([savepath0,num2str(ii),'_',DAT{8}{DAT{14}(3,ii),7}],'dir')==7)
mkdir([savepath0,num2str(ii),'_',DAT{8}{DAT{14}(3,ii),7}])
end
densitymap_filename{1,ii}=[savepath0,num2str(ii),'_',DAT{8}{DAT{14}(3,ii),7},filesep,filename_com];
save(densitymap_filename{1,ii},'FR','lsf','fr_nat','lsOg','MAT_filtered_neg','Ogive_end_prc','datalength','int2','ixmultLOG','Ogive30','yax','datalength','Flux_lindet','BEA_STATS_LOG');
end
catch me
txtlog=errordisp(SCRIPT_LEVEL,fname,txtlog,me,id);
runok=0;
WS=ws2struct();save([errordir,filesep,'workspace_',datestr(now,'yyyymmdd_HHMMSS_FFF'),'_Densitymapsub_error.mat'],'WS')
end
end