PYTHON 7
Plot_eta.py Guest on 20th November 2020 06:00:43 PM
  1. import matplotlib, sys, os
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. import scipy.optimize as optimization
  5.  
  6.  
  7. def t(g, eta, N=1):
  8.         i = N/2.
  9.         nuc = (i*(i+2)/(5/4.))**0.5
  10.         if g==0:
  11.                 return (1. + eta**2 + 4 /nuc*eta**4)/(2 *nuc*eta + 4 *eta**3)
  12.         return (1 + 4 *g**4 + eta**2 + 4 *eta**4 +  g**2 *(5 + 8 *eta**2))/(2 *(eta + 2 *g**2 *eta + 2 *eta**3))
  13. #(4 + g**4 + eta**2 + eta**4 + g**2 *(5 + 2 *eta**2))/(2*eta* (2 + g**2 + eta**2))
  14.         #return (1 + 4 g^4 + \[Eta]^2 + 4 \[Eta]^4 +  g^2 (5 + 8 \[Eta]^2))/(2 (\[Eta] + 2 g^2 \[Eta] + 2 \[Eta]^3))
  15. #(1 + 4 g^4 + \[Eta]^2 + 4 \[Eta]^4 +
  16.  #g^2 (5 + 8 \[Eta]^2))/(2 (\[Eta] + 2 g^2 \[Eta] + 2 \[Eta]^3))
  17. try:
  18.         option = sys.argv[1]
  19. except IndexError:
  20.         option = '-t'
  21. COLORS = ['b','r','y','g','black','white','magenta','purple', 'cyan','grey','brown','pink','b','r','y','g','black','white','magenta','purple', 'cyan','grey','brown','pink']
  22. t_scale = []
  23. eta_scale = []
  24. short_times = []
  25. long_times = []
  26. MAX = 50
  27. for n in xrange(1, MAX):
  28.         x, y = [], []
  29.         try:
  30.                 with open('archive/mean/mean_times_N=%d'%n, 'r+') as f:
  31.                         for l in f.readlines():
  32.                                 y.append(float(l.split()[0]))
  33.                                 x.append(float(l.split()[2]))
  34.                         if option == '-t':             
  35.                                 plt.plot(x, y, label='N=%d'%n)
  36.                         t_scale.append((float(min(y)), n))
  37.                         eta_scale.append((float(x[y.index(min(y))]), n))
  38.                         short_times.append((y[2], n))
  39.                         try:
  40.                                 long_times.append((y[90], n))
  41.                         except:
  42.                                 continue
  43.         except IOError:
  44.                 continue
  45.  
  46.  
  47. tt = zip(*t_scale)
  48. ee = zip(*eta_scale)
  49. st = zip(*short_times)
  50. lt = zip(*long_times)
  51.  
  52.  
  53. if option == '-t':     
  54.         plt.ylabel(r'$\tau$', fontsize=24)
  55.         plt.xlabel(r'$\eta$', fontsize=24)
  56.         plt.xlim([0,50])
  57.         plt.ylim([0,20])
  58.         plt.legend()
  59.         #g = 1
  60. if option == '-an':    
  61.         x, y = [], []
  62.         with open('archive/mean_times_many_N=%d'%n, 'r+') as f:
  63.                 for l in f.readlines():
  64.                         y.append(float(l.split()[0]))
  65.                         x.append(float(l.split()[2]))
  66.                 plt.scatter(x, y, label='N=%d'%n, color=COLORS[n])
  67.                
  68.         for n in xrange(1,2):
  69.                 plt.plot(x, [t(0, eta/4, n) for eta in x], label='analytical N=%d'%n)
  70.         plt.ylabel(r'$\tau$', fontsize=24)
  71.         plt.xlabel(r'$\eta$', fontsize=24)
  72.         plt.xlim([0,5])
  73.         plt.ylim([0,20])
  74.         plt.legend()
  75.  
  76.  
  77. elif option == '-eta'
  78.         plt.scatter(ee[1], ee[0])
  79.         plt.ylabel(r'$\eta_{min}$', fontsize=24)
  80.         plt.xlabel(r'$N$', fontsize=24)
  81.         plt.xlim([0,MAX])
  82.         def func(x, a, b):
  83.                 return a + b * x
  84.         x0 = np.array([0.0, 0.0])
  85.         (a,b), corr = optimization.curve_fit(func, np.array(ee[1]), np.array(ee[0]), x0)
  86.         plt.annotate(r'$f(N) = %.3f  %+.3f\frac{N}{2}$'%(a,b), xy=(15, 1))
  87.         plt.plot(tt[1], [func(n, a, b) for n in ee[1]], label='analytical N=%d'%n)
  88. elif option == '-tau'
  89.         plt.scatter(tt[1], tt[0])
  90.         plt.ylabel(r'$\tau_{min}$', fontsize=24)
  91.         plt.xlabel(r'$N$', fontsize=24)
  92.         plt.xlim([0,MAX])
  93.         def func(x, a, b, c):
  94.                 return a + (x/2.*(b*x/2.+c))**0.5
  95.         x0 = np.array([0.0, 0.0, 0.0])
  96.         (a,b,c), corr = optimization.curve_fit(func, np.array(tt[1]), np.array(tt[0]), x0)
  97.        
  98.         plt.annotate(r'$f(N) = %.3f + \sqrt{\frac{N}{2}(%.3f\frac{N}{2}%+.3f)}$'%(a,b,c), xy=(15, 1))
  99.         plt.plot(tt[1], [func(n, a,  b,  c) for n in tt[1]], label='analytical N=%d'%n)
  100.        
  101. elif option == '-small':       
  102.         plt.scatter(st[1], st[0],label=r'$\eta \ll 1$')
  103.         plt.ylabel(r'$\tau$', fontsize=24)
  104.         plt.xlabel(r'$N$', fontsize=24)
  105.         plt.legend()
  106.         plt.xlim([0,MAX])
  107.         def func(x, a, b, c):
  108.                 return a + (x/2.*(b*x/2.+c))**0.5
  109.         x0 = np.array([0.0, 0.0, 0.0])
  110.         (a,b,c), corr = optimization.curve_fit(func, np.array(st[1]), np.array(st[0]), x0)
  111.         plt.annotate(r'$f(N) = %.3f + \sqrt{\frac{N}{2}(%.3f\frac{N}{2}%+.3f)}$'%(a,b,c), xy=(15, 1))
  112.         plt.plot(tt[1], [func(n, a,  b,  c) for n in st[1]], label='analytical N=%d'%n)
  113. elif option == '-large':       
  114.         plt.scatter(lt[1], lt[0],label=r'$\eta \gg 1$')
  115.         plt.ylabel(r'$\tau$', fontsize=24)
  116.         plt.xlabel(r'$N$', fontsize=24)
  117.         plt.legend()
  118.         plt.xlim([0,MAX])
  119.         def func(x, a, b, c):
  120.                 return a + (x/2.*(b*x/2.+c))**-0.5
  121.         x0 = np.array([0.0, 0.0, 0.0])
  122.         (a,b,c), corr = optimization.curve_fit(func, np.array(lt[1]), np.array(lt[0]), x0)
  123.         plt.annotate(r'$f(N) = %.3f + \frac{1}{\sqrt{\frac{N}{2}(%.3f\frac{N}{2}%+.3f)}}$'%(a,b,c), xy=(15, 1))
  124.         plt.plot(tt[1], [func(n, a,  b,  c) for n in lt[1]], label='analytical N=%d'%n)
  125.  
  126.  
  127.  
  128. #plt.ylim([0,20])
  129. #plt.ylim([1,20])
  130.  
  131. #
  132.  
  133.  
  134. plt.show()

Paste is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Raw Paste

Login or Register to edit or fork this paste. It's free.