# -*- coding: utf-8 -*- """ Created on Fri Dec 9 11:59:18 2022 @author: regin """ import numpy as np import scipy.linalg as splin from scipy.linalg import eig import matplotlib.pyplot as plt v=1. #------------------------------------------------------------------------------------- #for reading the Data from the Paper "Resonant diffusion of a gravitactic circle swimmer" def positions(data): f1 = open(data,'r') lines1 = f1.readlines() f1.close() t=np.array([]) x1=np.array([]) for line in lines1[2:]: p=line.split() t=np.append(t,float(p[0])) x1=np.append(x1,float(p[1])) return t,x1 #------------------------------------------------------------------------------------- num =101; #define size of matrix #L_0 in matrix form def L_mn(Omega,gamma,Drot): L = np.zeros((num,num),dtype=complex); for i in range(0,num): n = (i+1)-round(num/2); #could also be n = (i)-round(num/2); in the limit it doesnt matter L[i,i]=-(1.0j)*Omega*float(n) - Drot*float(n)*float(n); if (i0): L[i-1,i] += -gamma*n/2.; D, W0,V0 = eig(L, left=True) #eigenvalue, left eigenvector, right eigenvector V=V0.copy() W=W0.copy() for i in range(0,num): norm=np.sum(np.conj(W[:,i])*V[:,i]) V[:,i]=V[:,i]/np.sqrt(norm); W[:,i]=W[:,i]/np.conj(np.sqrt(norm)); return (-1.)*D, W, V #(-1.) because of our definition: check np.dot(L,V[:,50])+D[50]*V[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation) #D, W, V=L_mn(Omega,gamma,Drot) #Just to check the relations #idenity=np.zeros((num,num)) #for i in range(num): # idenity=idenity+np.dot(np.transpose(np.conj([W[:,i]])),np.array([V[:,i]])) #print(idenity) #define the Delta L matrices. One for Cosine, Sinus and Total #--------------------------------------------------------------------------------------------- c=np.zeros((num,num),dtype=complex); #for sinus s= np.zeros((num,num),dtype=complex); #for cosinus tot= np.zeros((num,num),dtype=complex); #for totat para= np.zeros((num,num),dtype=complex); #for totat perp= np.zeros((num,num),dtype=complex); #for totat #ACHTUNG: hier ein minus??? for i in range(0,num): if (i0): c[i-1,i]=-v*(0.5j) for i in range(0,num): if (i0): s[i-1,i]=v*(0.5) for i in range(0,num): if (i0): tot[i-1,i]=(0.5-0.5j) def bralLketr(l,n,M,r,m): bra_l=np.conj([l[:,n]]) ket_r=np.transpose([r[:,m]]) return np.matmul(bra_l,np.matmul(M,ket_r))[0,0] def braLket_fill_x(Omega,gamma,Drot): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bralLketr(W,index_l,c,V,index_r) return braLket def braLket_fill_y(Omega,gamma,Drot): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bralLketr(W,index_l,s,V,index_r) return braLket def braLket_fill_para(Omega,gamma,Drot): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) para=delta_L_mn_para(Omega,gamma,Drot) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bralLketr(W,index_l,para,V,index_r) return braLket def braLket_fill_perp(Omega,gamma,Drot): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) perp=delta_L_mn_perp(Omega,gamma,Drot) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bralLketr(W,index_l,perp,V,index_r) return braLket def mean_drift_x_const(Omega,gamma,Drot): braLket_x=braLket_fill_x(Omega,gamma,Drot) return braLket_x[0,0]*1.0j def mean_drift_y_const(Omega,gamma,Drot): braLket_y=braLket_fill_y(Omega,gamma,Drot) return braLket_y[0,0]*1.0j def delta_L_mn_para(Omega,gamma,Drot): #needed if the matrix depends on other variables para= np.zeros((num,num),dtype=complex); #for totat cos=mean_drift_x_const(Omega,gamma,Drot) sin=mean_drift_y_const(Omega,gamma,Drot) for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): para[i-1,i]=(sin*0.5-cos*0.5j)/np.sqrt(cos**2+sin**2) #gravitaxis return para def delta_L_mn_perp(Omega,gamma,Drot): #needed if the matrix depends on other variables perp= np.zeros((num,num),dtype=complex); #for totat cos=mean_drift_x_const(Omega,gamma,Drot) sin=mean_drift_y_const(Omega,gamma,Drot) for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): perp[i-1,i]+=(+cos*0.5+sin*0.5j)/np.sqrt(cos**2+sin**2) #gravitaxis return perp def L_plus_delta_L_mn_para(Omega,gamma,Drot,k): #needed if the matrix depends on other variables L = np.zeros((num,num),dtype=complex); cos=mean_drift_x_const(Omega,gamma,Drot) sin=mean_drift_y_const(Omega,gamma,Drot) for i in range(0,num): n = (i+1)-round(num/2); #could also be n = (i)-round(num/2); in the limit it doesnt matter L[i,i]+=-(1.0j)*Omega*float(n) - Drot*float(n)*float(n); if (i0): L[i-1,i] += -gamma*n/2.; #delta for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): sdgf=1 L[i-1,i] += k*(-0.5j*cos+0.5*sin)/np.sqrt(cos**2+sin**2)#cos #L[i-1,i]+=k*(-0.5) #s D, W0,V0 = eig(L, left=True) #eigenvalue, left eigenvector, right eigenvector V=V0.copy() W=W0.copy() for i in range(0,num): norm=np.sum(np.conj(W[:,i])*V[:,i]) V[:,i]=V[:,i]/np.sqrt(norm); W[:,i]=W[:,i]/np.conj(np.sqrt(norm)); return (-1.)*D, W, V #(-1.) because of our definition: check np.dot(L,V[:,50])+D[50]*V[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation) def L_plus_delta_L_mn_perp(Omega,gamma,Drot,k): #needed if the matrix depends on other variables L = np.zeros((num,num),dtype=complex); cos=mean_drift_x_const(Omega,gamma,Drot) sin=mean_drift_y_const(Omega,gamma,Drot) for i in range(0,num): n = (i+1)-round(num/2); #could also be n = (i)-round(num/2); in the limit it doesnt matter L[i,i]+=-(1.0j)*Omega*float(n) - Drot*float(n)*float(n); #Achtung minus? if (i0): L[i-1,i] += -gamma*n/2.; #delta for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): L[i-1,i] += k*(0.5j*sin+0.5*cos)/np.sqrt(cos**2+sin**2)#cos #L[i-1,i]+=k*(-0.5) #s D, W0,V0 = eig(L, left=True) #eigenvalue, left eigenvector, right eigenvector V=V0.copy() W=W0.copy() for i in range(0,num): norm=np.sum(np.conj(W[:,i])*V[:,i]) V[:,i]=V[:,i]/np.sqrt(norm); W[:,i]=W[:,i]/np.conj(np.sqrt(norm)); return (-1.)*D, W, V #(-1.) because of our definition: check np.dot(L,V[:,50])+D[50]*V[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation) def mean_displacement_x(t,Omega,gamma,Drot): braLket_x=braLket_fill_x(Omega,gamma,Drot) return t*braLket_x[0,0]*1.0j def mean_displacement_y(t,Omega,gamma,Drot): braLket_y=braLket_fill_y(Omega,gamma,Drot) return t*braLket_y[0,0]*1.0j def mean_displacement_perp(t,Omega,gamma,Drot): #D, W, V=L_mn(Omega,gamma,Drot) braLket=braLket_fill_perp(Omega,gamma,Drot) return t*braLket[0,0]*1.0j def mean_displacement_para(t,Omega,gamma,Drot): D, W, V=L_mn(Omega,gamma,Drot) braLket=braLket_fill_para(Omega,gamma,Drot) return t*braLket[0,0]*1.0j def L_plus_delta_L_cos(Omega,gamma,Drot,k): L = np.zeros((num,num),dtype=complex); for i in range(0,num): n = (i+1)-round(num/2); #could also be n = (i)-round(num/2); in the limit it doesnt matter L[i,i]+=-(1.0j)*Omega*float(n) - Drot*float(n)*float(n); if (i0): L[i-1,i]+= -gamma*n/2.; #delta for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): L[i-1,i]+=-k*(0.5j) D, W0,V0 = eig(L, left=True) #eigenvalue, left eigenvector, right eigenvector V=V0.copy() W=W0.copy() for i in range(0,num): norm=np.sum(np.conj(W[:,i])*V[:,i]) V[:,i]=V[:,i]/np.sqrt(norm); W[:,i]=W[:,i]/np.conj(np.sqrt(norm)); return (-1.)*D, W, V #(-1.) because of our definition: check np.dot(L,V[:,50])+D[50]*V[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation) #achtung minus def L_plus_delta_L_sin(Omega,gamma,Drot,k): L = np.zeros((num,num),dtype=complex); for i in range(0,num): n = (i+1)-round(num/2); #could also be n = (i)-round(num/2); in the limit it doesnt matter L[i,i]+=-(1.0j)*Omega*float(n) - Drot*float(n)*float(n); if (i0): L[i-1,i]+= -gamma*n/2.; #delta for i in range(0,num): n = (i+1)-round(num/2); #tot[i,i]=-(1.0j)*Omega*k -2.*Drot*n*k - Drot*k**2 #washborad potential if (i0): L[i-1,i]+=k*(0.5) #s D, W0,V0 = eig(L, left=True) #eigenvalue, left eigenvector, right eigenvector V=V0.copy() W=W0.copy() for i in range(0,num): norm=np.sum(np.conj(W[:,i])*V[:,i]) V[:,i]=V[:,i]/np.sqrt(norm); W[:,i]=W[:,i]/np.conj(np.sqrt(norm)); return (-1.)*D, W, V #(-1.) because of our definition: check np.dot(L,V[:,50])+D[50]*V[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation) def read_gravitaxis_2_2(path): data=np.loadtxt(path,skiprows=1) i, j, count_block, time , md , msd ,ISFx, ISFy, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = [data[:,i] for i in range(0,24)] return [i, j, count_block, time , md , msd ,ISFx, ISFy, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2] def read_gravitaxis_2_isf(path): data=np.loadtxt(path,skiprows=1) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = [data[:,i] for i in range(0,27)] return [i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2] def read_gravitaxis_2_isf_comoving(path): data=np.loadtxt(path,skiprows=1) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = [data[:,i] for i in range(0,39)] return [i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2] ############################################ #OBSERVABLES ############################################ def bra_lketr(l,n,r,m): bra_l=np.conj([l[:,n]]) ket_r=np.transpose([r[:,m]]) return np.matmul(bra_l,ket_r)[0,0] def braLketr_isf_x(Omega,gamma,Drot,k): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) Dnk, Wnk, Vnk =L_plus_delta_L_cos(Omega,gamma,Drot,k) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bra_lketr(W,0,Vnk,index_l)*bra_lketr(Wnk,index_r,V,0) return braLket def braLketr_isf_y(Omega,gamma,Drot,k): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) Dnk, Wnk, Vnk =L_plus_delta_L_sin(Omega,gamma,Drot,k) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bra_lketr(W,0,Vnk,index_l)*bra_lketr(Wnk,index_r,V,0) return braLket def braLketr_isf_para(Omega,gamma,Drot,k): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) Dnk, Wnk, Vnk =L_plus_delta_L_mn_para(Omega,gamma,Drot,k) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bra_lketr(W,0,Vnk,index_l)*bra_lketr(Wnk,index_r,V,0) return braLket def braLketr_isf_perp(Omega,gamma,Drot,k): braLket=np.zeros((num,num),dtype=np.complex_) D, W, V=L_mn(Omega,gamma,Drot) Dnk, Wnk, Vnk =L_plus_delta_L_mn_perp(Omega,gamma,Drot,k) for index_l in range(num): for index_r in range(num): braLket[index_l][index_r]=bra_lketr(W,0,Vnk,index_l)*bra_lketr(Wnk,index_r,V,0) return braLket def ISF_x(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_cos(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_x(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.real(isf) def ISF_y(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_sin(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_y(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.real(isf) def ISF_para(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_mn_para(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_para(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.real(isf) def ISF_perp(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_mn_perp(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_perp(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.real(isf) def ISF_x_imag(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_cos(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_x(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.imag(isf) def ISF_y_imag(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_sin(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_y(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.imag(isf) def ISF_para_imag(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_mn_para(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_para(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.imag(isf) def ISF_perp_imag(t,Omega,gamma,Drot,k): Dnk, Wnk, Vnk =L_plus_delta_L_mn_perp(Omega,gamma,Drot,k) isf=0+0j for index_l,lamb in enumerate(Dnk): isf+=braLketr_isf_perp(Omega,gamma,Drot,k)[index_l,index_l]*np.exp((-lamb)*t) return np.imag(isf) def save_x(gamma,Drot,t,k): Omega=1.0 savearray=ISF_x(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_y(gamma,Drot,t,k): Omega=1.0 savearray=ISF_y(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_y_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_y_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_para(gamma,Drot,t,k): Omega=1.0 savearray=ISF_para(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_perp(gamma,Drot,t,k): Omega=1.0 savearray=ISF_perp(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_x_omega(Omega,gamma,Drot,t,k): savearray=ISF_x(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_y_omega(Omega,gamma,Drot,t,k): savearray=ISF_y(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_y_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_y_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_para_omega(Omega,gamma,Drot,t,k): savearray=ISF_para(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_para_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_para_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_perp_omega(Omega,gamma,Drot,t,k): savearray=ISF_perp(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_perp_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_perp_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_x_omega_imag(Omega,gamma,Drot,t,k): savearray=ISF_x_imag(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_y_omega_imag(Omega,gamma,Drot,t,k): savearray=ISF_y_imag(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_y_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_imag_y_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_para_omega_imag(Omega,gamma,Drot,t,k): savearray=ISF_para_imag(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) def save_perp_omega_imag(Omega,gamma,Drot,t,k): savearray=ISF_perp_imag(t,Omega,gamma,Drot,k) np.save("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".npy", savearray,allow_pickle=True, fix_imports=True) np.savetxt("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(num)+".txt", savearray) t=np.logspace(-3,3, 1000) for gamma in [0.001,0.25, 0.5, 0.0, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3,0.75,1.25,1.5]: for k in [0.25, 1.0, 5.0]: for D in [0.25]: refina=1 #save_y(gamma,D,t,k) # for gamma in [0.001,0.25,0.5,0.75,1.0, 1.25,1.5]: for k in [0.25]: for D in [0.025]: mega=1 #save_y(gamma,D,t,k) #save_x(gamma,D,t,k) #save_para(gamma,D,t,k) #save_perp(gamma,D,t,k) for gamma in [1.5]: for k in [0.01,0.05,0.1,0.25,0.5,1.0,2.0,5.0,11.0]: for D in [0.025]: mega=1 #save_y(gamma,D,t,k) #save_x(gamma,D,t,k) #save_para(gamma,D,t,k) #save_perp(gamma,D,t,k) for gamma in [1.5]: for k in [0.25]: for D in [0.05]: mega=1 #save_x(gamma,D,t,k) for gamma in [ 1.25]: for k in [5.0,1.0,0.25]: for D in [0.025]: Omega=1 #save_x_omega_imag(Omega,gamma,D,t,k) #save_y_omega_imag(Omega,gamma,D,t,k) #save_para_omega_imag(Omega,gamma,D,t,k) #save_perp_omega_imag(Omega,gamma,D,t,k) for gamma in [1.5]: for k in [0.01,0.05,0.1,0.25,0.5,1.0,2.0,5.0,11.0]: for D in [0.025]: Omega=1 #save_x_omega_imag(Omega,gamma,D,t,k) #save_y_omega_imag(Omega,gamma,D,t,k) #save_para_omega_imag(Omega,gamma,D,t,k) #save_perp_omega_imag(Omega,gamma,D,t,k) for gamma in [0.001 ]: for k in [0.25, 1.0, 5.0,112.0]: for Omega in [0.0, 5.0, 50.3]: D=1 #save_x_omega(Omega,gamma,D,t,k) #save_y_omega(Omega,gamma,D,t,k) #save_para_omega(Omega,gamma,D,t,k) #save_perp_omega(Omega,gamma,D,t,k) for gamma in [1.5, 1.25]: for k in [0.25]: for D in [0.025]: Omega=1 #save_x_omega_imag(Omega,gamma,D,t,k) #save_x(gamma,D,t,k) #save_y_omega_imag(Omega,gamma,D,t,k) #save_y(gamma,D,t,k) #save_para_omega(Omega,gamma,D,t,k) #save_perp_omega(Omega,gamma,D,t,k) #save_para_omega_imag(Omega,gamma,D,t,k) #save_perp_omega_imag(Omega,gamma,D,t,k) ''' save_y(0.8,0.25,t,1.0) save_y(0.9,0.25,t,0.25) save_y(1.0,0.25,t,0.25) save_y(1.1,0.25,t,0.25) save_y(1.2,0.25,t,0.25) save_y(1.3,0.25,t,0.25) save_y(0.25,0.25,t,0.25) ''' #save_x(0.0,0.005,t,0.25) def harmonic_approximation_2_x(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=k ny=0.0 v=1.0 tau=1/(gamma*np.cos(vstar) ) #Dn=(v*tau)**2*Drot*(nx*Omega/gamma-ny*1/gamma/tau)**2 Dn=(v*tau)**2*Drot*(nx*np.sin(vstar)-ny*np.cos(vstar))**2 return 2.0*Dn*tau*(t/tau-1+np.exp(-t/tau)) def harmonic_approximation_2_y(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=0.0 ny=k v=1.0 tau=1/(gamma*np.cos(vstar) ) #Dn=(v*tau)**2*Drot*(nx*Omega/gamma+ny*1/gamma/tau)**2 Dn=(v*tau)**2*Drot*(nx*np.sin(vstar)-ny*np.cos(vstar))**2 return 2.0*Dn*tau*(t/tau-1+np.exp(-t/tau)) def harmonic_approximation_2_para(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=k*np.cos(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) ny=k*np.sin(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) v=1.0 tau=1/(gamma*np.cos(vstar) ) #tau=1/np.sqrt(gamma**2-Omega**2) #Dn=(v*tau)**2*Drot*(nx*Omega/gamma-ny*1/gamma/tau)**2 Dn=(v*tau)**2*Drot*(nx*np.sin(vstar)-ny*np.cos(vstar))**2 return 2.0*Dn*tau*(t/tau-1+np.exp(-t/tau)) def harmonic_approximation_2_perp(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=-k*np.sin(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) ny=k*np.cos(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) v=1.0 tau=1/(gamma*np.cos(vstar) ) #tau=1/np.sqrt(gamma**2-Omega**2) #Dn=(v*tau)**2*Drot*(nx*Omega/gamma-ny*1/gamma/tau)**2 Dn=(v*tau)**2*Drot*(nx*np.sin(vstar)-ny*np.cos(vstar))**2 return 2.0*Dn*tau*(t/tau-1+np.exp(-t/tau)) def HA_mean_x(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=k ny=0.0 v=1.0 x=np.cos(vstar)*t*v y=np.sin(vstar)*t*v return x*nx+y*ny def HA_mean_y(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=0.0 ny=k v=1.0 x=np.cos(vstar)*t*v y=np.sin(vstar)*t*v return x*nx+y*ny def HA_mean_para(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=k*np.cos(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) ny=k*np.sin(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) v=1.0 x=np.cos(vstar)*t*v y=np.sin(vstar)*t*v return x*nx+y*ny def HA_mean_perp(t,Omega,gamma,Drot,k): vstar=np.arcsin(Omega/gamma) nx=-k*np.sin(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) ny=k*np.cos(vstar)/np.sqrt(np.cos(vstar)**2+np.sin(vstar)**2) v=1.0 x=np.cos(vstar)*t*v y=np.sin(vstar)*t*v return x*nx+y*ny def plot_isf_HA(Omega,Drot,gamma): b=100 e=1000 v=1 Deff= v*(Drot)/(2*((Drot)**2+v**2)) size=8 fig = plt.figure(figsize=(3.375, 5.8), layout="constrained") spec = fig.add_gridspec(3, 1) ax0 = fig.add_subplot(spec[0, :]) ax10 = fig.add_subplot(spec[1, :]) ax20 = fig.add_subplot(spec[2, :]) #ax11 = fig.add_subplot(spec[1, 1]) #ax21 = fig.add_subplot(spec[2, 1]) ax0.tick_params(axis='both', which='both', left='on', right='on', bottom='on',top='on') ax10.tick_params(axis='both', which='both', left='on', right='on', bottom='on',top='on') ax20.tick_params(axis='both', which='both', left='on', right='on', bottom='on',top='on') #ax11.tick_params(axis='both', which='both', left='on', right='on', bottom='on',top='on') #ax21.tick_params(axis='both', which='both', left='on', right='on', bottom='on',top='on') ax0.text(0.03, 0.88, '(a)', transform=ax0.transAxes, fontsize=size, verticalalignment='top') ax10.text(0.03, 0.88, '(b)', transform=ax10.transAxes, fontsize=size, verticalalignment='top') ax20.text(0.03, 0.88, '(c)', transform=ax20.transAxes, fontsize=size, verticalalignment='top') #ax11.text(0.03, 0.88, '(d)', transform=ax11.transAxes, fontsize=20, verticalalignment='top') #ax21.text(0.03, 0.88, '(e)', transform=ax21.transAxes, fontsize=20, verticalalignment='top') params=[[11.0,'s','#ff00aa',r'$11.0$'] ,[5.0,'o','#aa00ff',r'$5.0$'] ,[2.0,'^','#0000ff',r'$2.0$'],[1.0,'D','#00aaff',r'$1.0$'],[0.5,'v','#00ffd4',r'$0.5$'],[0.25,'d','#00ff2b',r'$0.25$'],[0.1,'*','#aaff00',r'$0.1$'],[0.05,'<','#ffaa00',r'$0.05$'],[0.01,'>','#ff0000',r'$0.01$']] #params=[[2.0,'^','#0000ff',r'$2.0$'],[1.0,'D','#00aaff',r'$1.0$']] #params=[[5.0,'d','#aa00ff',r'$5.0$'],[11.0,'d','#ff00aa',r'$11.0$']] ax20.plot(0, 0,color='black', linestyle='--', label='HA') ax10.plot(0, 0,color='black', linestyle='--', label='HA') ax0.plot(0, 0,color='black', linestyle='--', label='HA') for k, marker,color,label in params: var= np.load("Variance_analytics\Var_x_D_"+str(Drot)+"_g_"+str(gamma)+"_num_"+str(101)+".npy") isf= np.load("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_x_imag= np.load("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") #path_isf="..\data\Data_gravitaxis_2_isf_comoving_HA_D_0.025_g_1.5"+"\Data_gravitaxis_2_isf_comoving"+"_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" path_isf="..\data\Data_gravitaxis_2_isf_comoving_HA_D_0.025_g_1.5"+"\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(1.5)+"\data_gravitaxis_2_isf.dat" #path_isf="..\data\Data_gravitaxis_2_isf_comoving_self"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(1.5)+"_k_"+str(k)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2= data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23], data[24], data[25], data[26], data[27], data[28], data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] #i, j, count_block, count_block2, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s,ISFpara_s, ISFperp_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23], data[24], data[25], data[26] ax0.set_ylabel(r"$\textsf{Re}[F(t,\mathbf{k})]$") ax0.set_xscale('log') ax0.plot(t[b:e], np.real(np.exp(-(0+1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), color='black',linestyle='--') #ax0.plot(time_s_isf, ISFx_s ,color='red', marker='o', linestyle='none') ax0.plot(time, ISFx ,color=color, marker=marker, linestyle='none') ax0.plot(t[b:e],isf[b:e],color=color,linestyle='-') ax0.set_xlabel(r'$\omega t$') ax10.set_ylabel(r"$\textsf{Re}[\exp(ikv_x t) F(t,\mathbf{k})]$") #ax10.set_ylabel(r"$\exp(-k^2\textsf{Var} [\langle \Delta x (t) \rangle ] /2)$") ax10.plot(time, ISFx_comoving ,color=color, marker=marker, linestyle='none') ax10.set_xscale('log') ax10.plot(t[b:e], np.exp(-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k)), linestyle='--', color='black') comoving_factor=np.exp((0+1j)*k*mean_displacement_x(t[b:e],Omega,gamma,Drot)) ax10.plot(t[b:e],np.real((isf[b:e]+(1j)*isf_x_imag[b:e])*comoving_factor),color=color,linestyle='-') comoving_factor=np.exp((0+1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)) ax10.plot(t[b:e], np.real(np.exp(-(0+1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))*comoving_factor), color='black', linestyle='--') ax10.set_xlabel(r'$\omega t$') ax20.set_ylabel(r"$\textsf{Re}[\exp(ikv_x t) F(t,\mathbf{k})]$") #ax20.set_ylabel(r"$\exp(-k^2\textsf{Var} [\langle \Delta x (t) \rangle ] /2)$") ax20.set_xscale('log') ax20.plot(time*k**2, ISFx_comoving ,color=color, marker=marker, linestyle='none') ax20.plot(t[b:e]*k**2, np.real(np.exp(-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), linestyle='--', color='black') ax20.set_xlabel(r'$k^2 v^2 t/ \omega $') comoving_factor=np.exp((0+1j)*k*mean_displacement_x(t[b:e],Omega,gamma,Drot)) ax20.plot(t[b:e]*k**2,np.real((isf[b:e]+(1j)*isf_x_imag[b:e])*comoving_factor),color=color,linestyle='-') ax20.set_xlabel(r'$k^2 v^2 t/ \omega $') ax20.plot(0, 0,label=str(label),color=color, marker=marker) ax10.plot(0, 0,label=str(label),color=color, marker=marker) ax0.plot(0, 0,label=str(label),color=color, marker=marker) ax20.legend(title=r'$k v/\omega$',loc='lower left') ax10.legend(title=r'$k v/\omega$',loc='lower left') ax0.legend(title=r'$k v/\omega$',loc='lower left') def plot_reg_sim_isf_mit_comoving(Omega,Drot): #path2='..\code\data_gravitaxis_2.dat' #params=[[1.5,'s','#750787',r'$1.5$'],[1.25,'s','#004dff',r'$1.25$'],[1.0,'o','#008026',r'$1.0$'],[0.75,'v','#ffed00',r'$0.75$'],[0.5,'v','#ff8c00',r'$0.5$'],[0.25,'d','#e40303',r'$0.25$'],[0.001,'d','black',r'$0.0$']] #params=[[1.5,'s','#90ee90',r'$1.5$'],[1.25,'s','#82d682',r'$1.25$'],[1.0,'o','#65a765',r'$1.0$'],[0.75,'v','#487748',r'$0.75$'],[0.5,'v','#3a5f3a',r'$0.5$'],[0.25,'d','#1d301d',r'$0.25$'],[0.001,'d','#0e180e',r'$0.0$']] #params=[[0.001,'d','#009900',r'$0.0$'],[0.25,'d','#009900',r'$0.25$'],[0.5,'v','#bf00ff',r'$0.5$'],[0.7,'v','#bf00ff',r'$0.7$'],[0.8,'D','deepskyblue',r'$0.8$'],[0.9,'^', '#ffc40c',r'$0.9$'],[1.0,'o','r',r'$1.0$'],[1.2,'s','b',r'$1.2$']] #params=[[0.25,'d','#009900',r'$0.25$'],[0.7,'v','#bf00ff',r'$0.7$'],[0.8,'D','deepskyblue',r'$0.8$'],[0.9,'^', '#ffc40c',r'$0.9$'],[1.0,'o','r',r'$1.0$'],[1.1,'s','b',r'$1.1$']] # Sample data for your plots letters=[['(a)','(b)','(c)'],['(d)','(e)','(f)'], ['(g)','(h)','(i)'],['(j)','(k)','(l)']] fig, axes = plt.subplots(4, 3, figsize=(18, 16)) # Adjust the figure size as needed for i in range(4): for j in range(3): #axes[i,j].grid() if i != 3: axes[i,j].tick_params( axis='x', # changes apply to the x-axis labelbottom=False) if j != 0: axes[i,j].tick_params( axis='y', # changes apply to the x-axis labelleft=False) axes[i,j].tick_params( axis='both', # changes apply to the x-axis which='both', # both major and minor ticks are affected left='on', # ticks along the left edge are off right='on', # ticks along the right edge are off bottom='on', # ticks along the bottom edge are off top='on') # ticks along the top edge are off) ) axes[i,j].text(0.03, 0.88, letters[i][j], transform=axes[i,j].transAxes, fontsize=20, verticalalignment='top') # Plot on each subplot params=[[1.5,'s','#bae604',r'$1.5$'],[1.25,'s','#a6cc03',r'$1.25$'],[1.0,'o','#7c9902',r'$1.0$'],[0.75,'v','#536602',r'$0.75$'],[0.5,'v','#536602',r'$0.5$'],[0.25,'d','#293301',r'$0.25$'],[0.001,'d','#151900',r'$0.0$']] for gamma, marker,color,label in params: axes[3,0].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[3,0].legend(title=r'$\gamma/\omega$') b=350 e=10000 # Create a figure with a 2x2 grid of subplots #path2='..\code\data_gravitaxis_2.dat' #i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block,count_block2, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s, ISFpara_s, ISFperp_s,cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26] k=0.25 #if gamma==0.001: # gamma_s=0.0 # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(0.25)+"\data_gravitaxis_2.dat" #else: # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(0.25)+"\data_gravitaxis_2.dat" # data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] path_isf="..\data\Data_gravitaxis_2_isf_comoving_1\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] t=np.logspace(-3,3, 1000) #axes[0,0].plot(time_s,ISFx_s,color=color, marker=marker, linestyle='none') axes[0,0].set_ylabel(r"$\textsf{Re}[\exp(ikv_x t) F(t,\mathbf{k}_x)]$") axes[0,0].set_title(r"$kv/\omega=0.25$") #axes[0,0].set_yscale('symlog') axes[0,0].set_xscale('log') #if gamma > Omega: # axes[0,0].plot(t[b:e], np.real(np.exp(-(0+1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') isf_x_imag= np.load("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_x= np.load("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_x(t[b:e],Omega,gamma,Drot)) axes[0,0].plot(t[b:e],np.real((isf_x[b:e]+(0+1j)*isf_x_imag[b:e])*comoving_factor),color=color) Deff= (Drot)/(2*((Drot)**2+1)) axes[0,0].plot(t[b:e],np.exp(-k**2*t[b:e]*Deff), linestyle='dotted', color='red') axes[0,0].plot(time,ISFx_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_y(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_y(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[1,0].plot(time_s,ISFy_s,color=color, marker=marker, linestyle='none') axes[1,0].set_ylabel(r"$\textsf{Re}[\exp(ikv_yt )\textsf{ ISF}(t,\mathbf{k}_y)]$") axes[1,0].set_xscale('log') isf_y= np.load("ISF_analytics\ISF_y_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_y_imag= np.load("ISF_analytics\ISF_y_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_y(t[b:e],Omega,gamma,Drot)) axes[1,0].plot(t[b:e],np.real((isf_y[b:e]+(0+1j)*isf_y_imag[b:e])*comoving_factor),color=color) axes[1,0].plot(time,ISFy_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[2,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,0].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') axes[2,0].set_ylabel(r"$\textsf{Re}[\exp(ik v_\parallel t )F(t,\mathbf{k}_{\parallel})]$") axes[2,0].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,0].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,0].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,0].plot(time_s_isf,ISFperp_s,color=color, marker=marker, linestyle='none') axes[3,0].set_ylabel(r"$\textsf{Re}[F(t,\mathbf{k}_{\perp})]$") axes[3,0].set_xlabel(r'$\omega t$') axes[3,0].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,0].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[3,0].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') k=1.0 params=[[1.5,'s','#00ffff',r'$1.5$'],[1.25,'s','#00aeff',r'$1.25$'],[1.0,'o','#00b4d8',r'$1.0$'],[0.75,'v','#0096c7',r'$0.75$'],[0.5,'v','#0077b6',r'$0.5$'],[0.25,'d','#023e8a',r'$0.25$'],[0.001,'d','#03045e',r'$0.0$']] for gamma, marker,color,label in params: path_isf="..\data\Data_gravitaxis_2_isf_comoving_1\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] b=350 e=1000 axes[3,1].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[3,1].legend(title=r'$\gamma/\omega$') #path_isf="..\code\data_gravitaxis_2_isf.dat" #data=read_gravitaxis_2_isf(path_isf) #i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, count_block2, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s,ISFpara_s, ISFperp_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23], data[24], data[25], data[26] #i, j, count_block, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s,cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] #if gamma==0.001: # gamma_s=0.0 # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(1.0)+"\data_gravitaxis_2.dat" #else: # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(1.0)+"\data_gravitaxis_2.dat" #data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] #if gamma > Omega: # axes[0,1].plot(t[b:e], np.real(np.exp(-(0-1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[0,1].plot(time_s_isf,ISFx_s,color=color, marker=marker, linestyle='none') axes[0,1].set_title(r"$kv/\omega=1.0$") #axes[0,1].set_yscale('symlog') axes[0,1].set_xscale('log') #axes[0,1].set_ylim(10**(-7),10**(3)) isf_x_imag= np.load("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_x= np.load("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*mean_displacement_x(t[b:e],Omega,gamma,Drot)) axes[0,1].plot(t[b:e],np.real((isf_x[b:e]+(0+1j)*isf_x_imag[b:e])*comoving_factor),color=color) axes[0,1].plot(time,ISFx_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_y(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_y(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[1,1].plot(time_s_isf,ISFy_s,color=color, marker=marker, linestyle='none') #axes[1,0].set_ylim(10**(-9),10**(3)) isf_y= np.load("ISF_analytics\ISF_y_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_y_imag= np.load("ISF_analytics\ISF_y_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_y(t[b:e],Omega,gamma,Drot)) axes[1,1].plot(t[b:e],np.real((isf_y[b:e]+(0+1j)*isf_y_imag[b:e])*comoving_factor),color=color) #axes[1,1].set_ylim(10**(-9),10**(3)) # axes[1,1].set_yscale('symlog') axes[1,1].set_xscale('log') axes[1,1].plot(time,ISFy_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[2,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,1].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') #axes[2,1].plot(time_s_isf,ISFx_s*ISFy_s,color=color, marker=marker, linestyle='none') axes[2,1].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,1].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,1].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,1].plot(time_s_isf,ISFperp_s,color=color, marker=marker, linestyle='none') axes[3,1].set_xlabel(r'$\omega t$') axes[3,1].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,1].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[3,1].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') k=5.0 params=[[1.5,'s','#ffa9ab',r'$1.5$'],[1.25,'s','#ff575c',r'$1.25$'],[1.0,'o','#ff4d4d',r'$1.0$'],[0.75,'v','#ff060d',r'$0.75$'],[0.5,'v','#b30005',r'$0.5$'],[0.25,'d','#610002',r'$0.25$'],[0.001,'d','#100000',r'$0.0$']] for gamma, marker,color,label in params: b=250 e=1000 axes[3,2].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[3,2].legend(title=r'$\gamma/\omega$') ''' if gamma==0.001: gamma_s=0.0 path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(5.0)+"\data_gravitaxis_2.dat" else: path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(5.0)+"\data_gravitaxis_2.dat" data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] ''' path_isf="..\data\Data_gravitaxis_2_isf_comoving_1\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] #if gamma > Omega: # axes[0,2].plot(t[b:e], np.real(np.exp(-(0-1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') axes[0,2].set_title(r"$kv/\omega=5.0$") #axes[0,2].set_yscale('symlog') axes[0,2].set_xscale('log') #axes[0,2].set_ylim(10**(-7),10**(3)) isf_x_imag= np.load("ISF_analytics\ISF_x_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_x= np.load("ISF_analytics\ISF_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_x(t[b:e],Omega,gamma,Drot)) axes[0,2].plot(t[b:e],np.real((isf_x[b:e]+(0+1j)*isf_x_imag[b:e])*comoving_factor),color=color) #axes[0,2].plot(time_s,ISFx_s,color=color, marker=marker, linestyle='none') axes[0,2].plot(time,ISFx_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_y(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_y(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') isf_y_imag= np.load("ISF_analytics\ISF_y_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_y= np.load("ISF_analytics\ISF_y_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_y(t[b:e],Omega,gamma,Drot)) axes[1,2].plot(t[b:e],np.real((isf_y[b:e]+(0+1j)*isf_y_imag[b:e])*comoving_factor),color=color) #axes[1,1].set_ylim(10**(-9),10**(3)) #axes[1,2].set_yscale('symlog') axes[1,2].set_xscale('log') axes[1,2].plot(time,ISFy_comoving, color=color, marker=marker, linestyle='none') #axes[1,2].plot(time_s,ISFy_s,color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_y(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_y(t[b:e],Omega,gamma,Drot,k))), linestyle=(0, (5, 10)), color='black') #if gamma > Omega: # axes[2,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,2].plot(time_s,ISFx_s,color=color, marker=marker, linestyle='none') axes[2,2].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,2].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,2].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,2].plot(time_s,ISFy_s,color=color, marker=marker, linestyle='none') axes[3,2].set_xlabel(r'$\omega t$') axes[3,2].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,2].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[3,2].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') axes[0,0].set_ylim(-0.5,1.1) axes[0,1].set_ylim(-0.5,1.1) axes[0,2].set_ylim(-0.5,1.1) axes[1,0].set_ylim(-0.75,1.1) axes[1,1].set_ylim(-0.75,1.1) axes[1,2].set_ylim(-0.75,1.1) axes[2,0].set_ylim(-0.75,1.1) axes[2,1].set_ylim(-0.75,1.1) axes[2,2].set_ylim(-0.75,1.1) axes[3,0].set_ylim(-0.48,1.1) axes[3,1].set_ylim(-0.48,1.1) axes[3,2].set_ylim(-0.48,1.1) def plot_reg_sim_isf_mit_comoving_imag(Omega,Drot): letters=[['(a)','(b)','(c)'],['(d)','(e)','(f)'], ['(g)','(h)','(i)'],['(j)','(k)','(l)']] fig, axes = plt.subplots(4, 3, figsize=(6.75, 6.9)) # Adjust the figure size as needed size=8 axes[0,0].set_title(r"$kv/\omega=0.25$",fontsize=size) axes[0,1].set_title(r"$kv/\omega=1.0$",fontsize=size) axes[0,2].set_title(r"$kv/\omega=5.0$",fontsize=size) for i in range(4): for j in range(3): #axes[i,j].grid() if i != 3: axes[i,j].tick_params( axis='x', # changes apply to the x-axis labelbottom=False) if j != 0: axes[i,j].tick_params( axis='y', # changes apply to the x-axis labelleft=False) axes[i,j].tick_params( axis='both', # changes apply to the x-axis which='both', # both major and minor ticks are affected left='on', # ticks along the left edge are off right='on', # ticks along the right edge are off bottom='on', # ticks along the bottom edge are off top='on') # ticks along the top edge are off) ) axes[i,j].text(0.85, 0.20, letters[i][j], transform=axes[i,j].transAxes, fontsize=size, verticalalignment='top') params=[[0.001,'d','#100000',r'$0.0$'],[0.25,'d','#610002',r'$0.25$'],[0.5,'v','#b30005',r'$0.5$'],[0.75,'v','#ff060d',r'$0.75$'],[1.0,'o','#ff4d4d',r'$1.0$'],[1.25,'s','#ff575c',r'$1.25$'],[1.5,'s','#ffa9ab',r'$1.5$']] for gamma, marker,color,label in params: axes[0,2].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[0,2].legend(title=r'$\gamma/\omega$') params=[[0.001,'d','#03045e',r'$0.0$'],[0.25,'d','#023e8a',r'$0.25$'],[0.5,'v','#0077b6',r'$0.5$'],[0.75,'v','#0096c7',r'$0.75$'],[1.0,'o','#00b4d8',r'$1.0$'],[1.25,'s','#00aeff',r'$1.25$'],[1.5,'s','#00ffff',r'$1.5$']] for gamma, marker,color,label in params: axes[0,1].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[0,1].legend(title=r'$\gamma/\omega$') params=[[0.001,'d','#151900',r'$0.0$'],[0.25,'d','#293301',r'$0.25$'],[0.5,'v','#536602',r'$0.5$'],[0.75,'v','#536602',r'$0.75$'],[1.0,'o','#7c9902',r'$1.0$'],[1.25,'s','#a6cc03',r'$1.25$'],[1.5,'s','#bae604',r'$1.5$']] for gamma, marker,color,label in params: axes[0,0].plot(0, 0,marker=marker,label=str(label),color=color) #for legend axes[0,0].legend(title=r'$\gamma/\omega$') # Plot on each subplot params=[[1.5,'s','#bae604',r'$1.5$'],[1.25,'s','#a6cc03',r'$1.25$'],[1.0,'o','#7c9902',r'$1.0$'],[0.75,'v','#536602',r'$0.75$'],[0.5,'v','#536602',r'$0.5$'],[0.25,'d','#293301',r'$0.25$'],[0.001,'d','#151900',r'$0.0$']] for gamma, marker,color,label in params: b=350 e=10000 # Create a figure with a 2x2 grid of subplots #path2='..\code\data_gravitaxis_2.dat' #i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block,count_block2, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s, ISFpara_s, ISFperp_s,cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26] k=0.25 axes[0,0].set_ylabel(r"$\textsf{Re}[\exp(ik \mathbf{n}_{\parallel} \cdot \mathbf{v} t )F(t,k\mathbf{n}_{\parallel})]$") axes[1,0].set_ylabel(r"$\textsf{Re}[F(t,k\mathbf{n}_{\perp})]$") axes[2,0].set_ylabel(r"$\textsf{Im}[\exp(ik \mathbf{n}_{\parallel} \cdot \mathbf{v} t )F(t,k\mathbf{n}_{\parallel})]$") axes[3,0].set_ylabel(r"$\textsf{Im}[F(t,k\mathbf{n}_{\perp})]$") axes[0,0].set_ylim(-0.65,1.1) axes[0,1].set_ylim(-0.65,1.1) axes[0,2].set_ylim(-0.65,1.1) axes[1,0].set_ylim(-0.5,1.1) axes[1,1].set_ylim(-0.5,1.1) axes[1,2].set_ylim(-0.5,1.1) axes[2,0].set_ylim(-0.75,0.51) axes[2,1].set_ylim(-0.75,0.51) axes[2,2].set_ylim(-0.75,0.51) axes[3,0].set_ylim(-0.15,0.025) axes[3,1].set_ylim(-0.15,0.025) axes[3,2].set_ylim(-0.15,0.025) #if gamma==0.001: # gamma_s=0.0 # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(0.25)+"\data_gravitaxis_2.dat" #else: # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(0.25)+"\data_gravitaxis_2.dat" # data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] #path_isf="..\data\Data_gravitaxis_2_isf_comoving_1\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" path_isf="..\data\Data_gravitaxis_2_isf_comoving_D_0.025\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] t=np.logspace(-3,3, 1000) #if gamma > Omega: # axes[2,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[0,0].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') axes[0,0].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[0,0].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[0,0].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') Deff= (Drot)/(2*((Drot)**2+1)) axes[0,0].plot(t[b:e],np.exp(-k**2*t[b:e]*Deff), linestyle='dashed', color='red') axes[1,0].plot(t[b:e],np.exp(-k**2*t[b:e]*Deff), linestyle='dashed', color='red') axes[1,0].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") #if gamma > 1.1: # isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(201)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[1,0].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[1,0].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[2,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,0].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') axes[2,0].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,0].plot(t[b:e],np.imag((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,0].plot(time,ISFpara_comoving_imag, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,0].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,0].plot(time_s_isf,ISFperp_s,color=color, marker=marker, linestyle='none') axes[3,0].set_xlabel(r'$\omega t$') axes[3,0].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,0].plot(t[b:e],np.imag((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[3,0].plot(time,ISFperp_comoving_imag, color=color, marker=marker, linestyle='none') k=1.0 params=[[1.5,'s','#00ffff',r'$1.5$'],[1.25,'s','#00aeff',r'$1.25$'],[1.0,'o','#00b4d8',r'$1.0$'],[0.75,'v','#0096c7',r'$0.75$'],[0.5,'v','#0077b6',r'$0.5$'],[0.25,'d','#023e8a',r'$0.25$'],[0.001,'d','#03045e',r'$0.0$']] for gamma, marker,color,label in params: path_isf="..\data\Data_gravitaxis_2_isf_comoving_D_0.025\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] b=350 e=1000 #path_isf="..\code\data_gravitaxis_2_isf.dat" #data=read_gravitaxis_2_isf(path_isf) #i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, count_block2, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s,ISFpara_s, ISFperp_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23], data[24], data[25], data[26] #i, j, count_block, time_s_isf , md_s , msd_s ,ISFx_s, ISFy_s,cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] #if gamma==0.001: # gamma_s=0.0 # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(1.0)+"\data_gravitaxis_2.dat" #else: # path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(1.0)+"\data_gravitaxis_2.dat" #data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] #i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] #if gamma > Omega: # axes[0,1].plot(t[b:e], np.real(np.exp(-(0-1j)*HA_mean_x(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_x(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #if gamma > Omega: # axes[2,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[0,1].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') #axes[0,1].plot(time_s_isf,ISFx_s*ISFy_s,color=color, marker=marker, linestyle='none') axes[0,1].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[0,1].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[0,1].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[1,1].plot(time_s_isf,ISFperp_s,color=color, marker=marker, linestyle='none') axes[1,1].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[1,1].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[1,1].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[2,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,1].plot(time_s_isf,ISFpara_s,color=color, marker=marker, linestyle='none') #axes[2,1].plot(time_s_isf,ISFx_s*ISFy_s,color=color, marker=marker, linestyle='none') axes[2,1].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,1].plot(t[b:e],np.imag((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,1].plot(time,ISFpara_comoving_imag, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,1].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,1].plot(time_s_isf,ISFperp_s,color=color, marker=marker, linestyle='none') axes[3,1].set_xlabel(r'$\omega t$') axes[3,1].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,1].plot(t[b:e],np.imag((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[3,1].plot(time,ISFperp_comoving_imag, color=color, marker=marker, linestyle='none') k=5.0 params=[[1.5,'s','#ffa9ab',r'$1.5$'],[1.25,'s','#ff575c',r'$1.25$'],[1.0,'o','#ff4d4d',r'$1.0$'],[0.75,'v','#ff060d',r'$0.75$'],[0.5,'v','#b30005',r'$0.5$'],[0.25,'d','#610002',r'$0.25$'],[0.001,'d','#100000',r'$0.0$']] for gamma, marker,color,label in params: b=250 e=1000 ''' if gamma==0.001: gamma_s=0.0 path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma_s)+"_k_"+str(5.0)+"\data_gravitaxis_2.dat" else: path2="..\data\Data_gravitaxis_2_6_isf"+"\Data_observables_data_gravitaxis_2_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(5.0)+"\data_gravitaxis_2.dat" data=read_gravitaxis_2_2(path2) ##i, j, count_block, time_s , md_s , msd_s ,msd_s, cos_s, sin_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17] i, j, count_block, time_s , md_s , msd_s ,ISFx_s, ISFy_s, cos_theta_v_s, sin_theta_v_s, rx_s ,ry_s ,rx2_s ,ry2_s ,rxry_s, rx2ry_s, rxry2_s,rx3_s ,ry3_s, rx4_s, ry4_s, rx3ry_s, ry3rx_s, rx2ry2_s = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23] ''' path_isf="..\data\Data_gravitaxis_2_isf_comoving_D_0.025\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" #path_isf="..\data\Data_gravitaxis_2_isf_comoving_1\Data_gravitaxis_2_isf_comoving_k_"+str(k)+"_D_"+str(Drot)+"_g_"+str(gamma)+"\data_gravitaxis_2_isf.dat" data=read_gravitaxis_2_isf_comoving(path_isf) i, j, count_block, count_block_2, time , md , msd ,ISFx, ISFy, ISFpara, ISFperp, ISFx_imag, ISFy_imag, ISFpara_imag, ISFperp_imag, ISFx_comoving, ISFy_comoving, ISFx_comoving_imag, ISFy_comoving_imag, ISFpara_comoving, ISFperp_comoving, ISFpara_comoving_imag, ISFperp_comoving_imag, cos_theta_v, sin_theta_v,md_rx ,md_ry ,rx2 ,ry2 ,rxry, rx2ry, rxry2,rx3 ,ry3, rx4, ry4, rx3ry, ry3rx, rx2ry2 = data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9], data[10], data[11], data[12], data[13], data[14], data[15], data[16], data[17], data[18], data[19], data[20], data[21], data[22], data[23],data[24], data[25], data[26], data[27], data[28],data[29], data[30], data[31], data[32], data[33], data[34], data[35], data[36], data[37], data[38] #if gamma > Omega: # axes[2,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[0,2].plot(time_s,ISFx_s,color=color, marker=marker, linestyle='none') axes[0,2].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[0,2].plot(t[b:e],np.real((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[0,2].plot(time,ISFpara_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[1,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') axes[1,2].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[1,2].plot(t[b:e],np.real((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])*comoving_factor),color=color) axes[1,2].plot(time,ISFperp_comoving, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[2,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_para(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_para(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[2,2].plot(time_s,ISFx_s,color=color, marker=marker, linestyle='none') axes[2,2].set_xscale('log') t=np.logspace(-3,3, 1000) isf_para_imag= np.load("ISF_analytics\ISF_para_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_para= np.load("ISF_analytics\ISF_para_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_para(t[b:e],Omega,gamma,Drot)) axes[2,2].plot(t[b:e],np.imag((isf_para[b:e]+(0+1j)*isf_para_imag[b:e])*comoving_factor),color=color) axes[2,2].plot(time,ISFpara_comoving_imag, color=color, marker=marker, linestyle='none') #if gamma > Omega: # axes[3,2].plot(t[b:e], np.real(np.exp((0-1j)*HA_mean_perp(t[b:e],Omega,gamma,Drot,k)-0.5*harmonic_approximation_2_perp(t[b:e],Omega,gamma,Drot,k))), linestyle='dotted', color='black') #axes[3,2].plot(time_s,ISFy_s,color=color, marker=marker, linestyle='none') axes[3,2].set_xlabel(r'$\omega t$') axes[3,2].set_xscale('log') isf_perp= np.load("ISF_analytics\ISF_perp_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") isf_perp_imag= np.load("ISF_analytics\ISF_perp_imag_Omega_"+str(Omega)+"_D_"+str(Drot)+"_g_"+str(gamma)+"_k_"+str(k)+"_num_"+str(101)+".npy") comoving_factor=np.exp((0+1j)*k*mean_displacement_perp(t[b:e],Omega,gamma,Drot)) axes[3,2].plot(t[b:e],np.imag((isf_perp[b:e]+(0+1j)*isf_perp_imag[b:e])),color=color) axes[3,2].plot(time,ISFperp_comoving_imag, color=color, marker=marker, linestyle='none') plot_reg_sim_isf_mit_comoving_imag(1,0.025) plt.tight_layout() #plt.savefig("../figures/Observables_gravitaxis_2/ISF_comoving_D_0.025.pdf",bbox_inches='tight') #plt.savefig("../figures/Observables_gravitaxis_2/ISF_D_0.25.pdf",bbox_inches='tight') plt.show() plot_isf_HA(1,0.025,1.5) plt.tight_layout() #plt.savefig("../figures/Observables_gravitaxis_2/ISF_HA_D_0.025_g_1_5.pdf",bbox_inches='tight') plt.show()