In [30]:

%matplotlib inline







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




20160620-134525




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

# 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]



## HDF5

This might also be done using hdf5



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