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):
	xs=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	ys=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	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.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	ys=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	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 make_figure3x2 (xs,ys,data1,data2,data3,data4,data5,data6,t1,t2,t3):
	# 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.57
	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.1
	
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data1[i])
		if (value>vmax): vmax=value
		value = max(data4[i])
		if (value>vmax): vmax=value
	vmax1=vmax	
		
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data2[i])
		if (value>vmax): vmax=value
		value = max(data5[i])
		if (value>vmax): vmax=value
	vmax2=vmax
	
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data3[i])
		if (value>vmax): vmax=value
		value = max(data6[i])
		if (value>vmax): vmax=value
	vmax3=vmax
	
	
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data1[i])
		if (value<vmin): vmin=value
		value = max(data4[i])
		if (value<vmin): vmin=value
	vmin1=vmin	
		
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data2[i])
		if (value<vmin): vmin=value
		value = max(data5[i])
		if (value<vmin): vmin=value
	vmin2=vmin	
	
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data3[i])
		if (value<vmin): vmin=value
		value = max(data6[i])
		if (value<vmin): vmin=value
	vmin3=vmin	
	
	vmin1 = vmin2 = vmin3 = 0
	
	
	with PdfPages('fig2.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',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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data1.T, shading = 'gouraud', cmap=cmap, vmin=vmin1, vmax=vmax1)
		cbar = plt.colorbar(pcm)
		panel.text(0.3,1.082,(r'$t=%3.2f \tau$' % t1),fontsize=10,transform=panel.transAxes)
		panel.text(-0.55, 0.7, 'simulations',rotation='vertical',transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		######################### panel 2 ##############################
		panel = fig.add_axes([x2, y1, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data2.T, shading = 'gouraud', cmap=cmap, vmin=vmin2, vmax=vmax2)
		panel.text(0.3,1.082,(r'$t=%3.2f \tau$' % t2),fontsize=10,transform=panel.transAxes)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		
		######################### panel 4 ##############################
		panel = fig.add_axes([x4, y1, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data3.T, shading = 'gouraud', cmap=cmap, vmin=vmin3, vmax=vmax3)
		cbar = plt.colorbar(pcm)
		panel.text(0.25,1.082,(r'$t=%3.2f \tau$' % t3),fontsize=10,transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		
		
		
		
		
		######################### panel 5 ##############################
		panel = fig.add_axes([x1, y5, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data4.T, shading = 'gouraud', cmap=cmap, vmin=vmin1, vmax=vmax1)
		cbar = plt.colorbar(pcm)

		panel.text(-0.55, 0.64, 'numerics',rotation='vertical',transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		######################### panel 6 ##############################
		panel = fig.add_axes([x2, y5, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data5.T, shading = 'gouraud', cmap=cmap, vmin=vmin2, vmax=vmax2)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		#~ ######################### panel 8 ##############################
		panel = fig.add_axes([x4, y5, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data6.T, shading = 'gouraud', cmap=cmap, vmin=vmin3, vmax=vmax3)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		

		pdf.savefig(fig)
	return

def make_figure3x3 (xs,ys,data1,data2,data3,data4,data5,data6,data7,data8,data9,t1,t2,t3):
	# setto alcune variabili comuni
	axisticslabelfontsize=8
	axisticslabelfontsizeinset=7
	axislabelfontsize=11 
	axislabelfontsizeinset=9
	

	xsize = 7.0
	ysize = 4.9
	
	xgap = 0.07
	x1 = xgap+0.03
	x1size = 0.17
	y1size = x1size *  xsize/ysize
	x4size = x1size *1.25
	
	x2 = x1 + 1.5*xgap + x4size
	x4 = x2 + 1.5*xgap + x4size
	
	y1 = 0.7
	y2 = 0.38
	y3 = 0.06
	
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data1[i])
		if (value>vmax): vmax=value
		value = max(data4[i])
		if (value>vmax): vmax=value
		value = max(data7[i])
		if (value>vmax): vmax=value
	vmax1=vmax	
		
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data2[i])
		if (value>vmax): vmax=value
		value = max(data5[i])
		if (value>vmax): vmax=value
		value = max(data8[i])
		if (value>vmax): vmax=value
	vmax2=vmax
	
	vmax = 0
	for i in range(len(data1[0])):
		value = max(data3[i])
		if (value>vmax): vmax=value
		value = max(data6[i])
		if (value>vmax): vmax=value
		value = max(data9[i])
		if (value>vmax): vmax=value
	vmax3=vmax
	
	
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data1[i])
		if (value<vmin): vmin=value
		value = max(data4[i])
		if (value<vmin): vmin=value
		value = max(data7[i])
		if (value<vmin): vmin=value
	vmin1=vmin	
		
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data2[i])
		if (value<vmin): vmin=value
		value = max(data5[i])
		if (value<vmin): vmin=value
		value = max(data8[i])
		if (value<vmin): vmin=value
	vmin2=vmin	
	
	vmin = 0
	for i in range(len(data1[0])):
		value = max(data3[i])
		if (value<vmin): vmin=value
		value = max(data6[i])
		if (value<vmin): vmin=value
		value = max(data9[i])
		if (value<vmin): vmin=value
	vmin3=vmin	
	
	vmin1 = vmin2 = vmin3 = 0
	
	
	with PdfPages('fig2.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')
		cmap = plt.get_cmap('Reds')
		
		######################### panel 1 ##############################
		panel = fig.add_axes([x1, y1, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data1.T, shading = 'gouraud', cmap=cmap, vmin=vmin1, vmax=vmax1)
		cbar = plt.colorbar(pcm)
		panel.text(0.3,1.082,(r'$t=%3.2f \tau$' % t1),fontsize=10,transform=panel.transAxes)
		panel.text(-0.55, 0.7, 'simulations',rotation='vertical',transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		######################### panel 2 ##############################
		panel = fig.add_axes([x2, y1, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data2.T, shading = 'gouraud', cmap=cmap, vmin=vmin2, vmax=vmax2)
		panel.text(0.3,1.082,(r'$t=%3.2f \tau$' % t2),fontsize=10,transform=panel.transAxes)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		######################### panel 3 ##############################
		panel = fig.add_axes([x4, y1, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		pcm = panel.pcolormesh(xs, ys, data3.T, shading = 'gouraud', cmap=cmap, vmin=vmin3, vmax=vmax3)
		cbar = plt.colorbar(pcm)
		panel.text(0.25,1.082,(r'$t=%3.2f \tau$' % t3),fontsize=10,transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		
		
		
		
		
		######################### panel 4 ##############################
		panel = fig.add_axes([x1, y2, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data4.T, shading = 'gouraud', cmap=cmap, vmin=vmin1, vmax=vmax1)
		cbar = plt.colorbar(pcm)

		panel.text(-0.55, 0.64, 'numerics',rotation='vertical',transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		######################### panel 5 ##############################
		panel = fig.add_axes([x2, y2, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data5.T, shading = 'gouraud', cmap=cmap, vmin=vmin2, vmax=vmax2)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		#~ ######################### panel 6 ##############################
		panel = fig.add_axes([x4, y2, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data6.T, shading = 'gouraud', cmap=cmap, vmin=vmin3, vmax=vmax3)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		
		
		
		
		
		######################### panel 7 ##############################
		panel = fig.add_axes([x1, y3, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data7.T, shading = 'gouraud', cmap=cmap, vmin=vmin1, vmax=vmax1)
		cbar = plt.colorbar(pcm)

		panel.text(-0.55, 0.82, 'analytics (q=2)',rotation='vertical',transform=panel.transAxes)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		######################### panel 8 ##############################
		panel = fig.add_axes([x2, y3, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data8.T, shading = 'gouraud', cmap=cmap, vmin=vmin2, vmax=vmax2)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		
		#~ ######################### panel 9 ##############################
		panel = fig.add_axes([x4, y3, x4size, 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'$x/R$',fontsize=axislabelfontsize,labelpad=-0.2)
		panel.set_ylabel(r'$y/R$',fontsize=axislabelfontsize,labelpad=-1)
		
		pcm = panel.pcolormesh(xs, ys, data9.T, shading = 'gouraud', cmap=cmap, vmin=vmin3, vmax=vmax3)
		cbar = plt.colorbar(pcm)
		
		panel.set_xlim(-1.15,1.15)
		panel.xaxis.set_major_locator(MultipleLocator(1))
		panel.xaxis.set_minor_locator(MultipleLocator(0.2))
		panel.set_ylim(-1.15,1.15)
		panel.yaxis.set_major_locator(MultipleLocator(1))
		panel.yaxis.set_minor_locator(MultipleLocator(0.2))
		
		circle = plt.Circle((0,0),1.0,facecolor='None',lw=0.8, edgecolor='k')
		panel.add_patch(circle)
		
		

		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.001		# time integration step
	time_duration = 0.4	# 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)
	
	# number of independent realization for the statistic in the simulations
	Nparticles = 10000000
	
	# time snapshots
	t1,t2,t3 = 0.05,0.15,time_duration
	
	# 2D hystogram data
	binning = 0.05
	Nbins = 49			# 1 in zero, Nbins/2 on the left and Nbins/2 on the right
	xs=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	ys=np.asarray([ -(Nbins-1)/2*binning+ binning*i  for i in range(Nbins)])
	
	
	
	#~ data1,data2,data3 = hw.spatial_distribution(D_rot,D,v,dt,time_duration,x0,y0,theta0,Nparticles,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,data5,data6 = hw.spatial_distribution_numerics(Gamm,Pe,r0,phi0,theta0,t1,t2,t3,binning,Nbins)
	#~ print_data (Nbins,binning,data4,data5,data6,'data_from_numerics.dat')
	data4,data5,data6 = read_data (Nbins,binning,'data_from_numerics.dat')
	
	#~ data7,data8,data9 = hw.spatial_distribution_analytics(Gamm,Pe,r0,phi0,theta0,t1,t2,t3,binning,Nbins)
	#~ print_data (Nbins,binning,data7,data8,data9,'data_from_analytics.dat')
	data7,data8,data9 = read_data (Nbins,binning,'data_from_analytics.dat')


	# plots
	#~ make_figure3x2 (xs,ys,data1,data2,data3,data4,data5,data6,t1,t2,t3)
	make_figure3x3 (xs,ys,data1,data2,data3,data4,data5,data6,data7,data8,data9,t1,t2,t3)



main()