Chemistry in GEAR
*****************

This page described quickly the implementation of the chemistry in ``GEAR``.
This gives the user an overview of the routines that are used to run ``GEAR``.

Setting the mass range
""""""""""""""""""""""

in chimie.c

For each stellar particle, we first find the min and max live time of a particle::

  /* minimum live time for a given metallicity */
  minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmax*SOLAR_MASS/All.UnitMass_in_g)*All.HubbleParam;

  /* maximum live time for a given metallicity */
  maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmin*SOLAR_MASS/All.UnitMass_in_g)*All.HubbleParam;

where the function :c:func:`star_lifetime` has been used.
Then, we compute ``t01`` and ``t02``, which bracket the time interval between the current age and
the time during which the stellar particule was created. 

Then, ``m1`` and ``m2`` are determined, corresponding to the mass of stars that end their life
during the current time step.

``m1`` is either::


  m1 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t02/All.HubbleParam)*All.HubbleParam;
  
or::

  m1 = Cp->Mmin*SOLAR_MASS/All.UnitMass_in_g*All.HubbleParam;    


Similarly, ``m2`` is either::


  m2 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t01/All.HubbleParam)*All.HubbleParam;
  
or::

  m2 = Cp->Mmax*SOLAR_MASS/All.UnitMass_in_g*All.HubbleParam; 

where we have used :c:func:`star_mass_from_age`.


Computing the number of dying stars
"""""""""""""""""""""""""""""""""""

Once we know that some stars may die during the
current timestep, we need to compute the number of dying stars

The number(float) of SNIa is given by (:c:func:`SNIa_rate`):: 

  NSNIa = SNIa_rate(m1/All.HubbleParam,m2/All.HubbleParam)*M0/All.HubbleParam;

The number(float) of SNII is given by (:c:func:`SNII_rate`):: 

  NSNII = SNII_rate(m1/All.HubbleParam,m2/All.HubbleParam)*M0/All.HubbleParam;

The number(float) of dyning stars other than supernova is given by (:c:func:`DYIN_rate`):: 

  NDYIN = DYIN_rate(m1/All.HubbleParam,m2/All.HubbleParam)*M0/All.HubbleParam;



Computing ejected mass and yields
"""""""""""""""""""""""""""""""""


The continuous case
-------------------

In the continuous case (default case,``CHIMIE_MC_SUPERNOVAE`` is not used), 
we assume that the number of dying stars of each type is given by 
``NSNII``,``NSNIa``,``NDYIN``, **even if these values are lower than 1**.
The ejected total mass, as well as the ejected mass of each element 
is computed using the function :c:func:`Total_mass_ejection`::

  Total_mass_ejection(m1/All.HubbleParam,m2/All.HubbleParam,M0/All.HubbleParam,metals);

This functions fill the array ``EjectedMass``. 
This array is then copied to the variables
containing the total gas mass and the mass ejected per element::

  StP[m].TotalEjectedGasMass = EjectedMass[0]*All.HubbleParam;        

  for (k=0;k<NELEMENTS;k++) 
    StP[m].TotalEjectedEltMass[k] = EjectedMass[k+2]*All.HubbleParam;             


The discrete case
-----------------


In the case where flag ``CHIMIE_MC_SUPERNOVAE`` is used, 
``NSNIa``, ``NSNII`` and  ``NDYIN`` are rounded 
to integer values, according to some probabilities (Monte Carlo technique) ::


  /* discretize SNIa */
  fNSNIa = NSNIa-floor(NSNIa);
  NSNIa  = floor(NSNIa);        

  if (get_Chimie_random_number(P[i].ID) < fNSNIa)
    NSNIa = NSNIa+1;
  
  /* discretize SNII */
  fNSNII = NSNII-floor(NSNII);
  NSNII  = floor(NSNII);
  
  if (get_Chimie_random_number(P[i].ID) < fNSNII)
    NSNII = NSNII+1;
  
  /* discretize DYIN */
  fNDYIN = NDYIN-floor(NDYIN);
  NDYIN  = floor(NDYIN);
  
  if (get_Chimie_random_number(P[i].ID) < fNDYIN)
    NDYIN = NDYIN+1;   


The total mass, as well as the ejected mass of each element 
is computed using the function :c:func:`Total_single_mass_ejection`::

  Total_single_mass_ejection(0.5*(m1+m2)/All.HubbleParam,metals,NSNII,NSNIa,NDYIN);


This functions fill the array ``SingleEjectedMass``.
This array is then copied to the variables containing the total gas mass and the mass ejected per element::

  StP[m].TotalEjectedGasMass = SingleEjectedMass[0]*All.HubbleParam;
  
  for (k=0;k<NELEMENTS;k++) 
    StP[m].TotalEjectedEltMass[k] = SingleEjectedMass[k+2]*All.HubbleParam;    




Computing ejected energy
""""""""""""""""""""""""

Finally, the total ejected energy per mass unit is obtained by multiplying the number of SNIa and SNII
by the assumed supernova energy::

  StP[m].TotalEjectedEgySpec = All.ChimieSupernovaEnergy*  (NSNIa + NSNII) /StP[m].TotalEjectedGasMass; 


Variables used for the ejection of mass and energy
""""""""""""""""""""""""""""""""""""""""""""""""""

The variables used for the ejection of mass and energy are::

  StP[m].TotalEjectedGasMass
  StP[m].TotalEjectedEltMass[k]
  StP[m].TotalEjectedEgySpec 
  StP[m].NumberOfSNIa 
  StP[m].NumberOfSNII