Merge pull request #97 from bmad-sim/species_unit_fixes

Species unit fixes

src/AtomicAndPhysicalConstants.jl

@@ -8,13 +8,14 @@ using Reexport
 @reexport using Unitful
 
 include("units_definition.jl")
-include("physical_constants.jl")
-include("atomic_isotopes.jl")
+include("constants.jl")
+include("types.jl")
+include("constructors.jl")
+include("isotopes.jl")
 include("subatomic_species.jl")
-include("species_initialize.jl")
 include("update_pion_mass.jl")
 include("update_constants.jl")
-include("update_iso_masses.jl")
+include("update_isos.jl")
 include("particle_functions.jl")
 include("set_units.jl")
 
@@ -25,7 +26,7 @@ export C_LIGHT, H_PLANCK, H_BAR_PLANCK, R_E, R_P, E_CHARGE, MU_0_VAC, EPS_0_VAC,
 export SubatomicSpecies
 export AtomicSpecies
 export SUBATOMIC_SPECIES
-export Atomic_Particles
+export ATOMIC_SPECIES
 export setCODATA
 export @u_str, NewUnits
 

src/constants.jl

@@ -0,0 +1,134 @@
+
+# Constants pulled from the NIST table of
+# the 2022 CODATA release
+
+
+
+#####################################################################
+# constants with dimension [mass]
+#####################################################################
+
+const __b_m_electron = .51099895069 * u"MeV/c^2"
+# Electron Mass [MeV]/c^2
+const __b_m_proton = 9.382720894300001e2 * u"eV/c^2"
+# Proton Mass [MeV]/c^2
+const __b_m_neutron = 9.395654219399999e2 * u"MeV/c^2"
+# Neutron Mass [MeV]/c^2
+const __b_m_muon = 1.056583755e2 * u"MeV/c^2"
+# Muon Mass [MeV]/c^2
+const __b_m_helion = 2.80839161112e3 * u"MeV/c^2"
+# Helion Mass He3 nucleus [MeV]/c^2
+const __b_m_deuteron = 1.875612945e3 * u"MeV/c^2"
+# Deuteron Mass [MeV]/c^2
+
+# constants mysteriously missing from the release
+# picked up from PDG
+const __b_m_pion_0 = 1.349768277676847e2 * u"MeV/c^2"
+# uncharged pion mass [eV]/c^2
+const __b_m_pion_charged = 1.3957039098368132e2 * u"MeV/c^2"
+# charged pion mass [eV]/c^2
+
+
+
+
+#####################################################################
+# constants with dimension [magnetic moment]
+#####################################################################
+
+const __b_mu_deuteron = 4.330735087e-27 * u"J/T"
+# deuteron magnetic moment in [J/T]
+const __b_mu_electron = -9.2847646917e-24 * u"J/T"
+# electron magnetic moment in [J/T]
+const __b_mu_helion = -1.07461755198e-26 * u"J/T"
+# helion magnetic moment in [J/T]
+const __b_mu_muon = -4.4904483e-26 * u"J/T"
+# muon magnetic moment in [J/T]
+const __b_mu_neutron = -9.6623653e-27 * u"J/T"
+# neutron magnetic moment in [J/T]
+const __b_mu_proton = 1.41060679545e-26 * u"J/T"
+# proton magnetic moment in [J/T]
+const __b_mu_triton = 1.5046095178e-26 * u"J/T"
+# triton magnetic moment in [J/T]
+
+
+
+#####################################################################
+# dimensionless constants
+#####################################################################
+
+const __b_N_avogadro = 6.02214076e23;
+# Avogadro's constant: Number / mole (exact)
+const __b_fine_structure = 0.0072973525643;
+# fine structure constant
+
+
+
+#####################################################################
+# unit conversion constants
+#####################################################################
+
+const __b_kg_per_amu = 1.66053906892e-27 * u"kg/amu";
+# kg per standard atomic mass unit (dalton)
+const __b_eV_per_amu = 9.3149410372e8 * u"(eV/c^2)/amu";
+# eV per standard atomic mass unit (dalton)
+const __b_J_per_eV = 1.602176634e-19 * u"J/eV";
+# Joules per eV
+
+
+
+#####################################################################
+# constants with miscelaneous dimension
+#####################################################################
+
+const __b_e_charge = 1.602176634e-19 * u"C";
+# elementary charge [C]
+const __b_r_e = 2.8179403205e-15 * u"m";
+# classical electron radius [m]
+const __b_r_p = __b_r_e * __b_m_electron / __b_m_proton *u"m";
+# classical proton radius [m]
+const __b_c_light = 2.99792458e8 * u"m/s";
+# speed of light [m/s]
+const __b_h_planck = 4.135667696e-15 * u"eV*s";
+# Planck's constant [eV*s]
+const __b_h_bar_planck = __b_h_planck / 2pi * u"eV*s";
+# h_planck/twopi [eV*s]
+const __b_classical_radius_factor = __b_r_e * __b_m_electron * u"MeV*m/c^2";
+# e^2 / (4 pi eps_0) = classical_radius * mass * c^2.
+# Is same for all particles of charge +/- 1.
+
+const __b_eps_0_vac = 8.8541878188e-12 * u"F/m";
+# Permittivity of free space in [F/m]
+const __b_mu_0_vac = 1.25663706127e-6 * u"N/A^2";
+# Vacuum permeability in [N/A^2] (newtons per ampere squared)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#---------------------------------------------------------------------------------------------------

