Update

plotter/detector.py

@@ -5,11 +5,15 @@
 import matplotlib.patches as patches
 import matplotlib.cm as colormap
 import scipy.stats as st
+import qrcode
 
 
-def detector(shotnumber, time, runnumber, is_debug=False):
+def detector(shotnumber, time, runnumber, title=None, is_debug=False):
     results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber)
 
+    if title is None:
+        r'COMPASS #' + shotnumber + '\n ($t = $' + time + ' ms)'
+
     try:
         cell_counts = np.loadtxt(os.path.join(results_folder, 'detector', 'cellcounter.dat'))
         det_x = np.loadtxt(os.path.join(results_folder, 'detector', 'detx'))
@@ -25,6 +29,8 @@ def detector(shotnumber, time, runnumber, is_debug=False):
     cmax = cell_counts.max() if cell_counts.max() > 0 else 1
 
     try:
+        plt.rc('font', family='serif')
+        plt.rc('text', usetex=True)
         fig, ax = plt.subplots()
         ax.set_xlim(xmin, xmax)
         ax.set_ylim(ymin, ymax)
@@ -47,13 +53,21 @@ def detector(shotnumber, time, runnumber, is_debug=False):
         if is_debug:
             print(cell_counts)
 
-        plt.xlabel(r"$x \mathrm{ [mm]}$")
-        plt.ylabel(r"$y \mathrm{ [mm]}$")
-        plt.title(f"COMPASS #{shotnumber} (t={time} s)")
+        plt.xlabel(r"toroidal detector position [mm]")
+        plt.ylabel(r"vertical detector position [mm]")
+        plt.title(title)
+
+        qr_content = f'raw.githubusercontent.com/taiga-project/refs/run/{runnumber}'
+        print(f'Add QR code: {qr_content}')
+        qr_image = qrcode.make(qr_content)
+        ax_qr = fig.add_axes([0.175, 0.495, 0.08, 0.5], anchor='NE', zorder=10)
+        ax_qr.imshow(qr_image, cmap='gray')
+        ax_qr.axis('off')
 
         try:
-            print(f'Save plot to {os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf")}')
-            plt.savefig(os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf"))
+            filename = f"detector_{runnumber}.pdf"
+            print(f'Save plot to {os.path.join(results_folder, filename)}')
+            plt.savefig(os.path.join(results_folder, filename), dpi=300, bbox_inches='tight')
             plt.clf()
         except Exception as e:
             print(f'Unable to save to: {results_folder} ({e})')

plotter/detector_plane.py

@@ -3,9 +3,10 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import scipy.stats as st
+import qrcode
 
 
-def detector_plane(shotnumber, time, runnumber, detector_par):
+def detector_plane(shotnumber, time, runnumber, detector_par, title=None):
     max_distance = 1e-4
     detector_par_list = detector_par.split(',')
 
@@ -36,6 +37,9 @@ def detector_plane(shotnumber, time, runnumber, detector_par):
 
     results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber)
 
+    if title is None:
+        r'COMPASS #' + shotnumber + '\n ($t = $' + time + ' ms)'
+
     try:
         R = np.loadtxt(os.path.join(results_folder, 'rad.dat'))
         Z = np.loadtxt(os.path.join(results_folder, 'z.dat'))
@@ -50,6 +54,9 @@ def detector_plane(shotnumber, time, runnumber, detector_par):
         y = Z[particle_on_detector]
         xmin, xmax = -4e-2, 4e-2
         ymin, ymax = 18e-2, 26e-2
+        #ymin, ymax = 24e-2, 32e-2
+        #xmin, xmax = -8e-2, 0e-2
+        #ymin, ymax = 14e-2, 22e-2
 
         # Perform the kernel density estimate
         xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
@@ -59,6 +66,8 @@ def detector_plane(shotnumber, time, runnumber, detector_par):
         f = np.reshape(kernel(positions).T, xx.shape)
 
         fig, ax = plt.subplots()
