DISPATCH
particle_utils.py
1 # -*- coding: utf-8 -*-
2 """
3  Particles object, with methods:
4 
5  part=particles(it=2,read=True)
6  part.read_shell(r=100,dr=5)
7 
8 Created on Wed Nov 29 09:25:50 2017
9 
10 @author: Aake
11 """
12 from __future__ import print_function
13 from scipy.io import FortranFile
14 import sys
15 if sys.version_info[0]!=2:
16  from importlib import reload
17 from numpy import size,sqrt,product,mean,mod
18 import dispatch_utils; reload(dispatch_utils)
19 import dispatch_data; reload(dispatch_data)
20 from dispatch_utils import patch_files
21 from dispatch_data import Patch
22 from scaling import cgs, scaling
23 from time import time
24 
25 class particles:
26  def __init__(self,it=0,run='',data='../data',read=0,r=0,dr=1,gas_mass=0,verbose=0):
27  ''' Initialize, by reading in the list of patch files '''
28  start=time()
29  self.files=patch_files(it,run,data)
30  if gas_mass>0:
31  self.gas_mass=gas_mass
32  if verbose:
33  print('{} files ({:.3f} sec)'.format(size(self.files),time()-start))
34  if read:
35  self.read_shell(r=r,dr=dr,verbose=verbose)
36 
37  def get_gas_mass(self,verbose=0):
38  ''' Sum up the gas mass '''
39  if verbose:
40  start=time()
41  print("computing gas mass ...")
42  mass=0.0
43  n=0
44  for file in self.files:
45  p=Patch(file)
46  l=p.li
47  u=p.ui
48  mass+=product(p.size)*mean(p.data[l[0]:u[0],l[1]:u[1],l[2]:u[2],0])
49  if verbose>1:
50  n+=1
51  if mod(n,200)==0: print(n,mass)
52  self.gas_mass=mass
53  if verbose:
54  print('mass: {:.4f} ({:.2f} sec)'.format(mass,time()-start))
55 
56  def read_all(self,verbose=0):
57  ''' Trigger reading of all patches (may take a long time) '''
58  self.read_shell(r=0,verbose=verbose)
59 
60  def read_patch(self,file,verbose=0):
61  ''' Read particles in a patch '''
62  name=file+'.peb'
63  f=FortranFile(name,'rb')
64  fmt=f.readInts()
65  if verbose>2:
66  print('ioformat:',fmt[0])
67  dim=f.readInts()
68  dim=dim[0:2]
69  a=f.readReals()
70  a=a.reshape(dim[1],dim[0]).transpose()
71  idv = f.readInts()
72  idx = idv[0:dim[0]]
73  f.close()
74  dict={}
75  for i in range(size(idx)):
76  id=idx[i]
77  dict[id]={'p':a[i,0:3],'v':a[i,0:3],'t':a[i,6],'s':a[i,7],'w':a[i,8]}
78  return idx,dict
79 
80  def read_shell(self,r=100,dr=1,dust2gas=0.01,save=True,verbose=0):
81  ''' Read particles in a shell, by default saving them '''
82  dict={}
83  n=0
84  npatch=0
85  rate=0.0
86  start=time()
87  t0=start
88  if r==0:
89  r2min=0.0
90  r2max=1e30
91  else:
92  r2min=(r-dr)**2
93  r2max=(r+dr)**2
94 
95  # open each patch file and check if relevant
96  for file in self.files:
97  p=Patch(file)
98  rc=p.corner_radii()
99 
100  # If relevant, open the .peb file, and read the data
101  if rc.min() < (r+dr) and rc.max() > (r-dr) or r==0:
102  npatch+=1
103  idx,dd=self.read_file(file)
104 
105  # If in the shell, sum up rate contribution, and add particle to dict
106  if (r>0):
107  for i in range(size(idx)):
108  id=idx[i]
109  d=dd[id]
110  p=d['p']
111  p2=sum(p**2)
112  if p2>r2min and p2<r2max:
113  v=d['v']
114  w=d['w']
115  rate -= sum(p*v)/sqrt(p2)*w
116  n=n+1
117  if save:
118  dict[id]=d
119  else:
120  for i in range(size(idx)):
121  id=idx[i]
122  dict[id]=dd[id]
123  if verbose>1:
124  now=time()
125  print('{:.3f} sec'.format(now-start))
126  start=now
127  if r==0:
128  self.particles=dict
129  print('{:.3f} sec'.format(time()-start))
130  return
131  if not hasattr(self,'gas_mass'):
132  self.get_gas_mass(verbose=verbose)
133  rate=rate/(2.0*dr)*self.gas_mass
134  self.rate=rate
135  if save:
136  self.particles=dict
137  if verbose>1:
138  print('{} self.particles saved'.format(size(dict.keys())))
139  if verbose:
140  units=scaling()
141  rate1=dust2gas*rate*cgs.yr/units.t*units.m/cgs.m_earth
142  s='rate:{:9.2e} ={:9.2e} M_E/yr'.format(rate,rate1)
143  s=s+', based on {} particles from {} patches'.format(n,npatch)
144  print(s+' ({:.1f} sec)'.format(time()-t0))
145 
146  def read_plane(self,x=None,y=None,z=None,ds=1,wmin=1e-9,dust2gas=0.01,save=True,verbose=0):
147  ''' Read particles in a plane, by default saving them '''
148  dict={}
149  start=time()
150 
151  def read_patch_plane(j,s,ds,c):
152  ''' read particles in patch belonging to a plane '''
153  if c[1] > (s-ds) and c[0] < (s+ds):
154  if verbose>1:
155  print(file)
156  idx,dd=self.read_patch(file)
157  for i in range(size(idx)):
158  id=idx[i]
159  d=dd[id]
160  p=d['p']
161  if d['w']>wmin and p[j]>(s-ds) and p[j]<(s+ds):
162  dict[id]=d
163  return 1
164  else:
165  return 0
166  # open each patch file and check if relevant
167  npatch=0
168  for file in self.files:
169  p=Patch(file)
170  c=p.active_corners
171  if type(x) != type(None):
172  npatch+=read_patch_plane(0,x,ds,c[:,0])
173  elif type(y) != type(None):
174  npatch+=read_patch_plane(1,y,ds,c[:,1])
175  elif type(z) != type(None):
176  npatch+=read_patch_plane(2,z,ds,c[:,2])
177  if verbose>1:
178  now=time()
179  print('{}: {:.3f} sec'.format(file,now-start))
180  start=now
181  self.particles=dict
182  if verbose>0:
183  print('{} particles from {} patches ({:.2f} sec)'.format(size(dict.keys()),npatch,time()-start))
184 
185  def shell(self,r=100,dr=1,dust2gas=0.01,verbose=0):
186  ''' evaluate rate of particles in a shell '''
187  n=0
188  rate=0.0
189  t0=time()
190  if r==0:
191  r2min=0.0
192  r2max=1e30
193  else:
194  r2min=(r-dr)**2
195  r2max=(r+dr)**2
196  for id in (self.particles).keys():
197  part=self.particles[id]
198  p=part['p']
199  p2=sum(p**2)
200  if p2>r2min and p2<r2max:
201  v=part['v']
202  w=part['w']
203  rate -= sum(p*v)/sqrt(p2)*w
204  n=n+1
205  if not hasattr(self,'gas_mass'):
206  self.get_gas_mass(verbose=verbose)
207  rate=rate/(2.*dr)*self.gas_mass
208  self.rate=rate
209  if verbose:
210  units=scaling()
211  rate1=dust2gas*rate*cgs.yr/units.t*units.m/cgs.m_earth
212  s='rate:{:9.2e} ={:9.2e} M_E/yr, based on {} particles.'.format(rate,rate1,n)
213  s=s+' Time used: {:.1f} sec.'.format(time()-t0)
214  print(s)
def read_shell(self, r=100, dr=1, dust2gas=0.01, save=True, verbose=0)
def get_gas_mass(self, verbose=0)
def read_plane(self, x=None, y=None, z=None, ds=1, wmin=1e-9, dust2gas=0.01, save=True, verbose=0)
def shell(self, r=100, dr=1, dust2gas=0.01, verbose=0)
def read_all(self, verbose=0)
def read_patch(self, file, verbose=0)
def __init__(self, it=0, run='', data='../data', read=0, r=0, dr=1, gas_mass=0, verbose=0)