Coloumb explosion.py

import numpy as np
from random import uniform
import math as m
import matplotlib.pyplot as plt

'''constants'''
delta = 1.758820149*10**11 #q/m
pi = m.pi
e0 = 8.854184817*10**-12
sigma = 1/10
mu = 0.5
pc = 10**-5


def bin_search(mas, x):
    i = 0
    j = len(mas)-1
    m = j // 2
    while i < j:
        if x > mas[m]:
            i = m+1
        else:
            j = m-1
        m = int((i+j)/2)
    return i
    

def F(x):
    return m.sqrt(x*(x-1)) + m.acosh(m.sqrt(x))

def Q(r) :
    C1 = 4*pc*m.sqrt(pi/2)
    C2 = m.sqrt(pi/2)*(sigma**2 + mu**2)
    C3 = m.erf( (r-mu)/(m.sqrt(2)*sigma) ) + m.erf( mu/(m.sqrt(2)*sigma) )
    C4 = sigma**2
    C5 = mu*m.exp(-mu**2/(2*sigma**2))
    C6 = (r+mu)*m.exp(-(r-mu)**2/(2*sigma**2))
    return C1*(C2*C3 + C4*(C5 - C6))

def gamma(r) :
    C1 = delta/(4*pi*e0)
    return C1*Q(r)

def lamda(r) :
    return m.sqrt(2*gamma(r))/r**(3/2)

def density_start(r):
    C1 = pc/(m.sqrt(2*pi)*sigma)
    return C1*m.exp(-(r-mu)**2/(2*sigma**2))


def density(R0, R, t):
    p0 = density_start(R0)
    C1 = R0**2 * p0 / R**2
    lamda0 = lamda(R0)
    C2 = R/R0 + t*m.sqrt(1 - R0/R)*(delta*p0/(e0*lamda0) - 3*lamda0/2)
    return C1 / C2


def characteristic(R0, R, lamda0) :
    return F(R/R0) / lamda0


def run_R(n, r_max) :
    dx = r_max/n
    global cord
    cord = [x for x in np.arange(dx, r_max+dx, dx)]
    data = []
    for i in range(1, len(cord)) :
        dR = 0.0001
        data.append([ [], [] ])
        for R in np.arange(cord[i], 2+dR, dR) :
            lamda0 = lamda(cord[i])
            data[i-1][0].append(characteristic(cord[i], R, lamda0))
            data[i-1][1].append(R)
    return data  # [ [time], [cord] ]


def run_p(data, t) :
    dens = []
    crd = []
    for i in range(1, len(cord)) :
        ind = bin_search(data[i-1][0], t)
        if ind > 0 :
            R_sr = (data[i-1][1][ind] + data[i-1][1][ind-1])/2
        else :
            R_sr = data[i-1][1][ind]
        ddd = density(cord[i], R_sr, t)
        if ddd < 0 :
            ddd = -ddd
        dens.append(ddd)
        crd.append(R_sr)
    cdr = crd.sort()
    return dens, crd


def main():
    data = run_R(520, 1)

    graph = plt.figure()
    ax = graph.add_subplot(111)
    for i in range(17, 52, 7) :
        ax.plot(data[i*10][1], data[i*10][0])
    plt.grid(True)

    graph = plt.figure()
    ay = graph.add_subplot(111)
    p, c = run_p(data, 0)
    ay.plot(c, p, "red")
    
    for i in np.arange(0, 3.6, 0.7) :
        p, c = run_p(data, i*10**-9)
        ay.plot(c, p)
    p, c = run_p(data, 3.77*10**-9)
    ay.plot(c, p)
    p, c = run_p(data, 3.96*10**-9)
    ay.plot(c, p)
    plt.grid(True)

    plt.show()
    

if __name__ == "__main__" :
    main()