+        plt.rc('font', family='serif')
+        plt.rc('text', usetex=True)
         ax.set_xlim(xmin, xmax)
         ax.set_ylim(ymin, ymax)
         ax.set_aspect('equal')
@@ -67,12 +76,20 @@ def detector_plane(shotnumber, time, runnumber, detector_par):
         ax.hist2d(x, y, bins=250, range=[[xmin, xmax], [ymin, ymax]], cmap='afmhot')
 
         # Set labels
-        plt.xlabel(r"$T \mathrm{ [m]}$")
-        plt.ylabel(r"$Z \mathrm{ [m]}$")
+        plt.xlabel(r"$T$ [m]")
+        plt.ylabel(r"$Z$ [m]")
+        plt.title(title)
+
+        qr_content = f'raw.githubusercontent.com/taiga-project/refs/run/{runnumber}'
+        print(f'Add QR code: {qr_content}')
+        qr_image = qrcode.make(qr_content)
+        ax_qr = fig.add_axes([0.24, 0.495, 0.08, 0.5], anchor='NE', zorder=10)
+        ax_qr.imshow(qr_image, cmap='gray')
+        ax_qr.axis('off')
 
         try:
             print(f'Save plot to {os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf")}')
-            plt.savefig(os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf"))
+            plt.savefig(os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf"), dpi=300, bbox_inches='tight')
             plt.clf()
         except Exception as e:
             print(f'Unable to save to: {results_folder} ({e})')

plotter/traj_plotter.py

@@ -5,10 +5,10 @@
 from scipy.interpolate import interp1d
 import h5py
 import matplotlib.pyplot as plt
+import qrcode
 
 
-
-def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, title=None):
+def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, title=None, R_min=0.65):
     root_folder = os.path.abspath(os.getcwd())
     result_folder = os.path.join(root_folder, 'results')
     field_folder = os.path.join(root_folder, 'input', 'fieldGrid')
@@ -25,7 +25,7 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit
     ionR = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'rad.dat'))
     ionY = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'ionyeald.dat'))
     ionX1 = np.gradient(ionY)
-    ionF = interp1d(ionR, -ionX1 / np.amax(np.absolute(ionX1)))
+    ionF = interp1d(ionR, -ionX1 / np.amax(np.absolute(ionX1[ionR > R_min])))
 
     plt.plot(ionR, ionF(ionR))
 
@@ -67,7 +67,7 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit
     R, Z = np.meshgrid(rcord, zcord)
 
     # Filter data
-    t_index = cellid >= 0
+    t_index = (cellid >= 0) & (t_rad[1, :] > R_min)
     t_rad = t_rad[:, t_index]
     t_z = t_z[:, t_index]
     t_tor = t_tor[:, t_index]
@@ -113,11 +113,14 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit
 
     # Create a 2D figure
     plt.clf()
+    plt.rc('font', family='serif')
+    plt.rc('text', usetex=True)
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
     # Plot filled contours
     plt.contourf(R, Z, PSI, levels=np.arange(-2, 2, 0.1), cmap='GnBu')
+    #plt.contourf(R, Z, PSI, levels=np.arange(0, 12, 0.3), cmap='GnBu')
 
     # Plot boundary line
     plt.plot(boundary_r, boundary_z, color='black')
@@ -135,7 +138,6 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit
     for i in intensity_indices:
         # Plot atom trajectory
         plt.plot([rcord[-1], t_rad[1, i]], [t_z[0, i], t_z[0, i]], color=(1, 1, 0), alpha=ion_factor[i] * alpha_factor)
-    print(number_of_ions)
     try:
         # Detector parameters
         detector_y_size = 12.2e-3
@@ -146,20 +148,27 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit
         detector_y_tilt = detector_y_size * np.cos(detector_phi / 180 * np.pi)
 
         # Plot detector lines
+        detector_colour = (0, 0.5, 0)
         plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)],
