In [6]:

    

import numpy as np
import matplotlib.pyplot as plt
import h5py
import time
import progressbar
import sys
from Qutils import *
import progressbar

def calc_dipol(psi, dx, n0, n1, nt, x,nstep):
    """
    Calculates the Dipol Moment for a given
    wave-function psi for nt timesteps. Since
    it is just the mean of the location operator x, 
    also a numpy array containing the corresponding 
    discrete spatial values is required. 
    """
    roh = get_prob_dens(psi, nt, nstep)
    print(roh.shape)
    print(x.shape)
    res  = np.zeros(psi.shape[0])
    print("Calculating Dipole moment")
    with progressbar.ProgressBar(max_value=int(psi.shape[0])) as bar:  
        for i in range(0, psi.shape[0]):
                res[i] = np.trapz(roh[i,:]*x, dx=dx)
                bar.update(i)
    return res

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 [15]:

    

filepath = "../../build/res.h5"
fig = plt.figure(figsize=(14,10))
ax1 = plt.subplot()
nx = np.int32(1e5)
nt = np.int32(1e5)
xmax = 30.0
xmin = -xmax
tmax = 100.0
tmin = 0

nstep = 100
t = np.linspace(tmin, tmax, int(nt/nstep))

h = 0.0006 
n0 = 50000
n1 = 66667
dx = 0.0006
psi = load_vals(filepath, nt, nx, nstep)
x = np.linspace(xmin, xmax, psi.shape[1])
p = calc_dipol(psi, dx, n0, n1, nt, x, nstep)

p *= 1/np.max(p)


ax1.plot(t, p, color="r",lw=2,label=r"$\lambda_1=29.98 \, nm , \; I = 0.816 \, keV$")

vals = calc_dist(x,t)
res = int_dist(vals, h)
res *= 1/np.max(res)
ax1.set_xlabel("t $(at.u.)$",size=20)
ax1.set_ylabel("Normed quantities $(arb. u.)$",size=20)

ax1.plot(t,res,"g",label=r"$\int \, V(x) \, dx$ normed")
plt.legend(loc='best',prop={'size':15})


plt.title(r"$\vec{P}=-e<x>\vec{e_x}$",size=20)
plt.show()


    

    




  2% (  2600 of 100000) |#                                                       | Elapsed Time: 0:00:00 ETA: 0:00:06





    




Loading file






    




100% (100000 of 100000) |#######################################################| Elapsed Time: 0:00:08 Time: 0:00:08
 18% ( 189 of 1000) |###########                                                 | Elapsed Time: 0:00:00 ETA: 0:00:00





    




(1000, 100000)
(100000,)
Calculating Dipole moment






    




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

In [ ]:

    




    

In [ ]:

    




    

In [ ]: