diff --git a/examples/animated_corrugated_nanowire.py b/examples/animated_corrugated_nanowire.py index fa13a0a..e159d24 100644 --- a/examples/animated_corrugated_nanowire.py +++ b/examples/animated_corrugated_nanowire.py @@ -16,7 +16,6 @@ OUTPUT_SCATTERING_MAP = False OUTPUT_RAW_THERMAL_MAP = True NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Animation parameters: diff --git a/examples/anisotropy_study_horizontal.py b/examples/anisotropy_study_horizontal.py new file mode 100644 index 0000000..bd93a66 --- /dev/null +++ b/examples/anisotropy_study_horizontal.py @@ -0,0 +1,74 @@ +"""Config file to simulate a membrane with a staggered lattice of rectangular slits""" + +import numpy as np + + +# General parameters: +OUTPUT_FOLDER_NAME = 'Anisotropy study horizontal' +NUMBER_OF_PHONONS = 1000 +NUMBER_OF_TIMESTEPS = 30000 +NUMBER_OF_NODES = 400 +TIMESTEP = 1e-12 +T = 4.0 + +# Map & profiles parameters: +NUMBER_OF_PIXELS_X = 100 +NUMBER_OF_PIXELS_Y = 100 +NUMBER_OF_TIMEFRAMES = 6 + +# Material parameters: +MEDIA = 'Si' +SPECIFIC_HEAT_CAPACITY = 0.0176 # [J/kg/K] for Si at 4 K + + +# System dimensions [m]: +THICKNESS = 150e-9 +WIDTH = 1300e-9 +LENGTH = 1300e-9 + +# STRUCTURE: + +# ----> +#------------------- +# +# H C +# +#------------------- + +INCLUDE_RIGHT_SIDEWALL = False +INCLUDE_LEFT_SIDEWALL = False +INCLUDE_TOP_SIDEWALL = True +INCLUDE_BOTTOM_SIDEWALL = True + +# Hot and cold sides [m]: +COLD_SIDE_POSITION = 'right' +HOT_SIDE_POSITION = 'left' +HOT_SIDE_X = -WIDTH/2 +HOT_SIDE_Y = LENGTH/2 +HOT_SIDE_WIDTH_X = 0 +HOT_SIDE_WIDTH_Y = LENGTH +HOT_SIDE_ANGLE_DISTRIBUTION = 'random_right' + + +# Hole array parameters [m]: +INCLUDE_HOLES = True +CIRCULAR_HOLE_DIAMETER = 0 +RECTANGULAR_HOLE_SIDE_X = 200e-9 +RECTANGULAR_HOLE_SIDE_Y = 100e-9 +PERIOD_X = 300e-9 +PERIOD_Y = 300e-9 + + +# Staggered attice of holes: +FIRST_HOLE_COORDINATE = 200e-9 +NUMBER_OF_PERIODS_X = 4 +NUMBER_OF_PERIODS_Y = 4 +HOLE_COORDINATES = np.zeros((NUMBER_OF_PERIODS_X * NUMBER_OF_PERIODS_Y, 3)) +HOLE_SHAPES = ['rectangle' for x in range(HOLE_COORDINATES.shape[0])] +hole_number = 0 +for i in range(NUMBER_OF_PERIODS_Y): + for j in range(NUMBER_OF_PERIODS_X): + HOLE_COORDINATES[hole_number, 0] = -(NUMBER_OF_PERIODS_X - 1) * PERIOD_X / 2 + j * PERIOD_X + HOLE_COORDINATES[hole_number, 1] = FIRST_HOLE_COORDINATE + i * PERIOD_Y + HOLE_COORDINATES[hole_number, 2] = 0 + hole_number += 1 diff --git a/examples/anisotropy_study_vertical.py b/examples/anisotropy_study_vertical.py new file mode 100644 index 0000000..1af7335 --- /dev/null +++ b/examples/anisotropy_study_vertical.py @@ -0,0 +1,74 @@ +"""Config file to simulate a membrane with a staggered lattice of rectangular slits +Here we impose thermal gradient in vertical direction""" + +import numpy as np + + +# General parameters: +OUTPUT_FOLDER_NAME = 'Anisotropy study vertical' +NUMBER_OF_PHONONS = 1000 +NUMBER_OF_TIMESTEPS = 30000 +NUMBER_OF_NODES = 400 +TIMESTEP = 1e-12 +T = 4.0 + +# Map & profiles parameters: +NUMBER_OF_PIXELS_X = 100 +NUMBER_OF_PIXELS_Y = 100 +NUMBER_OF_TIMEFRAMES = 6 + +# Material parameters: +MEDIA = 'Si' +SPECIFIC_HEAT_CAPACITY = 0.0176 # [J/kg/K] for Si at 4 K + + +# System dimensions [m]: +THICKNESS = 150e-9 +WIDTH = 1300e-9 +LENGTH = 1300e-9 + +# STRUCTURE: + +# | C | +# | | ^ +# | | | +# | | | +# | H | + +INCLUDE_RIGHT_SIDEWALL = True +INCLUDE_LEFT_SIDEWALL = True +INCLUDE_TOP_SIDEWALL = False +INCLUDE_BOTTOM_SIDEWALL = False + +# Hot and cold sides [m]: +COLD_SIDE_POSITION = 'top' +HOT_SIDE_POSITION = 'bottom' +HOT_SIDE_X = 0 +HOT_SIDE_Y = 0 +HOT_SIDE_WIDTH_X = WIDTH +HOT_SIDE_WIDTH_Y = 0 +HOT_SIDE_ANGLE_DISTRIBUTION = 'random_up' + +# Hole array parameters [m]: +INCLUDE_HOLES = True +CIRCULAR_HOLE_DIAMETER = 0 +RECTANGULAR_HOLE_SIDE_X = 200e-9 +RECTANGULAR_HOLE_SIDE_Y = 100e-9 +PERIOD_X = 300e-9 +PERIOD_Y = 300e-9 + + +# Staggered attice of holes: +FIRST_HOLE_COORDINATE = 200e-9 +NUMBER_OF_PERIODS_X = 4 +NUMBER_OF_PERIODS_Y = 4 +HOLE_COORDINATES = np.zeros((NUMBER_OF_PERIODS_X * NUMBER_OF_PERIODS_Y, 3)) +HOLE_SHAPES = ['rectangle' for x in range(HOLE_COORDINATES.shape[0])] +hole_number = 0 +for i in range(NUMBER_OF_PERIODS_Y): + for j in range(NUMBER_OF_PERIODS_X): + HOLE_COORDINATES[hole_number, 0] = -(NUMBER_OF_PERIODS_X - 1) * PERIOD_X / 2 + j * PERIOD_X + HOLE_COORDINATES[hole_number, 1] = FIRST_HOLE_COORDINATE + i * PERIOD_Y + HOLE_COORDINATES[hole_number, 2] = 0 + hole_number += 1 + diff --git a/examples/fishbone_nanowire.py b/examples/fishbone_nanowire.py index a192516..6c4453b 100644 --- a/examples/fishbone_nanowire.py +++ b/examples/fishbone_nanowire.py @@ -14,7 +14,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Map & profiles parameters: diff --git a/examples/parabolic_lens_focusing.py b/examples/parabolic_lens_focusing.py index 9184cd0..dfc267a 100644 --- a/examples/parabolic_lens_focusing.py +++ b/examples/parabolic_lens_focusing.py @@ -13,7 +13,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'directional' # Map & profiles parameters: @@ -43,6 +42,7 @@ COLD_SIDE_POSITION = 'top' HOT_SIDE_X = 0.0 HOT_SIDE_WIDTH_X = WIDTH +HOT_SIDE_ANGLE_DISTRIBUTION = 'directional' # Roughness [m]: SIDE_WALL_ROUGHNESS = 2e-9 diff --git a/examples/phononic_crystal.py b/examples/phononic_crystal.py index d74d671..7b009c4 100644 --- a/examples/phononic_crystal.py +++ b/examples/phononic_crystal.py @@ -15,7 +15,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Map & profiles parameters: diff --git a/examples/simple_nanowire.py b/examples/simple_nanowire.py index ff4fde7..b2a60a5 100644 --- a/examples/simple_nanowire.py +++ b/examples/simple_nanowire.py @@ -15,7 +15,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Map & profiles parameters: diff --git a/examples/slits_array.py b/examples/slits_array.py index 59d769c..cd481a7 100644 --- a/examples/slits_array.py +++ b/examples/slits_array.py @@ -15,7 +15,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Map & profiles parameters: diff --git a/examples/triangles_array.py b/examples/triangles_array.py index 5619d3e..2050f2d 100644 --- a/examples/triangles_array.py +++ b/examples/triangles_array.py @@ -15,7 +15,6 @@ OUTPUT_RAW_THERMAL_MAP = True OUTPUT_TRAJECTORIES_OF_FIRST = 30 NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = 'random' # Map & profiles parameters: diff --git a/freepaths/__main__.py b/freepaths/__main__.py index c8983ce..5068349 100644 --- a/freepaths/__main__.py +++ b/freepaths/__main__.py @@ -5,13 +5,13 @@ import freepaths.main_tracing import freepaths.main_mfp_sampling -__version__ = "1.2.2" +__version__ = "1.3" # Parse user arguments: parser = argparse.ArgumentParser( prog = 'FreePATHS', description = 'Monte Carlo simulator', - epilog = 'For more information, visit: https://github.com/anufrievroman/freepaths' + epilog = 'For more information, visit: https://anufrievroman.gitbook.io/freepaths' ) parser.add_argument('input_file', nargs='?', default=None, help='The input file') parser.add_argument("-s", "--sampling", help="Run in MFP sampling mode", action="store_true") diff --git a/freepaths/animation.py b/freepaths/animation.py index 7163b88..82c7384 100644 --- a/freepaths/animation.py +++ b/freepaths/animation.py @@ -80,9 +80,18 @@ def generate_animation_xy(): filenames = os.listdir("Frames/") for filename in filenames: images.append(imageio.imread(f"Frames/{filename}")) + + # Create a GIF file: imageio.mimsave("Animated paths XY.gif", images, fps=cf.output_animation_fps, subrectangles=True) + # Create an MP4 file: + # writer = imageio.get_writer('Animated paths XY.mp4', fps=cf.output_animation_fps) + + # for im in images: + # writer.append_data(imageio.imread(im)) + # writer.close() + def delete_frames(): """Delete frames after creating the animation""" diff --git a/freepaths/animation__.py b/freepaths/animation__.py new file mode 100644 index 0000000..7163b88 --- /dev/null +++ b/freepaths/animation__.py @@ -0,0 +1,97 @@ +"""Module that creates animations from recorded phonon paths""" + +import sys +import os +import shutil +import imageio +import numpy as np +import matplotlib.pyplot as plt + +from freepaths.config import cf +from freepaths.output_structure import draw_structure + + +def generate_frames_xy(): + """Generate animation frames with phonon paths""" + + data = np.genfromtxt("Data/Phonon paths.csv", + unpack=False, + delimiter=',', + skip_header=1, + encoding='utf-8') + + # Create XY plots for each timestep where all phonon shown at the same time: + number_of_steps = np.shape(data)[0] + number_of_phonons = np.shape(data)[1]//3 + + for step in range(1, number_of_steps): + fig, ax = plt.subplots() + + # Draw the structure: + patches = draw_structure(cf) + for patch in patches: + ax.add_patch(patch) + + # Draw the paths: + for phonon_num in range(number_of_phonons): + x_coords = np.trim_zeros(data[:, 3 * phonon_num], trim='b') + y_coords = np.trim_zeros(data[:, 3 * phonon_num + 1], trim='b') + ax.plot(x_coords[:step], y_coords[:step], linewidth=0.5) + + # Plot settings: + ax.set_xlim([-0.55*cf.width*1e6, 0.55*cf.width*1e6]) + ax.set_ylim([0, cf.length*1e6]) + ax.set_aspect('equal') + ax.axis('off') + for spine in ['top', 'right', 'left', 'bottom']: + ax.spines[spine].set_visible(False) + fig.savefig(f"Frames/frame_{step:0>6}.png", dpi=600, bbox_inches="tight") + plt.close(fig) + + # Progress: + sys.stdout.write(f"\rAnimation: {step}/{number_of_steps - 1} frames") + + # Create XY plots for each step for each phonon one by one: + # cmap = plt.get_cmap("tab10") + # frame_number = 0 + # number_of_phonons = np.shape(data)[1]//3 + # for phonon_num in range(number_of_phonons): + # x_coords = np.trim_zeros(data[:, 3 * phonon_num], trim='b') + # y_coords = np.trim_zeros(data[:, 3 * phonon_num + 1], trim='b') + # steps = np.shape(x_coords)[0] + # for step in range(1, steps): + # fig, ax = plt.subplots() + # ax.plot(x_coords[:step], y_coords[:step], color=cmap(phonon_num), linewidth=0.5) + # ax.axis('off') + # for spine in ['top', 'right', 'left', 'bottom']: + # ax.spines[spine].set_visible(False) + # ax.set_xlim([-0.55*cf.width*1e6, 0.55*cf.width*1e6]) + # ax.set_ylim([0, cf.length*1e6]) + # ax.set_aspect('equal') + # fig.savefig(f"Frames/frame_{frame_number:0>6}.png", dpi=600, bbox_inches="tight") + # plt.close(fig) + # frame_number += 1 + + +def generate_animation_xy(): + """Generate animation of phonon path in XY plane""" + sys.stdout.write(f"\rAnimation: creating animation file") + images = [] + filenames = os.listdir("Frames/") + for filename in filenames: + images.append(imageio.imread(f"Frames/{filename}")) + imageio.mimsave("Animated paths XY.gif", images, + fps=cf.output_animation_fps, subrectangles=True) + + +def delete_frames(): + """Delete frames after creating the animation""" + folder_path = "Frames/" + shutil.rmtree(folder_path) + + +def create_animation(): + """Main function that creates the animation""" + generate_frames_xy() + generate_animation_xy() + delete_frames() diff --git a/freepaths/config.py b/freepaths/config.py index d831fe4..b524026 100644 --- a/freepaths/config.py +++ b/freepaths/config.py @@ -1,4 +1,4 @@ -"""Module that reads the user input file, privides default values, and converts the variables into enums""" +"""Module that reads the user input file, provides default values, and converts the variables into enums""" import sys import argparse @@ -10,7 +10,7 @@ # Parse user arguments: -WEBSITE = 'https://github.com/anufrievroman/freepaths' +WEBSITE = 'https://anufrievroman.gitbook.io/freepaths' parser = argparse.ArgumentParser(prog='FreePATHS', description='Monte Carlo simulator', epilog=f'For more information, visit: {WEBSITE}') parser.add_argument('input_file', nargs='?', default=None, help='The input file') @@ -22,7 +22,7 @@ if args.input_file: exec(open(args.input_file, encoding='utf-8').read(), globals()) else: - print("You didn't provide any input file, so we run a demo simulation:\n") + print("You didn't provide any input file, so let's run a demo simulation!\n") class Config: @@ -68,10 +68,15 @@ def __init__(self): self.thickness = THICKNESS self.width = WIDTH self.length = LENGTH + self.include_right_sidewall = INCLUDE_RIGHT_SIDEWALL + self.include_left_sidewall = INCLUDE_LEFT_SIDEWALL + self.include_top_sidewall = INCLUDE_TOP_SIDEWALL + self.include_bottom_sidewall = INCLUDE_BOTTOM_SIDEWALL # Hot and cold sides: self.frequency_detector_size = FREQUENCY_DETECTOR_SIZE self.cold_side_position = COLD_SIDE_POSITION + self.hot_side_position = HOT_SIDE_POSITION self.hot_side_x = HOT_SIDE_X self.hot_side_y = HOT_SIDE_Y self.hot_side_width_x = HOT_SIDE_WIDTH_X @@ -111,8 +116,11 @@ def __init__(self): self.pillar_height = PILLAR_HEIGHT self.pillar_wall_angle = PILLAR_WALL_ANGLE + def convert_to_enums(self): """Convert some user generated parameters into enums""" + + # Distributions: valid_distributions =[member.name.lower() for member in Distributions] if self.hot_side_angle_distribution in valid_distributions: self.hot_side_angle_distribution = Distributions[self.hot_side_angle_distribution.upper()] @@ -122,6 +130,7 @@ def convert_to_enums(self): print(*valid_distributions, sep = ", ") sys.exit() + # Materials: valid_materials = [member.name for member in Materials] if self.media in valid_materials: self.media = Materials[self.media] @@ -131,6 +140,7 @@ def convert_to_enums(self): print(*valid_materials, sep = ", ") sys.exit() + # Positions: valid_positions = [member.name.lower() for member in Positions] if self.cold_side_position in valid_positions: self.cold_side_position = Positions[self.cold_side_position.upper()] @@ -140,8 +150,17 @@ def convert_to_enums(self): print(*valid_positions, sep = ", ") sys.exit() + if self.hot_side_position in valid_positions: + self.hot_side_position = Positions[self.hot_side_position.upper()] + else: + print("ERROR: Parameter HOT_SIDE_POSITION is not set correctly.") + print("HOT_SIDE_POSITION should be one of the following:") + print(*valid_positions, sep = ", ") + sys.exit() + + def check_parameter_validity(self): - """Check if some parameteres are valid""" + """Check if various parameters are valid""" if self.number_of_phonons < self.output_trajectories_of_first: self.output_trajectories_of_first = self.number_of_phonons print("WARNING: Parameter OUTPUT_TRAJECTORIES_OF_FIRST exceeded NUMBER_OF_PHONONS.\n") @@ -155,7 +174,7 @@ def check_parameter_validity(self): print("WARNING: Parameter HOT_SIDE_Y was negative.\n") if self.hot_side_y - self.hot_side_width_y / 2 < 0: - self.self.hot_side_width_y = self.hot_side_y * 2 + self.hot_side_width_y = self.hot_side_y * 2 print("WARNING: Parameter HOT_SIDE_WIDTH_Y was too large.\n") if self.hot_side_x > self.width/2: @@ -166,6 +185,10 @@ def check_parameter_validity(self): self.hot_side_width_x = self.width print("WARNING: Parameter HOT_SIDE_WIDTH_X exceeds WIDTH.\n") + if self.cold_side_position == self.hot_side_position: + print(f"ERROR: Hot and cold sides are set at {self.cold_side_position}.") + sys.exit() + if self.output_path_animation and self.number_of_timesteps > 5000: print("WARNING: NUMBER_OF_TIMESTEPS is rather large for animation.\n") diff --git a/freepaths/data.py b/freepaths/data.py index 5894c3d..e9c24b1 100644 --- a/freepaths/data.py +++ b/freepaths/data.py @@ -78,17 +78,17 @@ class ScatteringData: def __init__(self): """Initialize arrays according to the number of segments""" - self.wall_diffuse = np.zeros(cf.number_of_length_segments) - self.wall_specular = np.zeros(cf.number_of_length_segments) - self.top_diffuse = np.zeros(cf.number_of_length_segments) - self.top_specular = np.zeros(cf.number_of_length_segments) - self.hole_diffuse = np.zeros(cf.number_of_length_segments) - self.hole_specular = np.zeros(cf.number_of_length_segments) - self.pillar_diffuse = np.zeros(cf.number_of_length_segments) - self.pillar_specular = np.zeros(cf.number_of_length_segments) - self.hot_side = np.zeros(cf.number_of_length_segments) - self.internal = np.zeros(cf.number_of_length_segments) - self.total = np.zeros(cf.number_of_length_segments) + self.wall_diffuse = np.zeros(cf.number_of_length_segments+1) + self.wall_specular = np.zeros(cf.number_of_length_segments+1) + self.top_diffuse = np.zeros(cf.number_of_length_segments+1) + self.top_specular = np.zeros(cf.number_of_length_segments+1) + self.hole_diffuse = np.zeros(cf.number_of_length_segments+1) + self.hole_specular = np.zeros(cf.number_of_length_segments+1) + self.pillar_diffuse = np.zeros(cf.number_of_length_segments+1) + self.pillar_specular = np.zeros(cf.number_of_length_segments+1) + self.hot_side = np.zeros(cf.number_of_length_segments+1) + self.internal = np.zeros(cf.number_of_length_segments+1) + self.total = np.zeros(cf.number_of_length_segments+1) def save_scattering_events(self, y, scattering_types): """Analyze types of scattering at the current timestep and add it to the statistics""" diff --git a/freepaths/default_config.py b/freepaths/default_config.py index 129696a..adc1cea 100644 --- a/freepaths/default_config.py +++ b/freepaths/default_config.py @@ -15,7 +15,7 @@ OUTPUT_TRAJECTORIES_OF_FIRST = 10 OUTPUT_STRUCTURE_COLOR = "#F0F0F0" NUMBER_OF_LENGTH_SEGMENTS = 10 -HOT_SIDE_ANGLE_DISTRIBUTION = "random" +HOT_SIDE_ANGLE_DISTRIBUTION = "random_up" # Animation: OUTPUT_PATH_ANIMATION = False @@ -39,10 +39,15 @@ THICKNESS = 150e-9 WIDTH = 500e-9 LENGTH = 2200e-9 +INCLUDE_RIGHT_SIDEWALL = True +INCLUDE_LEFT_SIDEWALL = True +INCLUDE_TOP_SIDEWALL = False +INCLUDE_BOTTOM_SIDEWALL = False # Hot and cold sides [m]: FREQUENCY_DETECTOR_SIZE = WIDTH COLD_SIDE_POSITION = 'top' +HOT_SIDE_POSITION = 'bottom' HOT_SIDE_X = 0 HOT_SIDE_WIDTH_X = WIDTH HOT_SIDE_Y = 0 diff --git a/freepaths/options.py b/freepaths/options.py index 2608ca1..89df77d 100644 --- a/freepaths/options.py +++ b/freepaths/options.py @@ -13,18 +13,23 @@ class Materials(enum.Enum): class Distributions(enum.Enum): """Possible distributions of angles""" - RANDOM = 1 - LAMBERT = 2 - DIRECTIONAL = 3 - UNIFORM = 4 + RANDOM_UP = 1 + RANDOM_DOWN = 2 + RANDOM_RIGHT = 3 + RANDOM_LEFT = 4 + LAMBERT = 5 + DIRECTIONAL = 6 + UNIFORM = 7 class Positions(enum.Enum): - """Possible positions of cold side""" + """Possible positions of cold and hot sides""" TOP = 1 - RIGHT = 2 - TOP_AND_RIGHT = 3 - TOP_AND_BOTTOM = 4 + BOTTOM = 2 + RIGHT = 3 + LEFT = 4 + # TOP_AND_RIGHT = 5 + # TOP_AND_BOTTOM = 6 class Polarizations(enum.Enum): diff --git a/freepaths/output_plots.py b/freepaths/output_plots.py index b1c1fff..64cc36f 100644 --- a/freepaths/output_plots.py +++ b/freepaths/output_plots.py @@ -308,6 +308,7 @@ def plot_scattering_statistics(): def plot_data(): """Create plots of various distributions""" + plot_trajectories() plot_angle_distribution() plot_free_path_distribution() plot_frequency_distribution() @@ -320,7 +321,6 @@ def plot_data(): plot_temperature_profile() plot_heat_flux_profile() plot_thermal_map() - plot_trajectories() plot_scattering_statistics() if cf.output_scattering_map: plot_scattering_map() diff --git a/freepaths/phonon.py b/freepaths/phonon.py index 5207197..bd17aca 100644 --- a/freepaths/phonon.py +++ b/freepaths/phonon.py @@ -50,13 +50,19 @@ def is_in_system(self): Depending on where we set the cold side, we check if phonon crossed that line""" small_offset = 10e-9 if cf.cold_side_position == Positions.TOP: - return cf.length > self.y + return self.y < cf.length + if cf.cold_side_position == Positions.BOTTOM: + return self.y > 0 if cf.cold_side_position == Positions.RIGHT: - return (self.y < cf.length - 1.1e-6) or (self.y > cf.length - 1.1e-6 and self.x < cf.width / 2.0 - small_offset) - if cf.cold_side_position == Positions.TOP_AND_RIGHT: - return (self.y < 1.0e-6) or (1.0e-6 < self.y < cf.length and self.x < cf.width / 2.0 - small_offset) - if cf.cold_side_position == Positions.TOP_AND_BOTTOM: - return cf.length > self.y > 0 + # return (self.y < cf.length - 1.1e-6) or (self.y > cf.length - 1.1e-6 and self.x < cf.width / 2.0 - small_offset) + return self.x < cf.width / 2.0 + if cf.cold_side_position == Positions.LEFT: + # return (self.y < cf.length - 1.1e-6) or (self.y > cf.length - 1.1e-6 and self.x < cf.width / 2.0 - small_offset) + return self.x > - cf.width / 2.0 + # if cf.cold_side_position == Positions.TOP_AND_RIGHT: + # return (self.y < 1.0e-6) or (1.0e-6 < self.y < cf.length and self.x < cf.width / 2.0 - small_offset) + # if cf.cold_side_position == Positions.TOP_AND_BOTTOM: + # return cf.length > self.y > 0 raise ValueError('Specified "cold_side" is not valid. Only TOP, RIGHT, TOP_AND_RIGH, TOP_AND_BOTTOM') def assign_polarization(self): @@ -71,9 +77,19 @@ def assign_coordinates(self): def assign_angles(self): """Depending on angle distribution, assign angles""" - if cf.hot_side_angle_distribution == Distributions.RANDOM: + if cf.hot_side_angle_distribution == Distributions.RANDOM_UP: self.theta = -pi/2 + pi*random() self.phi = asin(2*random() - 1) + if cf.hot_side_angle_distribution == Distributions.RANDOM_DOWN: + rand_sign = sign((2*random() - 1)) + self.theta = rand_sign*(pi/2 + pi/2*random()) + self.phi = asin(2*random() - 1) + if cf.hot_side_angle_distribution == Distributions.RANDOM_RIGHT: + self.theta = pi*random() + self.phi = asin(2*random() - 1) + if cf.hot_side_angle_distribution == Distributions.RANDOM_LEFT: + self.theta = - pi*random() + self.phi = asin(2*random() - 1) if cf.hot_side_angle_distribution == Distributions.DIRECTIONAL: self.theta = 0 self.phi = -pi/2 + pi*random() diff --git a/freepaths/scattering.py b/freepaths/scattering.py index 573a7b3..99fce21 100644 --- a/freepaths/scattering.py +++ b/freepaths/scattering.py @@ -7,6 +7,7 @@ from freepaths.config import cf from freepaths.move import move from freepaths.scattering_types import Scattering +from freepaths.options import Positions def specularity(angle, roughness, wavelength): @@ -25,10 +26,15 @@ def internal_scattering(ph, flight, scattering_types): def reinitialization(ph, scattering_types): """Re-thermalize phonon if it comes back to the hot side""" - _, y, _ = move(ph, cf.timestep) - - # If phonon returns to the staring line y = 0, generate it again: - if y < 0: + x, y, _ = move(ph, cf.timestep) + + # If phonon returns to the hot side, generate it again: + if ( + (cf.hot_side_position == Positions.BOTTOM and y < 0) or + (cf.hot_side_position == Positions.TOP and y > cf.length) or + (cf.hot_side_position == Positions.RIGHT and x > cf.width/2) or + (cf.hot_side_position == Positions.LEFT and x < -cf.width/2) + ): ph.assign_angles() scattering_types.hot_side = Scattering.DIFFUSE @@ -364,12 +370,12 @@ def no_new_scattering(ph): return True if (abs(z) < cf.thickness / 2 and abs(x) < cf.width / 2 and y > 0) else False -def side_wall_scattering(ph, scattering_types): - """Check if the phonon hits a side wall and output new vector""" +def scattering_on_right_sidewall(ph, scattering_types): + """Check if the phonon hits right side wall and output new vector""" x, y, z = move(ph, cf.timestep) # If phonon is beyond the side wall: - if abs(x) > cf.width/2: + if x > cf.width/2: # Calculate angle to the surface and specular scattering probability: a = acos(cos(ph.phi)*sin(abs(ph.theta))) # Angle to the surface @@ -400,6 +406,107 @@ def side_wall_scattering(ph, scattering_types): break +def scattering_on_left_sidewall(ph, scattering_types): + """Check if the phonon hits left side wall and output new vector""" + x, y, z = move(ph, cf.timestep) + + # If phonon is beyond the side wall: + if x < -cf.width/2: + + # Calculate angle to the surface and specular scattering probability: + a = acos(cos(ph.phi)*sin(abs(ph.theta))) # Angle to the surface + p = specularity(a, cf.side_wall_roughness, ph.wavelength) + + # Specular scattering: + if random() < p: + scattering_types.walls = Scattering.SPECULAR + ph.theta = - ph.theta + + # Diffuse scattering: + else: + scattering_types.walls = Scattering.DIFFUSE + attempt = 0 + while attempt < 10: + attempt += 1 + + # Random distribution: + # theta = -sign(x)*pi+sign(x)*pi*random() + # phi = asin(2*random() - 1) + + # Lambert cosine distribution: + ph.theta = -sign(x)*pi/2 + asin(2*random() - 1) + ph.phi = asin((asin(2*random() - 1))/(pi/2)) + + # Accept the angles if they do not cause new scattering: + if no_new_scattering(ph): + break + + +def scattering_on_top_sidewall(ph, scattering_types): + """Check if the phonon hits top side wall and output new vector""" + x, y, z = move(ph, cf.timestep) + + # If phonon is beyond the side wall: + if y > cf.length: + + # Calculate angle to the surface and specular scattering probability: + a = acos(cos(ph.phi)*cos(ph.theta)) + p = specularity(a, cf.hole_roughness, ph.wavelength) + + # Specular scattering: + if random() < p: + scattering_types.holes = Scattering.SPECULAR + ph.theta = sign(ph.theta)*pi - ph.theta + + # Diffuse scattering: + else: + scattering_types.holes = Scattering.DIFFUSE + attempt = 0 + while attempt < 10: + attempt += 1 + + # Lambert cosine distribution: + rand_sign = sign((2*random() - 1)) + ph.theta = rand_sign*pi/2 + rand_sign*acos(random()) + ph.phi = asin((asin(2*random() - 1))/(pi/2)) + + # Accept the angles only if they do not immediately cause new scattering: + if no_new_scattering(ph): + break + + +def scattering_on_bottom_sidewall(ph, scattering_types): + """Check if the phonon hits bottom side wall and output new vector""" + x, y, z = move(ph, cf.timestep) + + # If phonon is beyond the side wall: + if y < 0.0: + + # Calculate angle to the surface and specular scattering probability: + a = acos(cos(ph.phi)*cos(ph.theta)) + p = specularity(a, cf.hole_roughness, ph.wavelength) + + # Specular scattering: + if random() < p: + scattering_types.holes = Scattering.SPECULAR + ph.theta = sign(ph.theta)*pi - ph.theta + + # Diffuse scattering: + else: + scattering_types.holes = Scattering.DIFFUSE + attempt = 0 + while attempt < 10: + attempt += 1 + + # Lambert cosine distribution: + ph.theta = asin(2*random() - 1) + ph.phi = asin((asin(2*random() - 1))/(pi/2)) + + # Accept the angles only if they do not immediately cause new scattering: + if no_new_scattering(ph): + break + + def top_scattering(ph, scattering_types): """Check if the phonon hits the top surface and output new vector""" x, y, z = move(ph, cf.timestep) @@ -526,7 +633,14 @@ def surface_scattering(ph, scattering_types): bottom_scattering(ph, scattering_types) # Scattering on sidewalls: - side_wall_scattering(ph, scattering_types) + if cf.include_right_sidewall: + scattering_on_right_sidewall(ph, scattering_types) + if cf.include_left_sidewall: + scattering_on_left_sidewall(ph, scattering_types) + if cf.include_top_sidewall: + scattering_on_top_sidewall(ph, scattering_types) + if cf.include_bottom_sidewall: + scattering_on_bottom_sidewall(ph, scattering_types) # Scattering on parabolic walls: if cf.include_top_parabola: @@ -565,6 +679,8 @@ def surface_scattering(ph, scattering_types): Lx = cf.rectangular_hole_side_x * (cf.hole_coordinates[i, 2] + 1) Ly = cf.rectangular_hole_side_y * (cf.hole_coordinates[i, 2] + 1) scattering_on_triangle_down_holes(ph, x0, y0, Lx, Ly, scattering_types, x, y, z) + else: + pass # If there was any scattering, then no need to check other holes: if scattering_types.holes is not None: diff --git a/schemes/angle_distributions.svg b/schemes/angle_distributions.svg new file mode 100644 index 0000000..714ac85 --- /dev/null +++ b/schemes/angle_distributions.svg @@ -0,0 +1,425 @@ + + + +random_up(default)random_downrandom_rightrandom_leftlambertdirectionaluniformHOT_SIDE_ANGLE_DISTRIBUTION diff --git a/schemes/gradients.svg b/schemes/gradients.svg new file mode 100644 index 0000000..0a879b0 --- /dev/null +++ b/schemes/gradients.svg @@ -0,0 +1,2096 @@ + + + +Horizontal gradientVertical gradientHot sideCold sideTop side wallBottom side wallHot sideCold sideRight side wallLeft side wall diff --git a/schemes/scheme.svg b/schemes/scheme.svg new file mode 100644 index 0000000..f72132c --- /dev/null +++ b/schemes/scheme.svg @@ -0,0 +1,503 @@ + + + +Hot side (bottom)Cold side (top)Right side wall (on)Left side wall (on)LengthWidthBottom side wall (off)Top side wall (off)Lengh segments (4)StepExit angleInitial angleTimeTime stepTime segment Total time = time step × number of time steps