src/species_initialize.jl → src/constructors.jl

@@ -1,19 +1,13 @@
+# Package: AtomicAndPhysicalConstants
+# file: src/species_initialize.jl
+# purpose: define constructors
 
 
-struct Species
-    name::String # name of the particle to track
-    charge::Int32 # charge of the particle (important to consider ionized atoms) in [e]
-    mass_in_eV::Float64 # mass of the particle in [eV/c^2]
-    spin::Float64 # spin of the particle in [ħ]
-    mu::Float64 # magnetic moment of the particle (for now it's 0 unless we have a recorded value)
-    iso::Int # if the particle is an atomic isotope, this is the mass number, otherwise 0
-end;
-export Species
-
 
 
+#####################################################################
+#####################################################################
 
-# ------------------------------------------------------------------------------------------------------------
 
 """
 		subatomic_particle(name::String)
@@ -26,13 +20,15 @@ subatomic_particle
 function subatomic_particle(name::String)
     # write the particle out directly
     return Species(name, SUBATOMIC_SPECIES[name].charge,
-        SUBATOMIC_SPECIES[name].mass_in_eV,
-        SUBATOMIC_SPECIES[name].spin,
+        SUBATOMIC_SPECIES[name].mass,
+        SUBATOMIC_SPECIES[name].planck_spin,
         SUBATOMIC_SPECIES[name].mu,
         0)
 end
 
-# -----------------------------------------------------------------------------------------------
+
+#####################################################################
+#####################################################################
 
 
 """
@@ -44,22 +40,22 @@ of information specifice to the chosen particle.
 # Fields:
 1. `name::String': 				the name of the particle 
 
-2. `charge::Int32': 			the net charge of the particle in units of |e|
-													- bookkeeping only, thus in internal units
-													- use the 'charge()' function to get the charge 
-													- in the desired units
+2. `int_charge::typeof(1u"q")': 				 the net charge of the particle in units of |e|
+																		 	 - bookkeeping only, thus in internal units
+																			 - use the 'charge()' function to get the charge 
+																			 - in the desired units
+
+3. `mass::typeof(1.0u"eV/c^2")': 				 the mass of the particle in eV/c^2
+																			 - bookkeeping only, thus in internal units
+																		 	 - use the 'mass()' function to get the mass 
+																			 - in the desired units
 
