# -*- coding: utf-8 -*- """ Created on Wed Dec 4 15:26:29 2024 @author: Atsushi Kobayashi This program was developed with the assistance of ChatGPT 4o. """ import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Point, LineString from collections import defaultdict # 六方格子の格子点を定義 a = 2 # 格子定数 lattice_points = [] for i in range(-5, 6): # x方向 for j in range(-5, 6): # y方向 x = a * (i + 0.5 * j) y = a * (np.sqrt(3) / 2) * j lattice_points.append(Point(x, y)) # 原点を定義 origin = Point(0, 0) # 原点から指定された距離以内の格子点を取得 def points_within_distance(lattice_points, origin, max_distance): return [p for p in lattice_points if origin.distance(p) <= max_distance] # 2a以内の格子点を取得 points_within_2a = points_within_distance(lattice_points, origin, 2 * a) # 垂直二等分線を計算 def calculate_perpendicular_bisectors(points, origin): bisectors = [] for p1 in points: if np.isclose(origin.distance(p1), 0): # 原点を除外 continue midpoint = ((origin.x + p1.x) / 2, (origin.y + p1.y) / 2) if np.isclose(p1.x, origin.x): # 垂直線の場合 line = LineString([(midpoint[0] - 10, midpoint[1]), (midpoint[0] + 10, midpoint[1])]) elif np.isclose(p1.y, origin.y): # 水平線の場合 line = LineString([(midpoint[0], midpoint[1] - 10), (midpoint[0], midpoint[1] + 10)]) else: slope = -1 / ((p1.y - origin.y) / (p1.x - origin.x)) intercept = midpoint[1] - slope * midpoint[0] line = LineString([(x, slope * x + intercept) for x in np.linspace(-10, 10, 500)]) bisectors.append(line) return bisectors # 垂直二等分線の交点を計算 def find_intersections(bisectors): intersections = [] for i, line1 in enumerate(bisectors): for line2 in bisectors[i + 1:]: if line1.intersects(line2): intersection = line1.intersection(line2) if intersection.geom_type == 'Point': intersections.append(intersection) return intersections # 重複削除と分類を同時に行う def classify_and_remove_duplicates(intersections, center): unique_points = {} for inter in intersections: rounded_point = (round(inter.x, 4), round(inter.y, 4)) # 座標を丸めて重複排除 if rounded_point not in unique_points: unique_points[rounded_point] = center.distance(inter) # 距離を計算して保存 # 距離ごとに分類 distance_dict = defaultdict(list) for point, distance in unique_points.items(): distance_dict[round(distance, 4)].append(Point(*point)) return distance_dict # 第一ブリルアンゾーンを描画 def plot_first_brillouin_zone(c1_points, origin, color="yellow", edgecolor="black"): c1_sorted = sorted(c1_points, key=lambda p: np.arctan2(p.y - origin.y, p.x - origin.x)) x_coords = [p.x for p in c1_sorted] + [c1_sorted[0].x] y_coords = [p.y for p in c1_sorted] + [c1_sorted[0].y] plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5, label="First Brillouin Zone") # 第二ブリルアンゾーンを描画 def plot_second_brillouin_zone(c1_points, c2_points, origin, color="green", edgecolor="black"): second_zone_triangles = [] c1_sorted = sorted(c1_points, key=lambda p: np.arctan2(p.y - origin.y, p.x - origin.x)) c1_pairs = [(c1_sorted[i], c1_sorted[(i + 1) % len(c1_sorted)]) for i in range(len(c1_sorted))] # ラベルがすでに作成されているか確認 label_created = False for c1_1, c1_2 in c1_pairs: line_midpoint = Point((c1_1.x + c1_2.x) / 2, (c1_1.y + c1_2.y) / 2) closest_c2 = min(c2_points, key=lambda p: p.distance(line_midpoint)) second_zone_triangles.append([c1_1, c1_2, closest_c2]) for triangle in second_zone_triangles: x_coords = [p.x for p in triangle] + [triangle[0].x] y_coords = [p.y for p in triangle] + [triangle[0].y] if not label_created: # 初回のみ凡例を追加 plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5, label="Second Brillouin Zone") label_created = True else: plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5) # 第三ブリルアンゾーンを描画 def plot_third_brillouin_zone(c1_points, c2_points, origin, color="cyan", edgecolor="black"): third_zone_triangles = [] label_created = False # ラベルの作成確認用 for c1 in c1_points: closest_c2_points = sorted(c2_points, key=lambda p: p.distance(c1))[:2] if len(closest_c2_points) == 2: third_zone_triangles.append([c1, closest_c2_points[0], closest_c2_points[1]]) for triangle in third_zone_triangles: x_coords = [p.x for p in triangle] + [triangle[0].x] y_coords = [p.y for p in triangle] + [triangle[0].y] if not label_created: # 初回のみ凡例を追加 plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5, label="Third Brillouin Zone") label_created = True else: plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5) # 第四ブリルアンゾーンを描画 def plot_fourth_brillouin_zone(c2_points, c3_points, origin, color="magenta", edgecolor="black"): fourth_zone_triangles = [] c2_sorted = sorted(c2_points, key=lambda p: np.arctan2(p.y - origin.y, p.x - origin.x)) c2_pairs = [(c2_sorted[i], c2_sorted[(i + 1) % len(c2_sorted)]) for i in range(len(c2_sorted))] label_created = False # ラベルの作成確認用 for c2_1, c2_2 in c2_pairs: line_midpoint = Point((c2_1.x + c2_2.x) / 2, (c2_1.y + c2_2.y) / 2) closest_c3 = min(c3_points, key=lambda p: p.distance(line_midpoint)) fourth_zone_triangles.append([c2_1, c2_2, closest_c3]) for triangle in fourth_zone_triangles: x_coords = [p.x for p in triangle] + [triangle[0].x] y_coords = [p.y for p in triangle] + [triangle[0].y] if not label_created: # 初回のみ凡例を追加 plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5, label="Fourth Brillouin Zone") label_created = True else: plt.fill(x_coords, y_coords, color=color, edgecolor=edgecolor, alpha=0.5) # 垂直二等分線を計算 bisectors = calculate_perpendicular_bisectors(points_within_2a, origin) # 交点を分類しつつ重複削除 intersections = find_intersections(bisectors) classified_intersections = classify_and_remove_duplicates(intersections, origin) # 距離順に分類 sorted_distances = sorted(classified_intersections.keys()) intersection_groups = {f"c{i+1}": classified_intersections[dist] for i, dist in enumerate(sorted_distances)} # 垂直二等分線の交点c1, c2, c3を取得 c1_points = intersection_groups["c1"] c2_points = intersection_groups["c2"] c3_points = intersection_groups["c3"] # プロット開始 plt.figure(figsize=(8, 8)) # 格子点をプロット plt.scatter([p.x for p in lattice_points], [p.y for p in lattice_points], color="black", s=50, marker='o', label="Lattice Points") # 第一から第四ブリルアンゾーンをプロット plot_first_brillouin_zone(c1_points, origin) plot_second_brillouin_zone(c1_points, c2_points, origin) plot_third_brillouin_zone(c1_points, c2_points, origin) plot_fourth_brillouin_zone(c2_points, c3_points, origin) # 垂直二等分線を最後に描画 for line in bisectors: x, y = line.xy plt.plot(x, y, linestyle="--", color="red", alpha=0.6) # プロットの仕上げ plt.title("Brillouin zones for hexagonal reciprocal lattice") plt.xlabel("kx") plt.ylabel("ky") plt.xlim(-5, 5) plt.ylim(-5, 5) plt.legend(loc="upper right", fontsize="small", ncol=1) plt.gca().set_aspect("equal", adjustable="box") plt.grid(False) plt.show()