-                 color='blue', linewidth=1)
+                 color=detector_colour, linewidth=1)
         plt.plot([detector_x0 - detector_x_tilt, detector_x0 + detector_x_tilt],
                  [detector_y0 + detector_y_tilt, detector_y0 - detector_y_tilt],
-                 color='blue', linewidth=2)
+                 color=detector_colour, linewidth=2)
     except ValueError:
         print('Invalid detector coordinates.')
 
     plt.xlabel(r'$R$ [m]')
     plt.ylabel(r'$Z$ [m]')
     plt.title(title)
     ax.set_aspect('equal', 'box')
+    qr_content = f'https://taiga-project.github.io?{runnumber}'
+    print(f'Add QR code: {qr_content}')
+    qr_image = qrcode.make(qr_content)
+    ax_qr = fig.add_axes([0.24, 0.495, 0.08, 0.5], anchor='NE', zorder=10)
+    ax_qr.imshow(qr_image, cmap='gray')
+    ax_qr.axis('off')
     try:
-        plot_filename = os.path.join(working_folder, f'traj_{shot_and_time}')
+        plot_filename = os.path.join(working_folder, f'traj_{runnumber}')
         plt.savefig(plot_filename+'.pdf', dpi=300, bbox_inches='tight')
         print(f'Save plot to {plot_filename}.pdf')
         plt.savefig(plot_filename+'.png', dpi=300, bbox_inches='tight')

preproc/renate_od/shake_flux_surface.py

@@ -1,5 +1,8 @@
 import matplotlib.pyplot
 import pandas
+import sys
+import os
+sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..'))
 
 from preproc.renate_od.interface import get_lcfs_radial_coordinate
 from preproc.renate_od.manager import RenateODManager
