PyRate.py

import numpy as np
import math
import FreeSimpleGUI as sg
import PyRate_Graphs as graphs
import PyRate_BeamParameters as beam
import PyRate_FieldSimulation as field
import PyRate_AnalyseData as anal
import PyRate_CustomScanPattern as pattern

 ###############
# DEFINITIONS #
###############
#
# x = i = horizontal = columns
# y = j = vertical = rows
#
# Beam Parameters
scan_speed_horizontal = 10640 #8929 #26400  # mm/s
scan_speed_vertical = 2750 #2907 #7600  # mm/s
energy = 245  # MeV
cyclotron_current = 800  # nA

# Target Parameters
shape = 'rectangle'
target_horizontal_size = 23  # mm 
target_vertical_size = 100  # mm

# Simulation Parameters
resolution = 0.25  # mm

# Output Parameters
display_graphs = True
analyse_data = True

# Dose Map Parameters
dose_map_dose_target = 10  # Gy
overscan_fraction = 5/3 
spot_spacing = 2.5  # mm
define_spot_order_raster = True
define_spot_order_from_file = False
define_spot_order_line = False
define_spot_order_spiral_in = False
define_spot_order_spiral_out = False
define_spot_order_circular_raster = False

## Simulating Custom Set Ups
clinical_spots = False
varian_limit_nozzle_current = False
clinical_transmission = False
custom_spot_size = False
horizontal_spot_size = 0.0
vertical_spot_size = 0.0



