density

title.rb

import pandas as pd
import numpy as np

r_res = 0.01 # 距離方向の分解能（m）

# DataFrameの作成
df = pd.DataFrame(columns=['r', 'theta', 'phi'])

# r, theta, phiの値を生成するリスト内包表記
data = [{'r': r, 'theta': theta, 'phi': phi} 
        for r in [num * r_res for num in range(0, 1000)]
        for theta in range(-30, 32, 2)
        for phi in range(-30, 32, 2)]

sin_th = np.sin(np.radians(df['theta'].astype(float)))
cos_th = np.cos(np.radians(df['theta'].astype(float)))
sin_ph = np.sin(np.radians(df['phi'].astype(float)))
cos_ph = np.cos(np.radians(df['phi'].astype(float)))

# DataFrameにデータを追加
df = df.append(data, ignore_index=True)

df['x'] = df['r'] * np.sin(np.radians(df['theta'].astype(float))) * np.cos(np.radians(df['phi'].astype(float)))
df['y'] = df['r'] * np.sin(np.radians(df['phi'].astype(float))) * np.cos(np.radians(df['theta'].astype(float)))
df['z'] = df['r'] * ((1-(np.sin(np.radians(df['theta'].astype(float)))**2)*(np.cos(np.radians(df['phi'].astype(float)))**2)-(np.sin(np.radians(df['phi'].astype(float)))**2)*(np.cos(np.radians(df['theta'].astype(float)))**2))**0.5)
df['r_check'] = (df['x']**2 + df['y']**2 + df['z']**2)**0.5

# データフレームを表示
df.tail()




import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 半球のパラメータ
r = 10  # 球の半径
phi, theta = np.mgrid[0:np.pi/2:100j, 0:2*np.pi:100j]  # phiとthetaの範囲を設定

# 極座標から直交座標に変換
x = r * np.sin(phi) * np.cos(theta)
y = r * np.sin(phi) * np.sin(theta)
z = r * np.cos(phi)

# プロット
fig = plt.figure(figsize=(16.0, 12.0))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(df['x'], df['y'], df['z'], alpha=0.01)
# 半球をプロット
ax.plot_surface(x, y, z, alpha=0.2)

# 軸の範囲とラベル設定
ax.set_xlim([-10, 10])
ax.set_ylim([-10, 10])
ax.set_zlim([0, 20])

# 軸の縮尺を等しくする
ax.set_box_aspect([1,1,1])

#ax.tick_params(labelbottom=False, labelleft=False, labelright=False, labeltop=False)

# グラフを表示
plt.show()


# 検証として1 m^3の立方体の密度を考える
立方体の座標範囲はx:-0.5～0.5, y:-0.5～0.5, z:n～n+1　（nは2から6まで1ずつ動かしてみる）
立方体内に含まれる点数をカウントするだけで密度は求まる
## まずは補正せずに計算


density = []
for n in range(1,8):
    df_d = df[(df["x"]<0.5) & (df["x"]>-0.5) & (df["y"]<0.5) & (df["y"]>-0.5) & (df["z"]<n+1) & (df["z"]>n)]
    density.append(len(df_d))
    print("n = {} のときの立方体内の点密度{}".format(n,len(df_d)))
df_d.head()


## カウントする際に半径の2乗を重みにかけて補正計算



adjusted_density = []
for n in range(1,5):
    df_d = df[(df["x"]<0.5) & (df["x"]>-0.5)&(df["y"]<0.5) & (df["y"]>-0.5)&(df["z"]<n+1) & (df["z"]>n)]
    df_d["r^2"] = df_d["r"]**2
    pts_sum = df_d["r^2"].sum()
    adjusted_density.append(pts_sum)
    print("n = {} のときの立方体内の補正点密度{}".format(n,pts_sum))