#!/usr/bin/env python # coding: utf-8 # In[29]: #縦孔内輝度値反射シュミレーション # 均等拡散反射 # 円筒形空洞 # 影あり # カメラの捉える輝度 #　　計算 #ライブラリのインポート======================== 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[30]: #パラメ-タの設定************************************************************************ #縦孔の形状 r = 50.0 #縦孔の半径[m] d1 = 50.0 #縦孔の壁の高さ（厚さ）[m] d2 = 50.0 #縦孔の空洞の高さ[m] R = 0.1 #壁の反射率[] cr = 2.0 #空洞壁の位置[]（rの何倍の位置かどうか） #太陽 J = 1.0 #太陽の入射エネルギー][J/s/m^2] theta = 45.0 #太陽高度[°] #その他 PI = math.pi #円周率[] #セルの幅 dx = 0.5 #x(太陽方位に平行)方向の分解能[m] dy = 0.5 #y(太陽方位に垂直)方向の分解能[m] dz = 0.5 #z(鉛直)方向の分解能[m] dr = 0.5 #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[31]: #カメラのパラメ-タの設定************************************************************************ i_CAM = np.array([0.0,0.0,2.0]) * r #カメラの位置座標 s_CAM_i = np.array([0.0,1.0,0.0]) #カメラの向きの初期値(単位ベクトル）y軸＋方向（北向き） elevation_angle = - 45.0 #カメラの仰角[°] azimuth = 270.0 #カメラの方位[°] h_CAM = 100 #縦のピクセル数 w_CAM = 100 #横のピクセル数 pix_size_CAM = 130.0 * 1.0e-6 #ピクセルサイズ［μm］() f_CAM = 10.22 * 1.0e-3 #焦点距離［mm］ pix_CAM = np.zeros ((h_CAM, w_CAM)) #ピクセル数の配列 #*********************************************************************************************** # In[32]: #カメラの計算==================================== pix = np.zeros((h_CAM,w_CAM)) IFOV_rad = 2.0 * math.atan((pix_size_CAM)/ (2.0 * f_CAM)) IFOV_deg = math.degrees(IFOV_rad) FOV_h_rad = 2.0 * math.atan((pix_size_CAM * h_CAM) / (2.0 * f_CAM)) FOV_h_deg = math.degrees(FOV_h_rad) FOV_w_rad = 2.0 * math.atan((pix_size_CAM * w_CAM) / (2.0 * f_CAM)) FOV_w_deg = math.degrees(FOV_w_rad) #============================================ # In[33]: #関数の定義==================================== #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 #ベクトルの角度を比較する関数(radで入力) def CompVec (A,B,c): a = math.acos(DCOS(A,B)) if a <= c: return 1 else: return 0 #ある点（位置ベクトル）と最も近い点を参照する関数 def SerchNearPoint(P,F,E,n): #P(x,y,z)点,F[i,j](x,y,z)(検索配列)，E[i,j]（輝度）,n[i,j] L_max = np.finfo(np.float64).max E_max = np.zeros(4) for i in range(np.shape(F)[0]): for j in range(np.shape(F)[1]): if L_max > LL(P,F[i,j]): L_max = LL(P,F[i,j]) E_max[0] = E[i,j] E_max[1] = n[i,j,0] E_max[2] = n[i,j,1] E_max[3] = n[i,j,2] return E_max #円筒（縦孔壁面）との交差点を求める関数 def HitCyl_2(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 t1 else: return 0.0 elif t1 <= 0 and t2 > 0: z = s[2] + (t2 * d[2]) if z >= hmin and z <= hmax: return t2 else: return 0.0 elif t1 >= t2 : z = s[2] + (t1 * d[2]) if z >= hmin and z <= hmax: return t1 else: return 0.0 else: z = s[2] + (t2 * d[2]) if z >= hmin and z <= hmax: return t2 else: return 0.0 #一つの位置ベクトルから4隅の位置ベクトルを求める関数 def VecOneFour(F,u,v,h,w): A = np.zeros([4,3]) A[0] = F + ((h**2.0 + w**2.0)**0.5 *0.5 * ((u + v)/np.linalg.norm(u + v))) A[1] = F + ((h**2.0 + w**2.0)**0.5 *0.5 * ((u - v)/ np.linalg.norm(u - v))) A[2] = F + ((h**2.0 + w**2.0)**0.5 *0.5* ((-u + v)/np.linalg.norm(-u + v))) A[3] = F + ((h**2.0 + w**2.0)**0.5 *0.5* ((-u - v)/np.linalg.norm(-u - v))) return A #3点から面積を求める関数 def CalArea(A,B,C): #3点の位置ベクトルを入力（A(x,y,z),B(x,y,z),C(x,y,z)） return (1/2.0) * np.linalg.norm(np.cross((B-A),(C-A))) #ベクトルをx,y,z軸回りに回転させる関数(x→y→zの順) def LotVec(A,t,p,o): Rx = np.array([[1, 0, 0], [0, np.cos(t), -np.sin(t)], [0, np.sin(t), np.cos(t)]]) Ry = np.array([[np.cos(p), 0,np.sin(p)], [0, 1, 0], [-np.sin(p), 0, np.cos(p)]]) Rz = np.array([[np.cos(o), -np.sin(o), 0], [np.sin(o), np.cos(o), 0], [0, 0, 1]]) R = Rz.dot(Ry).dot(Rx) B = np.dot(R,A) return B #底面とのあたり判定 def HitFL(s,d,n,x_max,x_min,y_max,y_min): #s(始点の位置)，d(方向ベクトル)，n（底面の法線ベクトル） t = -(np.dot(s,n))/(np.dot(d,n)) x = s[0] + (t * d[0]) y = s[1] + (t * d[1]) z = s[2] + (t * d[2]) if t < 0.0: return 0 elif(x_min <= x <= x_max) and (y_min <= y <= y_max): return 1 else: return 0 def HitFL2(s,d,n,x_max,x_min,y_max,y_min): #s(始点の位置)，d(方向ベクトル)，n（底面の法線ベクトル） t = -(np.dot(s,n))/(np.dot(d,n)) x = s[0] + (t * d[0]) y = s[1] + (t * d[1]) z = s[2] + (t * d[2]) if t < 0.0: return 0 elif(x_min <= x <= x_max) and (y_min <= y <= y_max): return t else: return 0 #空洞天井とのあたり判定 def HitCC(s,d,n,p,x_max,x_min,y_max,y_min): #s(始点の位置)，d(方向ベクトル)，n（底面の法線ベクトル） t = ((np.dot(p,n))-(np.dot(s,n)))/(np.dot(d,n)) x = s[0] + (t * d[0]) y = s[1] + (t * d[1]) z = s[2] + (t * d[2]) if t < 0.0: return 0 elif(x_min <= x <= x_max) and (y_min <= y <= y_max): return 1 else: return 0 def HitCC2(s,d,n,p,x_max,x_min,y_max,y_min): #s(始点の位置)，d(方向ベクトル)，n（底面の法線ベクトル） t = ((np.dot(p,n))-(np.dot(s,n)))/(np.dot(d,n)) x = s[0] + (t * d[0]) y = s[1] + (t * d[1]) z = s[2] + (t * d[2]) if t < 0.0: return 0 elif(x_min <= x <= x_max) and (y_min <= y <= y_max): return t else: return 0 #============================================ # In[34]: #カメラのグリッド設定=========================(カメラの回転に注意) SUR_CCD = "CCD" a_CCD = np.linspace((((-1*w_CAM/2)* pix_size_CAM) + (pix_size_CAM/2)),(((w_CAM/2)* pix_size_CAM) - (pix_size_CAM/2)),num = w_CAM) b_CCD = np.linspace((((-1*h_CAM/2)* pix_size_CAM) + (pix_size_CAM/2)),(((h_CAM/2)* pix_size_CAM) - (pix_size_CAM/2)),num = h_CAM) #=========================================== #縦孔底面のグリッド設定========================= 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(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_CAM = len(a_CCD) #観測面のセル数 N_CAM = len(b_CCD) #観測面のセル数 #縦孔底 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[35]: #ベクトルの設定***************************************** #太陽光ベクトル(i,j,k) I = np.array([-COS(theta), 0.0, -SIN(theta)]) #カメラのベクトル設定=========================================================================== #カメラの回転 s_CAM = LotVec(s_CAM_i,np.radians(elevation_angle),np.radians(0.0),np.radians(- azimuth)) #カメラの基底直行ベクトルを計算 X = s_CAM Y = (np.cross(np.array([0.0,0.0,1.0]),X)/np.linalg.norm(np.cross(np.array([0.0,0.0,1.0]),X))) Z = np.cross(X,Y) #============================================================================================ #カメラ(CCD)の法線ベクトルs============= s_CAM = s_CAM #カメラ(CCD)の位置ベクトルi i_CAM = i_CAM #カメラ(CCD)の微小面積dS dS = pix_size_CAM ** 2 #仮想スクリーンの位置ベクトル座標 i_SCR_i = np.zeros((M_CAM, N_CAM, 3)) i_SCR = np.zeros((M_CAM, N_CAM, 3)) for i in range(M_CAM): for j in range(N_CAM): i_SCR_i[i,j,0] = a_CCD[i] i_SCR_i[i,j,1] = f_CAM i_SCR_i[i,j,2] = b_CCD[j] i_SCR[i,j]= LotVec(i_SCR_i[i,j], np.radians(elevation_angle), np.radians(0.0), np.radians(-azimuth)) i_SCR[i,j]= np.add(i_SCR[i,j],i_CAM) #カメラの位置ベクトルからのスクリーン面に向かうベクトルH H = np.zeros((M_CAM, N_CAM, 3)) for i in range(M_CAM): for j in range(N_CAM): H[i,j] = np.array(VEC2(i_SCR[i,j],i_CAM)) #=================================== #縦孔底面の法線ベクトルn================ n_FL = np.zeros((M_FL, N_FL, 3)) for i in range(M_FL): for j in range(N_FL): n_FL[i,j,0] = 0.0 n_FL[i,j,1] = 0.0 n_FL[i,j,2] = 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.zeros((M_CC, N_CC, 3)) # for i in range(M_CC): # for j in range(N_CC): # n_CC[i,j,0] = 0.0 # n_CC[i,j,1] = 0.0 # n_CC[i,j,2] = 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 #空洞天井の法線ベクトルn================ n_CC = np.zeros((M_CC, N_CC, 3)) for i in range(M_CC): for j in range(N_CC): n_CC[i,j,0] = 0.0 n_CC[i,j,1] = 0.0 n_CC[i,j,2] = -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 #==================================== # #空洞天井の微小面積dA # dA_CC = np.zeros((M_CC)) # for i in range(M_CC): # dA_CC[i] = dr * a_CC[i] * dphi # #===================================== #************************************************* # In[36]: #放射照度データの入力（CSV）入力**************************** #太陽直達光 E_StoFL = np.loadtxt('PPP_StoFL_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_StoWA = np.loadtxt('PPP_StoWA_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_StoCW = np.loadtxt('PPP_StoCW_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') #底からの反射光 E_FLtoWA = np.loadtxt('PPP_FLtoWA_t45.0h1.0h1.0_cr2.0_2023-03-08.csv',delimiter=',') E_FLtoCW = np.loadtxt('PPP_FLtoCW_t45.0h1.0h1.0_cr2.0_2023-03-08.csv',delimiter=',') E_FLtoCC = np.loadtxt('PPP_FLtoCC_t45.0h1.0h1.0_cr2.0_2023-03-09.csv',delimiter=',') #縦孔壁からの反射光 E_WAtoFL = np.loadtxt('PPP_WAtoFL_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_WAtoWA = np.loadtxt('PPP_WAtoWA_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_WAtoCW = np.loadtxt('PPP_WAtoCW_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') #空洞壁からの反射光 E_CWtoFL = np.loadtxt('PPP_CWtoFL_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_CWtoWA = np.loadtxt('PPP_CWtoWA_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_CWtoCW = np.loadtxt('PPP_CWtoCW_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') E_CWtoCC = np.loadtxt('PPP_CWtoCC_t45.0h1.0h1.0_cr2.0_2023-03-10.csv',delimiter=',') #****************************************************** # In[37]: #出力結果の結合 R_f = R E_FL_all = (E_StoFL) + (R_f * E_WAtoFL) + (R_f * E_CWtoFL) E_WA_all = (E_StoWA) + (R_f * E_FLtoWA) + (R_f * E_WAtoWA) + (R_f * E_CWtoWA) E_CW_all = (E_StoCW) + (R_f * E_FLtoCW) + (R_f * E_WAtoCW) + (R_f * E_CWtoCW) E_CC_all = (R_f * E_FLtoCC) + (R_f * E_CWtoCC) # In[38]: #解析1======================================================================================== #=========================================================================================== #計算開始の時刻=========== start_time = time.time() print("--start--") #======================= #空洞面（CW）からのカメラ（CAM）に入射する放射照度E_CAM_FL（放射照度） E_CAM_CW = np.zeros((M_CAM, N_CAM)) #観測面のセルの照度の初期化 MM = M_CAM NN = N_CAM M = M_CW N = N_CW F = F_CW n = n_CW E = E_CW_all H4 = np.zeros([4,3]) i_POR_4 = np.zeros([4,3]) t4 = np.zeros([4]) for ii in range(MM): for jj in range(NN): #対象の空洞壁に当たるか判定 if HitCyl(i_CAM,H[ii,jj],cr*r,d2,0.0)==0: E_CAM_CW[ii, jj] = 0.0 elif HitCyl(i_CAM,H[ii,jj],r,d2+d1,d2)==1: E_CAM_CW[ii, jj] = 0.0 else:#当たるとき i_SCR_4_iijj = VecOneFour(i_SCR[ii,jj],Z,Y,pix_size_CAM,pix_size_CAM)#ピクセルの4隅の位置ベクトルを取得 #カメラの点からピクセルの4隅の位置に向かうベクトルを計算 H4[0] = VEC2(i_SCR_4_iijj[0],i_CAM) H4[1] = VEC2(i_SCR_4_iijj[1],i_CAM) H4[2] = VEC2(i_SCR_4_iijj[2],i_CAM) H4[3] = VEC2(i_SCR_4_iijj[3],i_CAM) #4つのベクトルがすべて壁面に当たる場合 if HitCyl(i_CAM,H4[0],(cr*r),d2,0.0)==1 and HitCyl(i_CAM,H4[1],cr*r,d2,0.0)==1 and HitCyl(i_CAM,H4[2],(cr*r),d2,0.0)==1 and HitCyl(i_CAM,H4[3],(cr*r),d2,0.0)==1: #4点から面積を計算 t4[0] = HitCyl_2(i_CAM,H4[0],cr*r,d2,0.0) i_POR_4[0] = np.add(i_CAM, H4[0]*t4[0]) t4[1] = HitCyl_2(i_CAM,H4[1],cr*r,d2,0.0) i_POR_4[1] = np.add(i_CAM, H4[1]*t4[1]) t4[2] = HitCyl_2(i_CAM,H4[2],cr*r,d2,0.0) i_POR_4[2] = np.add(i_CAM, H4[2]*t4[2]) t4[3] = HitCyl_2(i_CAM,H4[3],cr*r,d2,0.0) i_POR_4[3] = np.add(i_CAM, H4[3]*t4[3]) dA = CalArea(i_POR_4[0],i_POR_4[1],i_POR_4[2]) + CalArea(i_POR_4[1],i_POR_4[2],i_POR_4[3]) #中心から最も近い点（放射照度がある）を参照 t = HitCyl_2(i_CAM,H[ii,jj],cr*r,d2,0.0) i_PRO = np.add(i_CAM,H[ii,jj] * t) E_dA = SerchNearPoint(i_PRO,F,E,n) L2 = LL(i_CAM,i_PRO) E_CAM_CW[ii, jj] = Ref("Lambert", I, -H[ii,jj], np.array([E_dA[1], E_dA[2], E_dA[3]])) * (E_dA[0] * dA * (DCOS(H[ii,jj], s_CAM))) / (L2) else: E_CAM_CW[ii, jj] = 0.0 if E_CAM_CW[ii, jj] != E_CAM_CW[ii, jj]: E_CAM_CW[ii, jj] = 0.0 print("CW:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) E_CW_all = np.sum(E_CAM_CW) print("CW:100.00 %") print("E_CW_all:{}".format(E_CW_all)) #==================================================================================================== #==================================================================================================== date = datetime.date.today() np.savetxt("CAM_CW_AZ{}EL{}_PO{}th{}cr{}R{}_{}.csv".format(azimuth, elevation_angle, i_CAM[2] ,theta, cr,R, date), E_CAM_CW, delimiter=',') #計算終了時刻======================== end_time = time.time() - start_time end_time = end_time / 60.0 #================================= # In[39]: #解析2======================================================================================== #=========================================================================================== #計算開始の時刻=========== start_time = time.time() print("--start--") #======================= #縦孔壁（WA）からのカメラ（CAM）に入射する放射照度E_CAM_WA（放射照度） E_CAM_WA = np.zeros((M_CAM, N_CAM)) #観測面のセルの照度の初期化 MM = M_CAM NN = N_CAM M = M_WA N = N_WA F = F_WA n = n_WA E = E_WA_all H4 = np.zeros([4,3]) i_POR_4 = np.zeros([4,3]) t4 = np.zeros([4]) for ii in range(MM): for jj in range(NN): #対象の空洞壁に当たるか判定 if HitCyl(i_CAM,H[ii,jj],r,d2+d1,d2)==0: E_CAM_WA[ii, jj] = 0.0 else:#当たるとき i_SCR_4_iijj = VecOneFour(i_SCR[ii,jj],Z,Y,pix_size_CAM,pix_size_CAM)#ピクセルの4隅の位置ベクトルを取得 #カメラの点からピクセルの4隅の位置に向かうベクトルを計算 H4[0] = VEC2(i_SCR_4_iijj[0],i_CAM) H4[1] = VEC2(i_SCR_4_iijj[1],i_CAM) H4[2] = VEC2(i_SCR_4_iijj[2],i_CAM) H4[3] = VEC2(i_SCR_4_iijj[3],i_CAM) #4つのベクトルがすべて壁面に当たる場合 if HitCyl(i_CAM,H4[0],r,d2+d1,d2)==1 and HitCyl(i_CAM,H4[1],r,d2+d1,d2)==1 and HitCyl(i_CAM,H4[2],r,d2+d1,d2)==1 and HitCyl(i_CAM,H4[3],r,d2+d1,d2)==1: #4点から面積を計算 t4[0] = HitCyl_2(i_CAM,H4[0],r,d2+d1,d2) i_POR_4[0] = np.add(i_CAM, H4[0]*t4[0]) t4[1] = HitCyl_2(i_CAM,H4[1],r,d2+d1,d2) i_POR_4[1] = np.add(i_CAM, H4[1]*t4[1]) t4[2] = HitCyl_2(i_CAM,H4[2],r,d2+d1,d2) i_POR_4[2] = np.add(i_CAM, H4[2]*t4[2]) t4[3] = HitCyl_2(i_CAM,H4[3],r,d2+d1,d2) i_POR_4[3] = np.add(i_CAM, H4[3]*t4[3]) dA = CalArea(i_POR_4[0],i_POR_4[1],i_POR_4[2]) + CalArea(i_POR_4[1],i_POR_4[2],i_POR_4[3]) #中心から最も近い点（放射照度がある）を参照 t = HitCyl_2(i_CAM,H[ii,jj],r,d2+d1,d2) i_PRO = np.add(i_CAM,H[ii,jj] * t) E_dA = SerchNearPoint(i_PRO,F,E,n) L2 = LL(i_CAM,i_PRO) E_CAM_WA[ii, jj] = Ref("Lambert", I, -H[ii,jj], np.array([E_dA[1], E_dA[2], E_dA[3]])) * (E_dA[0] * dA * (DCOS(H[ii,jj], s_CAM))) / (L2) else: E_CAM_WA[ii, jj] = 0.0 if E_CAM_WA[ii, jj] != E_CAM_WA[ii, jj]: E_CAM_WA[ii, jj] = 0.0 print("CW:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) E_WA_all = np.sum(E_CAM_WA) print("WA:100.00 %") print("E_WA_all:{}".format(E_WA_all)) #==================================================================================================== #==================================================================================================== date = datetime.date.today() np.savetxt("CAM_WA_AZ{}EL{}_PO{}th{}cr{}R{}_{}.csv".format(azimuth, elevation_angle, i_CAM[2] ,theta, cr,R, date), E_CAM_WA, delimiter=',') #計算終了時刻======================== end_time = time.time() - start_time end_time = end_time / 60.0 #================================= # In[40]: #解析3======================================================================================== #=========================================================================================== #計算開始の時刻=========== start_time = time.time() print("--start--") #======================= #底面（FL）からのカメラ（CAM）に入射する放射照度E_CAM_FL（放射照度） E_CAM_FL = np.zeros((M_CAM, N_CAM)) #観測面のセルの照度の初期化 MM = M_CAM NN = N_CAM M = M_FL N = N_FL F = F_FL n = n_FL E = E_FL_all n1 = np.array([0.0,0.0,1.0]) H4 = np.zeros([4,3]) i_POR_4 = np.zeros([4,3]) t4 = np.zeros([4]) for ii in range(MM): for jj in range(NN): if HitFL(i_CAM,H[ii,jj],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min)==0: #対象の底に当たるか判定 E_CAM_FL[ii, jj] = 0.0 elif HitCyl(i_CAM,H[ii,jj],r,d2+d1,d2)==1: #縦孔の壁に当たるか判定 E_CAM_FL[ii, jj] = 0.0 elif HitCyl(i_CAM,H[ii,jj],cr*r,d2,0.0)==1: #空洞の壁に当たるか判定 E_CAM_FL[ii, jj] = 0.0 else: #対象の底面に当たるとき i_SCR_4_iijj = VecOneFour(i_SCR[ii,jj],Z,Y,pix_size_CAM,pix_size_CAM) #ピクセルの4隅の位置ベクトルを取得 #カメラの点からピクセルの4隅の位置に向かうベクトルを計算 H4[0] = VEC2(i_SCR_4_iijj[0],i_CAM) H4[1] = VEC2(i_SCR_4_iijj[1],i_CAM) H4[2] = VEC2(i_SCR_4_iijj[2],i_CAM) H4[3] = VEC2(i_SCR_4_iijj[3],i_CAM) #4つのベクトルがすべて壁面に当たる場合 if HitFL(i_CAM,H4[0],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min)==1 and HitFL(i_CAM,H4[1],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min)==1 and HitFL(i_CAM,H4[2],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min)==1 and HitFL(i_CAM,H4[3],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min)==1: #4点から面積を計算 t4[0] = HitFL2(i_CAM,H4[0],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_POR_4[0] = np.add(i_CAM, H4[0]*t4[0]) t4[1] = HitFL2(i_CAM,H4[1],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_POR_4[1] = np.add(i_CAM, H4[1]*t4[1]) t4[2] = HitFL2(i_CAM,H4[2],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_POR_4[2] = np.add(i_CAM, H4[2]*t4[2]) t4[3] = HitFL2(i_CAM,H4[3],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_POR_4[3] = np.add(i_CAM, H4[3]*t4[3]) dA = CalArea(i_POR_4[0],i_POR_4[1],i_POR_4[2]) + CalArea(i_POR_4[1],i_POR_4[2],i_POR_4[3]) #中心から最も近い点（放射照度がある）を参照 t = HitFL2(i_CAM,H[ii,jj],n1,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_PRO = np.add(i_CAM,H[ii,jj] * t) E_dA = SerchNearPoint(i_PRO,F,E,n) L2 = LL(i_CAM,i_PRO) #カメラに入射する光エネルギーを計算 E_CAM_FL[ii, jj] = Ref("Lambert", I, -H[ii,jj], np.array([E_dA[1], E_dA[2], E_dA[3]])) * (E_dA[0] * dA * (DCOS(H[ii,jj], s_CAM))) / (L2) else: E_CAM_FL[ii, jj] = 0.0 if E_CAM_FL[ii, jj] != E_CAM_FL[ii, jj]: E_CAM_FL[ii, jj] = 0.0 print("FL:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) #進行状況を表示 E_FL_all = np.sum(E_CAM_FL) print("FL:100.00 %") print("E_FL_all:{}".format(E_FL_all)) #==================================================================================================== #==================================================================================================== date = datetime.date.today() np.savetxt("CAM_FL_AZ{}EL{}_PO{}th{}cr{}R{}_{}.csv".format(azimuth, elevation_angle, i_CAM[2] ,theta, cr,R, date), E_CAM_FL, delimiter=',') #計算終了時刻======================== end_time = time.time() - start_time end_time = end_time / 60.0 #================================= # In[41]: #解析4======================================================================================== #=========================================================================================== #計算開始の時刻=========== start_time = time.time() print("--start--") #======================= #空洞天井（CC）からのカメラ（CAM）に入射する放射照度E_CAM_CC（放射照度） E_CAM_CC = np.zeros((M_CAM, N_CAM)) #観測面のセルの照度の初期化 MM = M_CAM NN = N_CAM M = M_CC N = N_CC F = F_CC n = n_CC E = E_CC_all n1 = np.array([0.0,0.0,-1.0]) p = np.array([0.0,0.0,1.0])*r H4 = np.zeros([4,3]) i_POR_4 = np.zeros([4,3]) t4 = np.zeros([4]) for ii in range(MM): for jj in range(NN): if HitCC(i_CAM,H[ii,jj],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min)==0: #対象の底に当たるか判定 E_CAM_CC[ii, jj] = 0.0 elif HitCyl(i_CAM,H[ii,jj],r,d2+d1,d2)==1: #縦孔の壁に当たるか判定 E_CAM_CC[ii, jj] = 0.0 elif HitCyl(i_CAM,H[ii,jj],cr*r,d2,0.0)==1: #空洞の壁に当たるか判定 E_CAM_CC[ii, jj] = 0.0 else: #対象の底面に当たるとき i_SCR_4_iijj = VecOneFour(i_SCR[ii,jj],Z,Y,pix_size_CAM,pix_size_CAM) #ピクセルの4隅の位置ベクトルを取得 #カメラの点からピクセルの4隅の位置に向かうベクトルを計算 H4[0] = VEC2(i_SCR_4_iijj[0],i_CAM) H4[1] = VEC2(i_SCR_4_iijj[1],i_CAM) H4[2] = VEC2(i_SCR_4_iijj[2],i_CAM) H4[3] = VEC2(i_SCR_4_iijj[3],i_CAM) #4つのベクトルがすべて壁面に当たる場合 if HitCC(i_CAM,H4[0],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min)==1 and HitCC(i_CAM,H4[1],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min)==1 and HitCC(i_CAM,H4[2],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min)==1 and HitCC(i_CAM,H4[3],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min)==1: #4点から面積を計算 t4[0] = HitCC2(i_CAM,H4[0],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min) i_POR_4[0] = np.add(i_CAM, H4[0]*t4[0]) t4[1] = HitCC2(i_CAM,H4[1],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min) i_POR_4[1] = np.add(i_CAM, H4[1]*t4[1]) t4[2] = HitCC2(i_CAM,H4[2],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min) i_POR_4[2] = np.add(i_CAM, H4[2]*t4[2]) t4[3] = HitCC2(i_CAM,H4[3],n1,p,CC_x_max,CC_x_min,CC_y_max,CC_y_min) i_POR_4[3] = np.add(i_CAM, H4[3]*t4[3]) dA = CalArea(i_POR_4[0],i_POR_4[1],i_POR_4[2]) + CalArea(i_POR_4[1],i_POR_4[2],i_POR_4[3]) #中心から最も近い点（放射照度がある）を参照 t = HitCC2(i_CAM,H[ii,jj],n1,p,FL_x_max,FL_x_min,FL_y_max,FL_y_min) i_PRO = np.add(i_CAM,H[ii,jj] * t) E_dA = SerchNearPoint(i_PRO,F,E,n) L2 = LL(i_CAM,i_PRO) #カメラに入射する光エネルギーを計算 E_CAM_CC[ii, jj] = Ref("Lambert", I, -H[ii,jj], np.array([E_dA[1], E_dA[2], E_dA[3]])) * (E_dA[0] * dA * (DCOS(H[ii,jj], s_CAM))) / (L2) else: E_CAM_CC[ii, jj] = 0.0 if E_CAM_CC[ii, jj] != E_CAM_CC[ii, jj]: E_CAM_CC[ii, jj] = 0.0 print("CC:{}.{} %".format(int(100 * ((ii /MM))), int(100 * (jj/NN)) )) #進行状況を表示 E_CC_all = np.sum(E_CAM_CC) print("FL:100.00 %") print("E_FL_all:{}".format(E_CC_all)) #==================================================================================================== #==================================================================================================== date = datetime.date.today() np.savetxt("CAM_CC_AZ{}EL{}_PO{}th{}cr{}R{}_{}.csv".format(azimuth, elevation_angle, i_CAM[2] ,theta, cr,R, date), E_CAM_CC, delimiter=',') #計算終了時刻======================== end_time = time.time() - start_time end_time = end_time / 60.0 #================================= # In[42]: E_CAM_ALL = np.flipud(np.add(np.add(np.add(E_CAM_WA.T,E_CAM_CW.T),E_CAM_FL.T),E_CAM_CC.T)) plt.imshow(E_CAM_ALL) date = datetime.date.today() plt.title("θ = {} , a = {}".format(theta,cr)) plt.colorbar (label="E/J []") print(np.max(E_CAM_ALL),np.min(E_CAM_ALL))