###############################################
# Start of Program
#
def run_simulation():
    print("-------------------")
    print("Start of Simulation")
    print("-------------------", end="\n\n")
    print(f"Resolution: {resolution} mm/pixel", end="\n")
    print(f'Using beamline commissioning data from {dataset}\n\n')
    print(">>> Simulating spot...", end="\n\n")

    print("Beam Details")
    print("-------------------")
    print(f'Beam Energy: {energy} MeV')
    print(f'Cyclotron Current: {cyclotron_current} nA')

    if clinical_spots:
        spot_sigma = (3.00, 3.00)  # clinical spot sizes
    elif custom_spot_size:
        spot_sigma = (float(horizontal_spot_size), float(vertical_spot_size))
    else:
        spot_sigma = beam.calculate_spot_sigma(energy, dataset)


    spot_pixels_max = beam.calculate_spot_extent(spot_sigma, resolution)
    spot_peak_dose = beam.gaussian_value(0, 0, (0, 0), spot_sigma[0], spot_sigma[1])
    spot_spacing_pixels = spot_spacing / resolution
    overscan_horizontal = math.ceil((spot_sigma[0] * overscan_fraction) / resolution)
    overscan_vertical = math.ceil((spot_sigma[1] * overscan_fraction) / resolution)
    overscan = (overscan_horizontal, overscan_vertical)
    spot_details = (spot_sigma, spot_pixels_max, spot_spacing_pixels)
    kernel = beam.create_spot_kernel(spot_details, resolution, 1)
    spot_peak_dose_rate, transmission, nozzle_current, reference_peak_dose_rate = (
        beam.calculate_spot_peak_dose_rate(kernel, energy, cyclotron_current, clinical_transmission,
                                           varian_limit_nozzle_current, dataset))
    print(f'Transmission: {transmission * 100:.3f} %')
    print(f'Nozzle Current: {nozzle_current:.2f} nA')
    print("", end="\n")
    print("Spot Details")
    print("-------------------")
    print(f'Reference Peak Dose Rate: {reference_peak_dose_rate:.2f} Gy/s')
    print(f'Sigma: ({spot_sigma[0]:.2f}, {spot_sigma[1]:.2f}) mm')
    print(f'Peak Signal: {spot_peak_dose:.5f}')
    print(f'Peak Dose Rate: {spot_peak_dose_rate:.2f} Gy/s')
    print(f'Size of Kernel: {2 * spot_pixels_max} pixels')
    print("", end="\n")

   
    test_field_details = field.create_dose_map(target_details, spot_details, overscan, resolution)

    # Determine spot map shape and calculate positions
    spot_map_order = []
    if define_spot_order_from_file:
        spot_map_order = np.loadtxt("spot_map_order.txt", dtype=int, delimiter=",")
    elif define_spot_order_line:
        spot_map_order = field.generate_line_spot_map_order(test_field_details)
    elif define_spot_order_spiral_in:
        spot_map_order = field.generate_spiral_spot_map_order(test_field_details, "in")
    elif define_spot_order_spiral_out:
        spot_map_order = field.generate_spiral_spot_map_order(test_field_details, "out")
    elif define_spot_order_raster:
        spot_map_order = field.generate_raster_dose_map(test_field_details)
    elif define_spot_order_circular_raster:
        spot_map_order = field.generate_circular_raster(test_field_details)

    spot_map = field.generate_spot_map(spot_map_order, spot_details)

    # Creating a point of interest for which we will calculate the average dose rate
    # data as we are progressing through the spots. Selected to be the middle of the
    # target.
    # point of interest is stored as [vertical, horizontal]
    mid_vertical = int((test_field_details[0].shape[0]) / 2)
    mid_horizontal = int((test_field_details[0].shape[1]) / 2)
    max_horizontal_spot = 0

    for x, y in spot_map:
        if x > max_horizontal_spot:
            max_horizontal_spot = x

    point_of_interest = []
    point_of_interest.append([mid_vertical, mid_horizontal])
    point_of_interest.append([mid_vertical, mid_horizontal + math.floor((target_horizontal_size / 2) / resolution / 2)])
    point_of_interest.append([mid_vertical, mid_horizontal + math.floor((target_horizontal_size / 2) / resolution)])
    point_of_interest.append([mid_vertical, max_horizontal_spot])

    print(">>> Simulating TEST field delivery...", end="\n\n")
    test_field_plateau_dose = field.simulate_field_with_defined_spots(test_field_details, kernel, spot_map, point_of_interest)

    field_to_spot_ratio = test_field_plateau_dose / spot_peak_dose
    spot_dose_target = dose_map_dose_target / field_to_spot_ratio
    spot_delivery_time = spot_dose_target / spot_peak_dose_rate
    dose_per_signal = dose_map_dose_target / test_field_plateau_dose
    kernel = beam.create_spot_kernel(spot_details, resolution, dose_per_signal)
    time_to_deliver = (spot_peak_dose * dose_per_signal) / spot_peak_dose_rate

    if display_graphs:
        graphs.plot_spot(kernel, resolution)
        graphs.plot_map_with_target(test_field_details[0], resolution, target_details, spot_pixels_max, overscan)

    print("Test Details")
    print("-------------------")
    print(f'Target Dose: {dose_map_dose_target} Gy')
    print(f'Field Peak Signal: {test_field_plateau_dose:.5f}')
    print(f'Ratio: {field_to_spot_ratio:.2f}')
    print(f'Dose Per Signal: {dose_per_signal:.2f} Gy')
    print(f'Spot Peak Dose: {spot_peak_dose*dose_per_signal:.3f} Gy')
    print(f'Time To Deliver Spot: {time_to_deliver*1000:.2f} ms')
    print("", end="\n")

    # Now we create the dose map we will actually use to simulate covering the specified
    # target in a homogeneous dose. The target size is padded with the maximum extent of
    # the spot calculated earlier and spots are placed to evenly cover this area. This
    # ensures that spots at the edge of this area are contributing ~0.1% of the dose to
    # the target volume (if using 3 fold sigma size) and so we can assume the target is
    # homogeneously covered. This size is padded with the same amount again for display
    # purposes.
    print("Field Details")
    print("-------------------")
    dose_map_details = field.create_dose_map(target_details, spot_details, overscan, resolution)

    print_details = dose_map_details[3]

    print(f'Target Size (HxV): {target_horizontal_size} x {target_vertical_size} mm')
    print(f'Target size (HxV): {print_details[0]} x {print_details[1]} pixels')
    print(f'Dose area (HxV): {print_details[2] * resolution} x {print_details[3] * resolution} mm')
    print(f'Dose area (HxV): {print_details[2]} x {print_details[3]} pixels')
    print(f'Spot spacing: {spot_spacing_pixels} pixels')
    print(f'Spot number (HxV): {print_details[4]} x {print_details[5]}')
    print(f'Display area (HxV): {print_details[6]} x {print_details[7]} pixel')
    print('', end='\n')

    # Creating a point of interest for which we will calculate the average dose rate
    # data as we are progressing through the spots. Selected to be the middle of the
    # target.
    # point of interest is stored as [vertical, horizontal]
    mid_vertical = int((dose_map_details[0].shape[0]) / 2)
    mid_horizontal = int((dose_map_details[0].shape[1]) / 2)
    max_horizontal_spot = 0

    for x, y in spot_map:
        if x > max_horizontal_spot:
            max_horizontal_spot = x


    # ------------- start looping over spots ----------------
    print(">>> Simulating Field Delivery...", end="\n\n")
    data, dose_map_history = field.simulate_field_and_time_with_defined_spots(dose_map_details,
                                                                              resolution,
                                                                              spot_delivery_time,
                                                                              scan_speed_horizontal,
                                                                              scan_speed_vertical,
                                                                              kernel,
                                                                              spot_map,
                                                                              point_of_interest)
    # -------------- done looping over spots ----------------

    dose_map = dose_map_history[len(dose_map_history) - 1][1]
    if display_graphs:
        graphs.plot_map_with_details(dose_map, resolution, point_of_interest, target_horizontal_size,
                                     target_vertical_size, spot_pixels_max, overscan, shape, spot_map)
        graphs.plot_simple_heat_map(dose_map, resolution, "Max Dose", "Dose (Gy)", target_details, spot_pixels_max, overscan)

    print("-------------------")
    print("Finished simulation")
    print("-------------------")
    print("", end="\n")
    print(">>> Running Analysis")
    print("> Writing Data for Points of Interest...")
    # Create file to output dose accumulation data into
    spot_t = f'{round(spot_sigma[0], 1)}x{round(spot_sigma[1], 1)}_spot'
    target_t = f'{target_horizontal_size}x{target_vertical_size}_target'
    dose_t = f'{dose_map_dose_target}Gy_{cyclotron_current}nA'
    pattern_t = ''
    if define_spot_order_from_file:
        pattern_t = 'FromFile'
    elif define_spot_order_line:
        pattern_t = 'Line'
    elif define_spot_order_spiral_in:
        pattern_t = 'SpiralIn'
    elif define_spot_order_spiral_out:
        pattern_t = 'SpiralOut'
    elif define_spot_order_raster:
        pattern_t = 'Raster'
    elif define_spot_order_circular_raster:
        pattern_t = 'RasterCircle'
    file_name = f'average_dose_rate_{spot_t}_{target_t}_{dose_t}_{pattern_t}.dat'
    out_file = open(file_name, 'w')
    # Pass file a dose accumulation data to write and calculate dose
    # rate for different thresholds
    time_to_target_data = anal.write_data(data, out_file, dose_map_dose_target)
    out_file.close()

    print("> Plotting Graph for Points of Interest...")
    graphs.make_graph(file_name,
                      spot_sigma,
                      time_to_target_data,
                      target_horizontal_size,
                      target_vertical_size,
                      dose_map_dose_target,
                      spot_peak_dose_rate,
                      cyclotron_current,
                      pattern_t)

    if analyse_data:
        print("> Calculating Time to Thresholds for Area...")
        time_to_dose, dose_rate_map = anal.dose_map_into_dose_rate(dose_map_history)

        if display_graphs:
            print("> Plotting Graph for Area...")
            graphs.plot_simple_heat_map(time_to_dose[7], resolution, "Time to Max Dose", "Time (s)", target_details, spot_pixels_max, overscan)
            graphs.plot_simple_heat_map(dose_rate_map[0], resolution, "Dose Rate 100%", "Dose Rate (Gy/s)", target_details, spot_pixels_max, overscan)
            graphs.plot_simple_heat_map(dose_rate_map[1], resolution, "Dose Rate 90%", "Dose Rate (Gy/s)", target_details, spot_pixels_max, overscan)
            graphs.plot_simple_heat_map(dose_rate_map[2], resolution, "Dose Rate 80%", "Dose Rate (Gy/s)", target_details, spot_pixels_max, overscan)
            graphs.plot_simple_heat_map(dose_rate_map[3], resolution, "Dose Rate 70%", "Dose Rate (Gy/s)", target_details, spot_pixels_max, overscan)

        averageD, stdevD, HID = anal.calculate_statistics(dose_map,
                                                     target_details,
                                                     spot_pixels_max,
                                                     overscan,
                                                     resolution)
        averageDR100, stdevDR100, HIDR100 = anal.calculate_statistics(dose_rate_map[0],
                                                             target_details,
                                                             spot_pixels_max,
                                                             overscan,
                                                             resolution)
        averageDR90, stdevDR90, HIDR90 = anal.calculate_statistics(dose_rate_map[1],
                                                           target_details,
                                                           spot_pixels_max,
                                                           overscan,
                                                           resolution)
        averageDR80, stdevDR80, HIDR80 = anal.calculate_statistics(dose_rate_map[2],
                                                           target_details,
                                                           spot_pixels_max,
                                                           overscan,
                                                           resolution)
        averageDR70, stdevDR70, HIDR70 = anal.calculate_statistics(dose_rate_map[3], target_details, spot_pixels_max, overscan,
                                                           resolution)

        print('', end="\n")
        print('-------------------')
        print('Summary')
        print('-------------------', end="\n\n")
        print(f'Dose to target: {averageD:.2f} +/- {stdevD:.2f} Gy (HI: {HID:.2f}%)')
        print(f'100% Dose Rate: {averageDR100:.2f} +/- {stdevDR100:.2f} Gy/s (HI: {HIDR100:.2f}%)')
        print(f'90% Dose Rate: {averageDR90:.2f} +/- {stdevDR90:.2f} Gy/s (HI: {HIDR90:.2f}%)')
        print(f'80% Dose Rate: {averageDR80:.2f} +/- {stdevDR80:.2f} Gy/s (HI: {HIDR80:.2f}%)')
        print(f'70% Dose Rate: {averageDR70:.2f} +/- {stdevDR70:.2f} Gy/s (HI: {HIDR70:.2f}%)')
        print('', end="\n")

        out_file = open(file_name, 'a')
        out_file.write('\n\n')
        out_file.write(f'Dose to target: {averageD:.2f} +/- {stdevD:.2f} Gy (HI: {HID:.2f}%)\n')
        out_file.write(f'100% Dose Rate: {averageDR100:.2f} +/- {stdevDR100:.2f} Gy/s (HI: {HIDR100:.2f}%)\n')
        out_file.write(f'90% Dose Rate: {averageDR90:.2f} +/- {stdevDR90:.2f} Gy/s (HI: {HIDR90:.2f}%)\n')
        out_file.write(f'80% Dose Rate: {averageDR80:.2f} +/- {stdevDR80:.2f} Gy/s (HI: {HIDR80:.2f}%)\n')
        out_file.write(f'70% Dose Rate: {averageDR70:.2f} +/- {stdevDR70:.2f} Gy/s (HI: {HIDR70:.2f}%)\n\n')

        out_file.write("-------------------\n")
        out_file.write("Simulation Details\n")
        out_file.write("-------------------\n\n")
        out_file.write(f'Resolution (mm/pixel): {resolution}\n\n')
        out_file.write("Beam Details\n")
        out_file.write("-------------------\n")
        out_file.write(f'Beam Energy: {energy} MeV\n')
        out_file.write(f'Cyclotron Current: {cyclotron_current} nA\n')
        out_file.write(f'Transmission: {transmission * 100:.3f} %\n')
        out_file.write(f'Nozzle Current: {nozzle_current:.2f} nA\n\n')
        out_file.write("Spot Details\n")
        out_file.write("-------------------\n")
        out_file.write(f'Sigma: ({spot_sigma[0]:.2f}, {spot_sigma[1]:.2f}) mm\n')
        out_file.write(f'Peak Signal: {spot_peak_dose:.5f}\n')
        out_file.write(f'Peak Dose Rate: {spot_peak_dose_rate:.2f} Gy/s\n')
        out_file.write(f'Size of Kernel: {2 * spot_pixels_max} pixels\n\n')
        out_file.write("Test Details\n")
        out_file.write("-------------------\n")
        out_file.write(f'Target Dose: {dose_map_dose_target} Gy\n')
        out_file.write(f'Field Peak Signal: {test_field_plateau_dose:.5f}\n')
        out_file.write(f'Ratio: {field_to_spot_ratio:.2f}\n')
        out_file.write(f'Dose Per Signal: {dose_per_signal:.2f} Gy\n')
        out_file.write(f'Spot Peak Dose: {spot_peak_dose*dose_per_signal:.3f} Gy\n')
        out_file.write(f'Time To Deliver Spot: {time_to_deliver*1000:.2f} ms\n\n')
        out_file.write("Field Details\n")
        out_file.write("-------------------\n")
        out_file.write(f'Target Size (HxV): {target_horizontal_size} x {target_vertical_size} mm\n')
        out_file.write(f'Target size (HxV): {print_details[0]} x {print_details[1]} pixels\n')
        out_file.write(f'Dose area (HxV): {print_details[2] * resolution} x {print_details[3] * resolution} mm\n')
        out_file.write(f'Dose area (HxV): {print_details[2]} x {print_details[3]} pixels\n')
        out_file.write(f'Spot spacing: {spot_spacing_pixels} pixels\n')
        out_file.write(f'Spot number (HxV): {print_details[4]} x {print_details[5]}\n')
        out_file.write(f'Display area (HxV): {print_details[6]} x {print_details[7]} pixel\n\n')

        out_file.close()

    print('Done.')

