In [1]:

    




    

In [2]:

    

def calc_dist(x, t):
    """
    Calculate the disturbance term 
    """
    a = 6.9314718055994524e-07
    b = 0.0069314718056
    t0 = 50.0
    w = 1.51939
    k = w/137
    I = 20.0
    res = np.zeros([t.size,x.size])
    for i in range(0, t.size):
        if t[i] < 50:
            g = t[i]/t0
        else:
            g = 1.0
        res[i] = I * np.sin(w*t[i]-k*x)*g
    return res

def int_dist(vals, h):
    """
    """
    res = np.zeros(vals.shape[0])
    for i in range(0, vals.shape[0]):
        res[i] = np.trapz(vals[i],dx=h)
    return res


    

In [3]:

    

import Qutils
filepath = "../../build/res.h5"
nx = np.int32(1e5)
nt = np.int32(1e5)
nstep = 100
h = 0.0006 
psi = load_vals(filepath, nt, nx, nstep)
res = integrate_prob_current(psi, 50000, 66667, 0.0006)
#res = integrate_prob(psi, 50000, 66667, 0.0006)
t = np.linspace(0,100, 1000)
x = np.linspace(0, 66667*0.0006-30.0, 5000)
res_2 = calc_dist(x,t)
res_2 = int_dist(res_2,h)
res_2 *= 1/np.max(res_2)
res *= 1/np.max(res)


    

    




  0% (   900 of 100000) |                                                        | Elapsed Time: 0:00:00 ETA: 0:00:13





    




Loading file






    




100% (100000 of 100000) |#######################################################| Elapsed Time: 0:00:09 Time: 0:00:09
  6% (  60 of 1000) |###                                                         | Elapsed Time: 0:00:00 ETA: 0:00:01





    




Calculating conjugate






    




100% (1000 of 1000) |###########################################################| Elapsed Time: 0:00:02 Time: 0:00:02






    




Calculating gradient...
Finished gradient!
Calculating probability current






    




100% (1000 of 1000) |###########################################################| Elapsed Time: 0:00:00 Time: 0:00:00






    




Finished calculating the integrated prob. current!

In [13]:

    

fig = plt.figure(figsize=(14,10))
plt.plot(t, res_2, color="red",label=r"$\int \, dx \, V_1(x,t)$")
plt.xlabel(r"$t \; (a.u)$", size=20)
plt.ylabel(r"Integrated quantities", size=20)

plt.plot(t, -res, color="blue", label=r"$-\int \, dx \, j(x)$")
plt.title("Integrated ")
plt.legend(loc='best')
plt.show()


    

In [ ]: