Source Code for Module chianti.ionrec

  1  #import types 
  2  import numpy as np 
  3  import chianti.core as ch 
  4  from chianti import util 
  5 -class ioneq(ch.ion): 
  6 -    def __init__(self,z, temperature, verbose=False): 
  7          '''calculates the ionization equilibrium for element z at a single temperatures''' 
  8  #        self.Defaults=defaults 
  9          self.Z=z 
 10          self.Temperature=temperature 
 11          # 
 12 -    def getIoneq(self): 
 13          ''' get the ionization equilibrium for this ion at the specified temperature 
 14          getIoneqFromEvolver does the same but is a bit less accurate''' 
 15          z = self.Z 
 16          temperature = self.Temperature 
 17          ionList=[] 
 18          chIons=[] 
 19          for ion in range(1, z+2): 
 20              ions=util.zion2name(z, ion) 
 21              ionList.append(ions) 
 22  #            print z, ion, ions 
 23              atom=ch.ion(ions, temperature=temperature, setup=0) 
 24              atom.ionizRate() 
 25              atom.recombRate() 
 26              chIons.append(atom) 
 27          ioneq=np.zeros((z+1), 'float64') 
 28          factor=[] 
 29          for anIon in chIons: 
 30              if type(anIon.IonizRate) != type(None) and type(anIon.RecombRate) != type(None): 
 31                  rat=anIon.IonizRate['rate']/anIon.RecombRate['rate'] 
 32                  factor.append(rat**2 + rat**(-2)) 
 33              else: 
 34                  factor.append(0.) 
 35          factor[0]=max(factor) 
 36          factor[-1]=max(factor) 
 37          ionmax=factor.index(min(factor)) 
 38          ioneq[ionmax]=1. 
 39          # 
 40          for iz in range(ionmax+1, z+1): 
 41              ionrate=chIons[iz-1].IonizRate['rate'] 
 42              recrate=chIons[iz].RecombRate['rate'] 
 43              ioneq[iz]=ionrate*ioneq[iz-1]/recrate 
 44  #            print ' iz, ioneq = ',iz,ioneq[iz] 
 45          # 
 46          for iz in range(ionmax-1, -1, -1): 
 47              ionrate=chIons[iz].IonizRate['rate'] 
 48              recrate=chIons[iz+1].RecombRate['rate'] 
 49              ioneq[iz]=recrate*ioneq[iz+1]/ionrate 
 50  #            print ' iz, ioneq = ',iz,ioneq[iz] 
 51          ionsum=ioneq.sum() 
 52  #        print ' ionsum = ', ionsum 
 53          ioneq=ioneq/ionsum 
 54          self.Ioneq=ioneq 
 55          # 
 56          # -------------------------------------- 
 57          # 
 58 -    def getIoneqFromStatic(self): 
 59          '''gets the ionization equilibrium for the specified ion at the specified temperature 
 60          not as accurate as getIoneq since this methods uses a matrix inversion which can result in 
 61          nonphysical values (usually small) that are zero'd out ''' 
 62          z = self.Z 
 63          self.getEvolver(self.Temperature) 
 64          b = np.zeros((z+1),'float64') 
 65          b[0] = 1. 
 66          static = self.Static 
 67          ioneq = np.linalg.solve(evolver,b) 
 68          ioneq = np.where(ioneq >= 0., ioneq, 0.) 
 69          self.Ioneq = ioneq 
 70          # 
 71          # -------------------------------------- 
 72          # 
 73   
 74 -    def getStatic(self): 
 75          ''' calculate the static matrix for the specified temperature''' 
 76          temperature = self.Temperature 
 77          z = self.Z 
 78          ionList=[] 
 79          chIons=[] 
 80          evolver = np.zeros((z+1,z+1),'float64') 
 81          for ion in range(1, z+2): 
 82              ions=util.zion2name(z, ion) 
 83              ionList.append(ions) 
 84  #            print z, ion, ions 
 85              atom=ch.ion(ions, temperature=temperature) 
 86              atom.ionizRate() 
 87              atom.recombRate() 
 88              chIons.append(atom) 
 89          for stage in range(0,z): 
 90              atom = chIons[stage] 
 91  #            if type(atom.IonizRate) != types.NoneType: 
 92              evolver[stage,stage] -= atom.IonizRate['rate'] 
 93              atom = chIons[stage+1] 
 94              evolver[stage,stage+1] += atom.RecombRate['rate'] 
 95  #                evolver[stage-2,stage] += atom.IonizRate['rate'] 
 96          for stage in range(1,z+1): 
 97              atom = chIons[stage] 
 98  #            if type(atom.RecombRate) != types.NoneType: 
 99              evolver[stage,stage] -= atom.RecombRate['rate'] 