sg.theme('DarkTeal1')
###########################
# Defining Layout Elements
#
data_output = [[sg.Text('Simulation Parameters', font=('Ubuntu',15))]]
data_output.append([sg.Text('Resolution (mm/pixel):\t'), sg.Input('0.25', key='-RES-', size=15)])
data_output.append([sg.Text('Beam Parameter Dataset:\t'), sg.Combo(['Pre-2023', 'January-2023', 'January-2024', 'January-2025'], key='-DATASET-', default_value='January-2025')])
data_output.append(sg.Checkbox('Display Graph', default=True, key='-GRAPH-'))
data_output.append(sg.Checkbox('Analyse Data', default=True, key='-ANAL-'))

beam_parameters = [[sg.Text('Beam Parameters', font=('Ubuntu',15))]]
beam_parameters.append([sg.Text('Energy (MeV):\t\t'), sg.Input('245', key='-ENERGY-', size=15)])
beam_parameters.append([sg.Text('Cyclotron Current (nA):\t'), sg.Input('800', key='-CURRENT-', size=15)])
beam_parameters.append([sg.Text('Horizontal Scan Speed (mm/s):'), sg.Input('10640', key='-HSCAN-', size=13)])
beam_parameters.append([sg.Text('Vertical Scan Speed (mm/s):'), sg.Input('2750', key='-VSCAN-', size=15)])

