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 cmath from scipy.linalg import eig import scipy.special as sc from scipy import linalg as sciLA 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): rs = np.asarray([binning*i for i in range(Nbins)]) f = open(filename,'w') for ix in range(Nbins): f.write('%2d %5.2f %16.10f %16.10f %16.10f \n' % (ix,rs[ix],data1[ix],data2[ix],data3[ix])) f.close() return def read_data (Nbins,binning,filename): rs = np.asarray([binning*i for i in range(Nbins)]) data1 = np.zeros((Nbins)) data2 = np.zeros((Nbins)) data3 = np.zeros((Nbins)) f = open(filename,'r') lines = f.readlines() for line in lines: p=line.split() ix = int(p[0]) data1[ix] = float(p[2]) data2[ix] = float(p[3]) data3[ix] = float(p[4]) return data1,data2,data3 def get_trajectory1 (D_rot,D,v,dt,time_duration,x0,y0,theta0,xf_target,yf_target): factor = sqrt(2*D*dt) factor2 = sqrt(2*D_rot*dt) factor4 = v * dt xf = x0 yf = y0 Nt = (int) (time_duration/dt) print 'dt=',dt,factor,factor4,factor2/(2*np.pi) print 'dt=',100*dt,10*factor,100*factor4,10*factor2/(2*np.pi) seed(7) i=0 while (abs(xf-xf_target)>0.03 or abs(yf-yf_target)>0.03): i += 1 if (i % 100 ==0): print i x = x0 y = y0 xs=[x0] ys=[y0] theta = theta0 for t in range(1,Nt+1): x = x + factor4 * cos(theta) + factor * gauss(0, 1) y = y + factor4 * sin(theta) + factor * gauss(0, 1) theta = theta + factor2 *gauss(0, 1) xs.append(x) ys.append(y) xf = x yf = y return xs,ys def make_figure (rs,data1,data2,data3,rs2,data4,data5,data6,data7,data8,data9,data10,data11,data12,t1,t2,t3,xs1,ys1,xs2,ys2): # setto alcune variabili comuni axisticslabelfontsize=8 axisticslabelfontsizeinset=7 axislabelfontsize=11 axislabelfontsizeinset=9 xsize = 3.5 ysize = 3.0 x1 = 0.14 x1size = 0.84 y1 = 0.55 y1size = 0.43 y2 = 0.12 x3 = 0.43 y3 = 0.2 x3size = 0.23 y3size = 0.26 with PdfPages('fig6.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())) ######################### panel 1 ############################## panel = fig.add_axes([x1, y1, x1size, y1size]) panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True) 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_ylabel(r'$P(r,t)$',fontsize=axislabelfontsize) panel.set_xticklabels([]) panel.plot(rs,data1,'o',color='orange',markersize=3,label=r'simulations $dt=10^{-4}\tau$') panel.plot(rs,data10,'go',markersize=2,label=r'simulations $dt=10^{-3}\tau$') panel.plot(rs,data7,'bo',markersize=1,label=r'simulations $dt=10^{-2}\tau$') panel.plot(rs2,data4,'r-',label='numerics',linewidth=1) panel.legend(loc='upper right', bbox_to_anchor=(0.9, 0.99),ncol=1,fontsize=7,handlelength=1.0,labelspacing=0.2) panel.text(0.05,0.86,(r'$t=%3.2f \tau$' % t1),fontsize=10,transform=panel.transAxes) panel.set_xlim(0,1.1) panel.xaxis.set_major_locator(MultipleLocator(0.2)) panel.xaxis.set_minor_locator(MultipleLocator(0.1)) panel.axvline(x = 1, color = 'k', linewidth=2.0) panel.set_ylim(0,3.7) panel.yaxis.set_major_locator(MultipleLocator(0.5)) panel.yaxis.set_minor_locator(MultipleLocator(0.25)) ######################### panel 2 ############################## panel = fig.add_axes([x1, y2, x1size, y1size]) panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True) 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'$r/R$',fontsize=axislabelfontsize) panel.set_ylabel(r'$P(r,t)$',fontsize=axislabelfontsize) panel.plot(rs,data3,'o',color='orange',markersize=3) panel.plot(rs,data9,'go',markersize=2) panel.plot(rs,data12,'bo',markersize=1) panel.plot(rs2,data6,'r-',label='numerics',linewidth=1) panel.text(0.05,0.86,(r'$t=%3.2f \tau$' % t3),fontsize=10,transform=panel.transAxes) panel.set_xlim(0,1.1) panel.xaxis.set_major_locator(MultipleLocator(0.2)) panel.xaxis.set_minor_locator(MultipleLocator(0.1)) panel.axvline(x = 1, color = 'k', linewidth=2.0) panel.set_ylim(0,1.2) panel.yaxis.set_major_locator(MultipleLocator(0.5)) panel.yaxis.set_minor_locator(MultipleLocator(0.25)) ######################### panel 3 ############################## panel = fig.add_axes([x3, y3, x3size, y3size]) panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True) for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsizeinset-1) for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsizeinset-1) panel.set_xlabel(r'$x/R$',fontsize=axislabelfontsizeinset-1,labelpad=-1) panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsizeinset-1) index = 468 panel.plot(xs1[:index],ys1[:index],'r-',linewidth=0.4) panel.plot(xs1[index:],ys1[index:],'-',color='orange',linewidth=0.4) panel.plot(xs2,ys2,'b-',linewidth=0.7) panel.plot(xs1[1],ys1[1],'ko',markersize=2) panel.plot(xs1[-1],ys1[-1],'ko',markersize=2) circle = plt.Circle((0,0),1.0,facecolor='None',lw=1, edgecolor='k') panel.add_patch(circle) panel.text(0.19,0.19,('A'),fontsize=8,transform=panel.transAxes) panel.text(0.65,0.1,('B'),fontsize=8,transform=panel.transAxes) panel.set_xlim(0.05,0.85) panel.xaxis.set_major_locator(MultipleLocator(0.2)) panel.xaxis.set_minor_locator(MultipleLocator(0.1)) panel.set_ylim(0.3,1.1) panel.yaxis.set_major_locator(MultipleLocator(0.2)) panel.yaxis.set_minor_locator(MultipleLocator(0.1)) pdf.savefig(fig) return def main(): # parameters Gamm=0.8 Pe = 4.0 # derived parameters for integrating the equation of motion D_rot = Gamm D = 1.0 v = Pe # other parameters dt = 0.0001 # time integration step time_duration = 0.1 # time duration of a given trajectory # initial position r0 = 0.0 phi0 = 0.0 # number of independent realization for the statistic in the simulations Nparticles = 10000000 # time snapshots t1,t2,t3 = 0.02,0.05,time_duration # 1D hystogram data binning = 0.02 Nbins = int(1.0/binning) rs = np.asarray([binning*i for i in range(Nbins)]) #~ data1,data2,data3 = hw.radial_distribution_avg(D_rot,D,v,dt,time_duration,r0,phi0,Nparticles,t1,t2,t3,binning,Nbins) #~ print_data (Nbins,binning,data1,data2,data3,'radialPDF_from_simulations.dat') data1,data2,data3 = read_data (Nbins,binning,'radialPDF_from_simulations.dat') dt=0.001 #~ data7,data8,data9 = hw.radial_distribution_avg(D_rot,D,v,dt,time_duration,r0,phi0,Nparticles,t1,t2,t3,binning,Nbins) #~ print_data (Nbins,binning,data7,data8,data9,'radialPDF_from_simulations_dt0.001.dat') data7,data8,data9 = read_data (Nbins,binning,'radialPDF_from_simulations_dt0.001.dat') dt=0.01 #~ data10,data11,data12 = hw.radial_distribution_avg(D_rot,D,v,dt,time_duration,r0,phi0,Nparticles,t1,t2,t3,binning,Nbins) #~ print_data (Nbins,binning,data10,data11,data12,'radialPDF_from_simulations_dt0.01.dat') data10,data11,data12 = read_data (Nbins,binning,'radialPDF_from_simulations_dt0.01.dat') binning = 0.01 Nbins = int(1.0/binning) rs2 = np.asarray([binning*i for i in range(Nbins)]) data4,data5,data6 = read_data (Nbins,binning,'radialPDF_from_numerics.dat') dt = 0.0001 time_duration = 0.06 x0,y0 = 0.3,0.6 xf_target,yf_target = 0.7,0.5 theta0 = np.pi * 0.2 xs1,ys1 = get_trajectory1 (D_rot,D,v,dt,time_duration,x0,y0,theta0,xf_target,yf_target) indexes = [1,150,190,310,450,560,-1] xs2=[] ys2=[] for i in indexes: xs2.append(xs1[i]) ys2.append(ys1[i]) # plots make_figure (rs,data1,data2,data3,rs2,data4,data5,data6,data7,data8,data9,data10,data11,data12,t1,t2,t3,xs1,ys1,xs2,ys2) main()