#!/usr/bin/env python # coding: utf-8 # In[33]: #縦孔内輝度値反射シュミレーション # 均等拡散反射 # 円筒形空洞 # 影あり # カメラの捉える輝度 #ライブラリのインポート======================== import math from tkinter import N from urllib.parse import DefragResult import numpy as np import matplotlib.pyplot as plt import time import datetime from mpl_toolkits.axes_grid1 import make_axes_locatable #========================================== # In[34]: #パラメ-タの設定************************************************************************ #縦孔の形状 r = 1.0 #縦孔の半径[m] d1 = 1.0 #縦孔の壁の高さ（厚さ）[m] d2 = 1.0 #縦孔の空洞の高さ[m] R = 1.0 #壁の反射率[] cr = 2.0 #空洞壁の位置[]（rの何倍の位置かどうか） #太陽 J = 1.0 #太陽の入射エネルギー][J/s/m^2] theta = 45.0 #太陽高度[°] #その他 PI = math.pi #円周率[] #セルの幅 dx = 0.01 #x(太陽方位に平行)方向の分解能[m] dy = 0.01 #y(太陽方位に垂直)方向の分解能[m] dz = 0.01 #z(鉛直)方向の分解能[m] dr = 0.01 #r方向の分解能[m] dphi = 1.0 #φ方向の分解能[°] #計算範囲 #縦孔底 FL_x_min = - 2.0 * r #x軸の最小値[m] FL_x_max = 2.0 * r #x軸の最大値[m] FL_y_min = - 2.0 * r #y軸の最小値[m] FL_y_max = 2.0 * r #y軸の最大値[m] FL_r_min = 0.0 * r #r方向の最小値[m] FL_r_max = 2.0 * r #r方向の最大値[m] FL_phi_min = -90.0 #φ方向の最小値[°] FL_phi_max = 270.0 #φ方向の最大値[°] #縦孔壁 WA_z_min = d2 #z軸の最小値[m] WA_z_max = d1 + d2 #z軸の最大値[m] WA_phi_min = -90.0 #φ方向の最小値[°] WA_phi_max = 270.0 #φ方向の最大値[°] #空洞壁 CW_z_min = 0 #z軸の最小値[m] CW_z_max = d2 #z軸の最大値[m] CW_phi_min = -90.0 #φ方向の最小値[°] CW_phi_max = 270.0 #φ方向の最大値[°] #空洞天井 CC_x_min = - 2.0 * r #x軸の最小値[m] CC_x_max = 2.0 * r #x軸の最大値[m] CC_y_min = - 2.0 * r #y軸の最小値[m] CC_y_max = 2.0 * r #y軸の最大値[m] CC_r_min = 1.0 * r #z軸の最小値[m] CC_r_max = 2.0 * r #z軸の最大値[m] CC_phi_min = -90.0 #φ方向の最小値[°] CC_phi_max = 270.0 #φ方向の最大値[°] #*********************************************************************************************** # In[35]: #関数の定義==================================== #rootを求める関数 def SQRT(x): return np.sqrt(x) #角度からsinを求める関数 def SIN(a): return np.sin(math.radians(a)) #角度からcosを求める関数 def COS(a): return np.cos(math.radians(a)) #2つベクトルよりcosθを求める関数 def DCOS(a, b): over = (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]) under = (SQRT(((a[0]) ** 2.0) + ((a[1]) ** 2.0) + ((a[2]) ** 2.0))) * (SQRT(((b[0]) ** 2.0) + ((b[1]) ** 2.0) + ((b[2]) ** 2.0))) if under == 0.0: return(0.0) else: return (over / under) #2点より単位ベクトルをもとめる関数 def VEC2(a,b): #bベクトルからaベクトル I = a[0] - b[0] J = a[1] - b[1] K = a[2] - b[2] under = SQRT((I ** 2.0) + (J ** 2.0) + (K ** 2.0)) if under == 0.0: return((0.0, 0.0, 0.0)) else: return((I/under), (J/under), (K/under)) #2点よりベクトルをもとめる関数 def VEC(a,b): #bベクトルからaベクトル I = a[0] - b[0] J = a[1] - b[1] K = a[2] - b[2] return(I, J, K) #2点間の距離の2乗を求める関数 def LL(a,b): return (((a[0]- b[0]) ** 2.0) + ((a[1]- b[1]) ** 2.0) + ((a[2]- b[2]) ** 2.0)) #グリッドの設定 def grid (d,max,min): l = int((max - min) / d) r = np.zeros(l+1) for i in range(l+1): r[i] = min + (d * float(i)) return r #反射パターン def Ref(pattern,Inc,Emi,Nor): if pattern == "Lambert": #均等拡散反射（ランバート反射） return DCOS(Nor,Emi) * R / PI elif pattern == "end": return DCOS(Emi,Inc) * R / PI else: print("stop") #2次関数の解の存否の関数 def dis(A,B,C):#A,B,Cの値 D = (B ** 2) - (4.0 * A * C) if(D < 0.0): return 0 #解がないとき0を返す else: return 1 #解があるとき1を返す #2次関数の解の公式の関数 def kai(A,B,C,s):#A,B,Cの値，sは−1or1 return ((- 1.0 * B )+ (s * SQRT(B ** 2 - (4.0 * A * C))))/(2.0 * A ) #円筒（縦孔壁面）とのあたり判定 def HitCyl(s,d,r,hmax,hmin): A = (d[0] ** 2) + (d[1] ** 2) B = 2.0 * ((s[0] * d[0]) + (s[1] * d[1])) C = (s[0] ** 2) + (s[1] ** 2) - (r ** 2) if dis(A,B,C) < 0.0: return 0 else: t1 = kai(A, B, C, 1.0) t2 = kai(A, B, C, -1.0 ) if t1 <= 0.0 and t2 <= 0.0 : return 0 elif t1 > 0 and t2 <=0: z = s[2] + (t1 * d[2]) if z >= hmin and z <= hmax: return 1 else: return 0 elif t1 <= 0 and t2 > 0: z = s[2] + (t2 * d[2]) if z >= hmin and z <= hmax: return 1 else: return 0 elif t1 >= t2 : z = s[2] + (t1 * d[2]) if z >= hmin and z <= hmax: return 1 else: return 0 else: z = s[2] + (t2 * d[2]) if z >= hmin and z <= hmax: return 1 else: return 0 #円盤（縦孔開口部）とのあたり判定 def HitDisk(s,d,r,h): if d[2]== 0.0: if s[2] == h: return 1 else: # print("yoko") return 0 elif (((s[0]+ ((h - s[2]) * d[0] / d[2]))**2)+(((s[1]+((h - s[2]) * d[1] / d[2]))**2))) <= (r ** 2): return 1 else: # print("soto") return 0 #============================================ # In[36]: #縦孔底面のグリッド設定========================= SUR_FL = "FL" a_FL = grid(dx, FL_x_max, FL_x_min) b_FL = grid(dy, FL_x_max, FL_x_min) #=========================================== #縦孔壁面のグリッド設定========================= SUR_WA = "WA" a_WA = grid(dz, WA_z_max, WA_z_min) b_WA = grid(dphi, WA_phi_max, WA_phi_min) #=========================================== #空洞壁面のグリッド設定========================= SUR_CW = "CW" a_CW = grid(dz, CW_z_max, CW_z_min) b_CW = grid(dphi, CW_phi_max, CW_phi_min) #=========================================== # #空洞天井のグリッド設定========================= # SUR_CC = "CC" # a_CC = grid(dr, CC_r_max, CC_r_min) # b_CC = grid(dphi, CC_phi_max, CC_phi_min) # #=========================================== #空洞天井のグリッド設定========================= SUR_CC = "CC" a_CC = grid(dx, CC_x_max, CC_x_min) b_CC = grid(dy, CC_y_max, CC_y_min) #=========================================== #計算ループの設定=============================== #縦孔底 M_FL = len(a_FL) #反射面のセル数 N_FL = len(b_FL) #反射面のセル数 #縦孔壁 M_WA = len(a_WA) #反射面のセル数 N_WA = len(b_WA) #反射面のセル数 #空洞壁 M_CW = len(a_CW) #反射面のセル数 N_CW = len(b_CW) #反射面のセル数 #空洞天井 M_CC = len(a_CC) #反射面のセル数 N_CC = len(b_CC) #反射面のセル数 #============================================ # In[37]: #ベクトルの設定***************************************** #太陽光ベクトル(i,j,k) I = np.array([-COS(theta), 0.0, -SIN(theta)]) SUR_S ="S" #縦孔底面の法線ベクトルn================ n_FL = np.array([0.0, 0.0, 1.0]) #縦孔底面の位置ベクトルF F_FL = np.zeros((M_FL, N_FL, 3)) for i in range(M_FL): for j in range(N_FL): F_FL[i,j,0] = a_FL[i] F_FL[i,j,1] = b_FL[j] F_FL[i,j,2] = 0.0 #縦孔底面の微小面積dA dA_FL = dx * dy #===================================== #縦孔壁面の法線ベクトルn================ n_WA = np.zeros((M_WA, N_WA, 3)) for i in range(M_WA): for j in range(N_WA): n_WA[i,j,0] = - COS(b_WA[j]) n_WA[i,j,1] = - SIN(b_WA[j]) n_WA[i,j,2] = 0.0 #縦孔壁面の位置ベクトルF F_WA = np.zeros((M_WA, N_WA, 3)) for i in range(M_WA): for j in range(N_WA): F_WA[i,j,0] = r * COS(b_WA[j]) F_WA[i,j,1] = r * SIN(b_WA[j]) F_WA[i,j,2] = a_WA[i] #縦孔底面の微小面積dA dA_WA = dz * r * dphi #===================================== #空洞壁面の法線ベクトルn================ n_CW = np.zeros((M_CW, N_CW, 3)) for i in range(M_CW): for j in range(N_CW): n_CW[i,j,0] = - COS(b_CW[j]) n_CW[i,j,1] = - SIN(b_CW[j]) n_CW[i,j,2] = 0.0 #空洞壁面の位置ベクトルF F_CW = np.zeros((M_CW, N_CW, 3)) for i in range(M_CW): for j in range(N_CW): F_CW[i,j,0] = r * cr * COS(b_CW[j]) F_CW[i,j,1] = r * cr * SIN(b_CW[j]) F_CW[i,j,2] = a_CW[i] #縦孔底面の微小面積dA dA_CW = dz * r * cr * dphi #===================================== #空洞天井の法線ベクトルn================ # n_CC = np.array([0.0, 0.0, -1.0]) # #空洞天井の位置ベクトルF # F_CC = np.zeros((M_CC, N_CC, 3)) # for i in range(M_CC): # for j in range(N_CC): # F_CC[i,j,0] = a_CC[i] * COS(b_CW[j]) # F_CC[i,j,1] = a_CC[i] * SIN(b_CW[j]) # F_CC[i,j,2] = d2 # #縦孔底面の微小面積dA # dA_CC = np.zeros((M_CC)) # for i in range(M_CC): # dA_CC[i] = dr * a_CC[i] * dphi #===================================== #空洞天井の法線ベクトルn================ n_CC = np.array([0.0, 0.0, -1.0]) #空洞天井の位置ベクトルF F_CC = np.zeros((M_CC, N_CC, 3)) for i in range(M_CC): for j in range(N_CC): F_CC[i,j,0] = a_CC[i] F_CC[i,j,1] = b_CC[j] F_CC[i,j,2] = d2 #縦孔底面の微小面積dA dA_CC = dx * dy #==================================== #************************************************* # In[38]: #計算開始の時刻=========== start_time = time.time() print("--start--") #======================= # In[39]: #解析1-1（太陽直達光に照らされる面（反射面の計算））================================================= #底_FL=========================================================================================== #放射照度E_StoFL（放射照度） E_StoFL = np.zeros((M_FL, N_FL)) #観測面のセルの照度の初期化 MM = M_FL NN = N_FL F = F_FL n = n_FL for ii in range(MM): for jj in range(NN): if DCOS(n,-1*I) < 0.0: # 壁面の法線方向と太陽光とのなす角が90°以上のとき E_StoFL[ii,jj] = 0.0 elif HitCyl(F[ii,jj], (-1.0 * I), r, (d1 + d2), d2) == 1: # 影のとき E_StoFL[i,j] = 0.0 elif HitDisk(F[ii,jj],(-1.0 * I), r, (d1 + d2)) == 1: #縦孔開口部からの光の時 if HitDisk(F[ii,jj],(-1.0 * I), r, d2) == 1: E_StoFL[ii,jj] = J * DCOS(n,-1*I) else: E_StoFL[ii,jj] = 0.0 else: E_StoFL[ii,jj] = 0.0 #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_S, SUR_FL, theta ,d1, d2, cr, date), E_StoFL, delimiter=',') #==================================================================================================== #==================================================================================================== # In[41]: #解析1-2（太陽直達光に照らされる面（反射面の計算））================================================= #縦孔壁_WA================================================================================== #放射照度E_StoWA（放射照度） E_StoWA = np.zeros((M_WA, N_WA)) #観測面のセルの照度の初期化 M = M_WA N = N_WA F = F_WA n = n_WA for i in range(M): for j in range(N): if HitCyl(F[i,j], (-1.0 * I), r, (d1 + d2), d2) == 1: # 影のとき E_StoWA[i,j] = 0.0 elif DCOS(n[i,j],-1*I) < 0.0: #反射面の裏方向の値は計算しない E_StoWA[i,j] = 0.0 else: E_StoWA[i,j] = J * DCOS(n[i,j],-1*I) #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_S, SUR_WA, theta ,d1, d2, cr, date), E_StoWA, delimiter=',') #==================================================================================================== #==================================================================================================== # In[42]: # In[43]: #解析1-3（太陽直達光に照らされる面（反射面の計算））================================================= #空洞壁_CW================================================================================== #放射照度E_StoCW（放射照度） E_StoCW = np.zeros((M_CW, N_CW)) #観測面のセルの照度の初期化 M = M_CW N = N_CW F = F_CW n = n_CW for i in range(M): for j in range(N): if DCOS(n[i,j],-1*I) < 0.0: # 壁面の法線方向と太陽光とのなす角が90°以上のとき E_StoCW[i,j] = 0.0 elif HitCyl(F[i,j], (-1.0 * I), r, (d1 + d2), d2) == 1: # 影のとき E_StoCW[i,j] = 0.0 elif HitDisk(F[i,j],(-1.0 * I), r, (d1 + d2)) == 1: #縦孔開口部からの光の時 if HitDisk(F[i,j],(-1.0 * I), r, d2) == 1: E_StoCW[i,j] = J * DCOS(n[i,j],-1*I) else: E_StoCW[i,j] = 0.0 else: E_StoCW[i,j] = 0.0 #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_S, SUR_CW, theta ,d1, d2, cr, date), E_StoCW, delimiter=',') #==================================================================================================== #==================================================================================================== # In[44]: # In[45]: #解析2-1（縦孔底に照らされる面）================================================= #縦孔底_FL→縦孔壁_WA========================================================= E_FLtoWA = np.zeros((M_WA, N_WA)) #観測面のセルの照度の初期化 MM = M_WA NN = N_WA M = M_FL N = N_FL FF = F_WA s = n_WA F = F_FL n = n_FL dS = dA_WA dA = dA_FL #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n,k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitDisk(F[i,j],-k[i,j], r, (d2-dz)) == 0: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n) * (E_StoFL[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_FLtoWA[ii, jj] = E_dA_all print("1/10_FLtoWA:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("1/10_FLtoWA:100.00 %") #==================================================================================================== #==================================================================================================== # In[46]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_FL, SUR_WA, theta ,d1, d2, cr, date), E_FLtoWA, delimiter=',') #==================================================================================================== # In[47]: #解析2-2（縦孔底に照らされる面）================================================= #縦孔底_FL→空洞壁_CW========================================================= E_FLtoCW = np.zeros((M_CW, N_CW)) #観測面のセルの照度の初期化 MM = M_CW NN = N_CW M = M_FL N = N_FL FF = F_CW s = n_CW F = F_FL n = n_FL dS = dA_CW dA = dA_FL #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n,k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n) * (E_StoFL[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_FLtoCW[ii, jj] = E_dA_all print("2/10_FLtoCW:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("2/10_FLtoCW:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_FL, SUR_CW, theta ,d1, d2, cr, date), E_FLtoCW, delimiter=',') #==================================================================================================== # In[ ]: #解析2-3（縦孔底に照らされる面）================================================= #縦孔底_FL→空洞天井_CC========================================================= E_FLtoCC = np.zeros((M_CC, N_CC)) #観測面のセルの照度の初期化 MM = M_CC NN = N_CC M = M_FL N = N_FL FF = F_CC s = n_CC F = F_FL n = n_FL dS = dA_CC dA = dA_FL #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n,k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s,-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n) * (E_StoFL[i,j] * dA * (DCOS(-k[i,j], s))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_FLtoCC[ii, jj] = E_dA_all print("3/10_FLtoCC:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("3/10_FLtoCC:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_FL, SUR_CC, theta ,d1, d2, cr, date), E_FLtoCC, delimiter=',') #==================================================================================================== In[ ]: #解析3-1（縦孔壁に照らされる面）================================================= #縦孔壁_WA→縦孔底_FL========================================================= E_WAtoFL = np.zeros((M_FL, N_FL)) #観測面のセルの照度の初期化 MM = M_FL NN = N_FL M = M_WA N = N_WA FF = F_FL s = n_FL F = F_WA n = n_WA dS = dA_FL dA = dA_WA #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s,-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoWA[i,j] * dA * (DCOS(-k[i,j], s))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_WAtoFL[ii, jj] = E_dA_all print("4/10_WAtoFL:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("4/10_WAtoFL:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_WA, SUR_FL, theta ,d1, d2, cr, date), E_WAtoFL, delimiter=',') #==================================================================================================== # In[ ]: #解析3-2（縦孔壁に照らされる面）================================================= #縦孔壁_WA→縦孔底_WA========================================================= E_WAtoWA = np.zeros((M_WA, N_WA)) #観測面のセルの照度の初期化 MM = M_WA NN = N_WA M = M_WA N = N_WA FF = F_WA s = n_WA F = F_WA n = n_WA dS = dA_WA dA = dA_WA #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 # elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 # E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoWA[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_WAtoWA[ii, jj] = E_dA_all print("5/10_WAtoWA:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("5/10_WAtoWA:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_WA, SUR_WA, theta ,d1, d2, cr, date), E_WAtoWA, delimiter=',') #==================================================================================================== # In[ ]: #解析3-3（縦孔壁に照らされる面）================================================= #縦孔壁_WA→空洞壁_CW========================================================= E_WAtoCW = np.zeros((M_CW, N_CW)) #観測面のセルの照度の初期化 MM = M_CW NN = N_CW M = M_WA N = N_WA FF = F_CW s = n_CW F = F_WA n = n_WA dS = dA_CW dA = dA_WA #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoWA[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_WAtoCW[ii, jj] = E_dA_all print("6/10_WAtoCW:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("6/10_WAtoCW:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: # #結果出力================================================== # #日付の取得 # date = datetime.date.today() # #csv出力 # np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_WA, SUR_CW, theta ,d1, d2, cr, date), E_WAtoCW, delimiter=',') # #==================================================================================================== # In[ ]: #解析4-1（縦孔壁に照らされる面）================================================= #空洞壁_CW→縦孔底_FL========================================================= E_CWtoFL = np.zeros((M_FL, N_FL)) #観測面のセルの照度の初期化 MM = M_FL NN = N_FL M = M_CW N = N_CW FF = F_FL s = n_FL F = F_CW n = n_CW dS = dA_FL dA = dA_CW #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s,-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoCW[i,j] * dA * (DCOS(-k[i,j], s))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_CWtoFL[ii, jj] = E_dA_all print("7/10_CWtoFL:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("7/10_CWtoFL:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_CW, SUR_FL, theta ,d1, d2, cr, date), E_CWtoFL, delimiter=',') #==================================================================================================== # In[ ]: #解析4-2（縦孔壁に照らされる面）================================================= #空洞壁_CW→縦孔壁_WA========================================================= E_CWtoWA = np.zeros((M_WA, N_WA)) #観測面のセルの照度の初期化 MM = M_WA NN = N_WA M = M_CW N = N_CW FF = F_WA s = n_WA F = F_CW n = n_CW dS = dA_WA dA = dA_CW #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitDisk(F[i,j],-k[i,j], r, (d2-dz)) == 0: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoCW[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_CWtoWA[ii, jj] = E_dA_all print("8/10_CWtoWA:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("8/10_CWtoWA:100.00 %") #==================================================================================================== #==================================================================================================== In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_CW, SUR_WA, theta ,d1, d2, cr, date), E_CWtoWA, delimiter=',') #==================================================================================================== # In[ ]: #解析4-3（縦孔壁に照らされる面）================================================= #空洞壁_CW→空洞壁_CW========================================================= E_CWtoCW = np.zeros((M_CW, N_CW)) #観測面のセルの照度の初期化 MM = M_CW NN = N_CW M = M_CW N = N_CW FF = F_CW s = n_CW F = F_CW n = n_CW dS = dA_CW dA = dA_CW #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s[ii,jj],-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoCW[i,j] * dA * (DCOS(-k[i,j], s[ii,jj]))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_CWtoCW[ii, jj] = E_dA_all print("9/10_CWtoCW:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("9/10_CWtoCW:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_CW, SUR_CW, theta ,d1, d2, cr, date), E_CWtoCW, delimiter=',') #==================================================================================================== # In[ ]: #解析4-4（縦孔壁に照らされる面）================================================= #空洞壁_CW→空洞天井_CC========================================================= E_CWtoCC = np.zeros((M_CC, N_CC)) #観測面のセルの照度の初期化 MM = M_CC NN = N_CC M = M_CW N = N_CW FF = F_CC s = n_CC F = F_CW n = n_CW dS = dA_CC dA = dA_CW #反射面の微小面積からの観測面に向かうベクトルk k = np.zeros((M, N, 3)) #壁面の微小面積からの観測点までの距離L2 L2 = np.zeros((M, N)) #各微小壁からの反射による観測点の微小面積が受けるエネルギーE_o（放射照度）[J/s/m^2] E_dA = np.zeros((M, N)) for ii in range(MM): for jj in range(NN): k = k * 0.0 L2 = L2 * 0.0 E_dA = E_dA * 0.0 for i in range(M): for j in range(N): k[i,j] = VEC2(FF[ii,jj], F[i,j]) # k[i,j] = VEC2(観測点の座標，反射面の座標) L2[i,j] = LL(F[i,j],FF[ii,jj]) if DCOS(n[i,j],k[i,j]) < 0.0: #反射面の裏方向の値は計算しない E_dA[i,j] = 0.0 if DCOS(s,-k[i,j]) < 0.0: #受光面の裏方向の値は計算しない E_dA[i,j] = 0.0 elif L2[i,j]== 0.0: #反射面が重なる場合の値は計算しない E_dA[i,j] = 0.0 elif HitCyl(F[i,j], k[i,j], r, (d1 + d2), d2) == 1: #円柱（縦孔壁）との交差判定 E_dA[i,j] = 0.0 else: E_dA[i,j] = Ref("Lambert", I, k[i,j], n[i,j]) * (E_StoCW[i,j] * dA * (DCOS(-k[i,j], s))) / (L2[i,j]) if E_dA[i,j] < 0.0: E_dA[i,j] = 0.0 #壁からの反射による観測点の微小面積が受けるエネルギーE_f_all（放射照度）[J/s/m^2] E_dA_all = (np.sum(E_dA)) E_CWtoCC[ii, jj] = E_dA_all print("10/10_CWtoCC:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) print("10/10_CWtoCC:100.00 %") #==================================================================================================== #==================================================================================================== # In[ ]: #結果出力================================================== #日付の取得 date = datetime.date.today() #csv出力 np.savetxt("PPP_{}to{}_t{}h{}h{}_cr{}_{}.csv".format(SUR_CW, SUR_CC, theta ,d1, d2, cr, date), E_CWtoCC, delimiter=',') #==================================================================================================== #計算終了時刻 end_time = time.time() - start_time end_time = end_time / 360.0 print("---") print("計算時間[h] : {}".format(end_time)) print("---") print("--end--")