1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 | import numpy as np
from scipy.integrate import odeint
from scipy.integrate import ode
import matplotlib.pylab as plt
import csv
##### Define parameters ########
endpoint = 10000000; # 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.j], 0 # initial condition
savestep = 100000; # save to file every savestep steps
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])
####### Define the functions required ########
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]] )
def jac(t, psi, deltalambda, k, thetam):
return -1.0j * np.array(hamiltonian(t, deltalambda, k, thetam))
############ integrate and save ##########
sol = ode(deripsi,jac).set_integrator('zvode', method='adams', atol=1e-8, with_jacobian=True)
sol.set_initial_value(psi0, x0).set_f_params(dellam,ks,thm).set_jac_params(dellam,ks,thm)
flag = 0
flagsave = 0
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] # save the probability
with open(r'assets/probtrans-test-1e7.csv', 'a') as f_handle:
np.savetxt(f_handle, probsave[flagsave])
flagsave = flagsave + 1
print np.absolute(sol.y[1])**2 + np.absolute(sol.y[0])**2 # Check if the solution conserves the total probability
print np.array(hamiltonian(sol.t,dellam,ks,thm)) # check if the Hamiltonian is Hermitian
flag = flag + 1
#### Notify at the end######
print "CONGRATS"
|
© 2016-2018, Lei Ma | Created with Sphinx and . | On GitHub | Physics Notebook Statistical Mechanics Notebook Neutrino Physics Notes Intelligence | Index | Page Source