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 (xs,ys,data1,filename): 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 \n' % (ix,iy,xs[ix],ys[iy],data1[ix,iy])) f.close() return def read_data (Nbins,Nbins2,filename): xs = np.zeros((Nbins)) ys = np.zeros((Nbins2)) data1 = np.zeros((Nbins,Nbins2)) f = open(filename,'r') lines = f.readlines() for line in lines: p=line.split() ix = int(p[0]) iy = int(p[1]) xs[ix] = float(p[2]) ys[iy] = float(p[3]) data1[ix,iy] = float(p[4]) return xs,ys,data1 def print_data_time (Nbins,binning,data1,filename): ts = 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 \n' % (ix,ts[ix],data1[ix])) f.close() return def read_data_time (Nbins,binning,filename): ts = 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]) return data1 def make_figure (ts,data1,data1_0,data1_2,data1_3,data1_4,data1_5,ts2,data2,data2_0,data2_2,data2_3,data2_4,data2_5,ts3,data3,data3_0,data3_2,data3_3,data3_4,data3_5): # setto alcune variabili comuni axisticslabelfontsize=8 axisticslabelfontsizeinset=7 axislabelfontsize=11 axislabelfontsizeinset=9 xsize = 3.5 ysize = 2.5 x1 = 0.14 y1 = 0.16 x1size = 0.84 y1size = 0.82 with PdfPages('fig5.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_xlabel(r'$t/\tau$',fontsize=axislabelfontsize) panel.set_ylabel(r'$F(t)$',fontsize=axislabelfontsize) x_dummy = [1.1,2.1] y_dummy = [100,110] panel.plot(x_dummy ,y_dummy,'o',color='gray',markersize=3,label='simulations') panel.plot(x_dummy ,y_dummy,'k-',label='numerics') panel.plot(x_dummy ,y_dummy,'k--',label=r'analytics ($q=2$)',linewidth=0.6) panel.plot(x_dummy ,y_dummy,'w--',label=r'$\,$',linewidth=0.6) panel.plot(ts,data1_0,'o',color='cyan',markersize=3,label='Pe=$0$') panel.plot(ts,data1,'^',color='orange',markersize=3,label='Pe=$4, \, \gamma=2$') panel.plot(ts,data1_5,'s',color='lightgreen',markersize=3,label='Pe=$4, \, \gamma=10$') panel.plot(ts,data1_3,'d',color='violet',markersize=3,label='Pe=$8, \, \gamma=0$') panel.plot(ts,data1_2,'v',color='yellow',markersize=3,label='Pe=$8, \, \gamma=2$') panel.plot(ts,data1_4,'*',color='peru',markersize=3,label='Pe=$8, \, \gamma=10$') panel.plot(ts2,data2_0,'b-') panel.plot(ts2,data2,'r-') panel.plot(ts2,data2_5,'g-') panel.plot(ts2,data2_3,'m-') panel.plot(ts2,data2_2,'y-') panel.plot(ts2,data2_4,'-',color='brown') panel.plot(ts3,data3_0,'k--',linewidth=0.6) panel.plot(ts3,data3,'r--',linewidth=0.6) panel.plot(ts3,data3_5,'g--',linewidth=0.6) panel.plot(ts3,data3_3,'m--',linewidth=0.6) panel.plot(ts3,data3_2,'y--',linewidth=0.6) panel.plot(ts3,data3_4,'--',color='brown',linewidth=0.6) panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=7,handlelength=1.0,labelspacing=0.2) panel.set_xlim(0,0.62) panel.xaxis.set_major_locator(MultipleLocator(0.2)) panel.xaxis.set_minor_locator(MultipleLocator(0.05)) panel.set_ylim(0,11) panel.yaxis.set_major_locator(MultipleLocator(2)) panel.yaxis.set_minor_locator(MultipleLocator(1)) pdf.savefig(fig) return def main(): # parameters Gamm=2.0 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.62 # time duration of a given trajectory # initial values r0 = 0.2 phi0 = 0.0 theta0 = np.pi/2. x0 = r0*np.cos(phi0) y0 = r0*np.sin(phi0) # SIMULATIONS # number of independent realization for the statistic in the simulations Nparticles = 10000000 # times binning = 0.01 Nbins = int(time_duration/binning) ts = np.asarray([binning*i for i in range(Nbins)]) #~ data1 = hw.fpt(D_rot,D,v,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1,'fpt_vs_time_from_simulations_Pe4_gamma2.dat') data1 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe4_gamma2.dat') #~ data1_0 = hw.fpt(D_rot,D,0.0,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1_0,'fpt_vs_time_from_simulations_Pe0.dat') data1_0 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe0.dat') #~ data1_2 = hw.fpt(D_rot,D,8.0,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1_2,'fpt_vs_time_from_simulations_Pe8_gamma2.dat') data1_2 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe8_gamma2.dat') #~ data1_3 = hw.fpt(0.0,D,8.0,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1_3,'fpt_vs_time_from_simulations_Pe8_gamma0.dat') data1_3 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe8_gamma0.dat') #~ data1_4 = hw.fpt(10.0,D,8.0,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1_4,'fpt_vs_time_from_simulations_Pe8_gamma10.dat') data1_4 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe8_gamma10.dat') #~ data1_5 = hw.fpt(10.0,D,Pe,dt,r0,phi0,theta0,Nparticles,binning,Nbins) #~ print_data_time (Nbins,binning,data1_5,'fpt_vs_time_from_simulations_Pe4_gamma10.dat') data1_5 = read_data_time (Nbins,binning,'fpt_vs_time_from_simulations_Pe4_gamma10.dat') # NUMERICS # times binning = 0.005 Nbins = int(time_duration/binning) ts2 = np.asarray([binning*i for i in range(Nbins)]) data2 = hw.fpt_numerics(Gamm,Pe,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2,'fpt_vs_time_from_numerics_Pe4_gamma2.dat') #~ data2 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe4_gamma2.dat') data2_0 = hw.fpt_numerics(Gamm,0.0,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2_0,'fpt_vs_time_from_numerics_Pe0.dat') #~ data2_0 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe0.dat') data2_2 = hw.fpt_numerics(Gamm,8.0,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2_2,'fpt_vs_time_from_numerics_Pe8_gamma2.dat') #~ data2_2 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe8_gamma2.dat') data2_3 = hw.fpt_numerics(0.0,8.0,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2_3,'fpt_vs_time_from_numerics_Pe8_gamma0.dat') #~ data2_3 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe8_gamma0.dat') data2_4 = hw.fpt_numerics(10.0,8.0,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2_4,'fpt_vs_time_from_numerics_Pe8_gamma10.dat') #~ data2_4 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe8_gamma10.dat') data2_5 = hw.fpt_numerics(10.0,Pe,r0,phi0,theta0,binning,Nbins) print_data_time (Nbins,binning,data2_5,'fpt_vs_time_from_numerics_Pe4_gamma10.dat') #~ data2_5 = read_data_time (Nbins,binning,'fpt_vs_time_from_numerics_Pe4_gamma10.dat') # ANALYTICS # times binning = 0.01 Nbins = int(time_duration/binning) ts3 = np.asarray([binning*i for i in range(Nbins)]) #~ data3 = hw.fpt_analytics(Gamm,Pe,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3,'fpt_vs_time_from_analytics_Pe4_gamma2.dat') data3 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe4_gamma2.dat') #~ data3_0 = hw.fpt_analytics(Gamm,0.0,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3_0,'fpt_vs_time_from_analytics_Pe0.dat') data3_0 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe0.dat') #~ data3_2 = hw.fpt_analytics(Gamm,8.0,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3_2,'fpt_vs_time_from_analytics_Pe8_gamma2.dat') data3_2 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe8_gamma2.dat') #~ data3_3 = hw.fpt_analytics(0.0,8.0,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3_3,'fpt_vs_time_from_analytics_Pe8_gamma0.dat') data3_3 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe8_gamma0.dat') #~ data3_4 = hw.fpt_analytics(10.0,8.0,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3_4,'fpt_vs_time_from_analytics_Pe8_gamma10.dat') data3_4 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe8_gamma10.dat') #~ data3_5 = hw.fpt_analytics(10.0,Pe,r0,phi0,theta0,binning,Nbins) #~ print_data_time (Nbins,binning,data3_5,'fpt_vs_time_from_analytics_Pe4_gamma10.dat') data3_5 = read_data_time (Nbins,binning,'fpt_vs_time_from_analytics_Pe4_gamma10.dat') # plots make_figure (ts,data1,data1_0,data1_2,data1_3,data1_4,data1_5,ts2,data2,data2_0,data2_2,data2_3,data2_4,data2_5,ts3,data3,data3_0,data3_2,data3_3,data3_4,data3_5) main()