@@ -47,7 +50,6 @@ def __init__(self, beamlet_geometry, shot_number, time, species, energy, shift=0
 
 
 def shake_flux_surface(species, shot_number, time, energy):
-
     export_directory='.'
     z = 0
     tor = 0
@@ -58,26 +60,31 @@ def shake_flux_surface(species, shot_number, time, energy):
     matplotlib.pyplot.minorticks_on()
     matplotlib.pyplot.grid(which='major')
     matplotlib.pyplot.xlabel('$R$ [m]', labelpad=-10.5, loc='right')
-    matplotlib.pyplot.ylabel('neutral beam attenuation')
-    matplotlib.pyplot.title('COMPASS #' + shot_number + ' (' + time + ' ms)')
-    R_LCFS = get_lcfs_radial_coordinate(shot_number, time)
+    matplotlib.pyplot.ylabel('normalized\nneutral beam density')
+    matplotlib.pyplot.title('reference discharge, ' + species + ' beam, ' + str(energy) + ' keV')
+    R_LCFS = get_lcfs_radial_coordinate(shot_number=shot_number, time=time)
     matplotlib.pyplot.axvline(R_LCFS, c='red', ls='--')
+    matplotlib.pyplot.text(R_LCFS+0.005, 0.56, 'reference', c='red', fontsize=6)
     matplotlib.pyplot.text(R_LCFS+0.005, 0.45, 'LCFS', c='red', fontsize=12)
 
     export_directory = get_home_directory() + '/input/ionProf/' + shot_number + '_' + time
     for shift in numpy.linspace(-0.02, 0.02, 5):
         r = RenateODManagerMock(beamlet_geometry, shot_number, time, species, energy, shift)
         radial_coordinate, relative_attenuation = r.get_attenuation_profile()
-        ax.plot(radial_coordinate, relative_attenuation, '-', linewidth=2, label=r'$\rho$->$\rho$+%.2f'%shift)
+        label = r'with %3.f%% $\rho$' % (100*(1+shift))
+        if shift == 0:
+            label = r'reference  $\rho$'
+        ax.plot(radial_coordinate, relative_attenuation, '-', linewidth=2, label=label,
+                color=(max(abs(shift*10), shift*50), abs(shift*30), max(abs(shift*10), -shift*50)))
         relative_attenuation.fillna(0).to_csv(export_directory + '/ionyeald' + str(int(100*(1+shift)))+'.dat', index=False, header=False)
-    ax.legend()
-    matplotlib.pyplot.savefig(export_directory+'/attenuation.pdf')
-    matplotlib.pyplot.savefig(export_directory+'/attenuation.svg')
+    ax.legend(labelspacing=0.3, borderpad=0)
+    matplotlib.pyplot.xlim(0.6, 0.75)
+    matplotlib.pyplot.subplots_adjust(left=0.18, right=0.95, top=0.85, bottom=0.15)
+    matplotlib.pyplot.savefig(export_directory+'/attenuation_shaked.pdf')
+    matplotlib.pyplot.savefig(export_directory+'/attenuation_shaked.svg')
 
 
 def shake_flux_surface_silent(species, shot_number, time, energy):
-
-
     beamlet_geometry = BeamletGeometry()
     beamlet_geometry.rad = numpy.linspace(0.72, 0.6, 1000)
     beamlet_geometry.set_with_value(0, 'z', 'rad')
@@ -95,7 +102,7 @@ def shake_flux_surface_silent(species, shot_number, time, energy):
     return numpy.max(diff)
 
 
-if __name__ == "__main__":
+def test_multi_energy():
     a_species = 'Li'
     a_shot_number = '17178'
     a_time = '1097'
@@ -118,9 +125,21 @@ def shake_flux_surface_silent(species, shot_number, time, energy):
     matplotlib.pyplot.savefig(export_directory+'/maxdiff.svg')
 
 
-# if __name__ == "__main__":
-#    a_species = 'Li'
-#    a_shot_number = '17178'
-#    a_time = '1097'
-#    an_energy = 80
-#    shake_flux_surface(a_species, a_shot_number, a_time, an_energy)
+def test_single_energy():
+    a_species = 'Li'
+    a_shot_number = '17178'
+    a_time = '1097'
+    an_energy = 80
+    shake_flux_surface(a_species, a_shot_number, a_time, an_energy)
+
+
+def test_single_energy2():
+    a_species = 'Na'
+    a_shot_number = '17178'
+    a_time = '1097'
+    an_energy = 50
+    shake_flux_surface(a_species, a_shot_number, a_time, an_energy)
+
+
+if __name__ == "__main__":
+    test_single_energy2()

sd.py

@@ -46,8 +46,8 @@ def run_core(self):
 
     def run_plotters(self):
         os.chdir(get_home_directory())
-        detector(self.shot_number, self.time, self.runnumber)
-        detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par)
+        detector(self.shot_number, self.time, self.runnumber, self.title)
+        detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par, self.title)
         if self.is_trajectory_detailed:
             traj_plotter(self.shot_number, self.time, self.runnumber, self.detector_par, self.species, self.energy, self.title)
 
@@ -66,7 +66,7 @@ def write_parameter_file(self):
         f.write("energy=" + self.energy + " #keV\n")
         f.write("toroidal_deflection=0 #degree; + 200 V deflection\n")
         f.write("diameter=5 #mm\n")
-        f.write("particles=10000\n")
+        f.write("particles=1000\n")
         f.write("detector='" + self.detector_par + "'\n")
         f.write("electric_field_module=0\n")
         f.write("detector_mask='final'\n")
@@ -83,9 +83,10 @@ def delete_parameter_file(self):
     a_time = 1097
     a_species = 'Li'
     an_energy = 70
-    z_det = 0.253
-    a_detector = f"0.6846,{z_det},0.0,38,0"
-    a_runnumber = '20240523214835'
+    z_det = 22
+    t_det = -3.0
+    a_detector = f"0.6846,{z_det/100},{t_det/100},38,0"
+    a_runnumber = None#'20240530015119'
     is_plot_trajectory = True
-    a_title = f"Reference discharge\n{a_species} beam, $E ={an_energy}$ keV\n$Z_D={z_det*100}$ cm"
+    a_title = f"reference magnetic field\n{a_species} beam, $E ={an_energy}$ keV\n$Z_D={z_det}$ cm, $T_D={t_det}$ cm"
     SD(a_shot_number, a_time, a_species, an_energy, a_detector, a_runnumber, is_trajectory_detailed=is_plot_trajectory, title=a_title)