From 6038dc392e1adcbb5998ec88e6454fa870a58776 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 3 Aug 2023 18:24:42 +0200
Subject: [PATCH] use of the new routine growth rate
---
matlab/compute/mode_growth_meter.m | 62 ++++++++++++++++---
wk/lin_run_script.m | 23 ++++----
wk/lin_scan_script.m | 58 ++++++++----------
wk/load_metadata_scan.m | 95 +++++++++++++++++++-----------
4 files changed, 152 insertions(+), 86 deletions(-)
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 921dd743..9993b265 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -30,6 +30,8 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
+[wkykx, ekykx] = compute_growth_rates(DATA.PHI(:,:,:,FRAMES),DATA.Ts3D(FRAMES));
+
FIGURE.fig = figure; %set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
@@ -41,14 +43,20 @@ for i = 1:2
if MODES_SELECTOR == 2
switch OPTIONS.ik
case 'sum'
- plt_ = @(x,ik) movmean(squeeze(sum(abs((x(:,ik,FRAMES))),1)),Nma);
+ plt_ = @(f,ik) movmean(squeeze(sum(abs((f(:,ik,FRAMES))),1)),Nma);
MODESTR = '$\sum_{k_y}$';
+ omega= @(ik) wkykx(1,:,end);
+ err = @(ik) ekykx(1,:);
case 'max'
- plt_ = @(x,ik) movmean(squeeze(max(abs((x(:,ik,FRAMES))),[],1)),Nma);
+ plt_ = @(f,ik) movmean(squeeze(max(abs((f(:,ik,FRAMES))),[],1)),Nma);
MODESTR = '$\max_{k_y}$';
+ omega= @(ik) wkykx(1,:,end);
+ err = @(ik) ekykx(1,:);
otherwise
- plt_ = @(x,ik) movmean(abs(squeeze(x(OPTIONS.ik,ik,FRAMES))),Nma);
+ plt_ = @(f,ik) movmean(abs(squeeze(f(OPTIONS.ik,ik,FRAMES))),Nma);
MODESTR = ['$k_y=$',num2str(DATA.grids.ky(OPTIONS.ik))];
+ omega= @(ik) wkykx(OPTIONS.ik,:,end);
+ err = @(ik) ekykx(OPTIONS.ik,:);
end
kstr = 'k_x';
% Number max of modes to plot is kx>0 (1/2) of the non filtered modes (2/3)
@@ -59,13 +67,19 @@ for i = 1:2
case 'sum'
plt_ = @(x,ik) movmean(squeeze(sum(abs((x(ik,:,FRAMES))),2)),Nma);
MODESTR = '$\sum_{k_x}$';
+ omega= @(ik) wkykx(:,1,end);
+ err = @(ik) ekykx(:,1);
case 'max'
plt_ = @(x,ik) movmean(squeeze(max(abs((x(ik,:,FRAMES))),[],2)),Nma);
MODESTR = '$\max_{k_x}$';
+ omega= @(ik) wkykx(:,1,end);
+ err = @(ik) ekykx(:,1);
otherwise
plt_ = @(x,ik) movmean(abs(squeeze(x(ik,OPTIONS.ik,FRAMES))),Nma);
MODESTR = ['$k_x=$',num2str(DATA.grids.kx(OPTIONS.ik))];
- end
+ omega= @(ik) wkykx(:,OPTIONS.ik,end);
+ err = @(ik) ekykx(:,OPTIONS.ik);
+ end
kstr = 'k_y';
% Number max of modes to plot is ky>0 (1/1) of the non filtered modes (2/3)
Nmax = ceil(DATA.grids.Nky*2/3);
@@ -81,6 +95,8 @@ for i = 1:2
gamma = MODES;
amp = MODES;
+ % w_ky = zeros(numel(MODES),numel(FRAMES)-1);
+ % ce = zeros(numel(MODES),numel(FRAMES));
i_=1;
for ik = MODES
@@ -88,9 +104,29 @@ for i = 1:2
gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
gamma(i_) = gr(1);
amp(i_) = gr(2);
+
+ % % Second method for measuring growth rate (measures frequ. too)
+ % if MODES_SELECTOR == 1
+ % to_measure = reshape(DATA.PHI(ik,1,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
+ % else
+ % to_measure = reshape(DATA.PHI(1,ik,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
+ % end
+ % for it = 2:numel(FRAMES)
+ % phi_n = to_measure(:,it);
+ % phi_nm1 = to_measure(:,it-1);
+ % dt = t(it)-t(it-1);
+ % ZS = sum(phi_nm1,1); %sum over z
+ %
+ % wl = log(phi_n./phi_nm1)/dt;
+ % w_ky(i_,it) = squeeze(sum(wl.*phi_nm1,1)./ZS);
+ %
+ % % for iky = 1:numel(w_ky(:,it))
+ % % ce(iky,it) = abs(sum(abs(w_ky(iky,it)-wl(iky,:)).^2.*phi_nm1(iky,:),2)./ZS(iky,:));
+ % % end
+ % end
i_=i_+1;
end
- mod2plot = [2:min(Nmax,OPTIONS.NMODES+1)];
+ mod2plot = 2:min(Nmax,OPTIONS.NMODES+1);
clr_ = jet(numel(mod2plot));
%plot
% subplot(2+d,3,1+3*(i-1))
@@ -120,15 +156,25 @@ for i = 1:2
% subplot(2+d,3,3+3*(i-1))
FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
- plot(k(MODES),gamma,...
- 'DisplayName',['(',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); hold on;
+ % yyaxis("left")
+ errorbar(k(MODES),real(omega(ik)),real(err(ik)),'-k',...
+ 'DisplayName',...
+ ['$\gamma$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
+ hold on;
for i_ = 1:numel(mod2plot)
plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
end
if MODES_SELECTOR == 2
plot(k(ikzf),gamma(ikzf),'ok');
end
- xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
+ % ylabel('$\gamma$');
+ % yyaxis("right")
+ errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--k',...
+ 'DisplayName',...
+ ['$\omega$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
+ hold on;
+ ylabel('$\gamma,\omega$');
+ xlabel(['$',kstr,'\rho_s$']);
title('Growth rates')
end
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 46e94362..acb62f80 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -6,8 +6,8 @@
%% Set up the paths for the necessary Matlab modules
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
-% mpirun = 'mpirun';
-mpirun = '/opt/homebrew/bin/mpirun'; % for macos
+mpirun = 'mpirun';
+% mpirun = '/opt/homebrew/bin/mpirun'; % for macos
addpath(genpath([gyacomodir,'matlab'])) % Add matlab folder
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot folder
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
@@ -24,20 +24,23 @@ EXECNAME = 'gyacomo23_sp'; % single precision
% run lin_DTT_HM_rho85
% run lin_DTT_HM_rho98
% run lin_DTT_LM_rho90
-run lin_DTT_LM_rho95
-% run lin_JET_rho97
+% run lin_DTT_LM_rho95
+run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
+% run lin_KBM
%% Change parameters
NY = 2;
-PMAX = 4;
+NX = 4;
+PMAX = 2;
JMAX = PMAX/2;
-ky = 1.0;
+ky = 0.5;
LY = 2*pi/ky;
-DT = 1e-5;
-% SIGMA_E = 0.04;
-TMAX = 10;
-DTSAVE0D = 200*DT;
+DT = 1e-3;
+TMAX = 10;
+% % SIGMA_E = 0.04;
+% TMAX = 10;
+% DTSAVE0D = 200*DT;
DTSAVE3D = TMAX/50;
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index a587fc23..fa4b3509 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -16,37 +16,41 @@ addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder
%% Setup run or load an executable
RUN = 1; % To run or just to load
-RERUN = 1; % rerun if the data does not exist
+RERUN = 0; % rerun if the data does not exist
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
-%% Setup basic parameters
-run lin_DTT_AB_rho85_PT
+%% Setup parameters
+% run lin_DTT_AB_rho85
+% run lin_DTT_AB_rho98
+run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
-%% Modify parameters
-% NZ = 1;
-NY = 2;
-DT = 0.5e-3;
-TMAX = 50;
-MU_X = 0.1; MU_Y = 0.1;
+%% Change parameters
+NY = 2;
+% SIGMA_E = 0.023;
%% Scan parameters
SIMID = [SIMID,'_scan'];
P_a = [2 4 6 8];
-ky_a = 0.05:0.05:1.5;
+ky_a = logspace(-1.5,1.5,30);
%% Scan loop
% arrays for the result
-g_ky = zeros(numel(ky_a),numel(P_a),2);
-g_avg= g_ky*0;
+g_ky = zeros(numel(ky_a),numel(P_a));
g_std= g_ky*0;
+w_ky = g_ky*0;
+w_std= g_ky*0;
j = 1;
for PMAX = P_a
JMAX = P/2;
i = 1;
for ky = ky_a
- LY = 2*pi/ky;
+ LY = 2*pi/ky;
+ DT = 1e-4;%min(1e-2,1e-3/ky);
+ TMAX = 10;%min(10,1.5/ky);
+ DTSAVE0D = 0.1;
+ DTSAVE3D = 0.1;
%% RUN
setup
% naming
@@ -67,8 +71,7 @@ for PMAX = P_a
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
- if numel(data_.Ts3D)>10
- if numel(data_.Ts3D)>5
+ if numel(data_.Ts0D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
@@ -83,23 +86,14 @@ for PMAX = P_a
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
- field = 0;
- field_t = 0;
- for ik = 2:NY/2+1
- field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
- field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
- to_measure = log(field_t(it1:it2));
- tw = double(data_.Ts3D(it1:it2));
- % gr = polyfit(tw,to_measure,1);
- gr = fit(tw,to_measure,'poly1');
- err= confint(gr);
- g_ky(i,j,ik) = gr.p1;
- g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
- end
+ [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
+ g_ky (i,j) = real(wkykx(2,1,end));
+ g_std(i,j) = real(ekykx(2,1));
+ w_ky (i,j) = imag(wkykx(2,1,end));
+ w_std(i,j) = imag(ekykx(2,1));
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
- end
end
i = i + 1;
end
@@ -107,11 +101,11 @@ for PMAX = P_a
end
%% take max growth rate among z coordinate
-y_ = g_ky(:,:,2);
-e_ = g_std(:,:,2);
+y_ = g_ky + 1i*w_ky;
+e_ = g_std+ 1i*w_std;
%% Save scan results (gamma)
-if(numel(ky_a)>1 && numel(P_a)>1)
+if(numel(ky_a)>1 || numel(P_a)>1)
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
filename = [num2str(NX),'x',num2str(NZ),...
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 600797d1..6615ee53 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -13,26 +13,32 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.5_P_2_8_DGGK_0.05_be_0.0034.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.5_P_2_8_DGGK_0.05_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_1_P_2_8_DGGK_0.05_be_0.0039.mat';
-datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
+% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/8x24_ky_0.05_1.5_P_2_8_kN_1.33_DGGK_0.1_be_0.0039.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_2_kN_10_DGGK_0.05_be_0.0031.mat';
+datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_2_2_kN_10_DGGK_0.1_be_0.0031.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 0;
%% Load data
fname = ['../results/',datafname];
d = load(fname);
+gamma = real(d.data); g_err = real(d.err);
+omega = imag(d.data); w_err = imag(d.err);
if FILTERGAMMA
- d.data = d.data.*(d.data>0.025);
- d.err = d.err.*(d.data>0.025);
+ d.data = real(gamma).*(real(gamma)>0.025) + imag(gamma);
end
if 0
%% Pcolor of the peak
figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
-pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
-% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
-% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
+pclr=imagesc_custom(XX_,YY_,gamma'.*(gamma>0)');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),gamma'.*(gamma>0)');
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),gamma'.*(gamma>0)');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
@@ -48,23 +54,40 @@ if 1
%% Scan along first dimension
figure
colors_ = jet(numel(d.s2));
+subplot(121)
for i = 1:numel(d.s2)
- % plot(d.s1,d.data(:,i),'s-',...
- plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
- 'color',colors_(i,:));
- % errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
- % 'LineWidth',2.0,...
- % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
- % 'color',colors_(i,:));
+ % plot(d.s1,gamma(:,i),'s-',...
+ % plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
+ % 'LineWidth',2.0,...
+ % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ % 'color',colors_(i,:));
+ errorbar(d.s1,gamma(:,i),g_err(:,i),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ 'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
+
+subplot(122)
+for i = 1:numel(d.s2)
+ % plot(d.s1,gamma(:,i),'s-',...
+ % plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
+ % 'LineWidth',2.0,...
+ % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ % 'color',colors_(i,:));
+ errorbar(d.s1,omega(:,i),w_err(:,i),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ 'color',colors_(i,:));
+ hold on;
+end
+xlabel(d.s1name); ylabel('$\omega c_s/R$');title(d.title);
+xlim([d.s1(1) d.s1(end)]);
+
colormap(colors_);
clb = colorbar;
-% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
clim([1 numel(d.s2)+1]);
clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
clb.Ticks =1.5:numel(d.s2)+1.5;
@@ -78,11 +101,11 @@ if 0
figure
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
- plot(d.s2,d.data(i,:),'s-',...
+ plot(d.s2,gamma(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
-% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
+% errorbar(d.s2,gamma(i,:),g_err(i,:),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
% 'color',colors_(i,:));
@@ -104,13 +127,13 @@ end
if 0
%% Convergence analysis
figure
-% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
+% target_ = 0.25*(gamma(end,end)+gamma(end-i_,end)+gamma(end,end-i_)+gamma(end-i_,end-i_));
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
-% target_ = d.data(i,end);
+% target_ = gamma(i,end);
target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
if target_ > 0
- eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
+ eps_ = abs(target_ - gamma(i,1:end-1))/abs(target_);
semilogy(d.s2(1:end-1),eps_,'-s',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
@@ -140,30 +163,30 @@ figure; i_ = 0;
% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
-% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
-% eps_ = log(abs(target_ - d.data)/abs(target_));
-eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
-sign_ = 1;%sign(d.data - target_);
-eps_ = d.data;
+% target_ = 0.25*(gamma(end,end)+gamma(end-i_,end)+gamma(end,end-i_)+gamma(end-i_,end-i_));
+% eps_ = log(abs(target_ - gamma)/abs(target_));
+eps_ = max(-10,log(abs(target_ - gamma)/abs(target_)));
+sign_ = 1;%sign(gamma - target_);
+eps_ = gamma;
for i = 1:numel(d.s1)
for j = 1:numel(d.s2)
- % target_ = d.data(i,end);
- % target_ = d.data(i,end);
- % target_ = d.data(end,j);
- eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
+ % target_ = gamma(i,end);
+ % target_ = gamma(i,end);
+ % target_ = gamma(end,j);
+ eps_(i,j) = log(abs(target_ - gamma(i,j))/target);
if target_ > 0
- % eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
- % eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
- % eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
+ % eps_(i,:) = max(-12,log(abs(target_ - gamma(i,1:end))/abs(target_)));
+ % eps_(i,:) = log(abs(target_ - gamma(i,1:end))/abs(target_));
+ % eps_(i,:) = min(100,100*abs(target_ - gamma(i,1:end))/abs(target_));
else
end
end
end
[XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
-pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
-% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
-% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
+pclr=imagesc_custom(XX_,YY_,eps_'.*(gamma>0)'.*sign_');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(gamma>0)'.*sign_',5);
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(gamma>0)'.*sign_');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
--
GitLab