""" @title: ray3cer - raytracing software for photon shot-noise simulations @author: Ugo Lo Cicero @email: ugo.locicero@inaf.it @licence: GNU General Public License v3.0 v.15 03/05/2022 - save area of various filter surfaces (membrane, frame, baffle) - save temperature of various filter surfaces (membrane, frame, baffle) - save Nph, number of simulated photons v.14 18/03/2022 - import matplotlib module only in "__main__" code v.13 30/11/2021 - modificato per uso con file di configurazione unico v.12 25/06/2021 - mappa con scala di colore logaritmica 08/10/2021 - mappa e grafico convergenza aggiornati durante la simulazione (il profilo non funziona) 08/10/2021 - possibilità di indicare la riflettività di pareti o frame in funzione della lunghezza d'onda v.11 09/06/2021 - corretto l'ordine di priorità tra la keyword 'flt_transrefl' e la cfg.bypass_transrefl - selezione wave_arr solo >0 per le analisi, in modo da avere i risultati anche in caso di Ctrl-C - modificata la funzione Circle3D per mostrare il tilt (solo tilt attorno a X) - verifica delle keyword e interruzione in caso di keyword non valida - quando viene verificata l'intersezione con il frame viene calcolata la distanza lungo il piano Z, piuttosto che lungo il piano del filtro - generazione dei fotoni dai baffle: viene calcolato il punto di generazione considerando il tilt dei filtri - è possibile specificare il tilt per ogni filtro con la keyword 'tilt' (15/06) - corretto bug per cui non veniva utilizzata la riflettività specifica del baffle durante l'interazione, ma soltanto per l'emissione (15/06) - corretto bug per cui non veniva selezionato correttamente la superficie del frame (up/down) da utilizzare (15/06) - il detector può avere adesso keywords (15/06) - evitato il runtime warning di divisione per zero nel calcolo dello spettro emesso con controllo di area o Pbb =0 v.10 27/04/2021 - possibilità di specificare la riflettività di singoli frame diversa tra sopra e sotto, con le keyword frame_refl_down e frame_refl_up v.09 19/04/2021 - possibilità di forzare trasmissività e riflettività di singoli filtri ('flt_transrefl') v.08 29/03/2021 - implementata simulazione fori v.07 08/03/2021 - aggiunta opzione per specificare la composizione di ogni filtro - aggiunto supporto per filtri "open": trasmissività 1 e i fotoni vengono seguiti anche se tornano indietro - disaccoppiate T_filtro e T_frame v.06 18/01/2021 - aggiunta emissione dai frame e dai baffle v.05 15/11/2020 - corretta la distribuzione angolare di emissione dei fotoni v.04 22/06/2020 - opzione per salvare informazioni aggregate per le traiettorie v.03 29/05/2020 - salva lo stato random iniziale, in modo che in seguito sia possibile replicare v.02 06/05/2020 - aggiunte le pareti laterali v.01 28/04/2020 - rispetto a ray3cer2 non usa una struttura ad albero per esplorare, ma uno stack, non avendo bisogno di mantenere tutte le informazioni sui nodi in memoria """ #%% import math import numpy as np import copy import time, sys, os, datetime from transmission import getTransRefl import scipy.integrate import scipy.constants as const from scipy.stats import rv_continuous # for random generation with given pdf import multiprocessing as mp import queue import gc # garbage collector #rc('text', usetex=True) # enable latex syntax for plots from r3_config_file import obj def init(): import r3_config_file return(r3_config_file.init()) #%% 3d plot functions def circle3d(ax,Z,r,tilt=0,center=(0,0)): # p = Circle(center,r, facecolor='None', edgecolor='b', lw=1) # ax.add_patch(p) # art3d.pathpatch_2d_to_3d(p, z=Z, zdir="z") θ = np.linspace(0, 2 * np.pi, 201) c1, c2, c3 = center[0], center[1], Z a1, a2, a3 = 0, 1, 0 b = rotate_vector(np.array([1, 0, 0]),(0, tilt, 0)) #b1, b2, b3 = 0, 1, 0 print(b) b1, b2, b3 = b[0], b[1], b[2] x = c1 + r*np.cos(θ)*a1 + r*np.sin(θ)*b1 y = c2 + r*np.cos(θ)*a2 + r*np.sin(θ)*b2 z = c3 + r*np.cos(θ)*a3 + r*np.sin(θ)*b3 ax.plot(x, y, z) def plot_cylinder(ax,r,z1,z2): x=np.linspace(-r, r, 100) z=np.linspace(z1, z2, 100) Xc, Zc=np.meshgrid(x, z) Yc = np.sqrt(r**2-Xc**2) rstride = 20 cstride = 10 ax.plot_surface(Xc, Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride) ax.plot_surface(Xc, -Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride) #%% Progress bar def progressbar(percent): step = percent/5 #print("\r[{:20s}] {:.2f}%".format('='*int(step), percent)) sys.stdout.write("\r[{:20s}] {:.2f}%".format('='*int(step), percent)) sys.stdout.flush() #%% def Rx(theta): return( np.array([ [1, 0, 0], [0, np.cos(theta), -np.sin(theta)], [0, np.sin(theta), np.cos(theta)]])) def Ry(theta): return( np.array([ [np.cos(theta), 0, np.sin(theta)], [0, 1, 0], [-np.sin(theta), 0, np.cos(theta)]])) def Rz(theta): return( np.array([ [np.cos(theta),-np.sin(theta),0], [np.sin(theta),np.cos(theta), 0], [0, 0, 1]])) def rotate_vector(vector,angles): # Rotate a (column) vector around origin # angles is [angle around x-axis, angle around y-axis,angle around z-axis] if len(vector.shape) == 1: # if 1d array it has to be transformed in a col vector vector = np.array([vector]).T rot_vector = Rz(angles[2]) @ Ry(angles[1]) @ Rx(angles[0]) @ vector return(rot_vector) class filter: def __init__(self,z,r,T,baffles,wave,trans,refl,tilt_dir=0,emissivity=1,surf_idx=None, cfg=None, **kwargs): ''' z: z-coordinate, r: filter radius, T: temperature, wave: wavelenghts for which transmissivity and reflectivity are given trans and refl: transmissivity and reflectivity array (function of lambda) tilt: filter tilt (rad), tilt_dir: tilting direction (rad) tilt_dir is an angle around z. tilt_dir = 0 means tilt around y-axis, tilt_dir = pi means around x ''' self.z = z self.r = r self.T = T self.wave = np.array(wave) self.trans = np.array(trans) self.refl = np.array(refl) self.tilt = kwargs.get('tilt',0) tilt = self.tilt a = np.sin(tilt)*np.cos(tilt_dir) b = np.sin(tilt)*np.sin(tilt_dir) c = np.cos(tilt) d = -c*z # passing on (0,0,z) self.plane = np.array([a,b,c,d]) # plane eq.: ax + by + cz + d = 0 self.surfaces = kwargs.get('surfaces',None) self.surf_idx = surf_idx next_surf = cfg.surfdef[surf_idx+1] if surf_idx0: # if not detector if baffles is None or baffles ==[]: baffles = [dict()] for i,baffle in enumerate(baffles): # currently only 1 baffle is supported if baffle.get('rad1') is None: baffle['rad1'] = r # rad1: top radius of the baffle segment (if absent assumed = filter radius) if baffle.get('rad2') is None: baffle['rad2'] = baffle['rad1'] # rad2: bottom radius of the baffle segment if it is a cone (if absent assumed = top radius) if baffle.get('zrad1') is None: baffle['zrad1'] = z # zrad1: z position of top radius of the baffle segment (if absent assumed = filter z position) if baffle.get('zrad2') is None: baffle['zrad2'] = next_surf[0] if next_surf else None # zrad2: z position of bottom radius of the baffle segment if it is a cone(if absent assumed = colder filter z position) if baffle.get('refl') is None: baffle['refl'] = cfg.refl_walls # refl: baffle segment reflectivity (if absent cfg.refl_walls is assumed) if baffle.get('T') is None: baffle['T'] = next_surf[2] if next_surf else None # refl: baffle segment temperature (if absent next filter T is assumed) if baffle['zrad1'] <= baffle['zrad2']: raise ValueError(f'Error for surface {surf_idx}, baffle {i}, name:{kwargs.get("name_baffle","")}: zrad1 has to be > zrad2.') self.baffle_refl = baffle['refl'] if cfg.only_frame_emission: area = np.pi*(baffle['rad1']**2 - self.r**2) # frame emission area in mm^2 self.frame_area = area self.Tframe = kwargs.get('Tframe') if not self.Tframe: self.Tframe = self.T self.emissivity = 1-self.get_refl('frame_refl_down', T=self.Tframe) if cfg.verbose: print(f"FRAME: baffle['rad1']={baffle['rad1']}; r={self.r}; area={area}; emiss.={self.emissivity}; T={self.Tframe}") self.bb,self.probab,self.Pbb = self.bbRadiation(T=self.Tframe,area=area) elif cfg.only_baffle_emission: # Assuming a single segment. Can be a cone. r1,r2,h = baffle['rad1'],baffle['rad2'],baffle['zrad1']-baffle['zrad2'] area = np.pi*(r1+r2)*math.sqrt((r1-r2)**2+h**2) # conical frustum; frame emission area in mm^2 self.baffle_area = area self.Tbaffle = baffle['T'] self.emissivity = 1-self.get_refl('baffle_refl',T=baffle['T']) self.bb,self.probab,self.Pbb = self.bbRadiation(T=baffle['T'],area=area) if cfg.verbose: print(f"BAFFLE: baffle['rad1']={r1}; baffle['rad2']={r2};h={h}; r={self.r}; area={area}; T={baffle['T']}; emiss.={self.emissivity}; Pbb={self.Pbb}") else: area = np.pi*(self.r**2) self.area = area self.emissivity = 1-self.trans-self.refl self.bb,self.probab,self.Pbb = self.bbRadiation(T=self.T, area=area) # filter emission area in mm^2 if cfg.verbose: print(f"FILTER: r={self.r}; area={area}; T={self.T}; emiss.={self.emissivity}; Pbb={self.Pbb}") self.baffles = baffles def get_refl(self, refl_spec, wave=None, T=None): re = getattr(self, refl_spec) if type(re) is np.ndarray: if T: wave = 2.8977685e-3/T * 1e9 # Wien displacement law if len(re.shape) == 2: w_arr = re[0] r_arr = re[1] else: w_arr = self.wave r_arr = re refl = np.interp(wave,w_arr,r_arr) return(refl) else: return(re) def get_probab(self,wavelenght): trans = np.double(np.interp(wavelenght,self.wave,self.trans)) # transmission probability refl = np.double(np.interp(wavelenght,self.wave,self.refl)) # reflection probability abs = 1-trans-refl # absorption probability return([trans,refl,abs]) def bbRadiation(self,wavelenghts=None,emiss=None,T=None,area=None): # Calculate black body radiation if wavelenghts is None: wavelenghts = self.wave if T is None: T = self.T if area == 0: return([0,0,0]) if area is None: area = np.pi*(self.r**2) lm = wavelenghts*1e-9 # convert wavelength in meters h = const.Planck # Planck Constant [J*s]=[m^2 Kg s^-1] k = const.Boltzmann # Boltzman Constant [m^2 Kg s^-2 K^-1] c = const.c # Light Speed [m/s] bb = bb_ph = np.zeros(lm.size) exponent = h*c/(k*lm*T) expOK = exponent<700 # check if exponent would cause overflow if emiss is None: # default filter emissivity is calculated as 1-trans-refl emiss = self.emissivity if type(emiss) is np.ndarray: # if emissivity is in function of wavelenght select values if emiss.shape[0] > 2: emiss_bb = emiss[expOK] else: emiss_bb = emiss[0] else: emiss_bb = emiss bb[expOK] = emiss_bb * 2*h*c**2/(lm[expOK]**5)/(np.exp(exponent[expOK])-1)*1e-9 # W/m^2/sr/nm bb = bb * 1e-6 * area # W/sr/nm from the emitting surface of the filter bb = bb * np.pi # W/nm toward the full hemishpere (2*pi steradiant)/2: the /2 is due to the cos(theta) effect (Lambert's cosine law) Pbb = scipy.integrate.simps(bb,wavelenghts) # W (for normalization) if Pbb == 0: return([bb,0,0]) bb_ph[expOK] = bb[expOK] * lm[expOK] / (h*c) # divide by photon energies to find photons/s/nm probab = scipy.integrate.cumtrapz(bb_ph,lm,initial=None)/scipy.integrate.trapz(bb_ph,lm) return([bb,probab,Pbb]) # return photon flux spectrum, probability and integrated power for normalization class photon: def __init__(self,pos=None,dir=None,wav=None,prob=np.double(1),cfg=None): self.cfg = cfg self.pos = pos self.dir = dir self.wav = wav self.prob = prob self.ntrans = 0 self.planePoint = np.zeros(3) # initialized here for speed def generate(self,Z=0,r=0,T=0,filter=None,epsilon=1e-6): a,b,c,d = 0,0,1,Z # default plane if filter: Z,r,T = filter.z,filter.r,filter.T a,b,c,d = filter.plane #---- the following line assigns a fixed starting position if self.cfg.fixEmissionPosition: self.pos = np.array(cfg.pos0) phi_ran = math.atan2(cfg.pos0[1], cfg.pos0[0]) # phi = atan(x/y) for the baffles if cfg.enableTracingPlot: print(f"Phi = {phi_ran}") r = cfg.pos0[0] / math.cos(phi_ran) # r=x/cos(phi) for baffles elif self.cfg.only_baffle_emission: # assuming single-segment cylindrical baffle r,z1,z2 = filter.baffles[0]['rad1'], filter.baffles[0]['zrad1'], filter.baffles[0]['zrad2'] filter_low = filter.surfaces[filter.surf_idx+1] plane_low = filter_low.plane if filter_low.z == 0: # detector z_min = 0 else: z_min = z2 - filter_low.baffles[0]['rad1'] * abs(math.cos(filter_low.tilt)) z_max = z1 + r * abs(math.cos(filter.tilt)) while(1): z, phi_ran = np.random.uniform(low=[z_min,0], high=[z_max,2*np.pi]) # generate random position on the baffle x = r*math.cos(phi_ran) y = r*math.sin(phi_ran) #plane eq.: ax + by + cz + d = 0; z = - (ax + by + d)/c z_high = - (a*x + b*y + d)/c # z di intersezione tra il piano del filtro ed il baffle alle coordinate x,y: z_high z_low = - (plane_low[0]*x + plane_low[1]*y + plane_low[3])/plane_low[2] # z di intersezione tra il piano del filtro basso ed il baffle alle coordinate x,y: z_high if (z_low < z) and (z_high > z): # se z è al di fuori di z_low - z_high scartare il punto e ripetere break self.pos = np.array([x,y,z]) elif self.cfg.only_frame_emission: # assign initial position accordig to a uniform distribution within the radius p_rad2 = 0 first_run = True # force to enter the while loop while(not ((p_rad2 > r**2) and (p_rad2 <= filter.baffles[0]['rad1']**2)) or first_run): # check if point is within the radius first_run = False x,y = np.random.uniform(low=[-filter.baffles[0]['rad1'],-filter.baffles[0]['rad1']], high=[filter.baffles[0]['rad1'],filter.baffles[0]['rad1']]) # generate random x,y on a square 2r*2r z = (-a*x-b*y-d)/c p_rad2 = x**2+y**2+(z-Z)**2 # point distance from (0,0,Z), squared; not using norm for speed self.pos = np.array([x,y,z]) else: # assign initial position accordig to a uniform distribution within the radius p_rad2 = 0 first_run = True # force to enter the while loop while(not p_rad2 <= r**2 or first_run): # check if point is within the radius first_run = False x,y = np.random.uniform(low=[-r,-r], high=[r,r]) # generate random x,y on a square 2r*2r z = (-a*x-b*y-d)/c p_rad2 = x**2+y**2+(z-Z)**2 # point distance from (0,0,Z), squared; not using norm for speed self.pos = np.array([x,y,z]) # PHOTONS GENERATION ANGLE: ## Select this to calculate angle to the far end of the detector (use only without multiple reflections) # also for tilt=0 there are cases in which photon can only enter trough multiple reflections! Select large angle option when using multiple reflections! #self.gen_angle = np.arctan((cfg.detector_diam/2+p_rad)/Z) ## Select this to force angle to the far end of the second surface (for tilt >0) #self.gen_angle = np.arctan((cfg.surfdef[1][1]+p_rad)/(Z-cfg.surfdef[1][0])) ## Select this to emit toward the whole hemisphere self.gen_angle = np.pi/2 if self.cfg.fixEmissionDirection: self.dir = np.array(self.cfg.dir0) # #---- generate a random direction vector taking a uniform sampling in z [-1,-cos(accept_angle)), and phi [0,2*pi) # z, phi = np.random.uniform(low=[-1,0], high=[-np.cos(self.gen_angle),2*np.pi]) # # non serve, è un caso speciale incluso nella precedente: z, phi = np.random.uniform(low=[-1,0], high=[0,2*np.pi]) # self.dir = np.array([math.sqrt(1-z**2)*np.cos(phi), math.sqrt(1-z**2)*np.sin(phi),z]) else: #---- generate a random direction vector taking a uniform sampling between sin(theta_min=0)=0 and sin(theta_max=gen_angle) in theta_ax [-1,-cos(accept_angle)), and phi [0,2*pi) sin_theta_2, phi = np.random.uniform(low=[0,0], high=[np.sin(self.gen_angle),2*np.pi]) sin_theta = math.sqrt(sin_theta_2) theta = math.asin(sin_theta) self.dir = np.array([sin_theta*math.cos(phi), sin_theta*math.sin(phi), -math.cos(theta)]) if self.cfg.only_baffle_emission: if self.cfg.enableTracingPlot: print(f"Original dir: {self.dir}") # rotate to emit from cylindrical baffle self.dir = rotate_vector(self.dir,[0,np.pi/2,phi_ran]).flatten() # angles is [angle around x-axis, angle around y-axis,angle around z-axis] if self.cfg.fixEmissionDirection: self.dir = np.array(self.cfg.dir0) self.wav = self.cfg.default_photon_wav self.pos0 = self.pos self.dir0 = self.dir def intersect(self, plane=None,roi=0, filter=None, epsilon=1e-6, update=True, baffles=None): tStart = time.time() if filter: plane = filter.plane roi = filter.r Z = filter.z else: Z = -plane[3]/plane[2] # filter center z planeNormal = plane[:-1] ndotu = planeNormal.dot(self.dir) if abs(ndotu) < epsilon: return(None) #a[0],a[1],a[2]=x,y,z is much faster than a[:]=[x,y,z] self.planePoint[0] = 0 self.planePoint[1] = 0 self.planePoint[2] = Z # point on plane for x=0, y=0 (set single element for speed) # Plain with normal n, passing for p0: (p-p0).n=0 # Line eq. passing for l0: p = l0+l*d # Combining and solving for d: d=((p0-l0).n)/(l.n) Put back in line eq. to find p w = self.pos - self.planePoint si = -planeNormal.dot(w) / ndotu pos = self.pos + si * self.dir p_rad2 = pos[0]**2 + pos[1]**2 # distance from pos projected on Z plane and (0,0,Z); not using numpy norm for speed self.intersect_surf = "filter" if p_rad2 + (pos[2]-Z)**2 > roi**2: # filter missed (walls or frame): distance along filter plane (filter diameter) if (not self.cfg.enableReflections) or (not self.cfg.enableWalls) or (baffles is None): # if walls disabled return None return(None) r_cyl = baffles[0]['rad1'] if p_rad2 size_max_all: size_max_all = size_max self.ax.set_xlim(-size_max_all/2, size_max_all/2) self.ax.set_ylim(-size_max_all/2, size_max_all/2) self.ax.set_zlim(0, size_max_all)# surfaces[0].z*1.2)#400) plt.show() def plot(self, positions): positions = np.array(positions) if positions.shape[0]>1: # if more than one positions self.ax.plot(positions[:,0],positions[:,1],positions[:,2],'.-') plt.show() #--------------------------------------------------------------- def tracer(cfg): ''' Vengono calcolate le seguenti probabilità: det_prob: rivelato lost_prob: perso lateralmente out_prob: perso attraverso il filtro più caldo abs_prob: assorbito da un filtro negl_prob: neglected - low probability residue: 1- le probabilità di cui sopra ''' out = obj() if type(cfg.rand_state) is tuple: np.random.set_state(cfg.rand_state) # if cfg.rand_state is a saved state it is set to start from that point else: np.random.seed(cfg.rand_state) # Set random seed. Use None to randomize between executions rand_state = np.random.get_state() # get initial state for later saving progressUpdateTime = 10 surfaces = list() verso = -1 if cfg.verbose: print('-'*60) time_now = datetime.datetime.today().strftime("%d/%m/%Y-%H:%M") print("{} - Calculating transmissions.".format(time_now)) bb_table = [None]*len(cfg.surfdef) # for each surface will contain photon emission probability as a function of energy Pbb_table = [None]*len(cfg.surfdef) # for each surface will contain the integrated power for normalization for i,s in enumerate(cfg.surfdef): if len(s) == 3: s.append(None) # if baffle is not defined pass a None to the filter class constructor if len(s) > 4: # parameter dictionary present kwargs = s[-1] for kw in kwargs: # check for allowed optional keywords if not kw in ['Tframe','frame_refl_up','frame_refl_down','name','name_frame','name_baffle','flt_comp','flt_transrefl', 'tilt', 'holes', 'skip_emiss']: print(f"ERROR - Surface {i}: '{kw}' is not an allowed keyword") sys.exit(0) else: kwargs = {} if s[0] == 0: # if detector surfaces.append(filter(*s[:4],cfg.wave,np.zeros_like(cfg.wave),cfg.detector_refl*np.ones_like(cfg.wave),0,0,surf_idx=i,cfg=cfg,surfaces=surfaces, **kwargs)) # detector (center in 0,0, const. reflectivity) else: try: composition = kwargs["flt_comp"] if cfg.verbose: print(f'Overriding filter composition: "{composition}"') except(KeyError, IndexError): composition = cfg.filter_composition if composition == "open": if cfg.verbose: print("Filter OPEN") trans = np.ones_like(cfg.wave) refl = np.zeros_like(cfg.wave) else: if kwargs.get('flt_transrefl'): trans, refl = kwargs.get('flt_transrefl') if cfg.verbose: print(f"Using fixed values for surface {i}: trans={trans}, refl={refl}, emiss={1-trans-refl}") trans *= np.ones_like(cfg.wave) refl *= np.ones_like(cfg.wave) elif getattr(cfg,"bypass_transrefl",None): # forced single transmission and reflection values trans, refl = cfg.bypass_transrefl trans *= np.ones_like(cfg.wave) refl *= np.ones_like(cfg.wave) if cfg.verbose: print(f"Using fixed values for surface {i}: trans={trans}, refl={refl}, emiss={1-trans-refl}") else: n_files = getattr(cfg,"n_files",None) #print(cfg.n_files) trans_calc, refl_calc = getTransRefl(cfg.wave[cfg.wave<=cfg.wave_sup],composition,1,1,n_files=n_files) trans = np.ones_like(cfg.wave) * trans_calc[-1] # extend the available transmissivity to higher wavelenghts refl = np.ones_like(cfg.wave) * refl_calc[-1] # extend the available reflectivity to higher wavelenghts trans[cfg.wave<=cfg.wave_sup] = trans_calc refl[cfg.wave<=cfg.wave_sup] = refl_calc verso *= -1 if kwargs.get('tilt', None) is None: kwargs['tilt'] = verso*cfg.tilt surfaces.append(filter(*s[:4],cfg.wave,trans,refl,cfg.tilt_dir,surf_idx=i,cfg=cfg,surfaces=surfaces, **kwargs)) bb_table[i], Pbb_table[i] = surfaces[-1].probab, surfaces[-1].Pbb transRefl_table = [None]*len(surfaces) detected_list = list() node_list = list() M = -1 lost_prob = np.double() out_prob = np.double() det_prob = np.double() abs_prob = np.double() negl_prob = np.double() residual_prob = np.double() direct_hits = 0 direct_frac_tot = np.double() volume_frac = np.double() progress_n = list() M_list = list() converg = list() det_power_arr = list() if cfg.save_dir_array: dir_array = list() map_data = np.zeros((100,100),dtype="float64") traj_info = list() # contains info for each trajectory [pos, dir, prob, wav] global trans_n, refl_wall_n trans_n = np.zeros(10,dtype="float64") refl_wall_n = np.zeros(int(1e6),dtype="float64") refl_wall_max = 0 n_det_traj = 0 det_prob_arr = np.zeros(int(cfg.N)) # array di output, contiene le probabilità di rivelazione per ogni fotone wave_arr = np.zeros(int(cfg.N)) # array di output, contiene le lunghezze d'onda per ogni fotone if cfg.save_nReflections: n_refl = np.zeros(int(cfg.N)) # array di output, contiene le lunghezze d'onda per ogni fotone if cfg.enableTracingPlot: tplot = tracingPlot(surfaces,cfg) startTime = time.time() if cfg.verbose: time_now = datetime.datetime.today().strftime("%d/%m/%Y-%H:%M") print("{} - Starting raytracing.".format(time_now)) progressTime = startTime +1 pht_stack = [photon(cfg=cfg) for ip in range(cfg.stack_size)] idx_max = 0 last_converg_idx = 0 if cfg.enablePlots: f_conv,ax_conv=plt.subplots();plt_convergence=ax_conv.plot(progress_n,det_power_arr)[0];plt.title("Power convergence");plt.grid(True);plt.xlabel("Photons");plt.ylabel("Power to detector [w]"); plt.ticklabel_format(axis='both', style='sci',scilimits=(0,0)) plmap = plot_map(det_r=surfaces[-1].r) for nPhoton in range(int(cfg.N)): try: if time.time() >= progressTime: # progress bar every progressUpdateTime seconds progressTime += progressUpdateTime M_list.append(M) tElaps = time.time() - startTime tRemain = tElaps*(1/(nPhoton +1)*cfg.N-1) ETA = (datetime.datetime.today()+datetime.timedelta(seconds=tRemain)).strftime("%d/%m/%Y-%H:%M") # calculate detected power Pbb = surfaces[cfg.select_surf].Pbb if wave_arr[0] == 0: det_power = np.nan else: gen_energy = const.Planck * const.c*1e9*np.sum(1/wave_arr[wave_arr>0]) gen_time = gen_energy/Pbb det_power = const.Planck * const.c*1e9/gen_time*np.sum(1/wave_arr[wave_arr>0]*det_prob_arr[wave_arr>0]) det_power_arr.append(det_power) if cfg.verbose: print("E:{:9.3e} s,R:{:9.3e}, ETA:{} - det:{:.3e}, lost:{:.3e}, out:{:.3e}, abs:{:.3e}, negl:{:.3e}, res:{:.3e}, P:{:.4e} W, N:{:.3e}, ".format(tElaps, tRemain, ETA, det_prob/nPhoton, lost_prob/nPhoton, out_prob/nPhoton,abs_prob/nPhoton, negl_prob/nPhoton, residual_prob/nPhoton, det_power, nPhoton),end='') progress_n.append(nPhoton) converg.append(det_prob/(nPhoton+1)) if cfg.enablePlots: plt_convergence.set_data(progress_n,det_power_arr) ax_conv.relim() # reset plot limits ax_conv.autoscale_view(True,True,True) # autoscale plot f_conv.canvas.draw() f_conv.canvas.flush_events() plmap.update(map_data) if cfg.enableConvergence and (nPhoton >= cfg.converg_window): for wind in range(len(progress_n)-1,0,-1): if progress_n[wind] <= nPhoton - cfg.converg_window: window = i break conv_val = (max(converg[wind:])-min(converg[wind:]))/np.mean(converg[wind:])*100 # convergence check if conv_val <= cfg.converg_thr: if cfg.verbose: print("\nConvergence stop: {:.1f}%".format(conv_val),flush=True) break else: if cfg.verbose: print("conv: {:.1f}%".format(conv_val),flush=True) else: if cfg.verbose: print("conv: NA",flush=True) if cfg.check_ph_prob: ph_rand_state = np.random.get_state() ph = photon(cfg=cfg) ph.generate(filter=surfaces[cfg.select_surf]) start_pos = ph.pos volume_frac += (1-np.cos(ph.gen_angle)-volume_frac)/(nPhoton+1) #volume_frac = 1 if not cfg.use_const_energy: ph.wav = np.interp(np.random.random(),bb_table[cfg.select_surf],cfg.wave[:-1]) # assign wavelength to photon for sidx,s in enumerate(surfaces): transRefl_table[sidx] = s.get_probab(ph.wav) # get transmission and reflection probabilities if cfg.enableTracingPlot: print("pos0 = {}, dir0 = {}, wav = {}".format(ph.pos0,ph.dir0,ph.wav)) det_surf = len(surfaces)-1 if cfg.save_dir_array: dir_array.append(ph.dir) if not ph.intersect(filter=surfaces[-1],update=False) is None: # check if can reach the detector as direct hit direct_frac = 1.0 for surface in range(1,len(surfaces)-1): # direct detector hit probability direct_frac *= transRefl_table[surface][0] direct_frac *= transRefl_table[-1][2] # absorbtion from detector direct_hits += 1 direct_frac_tot += direct_frac dir = 1 idx = 0 if ph.dir[2] > 0: #z-axis is pointed far from the detector ph.sense = -1 #1: towards the detector, -1: away from the detector else: ph.sense = 1 #1: towards the detector, -1: away from the detector ph.next_surface = cfg.select_surf+1 if ph.sense==1 else cfg.select_surf n_det_traj_ph = 0 lost_prob_ph = np.double(0) out_prob_ph = np.double(0) det_prob_ph = np.double(0) abs_prob_ph = np.double(0) negl_prob_ph = np.double(0) residual_prob_ph = np.double(0) pht_stack[0] = ph transmitted_stack = [0 for ip in range(cfg.stack_size)] reflected_wall_stack = [0 for ip in range(cfg.stack_size)] while idx>=0: curr_idx = idx curr_ph = pht_stack[idx] # pop photon if curr_ph.prob < cfg.prob_thr: # if probability too low, discharge idx -= 1 # decrement stack index negl_prob_ph += curr_ph.prob # increment lost probability continue # continue loop surface = curr_ph.next_surface start_pos = curr_ph.pos if curr_ph.sense == 1: # toward detector baffles = surfaces[surface-1].baffles # filter hotter-than-current baffle else: # sense = -1: toward 300 K baffles = surfaces[surface].baffles # current filter baffle pos = curr_ph.intersect(filter=surfaces[surface], baffles=baffles) if pos is None: # verify surface interaction idx -= 1 # decrement stack index lost_prob_ph += curr_ph.prob # increment lost probability continue # continue loop if cfg.enableTracingPlot: print("idx:{}; Photon probability: {:.3e}; {:6s}; pos: {}".format(idx, curr_ph.prob,curr_ph.intersect_surf,pos)) tplot.plot([start_pos,pos]) plt.draw() plt.pause(1) if curr_ph.intersect_surf == "frame": # Intersection with the FRAME: reflected_wall_stack[idx] += 1 # increase number of reflections on the walls/frames if reflected_wall_stack[idx] > refl_wall_max: refl_wall_max = reflected_wall_stack[idx] if np.random.rand() < cfg.ratio_frame: curr_ph.prob /= cfg.ratio_frame # re-introduce photon increasing its probability if curr_ph.sense == 1: #1: towards the detector, -1: away from the detector frame_refl = surfaces[surface].get_refl('frame_refl_up',wave=curr_ph.wav) else: frame_refl = surfaces[surface].get_refl('frame_refl_down',wave=curr_ph.wav) lost_prob_ph += curr_ph.prob*(1-frame_refl) # increment lost probability curr_ph.prob *= frame_refl curr_ph.reflect_dir(surfaces[surface].plane) # calculate reflection direction curr_ph.sense *= -1 curr_ph.next_surface = surface + curr_ph.sense continue # skip absorbtion and transmission checks else: idx -= 1 # decrement stack index continue # ignore photon and continue loop if curr_ph.intersect_surf == "walls": reflected_wall_stack[idx] += 1 # increase number of reflections on the walls/frames if reflected_wall_stack[idx] > refl_wall_max: refl_wall_max = reflected_wall_stack[idx] if np.random.rand() < cfg.ratio_walls: # Intersection with the WALLS: curr_ph.prob /= cfg.ratio_walls # re-introduce photon increasing its probability if curr_ph.sense == 1: # toward detector. Use baffle of the hotter filter refl_walls = surfaces[surface-1].get_refl('baffle_refl',wave=curr_ph.wav) else: # sense = -1: toward 300 K refl_walls = surfaces[surface].get_refl('baffle_refl',wave=curr_ph.wav) lost_prob_ph += curr_ph.prob*(1-refl_walls) # increment lost probability curr_ph.prob *= refl_walls curr_ph.reflect_dir(curr_ph.planePoint) # calculate reflection direction continue # skip absorbtion and transmission checks else: idx -= 1 # decrement stack index continue # ignore photon and continue loop # interaction: [pTrans,pRefl,pAbs] = transRefl_table[surface] # get transmission, reflection and absorption probability # trasmission-absorption: if surface == det_surf: # DETECTOR - if interaction surface is detector: n_det_traj_ph += 1 # increment detected trajectories det_prob_ph += curr_ph.prob * (pAbs+pTrans) #* (-curr_ph.dir[2]) # increment detected probability if cfg.enableTracingPlot: print(f"Detected with prob. {curr_ph.prob * (pAbs+pTrans)}") if cfg.save_traj: traj_info.append([curr_ph.pos, curr_ph.dir, curr_ph.prob, curr_ph.wav, curr_ph.prob * (pAbs+pTrans)]) mapX = int(np.interp(curr_ph.pos[0],[-surfaces[surface].r,surfaces[surface].r],[0,map_data.shape[0]])) mapY = int(np.interp(curr_ph.pos[1],[-surfaces[surface].r,surfaces[surface].r],[0,map_data.shape[1]])) map_data[mapX,mapY] += curr_ph.prob * (pAbs+pTrans) out.map = map_data if reflected_wall_stack[idx] < refl_wall_n.size: idx_n = int(reflected_wall_stack[idx]) else: idx_n = refl_wall_n.size-1 refl_wall_n[idx_n] += curr_ph.prob * (pAbs+pTrans) trans_n[int(curr_ph.ntrans)] += curr_ph.prob * (pAbs+pTrans) else: flt_open = False try: # Check for OPEN filter if cfg.surfdef[surface][4]["flt_comp"].lower() == "open": flt_open = True curr_ph.next_surface = surface + curr_ph.sense continue # continue loop if filter is open except: pass if surfaces[surface].holes: # check for HOLES hole_hit = False for hole in surfaces[surface].holes: dist_hole = (curr_ph.pos[0] - hole[0])**2 + (curr_ph.pos[1]-hole[1])**2 if dist_hole < hole[2]**2: hole_hit = True if cfg.enableTracingPlot: print(f"Hole hit") break if hole_hit: # continue main loop if photon hits a hole curr_ph.next_surface = surface + curr_ph.sense continue if curr_ph.sense == -1: # if transmitted trough hotter filter is OUT out_prob_ph += curr_ph.prob * pTrans # increment out probability abs_prob_ph += curr_ph.prob * pAbs # increment absorbed probability else: # if not detector: passing trough filter. increment absorbed probability abs_prob_ph += curr_ph.prob * pAbs # increment absorbed probability if not cfg.enableReflections: curr_ph.prob *= pTrans curr_ph.next_surface = surface + curr_ph.sense lost_prob_ph += pRefl # increment lost probability continue # continue loop idx += 1 # increment stack index if idx > idx_max: idx_max = idx tr_ph = pht_stack[idx] # put transmitted photon in a new position of the stack tr_ph.ntrans = curr_ph.ntrans + 1 # update its number of filter crossings tr_ph.pos = curr_ph.pos tr_ph.dir = curr_ph.dir tr_ph.wav = curr_ph.wav tr_ph.prob = curr_ph.prob * pTrans tr_ph.sense = curr_ph.sense tr_ph.next_surface = surface + curr_ph.sense # reflection: if cfg.save_nReflections: n_refl[nPhoton] += 1 curr_ph.reflect_dir(surfaces[surface].plane) # calculate reflection direction and update reflected photon curr_ph.prob *= pRefl curr_ph.sense = -curr_ph.sense # change sense to reflected photon curr_ph.next_surface = surface + curr_ph.sense residual_prob_ph = 1-(det_prob_ph+lost_prob_ph+out_prob_ph+abs_prob_ph+negl_prob_ph) det_prob += det_prob_ph lost_prob += lost_prob_ph out_prob += out_prob_ph abs_prob += abs_prob_ph negl_prob += negl_prob_ph residual_prob += residual_prob_ph wave_arr[nPhoton] = ph.wav det_prob_arr[nPhoton] = det_prob_ph n_det_traj += n_det_traj_ph if direct_frac_tot: M = det_prob/direct_frac_tot if cfg.check_ph_prob: if det_prob_ph >= 1e-9: out.stop = True print("Prob: {}".format(det_prob_ph)) print("Det. traj: {}".format(n_det_traj_ph)) print("Start pos: {}".format(ph.pos0)) print("Start dir: {}".format(ph.dir0)) print("Wavelength: {}".format(ph.wav)) import pickle with open("rand_state","wb") as f: pickle.dump(ph_rand_state,f) break else: out.stop = False except(KeyboardInterrupt): print("CTRL-C received. Stopping simulation.") break #------------- #global Pbb, weights # for debug only lost_prob /= nPhoton+1 out_prob /= nPhoton+1 det_prob /= nPhoton+1 abs_prob /= nPhoton+1 if cfg.verbose: print("\nMax size reached for stack: {}".format(idx_max)) print("\nN: {:.2e}".format(cfg.N)) print("Simulated photons: {:.2e}".format(nPhoton+1)) print("Elapsed time: {:.3e} s".format(time.time()-startTime)) print("Solid angle frac.: {:.3e}".format(volume_frac)) print("Tilt: {}".format(cfg.tilt)) print("Detector reflect.: {}".format(cfg.detector_refl)) print("Filter compos.: " + repr(cfg.filter_composition)) print("Surfaces: " + repr(cfg.surfdef)) if getattr(cfg,'bypass_transrefl',None): print(f"Bypass transmissivity/reflectivity: {cfg.bypass_transrefl}") print(f"Default baffle/frame reflectivity: {cfg.refl_walls}/{cfg.refl_frame}") print("Detected trajectories: {}".format(n_det_traj)) print("Max number of reflections on walls/frames: {}".format(refl_wall_max)) print("Direct hits: {}".format(direct_hits)) print("Direct hits ratio: {:.3e}".format(direct_hits/(nPhoton+1))) print("Direct hits ratio norm.: {:.3e}".format(direct_hits/(nPhoton+1)*volume_frac)) print("M: {}".format(M)) print("Prob. - det: {:.3e}, lost: {:.3e}, out: {:.3e}, abs: {:.3e}, residue: {:.3e}".format(det_prob,lost_prob,out_prob,abs_prob,1-(det_prob+lost_prob+out_prob+abs_prob))) if cfg.enablePlots: plt_convergence.set_data(progress_n,det_power_arr) ax_conv.relim() # reset plot limits ax_conv.autoscale_view(True,True,True) # autoscale plot plt.show() # Build output object out.rand_state = rand_state#np.random.get_state() out.t_elapsed = time.time()-startTime out.cfg = cfg if cfg.save_prob: out.wave_arr = wave_arr out.det_prob = det_prob out.det_prob_arr = det_prob_arr out.progress_n = progress_n out.converg = converg out.Nph = nPhoton+1 # number of simulated photons out.M_list = M_list out.surfaces = surfaces out.volume_frac = volume_frac out.map = map_data if cfg.save_dir_array: out.dir_array = dir_array if cfg.save_traj: out.traj_info = traj_info if cfg.save_nReflections: out.n_refl = n_refl return(out) class plot_map(): def __init__(self, map_data=None, det_r=None): self.f_map = plt.figure(figsize=(8, 8)) gs = matplotlib.gridspec.GridSpec(4, 1) ax_map = plt.subplot(gs[0:3,0]) ax_profile = plt.subplot(gs[3,0], sharex=ax_map) my_cmap = copy.copy(matplotlib.cm.get_cmap('inferno')) # copy the default cmap my_cmap.set_bad((0,0,0)) # set 'bad' pixel color to black if map_data: m = map_data mn =(m - np.min(m)) / (np.max(m) - np.min(m)) # normalized map else: mn = np.zeros((100,100)) if det_r is None: det_r = 0.5 im = ax_map.imshow(mn,cmap=my_cmap,extent=[-det_r,det_r,-det_r,det_r], origin='lower',norm=matplotlib.colors.LogNorm(vmin=1e-6, vmax=1)) # use normalized log colorscale. Set minimum to 1e-5 plt.colorbar(im,ax=ax_map) circle1 = plt.Circle((0, 0), det_r, facecolor='None', edgecolor='r', lw=1) ax_map.add_artist(circle1) self.prof = ax_profile.plot(np.linspace(-det_r,det_r,int(mn.shape[1])), mn[:,int(mn.shape[1]/2)]) pos = ax_map.get_position() pos2 = ax_profile.get_position() ax_profile.set_position([pos.x0,pos2.y0,pos.width,pos2.height]) ax_profile.set_ylim([0,1]) plt.grid() self.im = im plt.show() self.f_map.canvas.draw() self.f_map.canvas.flush_events() def update(self, map_data): m = map_data mn =(m - np.min(m)) / (np.max(m) - np.min(m)) # normalized map self.im.set_data(mn) self.prof[0].set_ydata(mn[:,int(mn.shape[1]/2)]) self.f_map.canvas.draw() self.f_map.canvas.flush_events() def analyze(data, bins_gen=np.linspace(0,2e6,10000),calc_spectrum_det=True): cfg = data.cfg surfaces = data.surfaces wave_arr = data.wave_arr det_prob_arr = data.det_prob_arr volume_frac = data.volume_frac out = obj() if cfg.enablePlots: calc_spectrum_det=True Pbb = surfaces[cfg.select_surf].Pbb gen_energy = np.sum(const.Planck * const.c*1e9/wave_arr[wave_arr>0]) gen_time = gen_energy/Pbb det_power = np.sum(const.Planck * const.c*1e9/gen_time/wave_arr[wave_arr>0]*det_prob_arr[wave_arr>0]) if calc_spectrum_det: wave_gen = (bins_gen[:-1]+bins_gen[1:])/2 # wavelenghts are computed as averages of the bin edges hni = const.Planck * const.c/(wave_gen/1e9) # energy of the photons: J*m/(nm/1e9) -> J hist_det = np.histogram(wave_arr[wave_arr>0], bins=bins_gen, weights=det_prob_arr[wave_arr>0])[0]/gen_time # histogram of detected photons, normalized spectrum_det = hist_det/np.diff(bins_gen)* hni / volume_frac # convert into a spectrum [J/nm] if cfg.enablePlots: # calculate histograms hist_gen, _ = np.histogram(wave_arr[wave_arr>0], bins=bins_gen) # histogram of generated photons spectrum_gen_ph = hist_gen/np.diff(bins_gen) # convert into a flux spectrum [photons/nm] spectrum_gen = spectrum_gen_ph * hni / volume_frac # convert into a energy spectrum [J/nm], consider the photons that would have been generated on the full hemisphere #gen_time = scipy.integrate.simps(spectrum_gen,wave_gen)/Pbb # time normalization relative to bb radiation [J/W]=[s], correspond to the time interval during which the photons would have been generated spectrum_gen_norm = spectrum_gen/gen_time # [W/nm] = [J/s/nm] # PLOT generated PSD plt.figure();plt.plot(wave_gen*1e-3,spectrum_gen_norm) # plot generated PSD in [W/nm] vs um plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) plt.title("Generated PSD") plt.xlabel("Wavelenght [um]"); plt.ylabel("PSD [W/nm]"); plt.xlim(0,250);plt.grid() # PLOT detected PSD plt.figure(); plt.plot(wave_gen*1e-3,spectrum_det) plt.title("Detected PSD") plt.xlabel("Wavelenght [um]") plt.ylabel("Detected PSD [W/nm]") plt.xlim(0,250);plt.grid() plt.show() out.wave_gen = wave_gen out.spectrum_gen_norm = spectrum_gen_norm out.spectrum_det = spectrum_det # PHOTON SHOT NOISE NEP2 = 2*np.sum(const.Planck**2*det_prob_arr[wave_arr>0]/gen_time*(const.c*1e9/wave_arr[wave_arr>0])**2/cfg.pixel_num) de_fwhm = math.sqrt(NEP2*cfg.det_tau)*2.35/const.eV # eV (sqrt(sum(NEP^2)) is a quadrature sum. ) if cfg.verbose: print("Generated power: {:.6e} W".format(Pbb)) print("Detected power: {:.6e} W".format(det_power)) print("NEP: {:.4e} W/Hz^0.5/pixel\nΔE FWHW: {:.4e} eV".format(math.sqrt(NEP2),de_fwhm)) out.gen_time = gen_time out.gen_power = Pbb out.det_power = det_power out.NEP2 = NEP2 out.de_fwhm = de_fwhm if calc_spectrum_det: out.spectrum_det = spectrum_det out.wave_gen = wave_gen return(out) if __name__ == '__main__': import matplotlib import matplotlib.pyplot as plt from matplotlib import rc from matplotlib.patches import Circle, PathPatch from matplotlib.text import TextPath from matplotlib.transforms import Affine2D # This import registers the 3D projection, but is otherwise unused. from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import import mpl_toolkits.mplot3d.art3d as art3d __spec__ = None cfg = init() if cfg.redirectStdio: print("Switching to file output.") log_file = os.path.join(cfg.outputDir,cfg.redirectStdioFile) sys.stdout = open(log_file, 'a+') results = tracer(cfg) o = analyze(results)