target_parameters = [[sg.Text('Target Parameters', font=('Ubuntu',15))]]
target_parameters.append([sg.Text('Target Shape:'),
                          sg.Radio('rectangle', "RADIO_SHAPE", default=True, key='-REC-'),
                          sg.Radio('circle', "RADIO_SHAPE", default=False, key='-CIR-')])
target_parameters.append([sg.Text('Horizontal Size (mm):\t'), sg.Input('23', key='-HSIZE-', size=15)])
target_parameters.append([sg.Text('Vertical Size (mm):\t'), sg.Input('100', key='-VSIZE-', size=15)])

dose_map_parameters = [[sg.Text('Dose Map Parameters', font=('Ubuntu',15))]]
dose_map_parameters.append([sg.Text('Target Dose (Gy):\t\t'), sg.Input('10', key='-DOSE-', size=15)])
dose_map_parameters.append([sg.Text('Spot Spacing (mm):\t'), sg.Input('2.5', key='-SPACE-', size=15),
                            sg.Push(), sg.Button('Custom Map'), sg.Push()])
dose_map_parameters.append([sg.Text('Overscan (sigma):\t\t'), sg.Input('1.66666666667', key='-OVERSCAN-', size=15)])
dose_map_parameters.append([sg.Text('Scan Pattern:'),
                            sg.Radio('Raster', "RADIO_PATTERN", default=True, key='-RASTER-'),
                            sg.Radio('Circular Raster', "RADIO_PATTERN", default=False, key='-CIRCRASTER-'),
                            sg.Radio('From File', "RADIO_PATTERN", default=False, key='-FROMFILE-'),
                            sg.Radio('Line', "RADIO_PATTERN", default=False, key='-LINE-'),
                            sg.Radio('Spiral In', "RADIO_PATTERN", default=False, key='-SPIRALIN-'),
                            sg.Radio('Spiral Out', "RADIO_PATTERN", default=False, key='-SPIRALOUT-')])

