import numpy as np
import matplotlib.pyplot as plt

# Constants and parameters
mu0 = 4 * np.pi * 1e-7      # permeability of free space [H/m]
amp_turn = 100              # amp-turn rating of calibration coil
R = 0.1                    # radius of calibration coil [m]
Np = 50                     # number of turns of the B-dot probe
Dp = 0.005                  # diameter of probe [m]
A = np.pi * (Dp/2)**2       # cross-sectional area of the probe [m^2]

# Define magnetic field along the z-axis for a circular loop
def B(z):
    return (mu0 * amp_turn * R**2) / (2 * (R**2 + z**2)**(3/2))

# Compute dB/dz analytically
def dBdz(z):
    # d/dz [ (R^2+z^2)^(-3/2) ] = -3z/(R^2+z^2)^(5/2)
    return - (mu0 * amp_turn * R**2) * (3 * z) / (2 * (R**2 + z**2)**(5/2))

# Time array [s]
t = np.linspace(0, 0.25, 1000)  # from 0 to 250 ms

# Define probe motion: initial position z0 = 0.2 m, moves to 0 m then back to 0.2 m.
z = np.zeros_like(t)
dz_dt = np.zeros_like(t)
for i, ti in enumerate(t):
    if ti < 0.1:  # 0 to 100 ms, moving in negative z at -2 m/s
        dz_dt[i] = -2
        z[i] = 0.2 - 2 * ti
    elif ti < 0.15:  # 100 to 150 ms, stationary at z = 0
        dz_dt[i] = 0
        z[i] = 0.0
    else:  # 150 to 250 ms, moving in positive z at +2 m/s
        dz_dt[i] = 2
        z[i] = 2 * (ti - 0.15)

# Calculate induced voltage using V(t) = - Np * A * (dB/dz) * (dz/dt)
V = -Np * A * dBdz(z) * dz_dt

# Plot the voltage trace
plt.figure(figsize=(8,4))
plt.plot(t*1e3, V*1e6)  # converting time to ms and voltage to microvolts
plt.xlabel('time, ms')
plt.ylabel('induced voltage, µV')
plt.title('Voltage trace')
plt.grid(True)
plt.show()


import numpy as np
import matplotlib.pyplot as plt

m_e = 9.1093837E-31 # electron mass
M_Xe = 2.1801714E-25 # Xenon mass
k = 1.38E-23 # Boltzmann constant
T_tr = 11600


T_e = [5, 50]   # electron temp
T_i = [0.5, 5]  # ion temp

def VDF(m,V,T):
    return (m/(2*np.pi*T*k*T_tr))**0.5 * np.exp(-m*V**2/(2*T*T_tr*k))

def SDF(m,C,T):
    return 4*np.pi*C**2*(m/(2*np.pi*T*k*T_tr))**1.5 * np.exp(-m*C**2/(2*T*T_tr*k))

V = np.linspace(-10000, 10000, 200000)

plt.figure(figsize=(8, 5))

for T, label in zip(T_i, ['5 eV', '50 eV']):
    y = VDF(M_Xe, V, T)
    plt.plot(V, y, label=f'Temperature = {label}')

# Customize your plot
plt.title(r'Speed distribution function for Xe')
plt.xlabel(r'C, m/s')
plt.ylabel(r'normalized distribution function')
plt.grid(True)
plt.legend()

# Display the plot
plt.show()

import numpy as np
import matplotlib.pyplot as plt



def func(x):
    return (5.8e-17*x**4 - 3.71e-12*x**3 - 8e-8*x**2 + 6.37e-4*x + 6.97)/(-1.08e-21*x**5 + 
                                                                          1.86e-16*x**4 - 6.49e-12*x**3 + 8.97e-8*x**2 - 5.42e-4*x + 2.02)


x = np.linspace(0, 15000, 100000)
y = func(x)

plt.figure(figsize=(8, 5))  # Optional: set figure size
plt.plot(x, y)

# Customize your plot
plt.title('ratio of the partition functions of $Xe^+$ and Xe  ')
plt.xlabel(r'$\theta$')
plt.ylabel('$f_+/f_a$')
plt.grid(True)
plt.legend()

# Display the plot
plt.show()

import numpy as np
import matplotlib.pyplot as plt

m = 9.1E-31
k = 1.38E-23
h = 6.62607015E-34
I = 12.13*1.6E-19


pressures_Torr = [1E-3, 1E-1, 1E1]  # Torr
pressures_Pa = [p * 133.322 for p in pressures_Torr]  # conversion from mTorr to Pa

def C(T,p):
    return 2*(2*np.pi*m)**1.5 *(k*T)**2.5 /(p*h**3)

def func(T):
    return (5.8e-17*T**4 - 3.71e-12*T**3 - 8e-8*T**2 + 6.37e-4*T + 6.97)/(-1.08e-21*T**5 + 
                                                                          1.86e-16*T**4 - 6.49e-12*T**3 + 8.97e-8*T**2 - 5.42e-4*T + 2.02)

def RHS(T,p):
    return (C(T,p)*func(T)*np.exp(-I/(k*T))/(1+C(T,p)*func(T)*np.exp(-I/(k*T))))**0.5

def A(T,p):
    return (-RHS(T, p) + (RHS(T, p)**2 + 4*RHS(T, p))**0.5)/2

T = np.linspace(0.1, 10000, 10000)

plt.figure(figsize=(8, 5))  # Optional: set figure size

for p, label in zip(pressures_Pa, ['1E-3 Torr', '1E-1 Torr', '10 Torr']):
    y = A(T,p)
    plt.plot(T, y, label=f'Pressure = {label}')

# Customize your plot
plt.title(r'ionization fraction $\alpha$ as a function of temperature for different pressures')
plt.xlabel(r'T,K')
plt.ylabel(r'$\alpha$')
plt.grid(True)
plt.legend()

# Display the plot
plt.show()