defangle.py

   1 import matplotlib as mpl
   2 import matplotlib.mlab as mlab
   3 import matplotlib.pyplot as plt
   4 import matplotlib.axes as axe
   5 import matplotlib.collections as collections
   6 import numpy as np
   7 import random
   8 from pylab import *
   9 from matplotlib.pyplot import *
  10
  11 ANGLE = "../data/rk_defangle.dat"
  12 VELO = "../data/rk_defangle_velo.dat"
  13 FOUT = "../defangle.png"
  14
  15 xlabel = r'$h$'
  16 ylabel = r'$\delta y$'
  17 res = (1440/80,900/80)
  18 pad1 = 1
  19 pad2 = 0.1
  20 tsize = 16
  21 lsize = 24
  22 max = 0.875
  23 acangle = 58
  24
  25 ######################################################################
  26
  27 angle = mlab.csv2rec(ANGLE, delimiter='\t', comments='#')
  28 velo = mlab.csv2rec(VELO, delimiter='\t', comments='#')
  29
  30 r = angle.r
  31 psi = angle.x
  32 vx = velo.x
  33
  34
  35
  36 def mean_angle(param,angle):
  37         unique_param = unique(param)
  38         list_angles = [ [] for DUMMYVAR in range(len(unique_param)) ]
  39         ret_avg = []
  40         for i in range(len(unique_param)):
  41                 for j in range(len(angle)):
  42                         if (param[j] == unique_param[i]):
  43                                 list_angles[i].append(angle[j])
  44                                 #print len(list_angles[i])
  45                 ret_avg.append(abs(float(sum(list_angles[i])) /
  46                         len(list_angles[i])))
  47                 #print len(ret_avg)
  48         return ret_avg
  49
  50
  51 left, width = 0.1, 0.8
  52 c = max/2
  53 rect1 = [left, c+0.1, width, c]
  54 rect2 = [left, 0.1, width, c]
  55
  56 rad_ac = acangle*(np.pi/180)
  57
  58 fig = plt.figure(figsize=res)
  59
  60 ax1 = fig.add_axes(rect1)
  61 ax2 = fig.add_axes(rect2)
  62
  63 avg_angles = mean_angle(r,psi)
  64 avg_velo = mean_angle(r,vx)
  65 unique_r = unique(r)
  66
  67 ax1.plot(unique_r, avg_angles, marker='.', color='black')
  68 ax2.plot(unique_r, avg_velo, marker='.', color='black')
  69
  70 xticks = np.arange(0,360+1,45)
  71 yticks1 = np.arange(135,315+1,45)
  72 yticks2 = np.arange(0,1.1,0.5)
  73
  74 ax1.set_xlim(r.min()-pad1,r.max()+pad1)
  75 ax1.set_ylim(yticks1.min(),yticks1.max())
  76 ax2.set_xlim(r.min()-pad2,r.max()+pad2)
  77 ax2.set_ylim(0,vx.max()+pad2)
  78 ax1.set_xticks([])
  79 ax2.set_xticks(xticks)
  80 ax1.set_yticks(yticks1)
  81 ax2.set_yticks(yticks2)
  82 ax1.axvline(acangle, color='black', ls=":")
  83 ax1.axvline(acangle-45, color='black', ls="-.")
  84 ax1.axvline(acangle+45, color='black', ls="-.")
  85 ax1.axvline(acangle+180, color='black', ls=":")
  86 ax1.axvline(acangle+180-45, color='black', ls="-.")
  87 ax1.axvline(acangle+180+45, color='black', ls="-.")
  88 ax2.axvline(acangle, color='black', ls=":")
  89 ax2.axvline(acangle+180, color='black', ls=":")
  90 ax2.axvline(acangle+180-45, color='black', ls="-.")
  91 ax2.axvline(acangle+180+45, color='black', ls="-.")
  92 ax2.axvline(acangle-45, color='black', ls="-.")
  93 ax2.axvline(acangle+45, color='black', ls="-.")
  94
  95 ax1.set_ylabel(r'$\psi$', size=lsize)
  96 ax2.set_ylabel(r'$v_x$', size=lsize)
  97 ax2.set_xlabel(r'$\phi$', size=lsize)
  98 ax2.set_xticklabels(ax2.get_xticks(), size=tsize)
  99 ax1.set_yticklabels(ax1.get_yticks(), size=tsize)
 100 ax2.set_yticklabels(ax2.get_yticks(), size=tsize)
 101
 102 plt.savefig(FOUT)
 103
 104 #plt.show()