Source Code for Module ComboCode.cc.plotting.PlotMeixner

 1  import numpy as np 
 2  from math import sin,cos,exp 
 3  import matplotlib.pyplot as plt 
 4  import pylab as pl 
 5   
 6  from cc.plotting import Plotting2 
 7   
 8   
 9   
10 -def expDensFunc(r,rSw,exponent): 
11   
12      return exp(-(r/rSw)**exponent) 
13   
14   
15   
16  @np.vectorize 
17 -def densityFunc(r,theta,A,B,C,D,E,F,rMin,rSw): 
18   
19      exponent = -B*(1+C*(sin(theta)**F)*(expDensFunc(r,rSw,D)/\ 
20                                          expDensFunc(rMin,rSw,D))) 
21      rho = (r/rMin)**exponent 
22      rho = rho * (1+A*((1-cos(theta))**F)*(expDensFunc(r,rSw,E)/\ 
23                                            expDensFunc(rMin,rSw,E))) 
24   
25      return rho 
26   
27   
28   
29 -def plotDens(A,B,C,D,E,F,rPlot,rMin,rSw,nRad,nTheta,filename=None,\ 
30               extension='pdf',landscape=0,show=0,figsize=(5.5, 10)): 
31       
32      #-- Define the radial and angular grids 
33      R = np.logspace(np.log10(rMin), np.log10(rPlot), num=nRad, endpoint=True) 
34      Theta = np.arange(0, np.pi/2.0+1.0/float(nTheta), (np.pi/2.0)/(nTheta)) 
35      R,Theta = np.meshgrid(R,Theta) 
36       
37      #-- Calculate the Meixner density distribution 
38      rho = densityFunc(R,Theta,A,B,C,D,E,F,rMin,rSw) 
39   
40      #-- Figure consists of two subplots, one with the density map, and one with 
41      #   profiles at 0 and 90 degrees. 
42      fig = plt.figure(figsize=figsize) 
43       
44      #-- Make the color scale density plot 
45      ax = fig.add_subplot(211) 
46      pl.pcolor(R*np.sin(Theta),R*np.cos(Theta),np.log10(rho/rho.max())) 
47      pl.pcolor(-R*np.sin(Theta),R*np.cos(Theta),np.log10(rho/rho.max())) 
48      pl.pcolor(-R*np.sin(Theta),-R*np.cos(Theta),np.log10(rho/rho.max())) 
49      pl.pcolor(R*np.sin(Theta),-R*np.cos(Theta),np.log10(rho/rho.max())) 
50      fig.subplots_adjust(right=0.8) 
51      fig.subplots_adjust(bottom=0.1) 
52      cbar_ax = fig.add_axes([0.85, 0.54, 0.05, 0.36]) 
53      pl.colorbar(cax=cbar_ax) 
54       
55      #-- Add density distribution at 0 and 90 degrees, and a r^-2 distribution 
56      ax2 = fig.add_subplot(212) 
57      ax2.plot(np.log10(R[0]),np.log10(rho[0]/rho[0].max()),'black',\ 
58               lw=2,marker='x',label=r'$\theta = 0^\circ$') 
59      ax2.plot(np.log10(R[-1]),np.log10(rho[-1]/rho[-1].max()),'magenta',\ 
60               lw=2,marker='|', label=r'$\theta = 90^\circ$') 
61      ax2.plot(np.log10(R[0]),np.log10(rMin*rMin/(R[0]*R[0])),'green',\ 
62               lw=2,label=r'$r^{-2}$') 
63      ax2.legend() 
64      label = r'$\rho_0[r_{\rm min}] / \rho_{90}[r_{\rm min}] = $'+\ 
65               '{:10.3e}'.format(rho[0][0]/rho[-1][0]) 
66      ax2.text(0.25,1.02,label,transform=ax2.transAxes) 
67   
68       
69      if filename: filename = Plotting2.saveFig(filename,extension,landscape) 
70      if show: pl.show() 
71       
72      return filename 
73   
74   
75   
76  if __name__ == '__main__': 
77   
78      A = 100.0 
79      B = 2.0 
80      C = 1.0 
81      D = 0.5 
82      E = 0.0 
83      F = 2.0 
84      rPlot = 1E16 
85      rMin = 1E13 
86      rSw = 5E13 
87      nTheta = 50 
88      nRad = 50 
89   
90      plotDens(A,B,C,D,E,F,rPlot,rMin,rSw,nRad,nTheta,show=1) 
91