COMET AND ASTEROID IMPACT HAZARD MODEL - PYTHON VERSION

1 INTRODUCTION

In 1999 John Lewis (University of Arizona) wrote a book "Comet and Asteroid Hazards on a Populated Earth" (Academic Press) which is a presentation and discussion of a simple computer model he developed to explore the subject. The book includes full computer code (both in print form and on a disc) that allows the reader to explore the model. The simulation is fully stochastic in that any single run produces a different result dependent upon several variables that are randomly chosen by the program. The original code was written in the GW-BASIC language. Because this language is now rarely used and access to a GW-BASIC interpreter that will run on current computers is also rare, we have rewritten the program in the widely used and freely available PYTHON language.

2 THE MODEL - A BRIEF DESCRIPTION

The factors applicable to the model are summarised below. See the book for more detail.

Mass: The mass is calculated based on logarithmic transformations of the time step (the time increment) and influenced by a uniform random factor.

Type: 7 different asteroid types that can appear in the simulation (iron, pallasite, mesosiderite, achondrite, ordinary chondrite, CM chondrite, CI chondrite) and 2 different comet types (long period, periodic). Each have a different likelihood which is influenced by the mass. Corresponding physical and chemical properties are defined upon allocation of type.

Energy: assigns velocity and kinetic energy based on type. Also randomly assigns an entry elevation angle. energy dissipation is determined by:

dE = 0.5 V^2 dm + m V dV

Blowoff: evaluates if explosive blowoff of the atmosphere takes place and establishes a flag, while estimating the proportion of incident mass that is expelled.

Fragment: analyses the atmospheric effects on the entering body using a loop. Models various changes over the entry process including mass, velocity and ablation using:

dV/dt = -Cd ρ V² / ( 2 A ) + g sin( θ ) / m

dm = 0.413 A Γ ρ V³ dt / ( Hvap )

Uses these values to arrive at a fate for the body (burnup, fragmentation, deceleration, impact, orbit, escape). The code uses an exponential approximation for the air density.

Crater: If the object is an impactor or releases high-mass fragments, the crater diameter and depth, area affected and ejected mass are calculated based on impactors’ energy.

Hazard: for impactors, assigns an impact location (oceanic or continental), this is used to influence fatalities caused from cratering or tidal waves. In addition, blast injuries from airbursts and the fire hazard from fireball heating are assessed.

3 MONTE CARLO SIMULATION

Monte Carlo methods transform uniform random numbers into specific probability distributions.

In GW-BASIC, this requires manual transformations for each variable. In Python, libraries like numpy and scipy.stats provide built-in distributions, this code uses random.random() to generate uniform random numbers in the range [0,1]. These are then transformed to match the probability distributions of different impact parameters. Specifically:

• θ (Impact Angle) – Converted to match empirical distributions.
• V (Velocity at Impact) – Mapped to a speed distribution based on asteroid/comet populations.
• D (Diameter of Impactor) – Drawn from a power-law distribution to reflect size-frequency relationships.
• ρ (Density of Impactor) – Assigned based on the material composition (asteroid vs. comet).
• H (Initial Altitude) – Set depending on atmospheric entry conditions.

These transformations ensure the generated samples reflect real-world impact scenarios rather than a uniform spread.

4 THE MODEL IN PYTHON CODE

Due to limited access of a BASIC interpreter/compiler that will run on a newer computer we have ported the program code to the popular PYTHON language and present it here.

The first program neohazards.py outputs the following simulation parameters to a CSV file:

• "Type": The type of object (e.g., asteroid or comet).
• "Mass (kg)": The entry mass of the object in kilograms (converted from original units).
• "M(H2O)": The mass of water delivered by the object in kilograms.
• "Velocity (km/s)": The entry velocity of the object in kilometers per second.
• "Energy (J)": The kinetic energy of the object in joules.
• "Iridium (µg/cm2)": The average iridium injected by the object over a global area.
• "F": An indication of the object’s fate
• "Elev (deg)": The elevation angle of the object’s entry.
• "Zmin (km)": the altitude at which an object that burns up or fragments last broke up
• "Vfinal (km/s)": The final velocity of the object after the event.
• "Mmax": peak brightness of the fireball in stellar magnitude as seen from the ground directly beneath.
• "FATAL": The number of fatalities
• "Density (kg/m^3)": The density of the object in kilograms per cubic meter.
• "HVAP (J/kg)": The heat of vaporization of the object in joules per kilogram,
• "Radius (m)": The radius of the object in meters.
• "Crater (m)": The radius of the impact crater in meters if the object impacted Earth.
• "RFRAGSTR": A measure of the fragmentation strength of the object.
• "MASS LEFT (%)": The percentage of the original mass of the object that remains after the impact event.

The second program neograph.py is used to produce graphs of NEO atmospheric entry events. This code will plot the velocity, mass and height for specified individual events from the CSV file.

5 VALIDATION OF THE PYTHON CODE

To validate the Python code against the GW-BASIC code, I replaced all instances of random number generation (RND in GW-BASIC and random.random() in Python) with a fixed value (e.g., 0.5). This removed randomness, allowing for a direct comparison of the results. By running both simulations with the same fixed values, I compared the outputs step by step. Since the outputs matched, it confirmed that the Python code replicated the functionality of the GW-BASIC code correctly.

6 MODEL OUTPUT

The output of the program "neohazards.py" produces analysis.csv, a CSV (comma separated variable)file that can be used by Excel or another program.

7 THE ATMOSPHERIC ENTRY PROCESS

Three specific examples were used for plotting velocity, mass and height during reentry (using the python code neograph.py) :

Comet entering at 60km/s

- comet / periodic
- density 1100 kg/m³
- mass: 698,349 kg

Iron Asteroid entering at 15km/s

- asteroid / iron
- density 7900 kg/m³
- mass: 207,758 kg

Ordinary Chondrite Asteroid entering at 30km/s

- asteroid / Ordinary Chondrite
- density 3610 kg/m³
- mass: 3,915,096 kg

The 3 graphs below plot the height, velocity and mass of the comet in the first example given above during the entry process into the Earth' atmosphere.

Australian Space Academy