diff --git a/wk/NT_study_DTT_rho85.m b/wk/NT_study_DTT_rho85.m
new file mode 100644
index 0000000000000000000000000000000000000000..20bbe9c790b17c472599deac2613211d17465a0a
--- /dev/null
+++ b/wk/NT_study_DTT_rho85.m
@@ -0,0 +1,55 @@
+
+[t_PT, Gi_PT, Qi_PT, Ge_PT, Qe_PT] = read_flux_out_XX('/Users/ahoffmann/gyacomo/results/paper_3/DTT_rho85/3x2x192x48x32');
+half = ceil(numel(t_PT)/2);
+Gi_avg_PT = mean(Gi_PT(half:end));
+Ge_avg_PT = mean(Ge_PT(half:end));
+Qi_avg_PT = mean(Qi_PT(half:end));
+Qe_avg_PT = mean(Qe_PT(half:end));
+
+
+[t_NT, Gi_NT, Qi_NT, Ge_NT, Qe_NT] = read_flux_out_XX('/Users/ahoffmann/gyacomo/results/paper_3/DTT_rho85/3x2x192x48x32_NT');
+half = ceil(numel(t_NT)/2);
+Gi_avg_NT = mean(Gi_NT(half:end));
+Ge_avg_NT = mean(Ge_NT(half:end));
+Qi_avg_NT = mean(Qi_NT(half:end));
+Qe_avg_NT = mean(Qe_NT(half:end));
+
+prct_change = @(x,y) abs(y-x)/abs(x)*100;
+nmvm = 10;
+plt_ = @(t,f,stl,nme) plot(movmean(t,nmvm),movmean(f,nmvm), stl,'DisplayName',nme);
+
+figure
+
+nr = 1;
+if 1
+    nr = 2;
+    subplot(223)
+    plt_(t_PT,Gi_PT, '-r','ions PT'); hold on;
+    plt_(t_NT,Gi_NT,'--r','ions NT'); hold on;
+    xlabel('$tc_s/R$'); ylabel('$\Gamma_{x}$');
+    legend('show')
+    title(['Ion particle flux, NT/NP=',sprintf('%3.1f',prct_change(Gi_avg_PT,Gi_avg_NT)),'$\%$'])
+    subplot(224)
+    plt_(t_PT,Ge_PT, '-b','electrons PT'); hold on;
+    plt_(t_NT,Ge_NT,'--b','electrons NT'); hold on;
+    xlabel('$tc_s/R$'); ylabel('$\Gamma_{x}$');
+    legend('show')
+    title(['Elec. particle flux, NT/NP=',sprintf('%3.1f',prct_change(Ge_avg_PT,Ge_avg_NT)),'$\%$'])
+end
+
+subplot(nr,2,1)
+plt_(t_PT,Qi_PT, '-r','ions PT'); hold on;
+plt_(t_NT,Qi_NT,'--r','ions NT'); hold on;
+xlabel('$tc_s/R$'); ylabel('$Q_{x}$');
+legend('show')
+title(['Ion heat flux, NT/NP=',sprintf('%3.1f',prct_change(Qi_avg_PT,Qi_avg_NT)),'$\%$'])
+
+subplot(nr,2,2)
+plt_(t_PT,Qe_PT, '-b','electrons PT'); hold on;
+plt_(t_NT,Qe_NT,'--b','electrons NT'); hold on;
+xlabel('$tc_s/R$'); ylabel('$Q_{x}$');
+legend('show')
+title(['Elec. heat flux, NT/NP=',sprintf('%3.1f',prct_change(Qe_avg_PT,Qe_avg_NT)),'$\%$'])
+
+disp(['PT total heat flux : ',num2str(Qi_avg_PT + Qe_avg_PT)]);
+disp(['NT total heat flux : ',num2str(Qi_avg_NT + Qe_avg_NT)]);