save script

wk/p3_geom_scan_analysis.m

 % Get the current directory (wk)
 curdir  = pwd;
+% parition= '/misc/gyacomo23_outputs/paper_3/';
+partition= '../results/paper_3/';
 % Get the scan directory
-if 0 % delta kappa scan
-    scandir = '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/P2_J1_delta_kappa_scan/';
-    % scandir = '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/P4_J2/';
-    % scandir = '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/P4_J2_delta_kappa_scan/';
-    pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
-    pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53;
-else % shear safety factor scan
-    scandir = '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/P2_J1_PT_sfact_shear_scan/';
-    % scandir = '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/P2_J1_NT_sfact_shear_scan/';
-    pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63;
-    pnam2 = '$q_0$';    attr2 = 'q0';    pref2 =-2.15;
-end
+switch 3
+    case 1 % delta kappa scan
+    % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_kappa_scan/';
+    scandir = 'DTT_rho85_geom_scan/P4_J2_delta_kappa_scan/';
+    nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
+    nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53;
+    case 2 % shear safety factor scan
+    % scandir = 'DTT_rho85_geom_scan/P2_J1_PT_sfact_shear_scan/';
+    scandir = 'DTT_rho85_geom_scan/P2_J1_NT_sfact_shear_scan/';
+    nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63;
+    nml2 = 'GEOMETRY'; pnam2 = '$q_0$';    attr2 = 'q0';    pref2 =-2.15;
+    case 3
+    % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuDGGK_scan/';
+    % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuDGGK_scan/';
+    % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuSGGK_scan/';
+    scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuSGGK_scan/';
+    % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuLDGK_scan/';
+    % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuLDGK_scan/';
+    nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
+    nml2 = 'MODEL';    pnam2 = '$\nu$';    attr2 = 'nu';    pref2 = 0.5;
+end 
+scandir = [partition,scandir]; 
 % Get a list of all items in the current directory
 contents = dir(scandir);
 
 % give ref value and param names
-REFVAL= 0;
-
+REFVAL= 1;
+% Get and plot the fluxsurface
+GETFLUXSURFACE = 0;
 
 % Iterate through the contents
-Qxavg = []; Qxerr = []; para1 = []; para2 = [];
+Qxavg = []; Qxerr = []; para1 = []; para2 = []; R = []; Z = [];
 for i = 1:length(contents)
     % Check if the item is a directory and not '.' or '..'
-    if contents(i).isdir && ~strcmp(contents(i).name, '.') && ~strcmp(contents(i).name, '..')
+    if contents(i).isdir && ~strcmp(contents(i).name, '.') ...
+            && ~strcmp(contents(i).name, '..')
         % Get and display the name of the subdirectory
-        sdubdir = [scandir,contents(i).name];
+        subdir = [scandir,contents(i).name];
         disp(['Subdirectory: ' contents(i).name]);
         % Get parameters
-        param = read_namelist([sdubdir,'/fort_00.90']);
-        para1 = [para1 param.GEOMETRY.(attr1)];
-        para2 = [para2 param.GEOMETRY.(attr2)];        
+        param = read_namelist([subdir,'/fort_00.90']);
+        para1 = [para1 param.(nml1).(attr1)];
+        para2 = [para2 param.(nml2).(attr2)];        
         % Now you are in the subdirectory. You can perform operations here.
-        [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(sdubdir);
-        if(t_all(end) > 40)
+        [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir);
+        if(t_all(end) > 200)
             [fullAvg,sliceAvg,sliceErr] = sliceAverage(Qxi_all+Qxe_all,3);
             Qxavg = [Qxavg fullAvg];
             Qxerr = [Qxerr mean(sliceErr)];
@@ -43,6 +57,14 @@ for i = 1:length(contents)
             Qxerr = [Qxerr nan];
         end
     end
+    if GETFLUXSURFACE
+        data = load([subdir,'/RZ.txt']);
+        R_ = data(:, 1);
+        Z_ = data(:, 2);
+        R_ = [R_;R_(1)]'; Z_  = [Z_;Z_(1)]';
+        R  = [R ; R_];    Z   = [Z ; Z_];
+    end
+
 end
 %% reshaping, sorting and plotting
 p1 = unique(para1);