hack_parameters = [[sg.Text('Alternative Beam Simulations', font=('Ubuntu',15))]]
hack_parameters.append([sg.Text('Spot Sizes:'),
                        sg.Radio('Calculated', "RADIO_SPOT", default=True, key='-SPOT_NONE-'),
                        sg.Radio('Clinical (3.00x3.00 mm)', "RADIO_SPOT", default=False, key='-SPOT_CLIN-')])
hack_parameters.append([sg.Text('\t   '),
                        sg.Radio('Custom:', 'RADIO_SPOT', default=False, key='-SPOT_CUSTOM-', enable_events=True),
                        sg.Text('Horizontal (mm):'), sg.Input('-', key='-H_SPOT_SIZE-', size=8, disabled=True),
                        sg.Text('Vertical (mm):'), sg.Input('-', key='-V_SPOT_SIZE-', size=8, disabled=True)])
hack_parameters.append([sg.Text('Transmission:'),
                        sg.Radio('Calculated', "RADIO_TRAN", default=True, key='-TRAN_NONE-'),
                        sg.Radio('Max Limit (215 nA)', "RADIO_TRAN", default=False, key='-TRAN_VAR-'),
                        sg.Radio('Clinical (10%)', "RADIO_TRAN", default=False, key='-TRAN_CLIN-')])

###################
# Creating Windows
#
layout = [[sg.Push(), sg.Text('PyRate', font=('Ubuntu',24, 'bold')), sg.Push()],
          [sg.Push(), sg.Text('Dose Rate Simulator For Spot Scanning Beams', font=('Ubuntu',14)), sg.Push()],
          [sg.Push(), sg.Text('Yarrrr...', font=('Ubuntu',8)), sg.Push()],
          [sg.HorizontalSeparator()],
          data_output,
          [sg.HorizontalSeparator()],
          beam_parameters,
          [sg.HorizontalSeparator()],
          target_parameters,
          [sg.HorizontalSeparator()],
          dose_map_parameters,
          [sg.HorizontalSeparator()],
          hack_parameters,
          [sg.HorizontalSeparator()],
          [sg.Button('Run'), sg.Button('Finish')]]

window = sg.Window('PyRate', layout)

while True:
    event, values = window.read()
    if event == sg.WIN_CLOSED or event == 'Finish':
        break
    elif event == 'Custom Map':
        pattern.open_window()
        window['-FROMFILE-'].update(True)
    elif event == '-SPOT_CUSTOM-':
        print('Selecting Custom Spot Size')
        window['-H_SPOT_SIZE-'].update(disabled=False)
        window['-V_SPOT_SIZE-'].update(disabled=False)
    elif event == 'Run':

        #####################
        # Setting Parameters
        #
        resolution = float(values['-RES-'])
        dataset = values['-DATASET-']
        print(dataset)
        display_graphs = values['-GRAPH-']
        analyse_data = values['-ANAL-']

        energy = float(values['-ENERGY-'])
        cyclotron_current = float(values['-CURRENT-'])
        scan_speed_horizontal = float(values['-HSCAN-'])
        scan_speed_vertical = float(values['-VSCAN-'])

        if values['-REC-']:
            shape = 'rectangle'
        if values['-CIR-']:
            shape = 'circle'
        target_horizontal_size = float(values['-HSIZE-'])
        target_vertical_size = float(values['-VSIZE-'])
        target_details = [target_horizontal_size, target_vertical_size, shape]

        dose_map_dose_target = float(values['-DOSE-'])
        spot_spacing = float(values['-SPACE-'])
        overscan_fraction = float(values['-OVERSCAN-'])
        define_spot_order_raster = values['-RASTER-']
        define_spot_order_from_file = values['-FROMFILE-']
        define_spot_order_line = values['-LINE-']
        define_spot_order_spiral_in = values['-SPIRALIN-']
        define_spot_order_spiral_out = values['-SPIRALOUT-']
        define_spot_order_circular_raster = values['-CIRCRASTER-']

        clinical_spots = values['-SPOT_CLIN-']
        varian_limit_nozzle_current = values['-TRAN_VAR-']
        clinical_transmission = values['-TRAN_CLIN-']
        custom_spot_size = values['-SPOT_CUSTOM-']
        horizontal_spot_size = values['-H_SPOT_SIZE-']
        vertical_spot_size = values['-V_SPOT_SIZE-']

        #####################
        # Running Simulation
        #
        run_simulation()

window.close()