In [30]:
%matplotlib inline
%load_ext snakeviz
In [31]:
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import ode
import matplotlib.pylab as plt
import csv
import time
In [32]:
endpoint = 1000; # integration range
dx = 1.0; # step size
lam0 = 0.845258; # in unit of omegam, omegam = 3.66619*10^-17
dellam = np.array([0.00003588645221954444, 0.06486364865874367]); # deltalambda/omegam
ks = [1.0,1.0/90]; # two k's
thm = 0.16212913985547778; # theta_m
psi0, x0 = [1.0+0.j, 0.0], 0 # initial condition
savestep = 1; # save to file every savestep steps
In [36]:
xlin = np.arange(dx,endpoint+1*dx, dx)
psi = np.zeros([len(xlin) , 2], dtype='complex_')
xlinsave = np.zeros(len(xlin)/savestep);
psisave = np.zeros([len(xlinsave) , 2], dtype='complex_')
probsave = np.zeros([len(xlinsave) , 3])
def hamiltonian(x, deltalambda, k, thetam):
return [[ 0, 0.5* np.sin(2*thetam) * ( deltalambda[0] * np.sin(k[0]*x) + deltalambda[1] * np.sin(k[1]*x) ) * np.exp( 1.0j * ( - x - np.cos(2*thetam) * ( ( deltalambda[0]/k[0] * np.cos(k[0]*x) + deltalambda[1]/k[1] * np.cos(k[1]*x) ) ) ) ) ], [ 0.5* np.sin(2*thetam) * ( deltalambda[0] * np.sin(k[0]*x) + deltalambda[1] * np.sin(k[1]*x) ) * np.exp( -1.0j * ( - x - np.cos(2*thetam) * ( deltalambda[0] /k[0] * np.cos(k[0]*x) + deltalambda[1] /k[1] * np.cos(k[1]*x) ) ) ), 0 ]] # Hamiltonian for double frequency
def deripsi(t, psi, deltalambda, k , thetam):
return -1.0j * np.dot( hamiltonian(t, deltalambda,k,thetam), [psi[0], psi[1]] )
sol = ode(deripsi).set_integrator('zvode', method='bdf', atol=1e-8, with_jacobian=False)
sol.set_initial_value(psi0, x0).set_f_params(dellam,ks,thm)
flag = 0
flagsave = 0
In [37]:
timestampstr = time.strftime("%Y%m%d-%H%M%S")
print timestampstr
In [38]:
while sol.successful() and sol.t < endpoint:
sol.integrate(xlin[flag])
if np.mod(flag,savestep)==0:
probsave[flagsave] = [sol.t, np.absolute(sol.y[1])**2, np.absolute(sol.y[0])**2]
with open(r'probtrans-test-'+timestampstr+'.csv', 'a') as f_handle:
np.savetxt(f_handle, probsave[flagsave])
# with open(r'probtrans-test-'+timestampstr+'.csv', 'a') as f:
# writer = csv.writer(f)
# writer.writerow(1)
flagsave = flagsave + 1
flag = flag + 1
# print sol.t, sol.y
In [39]:
# # ploting using probsave array inside file
# plt.figure(figsize=(18,13))
# plt.plot(probsave[:,0], probsave[:,1],'-')
# plt.title("Probabilities",fontsize=20)
# plt.xlabel("$\hat x$",fontsize=20)
# plt.ylabel("Probability",fontsize=20)
# plt.show()
In [46]:
# Template for reading the csv file
# Ploting using data file
probsavefromfile = np.loadtxt("probtrans-test-"+timestampstr+".csv")
# print test
# print len(test[1::2]), test[1::2], len(test[::2]), test[::2]
plt.figure(figsize=(18,13))
plt.plot(probsavefromfile[::3], probsavefromfile[1::3],'-')
plt.title("Probabilities",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.show()
plt.figure(figsize=(18,13))
plt.plot(probsavefromfile[::3], probsavefromfile[2::3],'-')
plt.title("Probabilities",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.show()
plt.figure(figsize=(18,13))
plt.plot(probsavefromfile[::3], probsavefromfile[2::3] + probsavefromfile[1::3],'-')
plt.title("Probabilities",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.ylim([0.999999,1.000001])
plt.show()
In [ ]:
probsavefromfile[0]
In [ ]:
# http://stackoverflow.com/questions/30376581/save-numpy-array-in-append-mode
import tables
filename = 'outarray.h5'
ROW_SIZE = 100
NUM_COLUMNS = 2
f = tables.open_file(filename, mode='w')
atom = tables.Float64Atom()
array_c = f.create_earray(f.root, 'data', atom, (0, ROW_SIZE))
for idx in range(NUM_COLUMNS):
x = np.random.rand(1, ROW_SIZE)
array_c.append(x)
f.close()
In [ ]:
f = tables.open_file(filename, mode='a')
f.root.data.append(x)
In [ ]:
f = tables.open_file(filename, mode='r')
print(f.root.data[1:10,2:20]) # e.g. read from disk only this part of the dataset