100              atom = chIons[stage-1] 
101              evolver[stage,stage-1] += atom.IonizRate['rate'] 
102  #                evolver[stage-1,stage-2] += atom.RecombRate['rate'] 
103          # include normalization of ionization stages 
104          for stage in range(0,z+1): 
105              evolver[0,stage] = 1. 
106          self.Static = evolver 
107          # 
108          # ----------------------------------------- 
109          # 
110 -    def getEvolver(self,temperature,density,time): 
111          ''' calculate the evolver matrix for the new temperature''' 
112          self.NewTemperature = temperature 
113          z = self.Z 
114          ionList=[] 
115          chIons=[] 
116          evolver = np.zeros((z+1,z+1),'float64') 
117          for ion in range(1, z+2): 
118              ions=util.zion2name(z, ion) 
119              ionList.append(ions) 
120  #            print z, ion, ions 
121              atom=ch.ion(ions, temperature=temperature, setup=0) 
122              atom.ionizRate() 
123              atom.recombRate() 
124              chIons.append(atom) 
125          for stage in range(0,z): 
126  #            atom = chIons[stage] 
127  #            if type(atom.IonizRate) != types.NoneType: 
128              evolver[stage,stage] -= chIons[stage].IonizRate['rate'] 
129  #            atom = chIons[stage+1] 
130              evolver[stage,stage+1] += chIons[stage+1].RecombRate['rate'] 
131  #                evolver[stage-2,stage] += atom.IonizRate['rate'] 
132          for stage in range(1,z+1): 
133  #            atom = chIons[stage] 
134  #            if type(atom.RecombRate) != types.NoneType: 
135              evolver[stage,stage] -= chIons[stage].RecombRate['rate'] 
136  #            atom = chIons[stage-1] 
137              evolver[stage,stage-1] += chIons[stage-1].IonizRate['rate'] 
138  #                evolver[stage-1,stage-2] += atom.RecombRate['rate'] 
139          # include normalization of ionization stages 
140  #        for stage in range(0,z+1): 
141  #            evolver[0,stage] = 1. 
142          # this is the first order Euler solution 
143  #        implicitEvolver = np.invert(np.identity(z+1,'float64') - evolver) 
144          implicitEvolver = np.dual.inv(np.identity(z+1,'float64') - density*time*evolver) 
145          self.Evolver = implicitEvolver 
146  #        imvolver = np.asmatrix(np.identity(z+1) - evolver) 
147  #        self.Evolver = np.dual.inv(imvolver) 
148  #        self.Evolver = evolver 
149          # 
150          # ----------------------------------------- 
151          # 
152 -    def evolve(self,temperature, density, haveEvolver = 0): 
153          ''' to evolve the plasma to the final temperature''' 
154          z = self.Z 
155          self.finalTemperature=temperature 
156          if not haveEvolver: 
157              self.getEvolver(temperature,density,1.) 
158  #        diff = np.abs(density*np.dot(self.Evolver,self.Ioneq)/self.Ioneq) 
159  #        diff = np.abs(self.Evolver).max() 
160  #        print ' diffs = ',diff 
161          maxDiff = 1. 
162  #        print ' max time  = ',1./diff 
163  #        dt = (maxDiff/diff.max()) 
164  #        diff = np.abs((np.ones(z+1,'float64') - self.Evolver.diagonal()).max()) 
165          new = np.dot(self.Evolver,self.Ioneq) 
166          diff = np.where(self.Ioneq > 0.,np.abs(new - self.Ioneq)/self.Ioneq,100.) 
167          print((' diff = %10.2e'%(diff))) 
168          dt = maxDiff/diff.min() 
169          print((' dt = %12.2e'%(dt))) 
170          self.getEvolver(temperature,density,dt) 
171          niter = int(10.*maxDiff/diff.min())+1 
172          print((' niter = %i'%(niter))) 
173          state = np.zeros((niter+1,z+1),'float64') 
174          stateTot = np.zeros(niter+1,'float64') 
175          stateDiff = np.zeros((niter,z+1,z+1),'float64') 
176          change = np.zeros((niter+1),'float64') 
177          deltaChange = np.zeros((niter),'float64') 
178          deltaChange[0] = 1. 
179          state[0] = self.Ioneq 
180          stateTot[0] = self.Ioneq.sum() 
181          previousChange = 1. 
182          dChange = 10. 
183          iter = 0 
184  #        for iter in range(niter): 
185          while (dChange > 1.e-7) and (iter < niter): 
186              present = state[iter] 
187  #            diff = np.dot(self.Evolver,present) 
188  #            imdiff = np.asmatrix(np.identity(self.Z+1,dtype='float64') - diff) 
189  #            evolver = np.dual.inv(imdiff) 
190  #            diff = diff - diff.sum() 
191  #            new = np.dual.solve(imdiff,present) 
192              new = np.dot(self.Evolver,present) 
193              thisdiff = np.where(present > 0.,np.abs(new - present)/present,0.) 
194  #            print iter, change.sum, change 
195              state[iter+1] = new 
196              stateTot[iter+1] = new.sum() 
197  #            stateDiff[iter] = imdiff 
198              change[iter] = thisdiff.sum() 
199              deltaChange[iter] = np.abs(change[iter] - previousChange) 
200              dChange = deltaChange[iter] 
201              previousChange = change[iter] 
202              iter+=1 
203          print((' # of iterations = %i'%(iter))) 
204          self.Temperature = self.finalTemperature 
205          self.getIoneq() 
206          state[iter] = self.Ioneq 
207          self.State = state[:iter].transpose() 
208          self.StateTot = stateTot[:iter] 
209          self.StateDiff = stateDiff[:iter] 
210          self.Change = change[:iter] 
211          self.DeltaChange = deltaChange[:iter] 
212          self.Iter = iter 
213   
214   
215   
216   
217   
218   
219   
220  #    def plot(self, ions=None, xr=None, yr=None, oplot=False, label=True, title=True,  bw=False): 
221  #        ''' self.plot(xr=None, yr=None, oplot=False) 
222  #        ions = sequence of ions to be plotted, in spectroscopic notation 
223  #        xr = temperature range, yr = ion fraction range 
224  #        for overplotting: 
225  #        oplot="ioneqfilename" such as 'mazzotta' 
226  #        or oplot=True and a widget will come up so that a file can be selected''' 
227  #        if bw: 
228  #            linestyle=['k-','k--', 'k-.', 'k:'] 
229  #        else: 
230  #            linestyle=['b-','r--', 'g-.', 'm:'] 
231  #        # 
232  #        if type(ions) == NoneType: 
233  #            ions=range(1, self.Z+2) 
234  #        elif min(ions) < 1 or max(ions) > self.Z+1: 
235  #            ions=range(1, self.Z+2)  #  spectroscopic notation 
236  #        if type(xr) == NoneType: 
237  #            xr=[self.Temperature.min(), self.Temperature.max()] 
238  #        if type(yr) == NoneType: 
239  #            yr=[0.01, 1.1] 
240  #        xyr=list(xr) 
241  #        xyr.extend(list(yr)) 
242  #        # 
243  #        iz=ions[0] 
244  #        pl.loglog(self.Temperature, self.Ioneq[iz-1]) 
245  #        if label: 
246  #            idx=self.Ioneq[iz-1] == self.Ioneq[iz-1].max() 
247  #            if idx.sum() > 1: 
248  #                jdx=np.arange(len(idx)) 
249  #                idx=jdx[idx].max() 
250  #            ann=Ionstage[iz-1] 
251  #            pl.annotate(ann, [self.Temperature[idx], 0.7*self.Ioneq[iz-1, idx]], ha='center') 
252  #        for iz in ions[1:]: 
253  #            pl.plot(self.Temperature, self.Ioneq[iz-1], linestyle[0]) 
254  #            if label: 
255  #                idx=self.Ioneq[iz-1] == self.Ioneq[iz-1].max() 
256  #                if idx.sum() > 1: 
257  #                    jdx=np.arange(len(idx)) 
258  #                    idx=jdx[idx].mean() 
259  #                ann=Ionstage[iz-1] 
260  #                pl.annotate(ann, [self.Temperature[idx], 0.7*self.Ioneq[iz-1, idx]], ha='center') 
261  #        pl.xlabel('Temperature (K)') 
262  #        pl.ylabel('Ion Fraction') 
263  #        atitle='Chianti Ionization Equilibrium for '+El[self.Z-1].capitalize() 
264  #        # 
265  #        if oplot != False: 
266  #            if type(oplot) == BooleanType: 
267  #                result=self.ioneqRead(ioneqname='',default=False) 
268  #                if result != False: 
269  #                    atitle+='  & '+result['ioneqname'].replace('.ioneq', '') 
270  #                    atitle+=' '+linestyle[0] 
271  #                    for iz in ions: 
272  #                        pl.plot(self.IoneqTemperature, self.IoneqAll[self.Z-1, iz-1],linestyle[0], linestyle[1]) 
273  #            elif type(oplot) == StringType: 
274  #                atitle+='  & '+oplot+' '+linestyle[0] 
275  #                atitle+=' '+linestyle[0] 
276  #                result=self.ioneqRead(ioneqname=oplot,default=False) 
277  #                if result != False: 
278  #                    for iz in ions: 
279  #                        pl.plot(self.IoneqTemperature, self.IoneqAll[self.Z-1, iz-1],linestyle[0], linestyle[1]) 
280  #            elif type(oplot) == ListType: 
281  #                for iplot in range(len(oplot)): 
282  #                    result=self.ioneqRead(ioneqname=oplot[iplot],default=False) 
283  #                    if result != False: 
284  #                        atitle+='  & '+oplot[iplot]+' '+linestyle[iplot%3] 
285  #                        for iz in ions: 
286  #                            pl.plot(self.IoneqTemperature, self.IoneqAll[self.Z-1, iz-1], linestyle[iplot%4]) 
287  #            else: 
288  #                print ' oplot not understood ', oplot 
289  #        if title: 
290  #            pl.title(atitle) 
291  #        pl.axis(xyr) 
292          # 
293