-3. `mass_in_eV::Float64': the mass of the particle in eV/c^2
-													- bookkeeping only, thus in internal units
-													- use the 'mass()' function to get the mass 
-													- in the desired units
+4. `planck_spin::typeof(1.0u"ħ")': 					 the spin of the particle in ħ
 
-4. `spin::Float64': 			the spin of the particle in eV⋅s 
-													- (half/integer multiplied with ħ)
+5. `moment::typeof(1.0u"eV/T")': 					 the magnetic moment of the particle in eV/T
 
-5. `mu::Float64': 				the magnetic moment of the particle in eV/T.
-5. `iso::Int': 						if the particle is an atomic isotope, this is the 
-													- mass number, otherwise -1
+6. `iso::Int': 												 if the particle is an atomic isotope, this is the 
+																			 - mass number, otherwise -1
 
 The Species Struct also has a constructor called Species, 
 documentation for which follows.
@@ -147,35 +143,38 @@ function Species(name::String, charge::Int=0, iso::Int=-1)
                 charge = tryparse(Int, chstr)
             end
         end
-
+				if count('+', name) != 0 && count('-', name) != 0
+					error(f"""You made a typo in "{name}". You have both a + and a - in the name. """)
+					return
+				end
         if haskey(ATOMIC_SPECIES, AS) # is the particle in the Atomic_Particles dictionary?
             if iso ∉ keys(ATOMIC_SPECIES[AS].mass_in_amu) # error handling if the isotope isn't available
                 error("The isotope you specified is not available: Isotopes are specified by the atomic symbol and integer mass number.")
                 return
             end
-            mass_in_eV = begin
+            mass = begin
                 if anti_atom == false
-                    nmass = uconvert(u"eV/c^2", ATOMIC_SPECIES[AS].mass_in_amu[iso]u"amu"); # mass of the positively charged isotope in eV/c^2
+                    nmass = uconvert(u"MeV/c^2", ATOMIC_SPECIES[AS].mass_in_amu[iso]u"amu"); # mass of the positively charged isotope in eV/c^2
                     nmass.val + __b_m_electron.val * (ATOMIC_SPECIES[AS].Z - charge) # put it in eV/c^2 and remove the electrons
                 elseif anti_atom == true
-                    nmass = uconvert(u"eV/c^2", ATOMIC_SPECIES[AS].mass_in_amu[iso]u"amu"); # mass of the positively charged isotope in amu
+                    nmass = uconvert(u"MeV/c^2", ATOMIC_SPECIES[AS].mass_in_amu[iso]u"amu"); # mass of the positively charged isotope in amu
                     nmass.val + __b_m_electron.val * (-ATOMIC_SPECIES[AS].Z + charge) # put it in eV/c^2 and remove the positrons
                 end
             end
             if iso == -1 # if it's the average, make an educated guess at the spin
                 partonum = round(ATOMIC_SPECIES[AS].mass_in_amu[iso])
                 if anti_atom == false
-                    spin = 0.5 * (partonum + (ATOMIC_SPECIES[AS].Z - charge))
+                    planck_spin = 0.5 * (partonum + (ATOMIC_SPECIES[AS].Z - charge))
                 elseif anti_atom == true
-                    spin = 0.5 * (partonum + (ATOMIC_SPECIES[AS].Z + charge))
+                    planck_spin = 0.5 * (partonum + (ATOMIC_SPECIES[AS].Z + charge))
                 end
             else # otherwise, use the sum of proton and neutron spins
-                spin = 0.5 * iso
+                planck_spin = 0.5 * iso
             end
             if anti_atom == false
-                return Species(AS, charge, mass_in_eV, spin, 0, iso) # return the object to track
+                return Species(AS, charge*u"q", mass*u"MeV/c^2", planck_spin*u"ħ", 0*u"J/T", iso) # return the object to track
             elseif anti_atom == true
-                return Species("anti-" * AS, charge, mass_in_eV, spin, 0, iso)
+                return Species("anti-" * AS, charge*u"q", mass*u"MeV/c^2", planck_spin*u"ħ", 0u"J/T", iso)
             end