In [60]:
%matplotlib inline
%load_ext snakeviz


The snakeviz extension is already loaded. To reload it, use:
  %reload_ext snakeviz

In [2]:
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import ode
import matplotlib.pylab as plt
import csv
import time

In [42]:
endpoint = 10000; # integration range
dx = 100.0; # step size
dx100 = 100000.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.0, 0.0, 0.0], 0 # initial condition
savestep = 1000; # save to file every savestep steps

In [43]:
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 hamiltonian4(x, deltalambda, k, thetam):
    hr = np.array(hamiltonian(x, deltalambda, k, thetam)).real;
    hi = np.array(hamiltonian(x, deltalambda, k, thetam)).imag;
    
    # Ordering of 
    
    return [[hi[0][0],hi[0][1],hr[0][0],hr[0][1]], [hi[1][0],hi[1][1],hr[1][0],hr[1][1]], [- hr[0][0], - hr[0][1], hi[0][0], hi[0][1]],  [- hr[1][0], - hr[1][1], hi[1][0], hi[1][1]] ]

def sysdpsidt(psi, x, deltalambda, k, thetam):
    
    return np.dot(hamiltonian4(x, deltalambda, k, thetam), [psi[0], psi[1], psi[2], psi[3]])

def sysjac(psi, x, deltalambda, k, thetam):
    
    return hamiltonian4(x, deltalambda, k, thetam)

In [58]:
def integral_tol(total_error_needed,totalrange, stepsize): # tolenrance of the integral that we require
    
    return total_error_needed*stepsize/totalrange

In [59]:
# Test the function
#hamiltonian4(10,dellam,ks,thm)
#hamiltonian(10,dellam,ks,thm)
integral_tol(1e-4,endpoint,dx)


Out[59]:
1e-06

In [45]:
xlin = np.linspace(0, endpoint, np.floor(endpoint/dx) )
# xlin100 = np.linspace(0, endpoint, np.floor(endpoint/dx100) )
#print xlin

In [46]:
solodeint = odeint(sysdpsidt, psi0, xlin, args = (dellam,ks,thm), full_output = 1)

In [47]:
# solodeint100 = odeint(sysdpsidt, psi0, xlin100, args = (dellam,ks,thm), full_output = 1)

In [48]:
solodeint[1]['message']


Out[48]:
'Integration successful.'

In [49]:
prob0=solodeint[0][:,0]**2+solodeint[0][:,2]**2
prob1=solodeint[0][:,1]**2+solodeint[0][:,3]**2

#prob0_100=solodeint100[0][:,0]**2+solodeint100[0][:,2]**2
#prob1_100=solodeint100[0][:,1]**2+solodeint100[0][:,3]**2

In [50]:
#print prob0, prob1, prob0+prob1

In [51]:
# plt.figure(figsize=(18,13))

