In [9]:

    

#Pkg.clone("https://github.com/jhlq/Equations.jl")
#Pkg.update()
using Equations

#=
    λD is the Debye length, 
    ε0 is the permittivity of free space, 
    kB is Boltzmann's constant,
    qe is the charge on an electron, 
    Te and Ti are the temperatures of the electrons and ions, respectively, multiplied by the Boltzmann's constant,
    ne is the density of electrons
=#

constants=[:me≖9.1094e-31,:qe≖1.6022e-19,:kB≖1.3807e-23,:ε0≖8.8542e-12]

# Solar wind
ne=:ne≖10.0^7
Te=:Te≖10.0*:qe
variables=[ne,Te]
# Debye length lD
lD=:λD≖sqrt(:ε0*:Te/(:ne*:qe^2)) #valid when Te>>Ti
# plasma frequency wp
wp=:ωp≖sqrt(:ne*:qe^2/(:ε0*:me))
# particles in a debye cube ND
ND=:ND≖:ne*:λD^3
println(lD&variables&constants)
println(wp&lD&variables&constants)
println(ND&lD&variables&constants)


    

    




λD ≖ 7.433892903446526
ωp ≖ 178400.95597166813
ND ≖ 4.108174668936227e9

In [10]:

    

println("Lightning")
ne.rhs=10.0^23
println(lD&variables&constants)
println(wp&lD&variables&constants)
println(ND&lD&variables&constants)
println("Magnetic fusion plasma")
ne.rhs=10.0^20
Te.rhs=10.0^4*:qe
println(lD&variables&constants)
println(wp&lD&variables&constants)
println(ND&lD&variables&constants)
println("Magnetic fusion plasma")
ne.rhs=10.0^17
Te.rhs=10.0^4*:kB
println(lD&variables&constants)
println(wp&lD&variables&constants)
println(ND&lD&variables&constants)


    

    




Lightning
λD ≖ 7.433892903446526e-8
ωp ≖ 1.7840095597166812e13
ND ≖ 41.08174668936228
Magnetic fusion plasma
λD ≖ 7.433892903446526e-5
ωp ≖ 5.641533576218887e11
ND ≖ 4.108174668936228e7
Magnetic fusion plasma
λD ≖ 2.1822655603609336e-5
ωp ≖ 1.784009559716681e10
ND ≖ 1039.2566127081407

In [11]:

    

#parameters for solar core: ne=10^32, Te=10^7 => λD ~ 10^−11
ne.rhs=10.0^32
Te.rhs=10.0^7*:kB
println(lD&variables&constants)
println(wp&lD&variables&constants)
println(ND&lD&variables&constants)


    

    




λD ≖ 2.1822655603609335e-11
ωp ≖ 5.641533576218887e17
ND ≖ 1.0392566127081406

In [12]:

    

# Drifts
ve=:ve≖Cross(:E⊥,:B)/norm(:B)^2
ex=[1,0,0];ey=[0,1,0];ez=[0,0,1]
B=:B≖:Bt*Vec(ez)
Bt=:Bt≖1.5
r=:r≖1
dBt=:dBt≖0.3
Eperp=:E⊥≖-:r/2*:dBt*Vec(ey)
vars=[B,Bt,Eperp,r,dBt]
println(ve&vars)
println(ve&Eperp&B)


    

    




ve ≖ (-0.09999999999999999) Vec([1,0,0])
ve ≖ (-0.5) dBt r ╱(Bt) Vec([1,0,0])

In [13]:

    

equ=Equation

# Compression factor
op1=equ(:rt,:r0*sqrt(:Bz0/:Bzt))
op3=equ(:rt,:r0*:Bzt/:Bz0)

Bpz=equ(:Ḃz,0.1)
Bz0=equ(:Bz0,5)
r0=equ(:r0,3)
rp=equ(:ṙ,-0.5*:r/:Bz*:Ḃz)
rp0=equ(:ṙ0,-0.5*:r0/:Bz0*:Ḃz)
vars=[rp0,Bpz,Bz0,r0]
rd=equ(:rΔ,:r0+:ṙ0)
Bzd=equ(:BzΔ,:Bz0+:Ḃz)
Bzds=Bzd&vars
rds=rd&rp0&vars
rdn=(op1&equ(:Bzt,:BzΔ)&vars&Bzds).rhs #n=numeric

nsteps=9
rq=rd
rqn=rds.rhs
op1s=rdn
Bqn=Bzds.rhs
for step in 1:nsteps #q=r-1
	rqp=-0.5/Bqn*rq.rhs*:Ḃz
	rr=equ(symbol("r$(step)Δ"),rqn+rqp)
	rrs=rr&vars
	print(rrs,"   \t")
	Br=equ(symbol("Bz$(step)Δ"),Bqn+:Ḃz)
	Brs=Br&vars
	println((op1&equ(:Bzt,symbol("Bz$(step)Δ"))&vars&Brs))
	Bqn=Brs.rhs
	rqn=rrs.rhs
end


    

    




r1Δ ≖ 2.940882352941176   	rt ≖ 2.9417420270727606
r2Δ ≖ 2.912324660633484   	rt ≖ 2.913857587071793
r3Δ ≖ 2.8843057927089557   	rt ≖ 2.886751345948129
r4Δ ≖ 2.856805792708956   	rt ≖ 2.8603877677367775
r5Δ ≖ 2.8298057927089557   	rt ≖ 2.8347335475692046
r6Δ ≖ 2.8032879355660985   	rt ≖ 2.8097574347450824
r7Δ ≖ 2.777235303987151   	rt ≖ 2.785430072655779
r8Δ ≖ 2.751631855711289   	rt ≖ 2.7617238536949706
r9Δ ≖ 2.726462364185865   	rt ≖ 2.7386127875258315

In [14]:

    

# Temperatures
Ekin=equ(:Ekin,0.5*:Tpa+:Tpe)
parfac=1
Tpa1=equ(:Tpa1,parfac*:T0) #is parallel temperature affected?
Tpe1=equ(:Tpe1,2*:T0)
println(Tpa1,", ",Tpe1)
Ek1=equ(:Ek1,0.5*:Tpa1+:Tpe1)
Ek2=equ(:Ek1,0.5*:T2+:T2)
Ek1s=Ek1&Tpa1&Tpe1 #s=solved
T2=(Ek2&Ek1s)/1.5
T2=equ(T2.rhs,T2.lhs)
println(T2)
Tpa3=equ(:Tpa3,1/parfac*:T2) #is it?
Tpe3=equ(:Tpe3,:T2/2)
Ek3=equ(:Ek3,0.5*:Tpa3+:Tpe3)
Ekf=[equ(:Ekf,0.5*:Tf+:Tf),equ(:Ekf,:Ek3)]
Tf=(Ekf[1]&Ekf[2]&Ek3&Tpa3&Tpe3&T2)/1.5
println(Tf) #nope


    

    




Tpa1 ≖ 1 T0, Tpe1 ≖ 2 T0
T2 ≖ 1.6666666666666665 T0
Tf ≖ 1.111111111111111 T0