import sys import os import random from random import seed from random import random from random import gauss import datetime import numpy as np import math from mpmath import * import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages from matplotlib.ticker import MultipleLocator from matplotlib.ticker import ScalarFormatter from matplotlib.ticker import NullFormatter import matplotlib.image as image import matplotlib.patches as patches import matplotlib.lines as lines import matplotlib.ticker as ticker from matplotlib.colors import BoundaryNorm from matplotlib.ticker import MaxNLocator import helloworld as hw mp.dps = 16; mp.pretty = False; def print_data (Nbins,binning,data1,data2,data3,filename): xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) f = open(filename,'w') for ix in range(len(xs)): for iy in range(len(ys)): f.write('%2d %2d %5.2f %5.2f %16.10f %16.10f %16.10f \n' % (ix,iy,xs[ix],ys[iy],data1[ix,iy],data2[ix,iy],data3[ix,iy])) f.close() return def read_data (Nbins,binning,filename): xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) data1 = np.zeros((Nbins,Nbins)) data2 = np.zeros((Nbins,Nbins)) data3 = np.zeros((Nbins,Nbins)) f = open(filename,'r') lines = f.readlines() for line in lines: p=line.split() ix = int(p[0]) iy = int(p[1]) data1[ix,iy] = float(p[4]) data2[ix,iy] = float(p[5]) data3[ix,iy] = float(p[6]) return data1,data2,data3 def getMax(M): maxVal = 0 for row in M: if max(row) > maxVal: maxVal = max(row) return maxVal def make_figure (xs,ys,d,data1,data2,data3,data4,data5,data6): # setto alcune variabili comuni axisticslabelfontsize=8 axisticslabelfontsizeinset=7 axislabelfontsize=11 axislabelfontsizeinset=9 xsize = 7.0 ysize = 3.3 xgap = 0.07 x1 = xgap+0.03 y1 = 0.56 x1size = 0.17 y1size = x1size * xsize/ysize x4size = x1size *1.25 x2 = x1 + 1.5*xgap + x4size x4 = x2 + 1.5*xgap + x4size y5 = 0.09 with PdfPages('fig1.pdf') as pdf: fig = plt.figure(figsize=(xsize,ysize)) plt.rc('text', usetex=True) plt.rc('text.latex', preamble = ','.join('''\usepackage{txfonts} \usepackage{lmodern}'''.split())) cmap = plt.get_cmap('bwr') ######################### panel 1 ############################## panel = fig.add_axes([x1, y1, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data1.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=VMAX) pcm = panel.pcolormesh(xs/d, ys/d, data1.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.4) cbar = plt.colorbar(pcm) panel.text(0.3,1.08,r'$t=0.2$',fontsize=10,transform=panel.transAxes) panel.text(-0.55, 0.7, 'simulations',rotation='vertical',transform=panel.transAxes) ######################### panel 2 ############################## panel = fig.add_axes([x2, y1, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data2.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=VMAX) pcm = panel.pcolormesh(xs/d, ys/d, data2.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.08) panel.text(0.3,1.08,r'$t=1.0$',fontsize=10,transform=panel.transAxes) cbar = plt.colorbar(pcm) ######################### panel 4 ############################## panel = fig.add_axes([x4, y1, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data4.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=VMAX) pcm = panel.pcolormesh(xs/d, ys/d, data3.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.02) cbar = plt.colorbar(pcm) panel.text(0.25,1.08,r'$t=100.0$',fontsize=10,transform=panel.transAxes) ######################### panel 5 ############################## panel = fig.add_axes([x1, y5, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data5.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=VMAX) pcm = panel.pcolormesh(xs/d, ys/d, data4.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.4) cbar = plt.colorbar(pcm) panel.text(-0.55, 0.64, 'numerics',rotation='vertical',transform=panel.transAxes) ######################### panel 6 ############################## panel = fig.add_axes([x2, y5, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data6.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=VMAX) pcm = panel.pcolormesh(xs/d, ys/d, data5.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.08) cbar = plt.colorbar(pcm) #~ ######################### panel 8 ############################## panel = fig.add_axes([x4, y5, x4size, y1size]) panel.tick_params(axis='both',which='both',direction='in') for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize) panel.set_xlabel(r'$x/d$',fontsize=axislabelfontsize,labelpad=-1) panel.set_ylabel(r'$y/d$',fontsize=axislabelfontsize,labelpad=-1) #~ pcm = panel.pcolormesh(xs, ys, data8.T, shading = 'gouraud', cmap=cmap) pcm = panel.pcolormesh(xs/d, ys/d, data6.T, shading = 'gouraud', cmap=cmap, vmin=0, vmax=0.02) cbar = plt.colorbar(pcm) pdf.savefig(fig) return def main(): # parameters D_rot = 0.8 D = 1.0 v = 4.0 k = 1.0 mu = 1.0 # derived adimensional parameters d = np.sqrt(D/(mu*k)) # better choice of the parameters is possible by removing the ones that fix the units. For the moment I keep everything equal to 1 # Eq.s of motion are # x = x + v * cos(theta) * dt - mu * k * x * dt + sqrt(2*D*dt) * eta # y = y + v * sin(theta) * dt - mu * k * y * dt + sqrt(2*D*dt) * eta # theta = theta + sqrt(2 * D_rot * dt) * eta print 'tau = ',1./(mu*k) print 'd = ',d print 'epsilon = ',v*d/D # other parameters dt = 0.01 # time integration step time_duration = 100.0 # time duration of a given trajectory # initial values x0 = 4.0 y0 = 0.0 theta0 = np.pi/2. # number of independent realization for the statistic in the simulations Nstories = 200000 # time snapshots t1,t2,t3 = 0.2,1.0,100.0 # 2D hystogram data binning = 0.5 Nbins = 29 # 1 in zero, Nbins/2 on the left and Nbins/2 on the right xs = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) ys = np.arange(-(Nbins-1)/2*binning,(Nbins-1)/2*binning + binning,binning) #~ data1,data2,data3 = hw.spatial_distribution(D_rot,D,v,k,mu,dt,time_duration,x0,y0,theta0,Nstories,t1,t2,t3,binning,Nbins) #~ print_data (Nbins,binning,data1,data2,data3,'data_from_simulations.dat') data1,data2,data3 = read_data (Nbins,binning,'data_from_simulations.dat') data4 = hw.spatial_distribution_numerics(D_rot, D, v, k, mu, x0, y0, theta0, t1, binning, Nbins, 40) data5 = hw.spatial_distribution_numerics(D_rot, D, v, k, mu, x0, y0, theta0, t2, binning, Nbins, 40) data6 = hw.spatial_distribution_numerics(D_rot, D, v, k, mu, x0, y0, theta0, t3, binning, Nbins, 40) print_data (Nbins,binning,data4,data5,data6,'data_from_numerics_q40.dat') #~ data4,data5,data6 = read_data (Nbins,binning,'data_from_numerics_q40.dat') #~ data4,data5,data6 = hw.spatial_distribution_numerics_L_diagonalization(D_rot, D, v, k, mu, x0, y0, theta0, t1, t2, t3, binning, Nbins, 20) #~ print_data (Nbins,binning,data4,data5,data6,'data_from_numerics_on_psi_q20.dat') #~ data4,data5,data6 = read_data (Nbins,binning,'data_from_numerics_on_psi_q10.dat') # plots make_figure (xs,ys,d,data1,data2,data3,data4,data5,data6) #~ make_figure_prova (xs,ys,data1,data2,data3,data4,data5,data6) main()