# plt.plot(xlin, prob0,'-')
# 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(xlin, prob1,'-')
# 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(xlin, prob0+prob1, '-')
# plt.plot(xlin, prob0+prob1, '-', xlin100, prob0_100 + prob1_100,'*')
# plt.plot(xlin100, prob0_100 + prob1_100,'*')
plt.title("Probabilities",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.ylim([0.999999,1.00001])
plt.show()



In [ ]:


In [3]:
import scipy
scipy.__version__


Out[3]:
'0.17.1'

Data Processing Unit


In [ ]:
xlin_supernova = np.load("assets/two-freq-real-ize-xlin-10000000-100.0.npy")
prob0_supernova = np.load("assets/two-freq-real-ize-prob0-10000000-100.0.npy")
prob1_supernova = np.load("assets/two-freq-real-ize-prob1-10000000-100.0.npy")

plt.figure(figsize=(18,13))

plt.plot(xlin_supernova, prob0_supernova,'-')
plt.title("Probability $P_{1\\to1}$ (step=100.0)",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.show()


plt.figure(figsize=(18,13))

plt.plot(xlin_supernova, prob1_supernova,'-')
plt.title("Probability $P_{1\\to2}$ (step=100.0)",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.show()

plt.figure(figsize=(18,13))

plt.plot(xlin_supernova, prob1_supernova+prob0_supernova,'-')
plt.title("Total Probability (step=100.0)",fontsize=20)
plt.xlabel("$\hat x$",fontsize=20)
plt.ylabel("Probability",fontsize=20)
plt.show()

Examples of odeint


In [33]:
def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

b = 0.25
c = 5.0

y0 = [np.pi - 0.1, 0.0]

t = np.linspace(0, 10, 101)

sol = odeint(pend, y0, t, args=(b, c),full_output=1)


plt.plot(t, sol[0][:, 0], 'b', label='theta(t)')
plt.plot(t, sol[0][:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()



In [36]:
sol[1]


Out[36]:
{'hu': array([ 0.02909913,  0.02909913,  0.02909913,  0.02461187,  0.04466197,
         0.04466197,  0.04975119,  0.04975119,  0.0394766 ,  0.0394766 ,
         0.0324055 ,  0.0324055 ,  0.0324055 ,  0.03697811,  0.03697811,
         0.03064065,  0.03064065,  0.03064065,  0.02819459,  0.02819459,
         0.02819459,  0.02819459,  0.02819459,  0.02819459,  0.03341736,
         0.03341736,  0.03341736,  0.03341736,  0.04176933,  0.03479012,
         0.03479012,  0.03479012,  0.03487164,  0.03487164,  0.03487164,
         0.03487164,  0.03487164,  0.02855564,  0.02855564,  0.02855564,
         0.03357137,  0.03357137,  0.03357137,  0.03754137,  0.03754137,
         0.03754137,  0.03754137,  0.05093141,  0.05093141,  0.05093141,
         0.05093141,  0.04016459,  0.04016459,  0.04016459,  0.04016459,
         0.04016459,  0.04016459,  0.04016459,  0.04016459,  0.04016459,
         0.04016459,  0.04974618,  0.0414047 ,  0.0414047 ,  0.0414047 ,
         0.05626441,  0.05626441,  0.05626441,  0.05626441,  0.04607684,
         0.04607684,  0.04607684,  0.04607684,  0.04607684,  0.05334296,
         0.05334296,  0.05334296,  0.05334296,  0.05334296,  0.05334296,
         0.05334296,  0.0451986 ,  0.0451986 ,  0.0451986 ,  0.0451986 ,
         0.0451986 ,  0.05230448,  0.05230448,  0.05230448,  0.05230448,
         0.05230448,  0.04399709,  0.04399709,  0.04399709,  0.04399709,
         0.04399709,  0.05666398,  0.05666398,  0.05666398,  0.05666398]),
 'imxer': -1,
 'leniw': 22,
 'lenrw': 52,
 'message': 'Integration successful.',
 'mused': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1], dtype=int32),
 'nfe': array([ 39,  45,  51,  61,  67,  71,  77,  81,  87,  93,  99, 105, 113,
        117, 123, 131, 137, 143, 153, 159, 167, 173, 181, 187, 193, 199,
        205, 211, 217, 223, 229, 235, 243, 247, 253, 259, 265, 273, 281,
        287, 295, 301, 307, 313, 317, 323, 329, 333, 337, 341, 345, 351,
        357, 361, 367, 371, 377, 381, 387, 391, 395, 399, 407, 411, 417,
        421, 423, 427, 431, 437, 441, 445, 451, 455, 459, 463, 467, 471,
        473, 477, 481, 485, 491, 495, 499, 503, 507, 511, 515, 519, 523,
        527, 531, 537, 541, 545, 549, 553, 557, 561], dtype=int32),
 'nje': array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0], dtype=int32),
 'nqu': array([5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
        7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8], dtype=int32),
 'nst': array([ 17,  20,  23,  28,  30,  32,  35,  37,  39,  42,  44,  47,  51,
         53,  56,  59,  62,  65,  69,  72,  76,  79,  83,  86,  89,  92,
         95,  98, 101, 103, 106, 109, 112, 114, 117, 120, 123, 126, 130,
        133, 137, 140, 143, 146, 148, 151, 154, 156, 158, 160, 162, 164,
        167, 169, 172, 174, 177, 179, 182, 184, 186, 188, 191, 193, 196,
        198, 199, 201, 203, 205, 207, 209, 212, 214, 216, 218, 220, 222,
        223, 225, 227, 229, 232, 234, 236, 238, 240, 242, 244, 246, 248,
        250, 252, 255, 257, 259, 261, 263, 265, 267], dtype=int32),
 'tcur': array([  0.12630625,   0.21360365,   0.30090104,   0.42396041,
          0.51328435,   0.60260828,   0.74168341,   0.84118578,
          0.92013898,   1.03856878,   1.10337979,   1.20059629,
          1.3302183 ,   1.40417452,   1.51510886,   1.61970573,
          1.71162767,   1.80354962,   1.92783714,   2.01242092,
          2.12519928,   2.20978306,   2.32256142,   2.4071452 ,
          2.50739729,   2.60764939,   2.70790148,   2.80815357,
          2.9251096 ,   3.00166905,   3.1060394 ,   3.21040976,
          3.33320553,   3.40294881,   3.50756374,   3.61217867,
          3.7167936 ,   3.80246052,   3.91668309,   4.00235002,
          4.12158831,   4.22230243,   4.32301655,   4.43167067,
          4.50675341,   4.61937753,   4.73200164,   4.82047442,
          4.92233724,   5.02420006,   5.12606288,   5.21715888,
          5.33765264,   5.41798181,   5.53847557,   5.61880475,
          5.73929851,   5.81962768,   5.94012144,   6.02045062,
          6.10077979,   6.20027215,   6.32448626,   6.40729567,
          6.53150978,   6.6440386 ,   6.70030301,   6.81283183,
          6.92536065,   7.01751432,   7.10966799,   7.20182167,
          7.34005218,   7.43220585,   7.53162565,   7.63831156,
          7.74499748,   7.8516834 ,   7.90502636,   8.01171228,
          8.1183982 ,   8.20879539,   8.34439118,   8.43478837,
          8.52518556,   8.61558275,   8.71308583,   8.81769479,
          8.92230375,   9.02691271,   9.13152167,   9.21951586,
          9.30751005,   9.43950133,   9.52749552,   9.61548971,
          9.71615078,   9.82947874,   9.94280669,  10.05613464]),
 'tolsf': array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]),
 'tsw': array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
         0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])}

In [ ]:


In [ ]: