from numpy import *

eps0 = 8.8541878128e-12
kB   = 1.380649e-23
qe   = 1.602176634e-19


T = 1.e5
n0 = 1.e7
r0 = 1.
r1 = 1./8
gamma = 5./3


def debye(T,n):
    return sqrt(eps0*kB*T/n/qe**2)

def plasma_par(n,l):
    return n*l**3


#Bestimmen Sie die Debyelaenge bei 1 AU
l_D = debye(T,n0)
print 'Die Debyelaenge be 1 AU betraegt {} km.'.format(l_D/1.e3)

#Bestimmen Sie den Plasmaparameter be 1 AU
pp = plasma_par(n0,l_D)
print 'Der Plasmaparameter be 1 AU betraegt {}'.format(pp)

#Dieselben Groessen bei 1/8 AU:

#spaerische und adiabatische Expansion:
#Dichte geht mit 1/r^2
#Temperatur:
"""
\begin{eqnarray*}
p &=& n k_B T, \\
p^{1 - \gamma} T^\gamma &=& {\rm const}.\\
[n k_B T|_{\rm 1 AU}]^{1 - \gamma} T_{\rm 1 AU}^\gamma = n_{\rm 1 AU}^{1 - \gamma} k_B T_{\rm 1 AU},\\
&=& n_{\rm 1/8 AU}^{1 - \gamma} k_B T_{\rm 1/8 AU},\\
T_{\rm 1/8 AU} = \left(\frac{n_{\rm a AU}}{n_{\rm 1/8 AU}}\right)^{1-\gamma} T_{\rm 1 AU}.\\
T{\rm 1/8 AU} &=& T_{\rm 1 AU} (r^2)^{1 - \gamma}
\end{eqnarray*}

"""
n1 = (r0/r1)**2
T1 = T*(r1/r0)**(2*(1-gamma))


l_D1 = debye(T1,n1)
pp1 = plasma_par(n1,l_D1)
print 'Die Debyelaenge bei 1/8 AU betraegt {} km'.format(l_D1/1.e3)
print 'Der Plasmaparameter be 1/8 AU betraegt {}'.format(pp1)