CasEx

Implementation of the mathematics behind JARUS WG-6 SORA Annex F.

Overview

This package was created to support the mathematical models used to derive the iGRC tables and supporting material found in JARUS WG-6 SORA Annex F [3]. It implements the models for computing critical areas, effect of obstacles, ballistic descent, effect of explosions and deflagration that are being used in Annex F. It also includes easy access to the specific choice of parameters that are used in Annex F.

There are a number of examples of how to use this package that operators and others can use as a starting point for doing their own computations as an alternative to using the iGRC table and supporting tables in Annex F and in the SORA.

Finally, this packages can be used for easily recreating some of the figures found in Annex F (see example 8).

The mathematical models implemented in the package are the results of the combined effort of the Quantitative Methods subgroup under Working Group 6.

Installation

To install CasEx, simply run the following command in your terminal.

pip install casex

Background

Since the summer 2019, a small group of people have been working under the auspice of JARUS on developing guidance material for use to determine the ground risk associated with operations of unmanned aircraft. This group is called the Quantitative Methods subgroup. The developed concept had to be applicable to most types of UAS and it has to be relatively simple to allow for everyday use by operators and authorities.

While the work of this group is still onging, the results so far of this development is described in Annex F [1] in the SORA set of document developed and published by JARUS. This Annex provide a elaborate explanation of how to determine the ground risk for a given operation, where a number of aircraft parameters and the density of people on the ground is known.

As part of the development of Annex F, the group has gone through quite a lot of programming in Matlab and Python to test out difference ideas, verifying assumptions, plotting graphs, etc. Some of this work might be useful to other people, and may also serve to document how the various models and math implementation have been used in the development of Annex F.

Therefor, the group decided to spent the effort to convert some of the code into a publicly available toolbox that can be used to apply the methods described in Annex F without having to implement it all from scratch. The toolbox also reproduces some of the figures in the Annex.

CasEx is the result of this implementation. It is intended only as a supplement to Annex F, and it is still the responsibility of the operator or authority using the toolbox to make sure that any computed value and use of such values is in compliance with Annex F and any other requirements.

Still, we hope that the toolbox can be useful. Comments and queries related to the implementation and use can be submitted to the main author Anders la Cour-Harbo at at alc(a)es.aau.dk.

1

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

Contributors

The mathematical models implemented in the package are the results of the combined effort of the Quantitative Methods subgroup under Working Group 6 in JARUS. The effort is a combination of gathering and unifying existing models and of developing models adapted to the specific needs of Annex F, including the JARUS critical area model and the obstacle model.

The group has the following core members:

• Terrence Martin (chair), Queensland University of Technology, Australia

• Henri Hohtari, AirMap, Switzerland

• Tony Nannini, Wing, USA

• Jörg Dittrich, DLR, Germany

• Anders la Cour-Harbo, Aalborg University, Denmark

• Thomas Putland, TASDCRC, Australia

The group has used various tools for developing the mathematical content of Annex F, including Matlab and Python. CasEx is a publicly available, streamlined, fully documented version of this work, and it can be installed as any other Python package.

The development of CasEx as a Python package was done by:

Lead programmer and math modelling

Anders la Cour-Harbo

Math modelling and Matlab implementations

Terrence Martin

Math support

Thomas Putland

Documentation layout and package integration

Christian Schmidt Godiksen

Examples

Example 1: Critical area

This first example demonstrates how to use the CasEx package for computing a critical area for a specific aircraft. This computation can be used for determining a more accurate iGRC value than what is possible for the iGRC table.

We start by setting two basic parameters required for computing the size of the critical area, namely the height of a standad person and the “radius” of a standard person. The latter is the radius of a circle circumscribing a person seen from above. It is possible to leave out these values when instantiating the critical area model, and it will revert to the standard value, which are the same as we use here (and also the values used throughout Annex F).

person_radius = 0.3
person_height = 1.8

We then instantiate first the class that can provide some friction coefficients, and then the class that holds the various models for computing the size of the critical area.

FC = FrictionCoefficients()
CA = CriticalAreaModels(person_radius, person_height)

We then need to decide the impact speed and the impact angle. The former is the speed the aircraft will have along the direction of travel at the time it impacts the ground. This is used for computing the glide distance after impact. The latter is the angle of the straight line that the aircraft follows just prior to impact. This is measured from horizontal, so that 10 degrees is a very shallow impact, while 90 is a vertical impact. For this example, we chose a speed of impact of 45 m/s and a 35 degree impact angle .

impact_speed = 45
impact_angle = 35

We also need to specify some properties of the aircraft. This is done using the AircraftSpecs class. But first, we set some of the parameters used in the instantiation of the class.

We have the wingspan (also known as width or characterstic dimension), the mass, and the coefficient of friction between the aircraft and the terrain. The latter is in this case captured from the FrictionCoefficients, but could also be set manually. Note that in Annex F this value is always 0.65.

aircraft_type = enums.AircraftType.GENERIC
width = 2.2
mass = 5
friction_coefficient = FC.get_coefficient(enums.AircraftMaterial.RUBBER,
                                          enums.GroundMaterial.CONCRETE)

Then we instantiate the aircraft class, and set the friction value.

aircraft = AircraftSpecs(aircraft_type, width, mass)
aircraft.set_friction_coefficient(friction_coefficient)

The deflagration model can be included by setting the type of fuel and the amount of fuel. Here, we set the amount to 0, which effectively deactivates the deflagration model. It can be activated by setting a non-zero amount of fuel.

aircraft.set_fuel_type(enums.FuelType.GASOLINE)
aircraft.set_fuel_quantity(0)

The coefficient of restitution is set as a function of the impact angle.

aircraft.set_coefficient_of_restitution(AnnexFParms.CoR_from_impact_angle(impact_angle))

When the deflagration model is used there will be an overlap between the critical area from the glide and slide and the critical area from the deflagration (since the deflagration happens to some extend in the same area). We therefore need to set the fraction of overlap. The total critical area (sum of glide/slide and deflagration) is reduced by the smaller of the two multiplied with the overlap factor.

Note that this value is not used when there is no fuel onboard, and therefore no deflagration.

critical_areas_overlap = 0.5

During the slide, the aircraft dissipates kinetic energy, which ends at 0 when the aircraft comes to a halt. It is possible to reduce the size of the critical area by assuming that the aircraft is no longer lethal at a higher kinetic energy than 0. The chosen value is set here. If the default value is to be used, the value is set to -1. The documentation for critical_area_models has more info on this.

lethal_kinetic_energy = -1

We are now ready to compute the size of the critical area. We do this by calling the critical_area method.

p = CA.critical_area(aircraft, impact_speed, impact_angle,
                     critical_areas_overlap, lethal_kinetic_energy)

The output from this method is a list with nine values. These are 1) total critical area, 2) glide area, 3) slide area, 4) the inert critical area, 5) the critical area from the deflagration, 6) the glide distance, 7) the slide distance until the kinetic energy has reached the non-lethal limit, 8) the speed at which the sliding aircraft has the kinetinc energy from 7), 9) the time it takes to reach the kinetic energy from 7). The first five are shown in the output below along with some of the aircraft parameters.

Finally, we can also compute the raw iGRC as described in Annex F. This is done using the AnnexFParms class. Here we need to input the population density, which is set to 3000 here.

raw_iGRC = AnnexFParms.iGRC(3000, p[0])

This iGRC value can now be used in the SORA process. The value should be rounded up to the nearest larger integer.

The output of the example is as follows. We see that the raw iGRC is 6.4, which must be rounded up to 7 for the actual iGRC value. This demonstrates how an aircraft, which belongs in the third column (due to the high impact speed above 35 m/s) can achieve a lower iGRC value than given in the table (which would be 8) by doing the actual calculations for the critical area.

Wingspan:                   2.2 m
Mass:                       5 kg
Horizontal impact speed:    36.9 m/s
Glide area:                 6.2 m^2
Slide area:                 78.1 m^2
Critical area inert:        84.2 m^2
Critical area deflagration: 0.0 m^2
Total critical area:        84 m^2
Raw iGRC:                   6.4

It is important to note that this example uses rubber against concrete as the friction coefficient. This obviously has to be adjusted depending on the aircraft and the overflown terrain. Alternatively, it is always acceptable to use the Annex F standard value of 0.65.

Example 2: Vector input

This example shows how to use the JARUS model with a vector input to generate a vector output. Using this feature it is possible to generate graphs for variations on one input parameter without resorting to a loop.

This example shows this for the following parameters:

• width

• impact angle

• impact speed

• lethal area overlap

• fuel quantity

• friction coefficient

We begin with setup, which was also done in example 1.

# Instantiate necessary classes.
FC = FrictionCoefficients()
CA = CriticalAreaModels()

# Set aircraft values.
aircraft_type = enums.AircraftType.FIXED_WING
width = 4
mass = 25
friction_coefficient = FC.get_coefficient(enums.AircraftMaterial.ALUMINUM,
                                          enums.GroundMaterial.CONCRETE)

# Instantiate and add data to AircraftSpecs class.
aircraft = AircraftSpecs(aircraft_type, width, mass)
aircraft.set_fuel_type(enums.FuelType.GASOLINE)
aircraft.set_fuel_quantity(5)
aircraft.set_friction_coefficient(friction_coefficient)

# Set impact speed and angle.
impact_speed = 50
impact_angle = 25

First, we set the width to an array instead of a scalar. Here we choose the width to vary from 1 to 5 meters over 100 steps. The input is set via the aircraft class. We collect the result in the p array.

v_width = np.linspace(1, 5, 100)
aircraft.set_width(v_width)
p.append(CA.critical_area(aircraft, impact_speed, impact_angle))

The we set the impact speed as a vector. Since only one parameters at a time can be a vector input, we reset the width to a scalar value.

v_impact_angle = np.linspace(10, 80, 100)
aircraft.set_width(width)
p.append(CA.critical_area(aircraft, impact_speed, v_impact_angle))

Since the impact angle is not set through the aircraft class, but directly as an input to critical_area, we do not need to reset the impact angle. We simply use the original variable impact_angle. And we set the impact speed to be an array.

v_impact_speed = np.linspace(5, 40, 100)
p.append(CA.critical_area(aircraft, v_impact_speed, impact_angle))

We do the same with the critical area overlap.

v_critical_areas_overlap = np.linspace(0, 1, 100)
p.append(CA.critical_area(aircraft, impact_speed, impact_angle))

The fuel quantity is set as an array.

v_fuel_quantity = np.linspace(0, 10, 100)
aircraft.set_fuel_quantity(v_fuel_quantity)
p.append(CA.critical_area(aircraft, impact_speed, impact_angle))

We can also set the friction coefficient as an array.

aircraft.set_fuel_quantity(5)
v_friction_coefficient = np.linspace(0.2, 0.9, 100)
aircraft.set_friction_coefficient(v_friction_coefficient)
p.append(CA.critical_area(aircraft, impact_speed, impact_angle))

Finally, we plot all the outputs.

Example 3: iGRC table

This example reproduces the nominal iGRC table in Annex F [1]. It also allows for seeing the table with different impact angles, for ballistic descent, as well as with and without both obstacles and the conservative assumption of -0.3 (see details in Annex F).

First, we call the AnnexFTables class (the called method is static).

console_output = AnnexFTables.iGRC_tables(show_with_obstacles = True,
                                          show_ballistic = False,
                                          show_with_conservative_compensation = False,
                                          show_glide_angle = False,
                                          show_additional_pop_density = True)

The returned value is a list of strings to print to the console, which we then do.

for s in console_output:
    print(s)

iGRC table (35 deg)
                          1 m       3 m       8 m      20 m      40 m
--------------------+-------------------------------------------------
Controlled (< 0.1)  |     1.1       2.3       3.3       4.3       5.1
      8 [ppl/km^2]  |     2.6       3.8       4.8       5.8       6.6
     25 [ppl/km^2]  |     3.1       4.3       5.3       6.3       7.1
     75 [ppl/km^2]  |     3.6       4.8       5.8       6.8       7.6
    250 [ppl/km^2]  |     4.1       5.3       6.3       7.3       8.1
    750 [ppl/km^2]  |     4.6       5.8       6.8       7.8       8.6
   2500 [ppl/km^2]  |     5.1       6.3       7.3       8.3       9.1
   7500 [ppl/km^2]  |     5.6       6.8       7.8       8.8       9.6
  25000 [ppl/km^2]  |     6.1       7.3       8.3       9.3      10.1
 250000 [ppl/km^2]  |     7.1       8.3       9.3      10.3      11.1
 750000 [ppl/km^2]  |     7.6       8.8       9.8      10.8      11.6
 750000 [ppl/km^2]  |     7.6       8.8       9.8      10.8      11.6
2500000 [ppl/km^2]  |     8.1       9.3      10.3      11.3      12.1
--------------------+-------------------------------------------------
     CA size [m^2]  |       5        85       878      8190     47875
Slide distance [m]  |       0        26       106       403      1178
--------------------+-------------------------------------------------

1

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

Example 4: iGRC tradeoff tables

This example produces the tradeoff tables in Annex F [1].

We want to see the T2 tradeoff table in rounded up integers. First, we set the tradeoff value

tradeoff_type = 2

and request integers

show_integer_iGRC = True

and not relative values

show_relative = False

If showing relative values, the table will show the difference to the nominal iGRC table. This feature is primarily for developers that want to see how big the difference is when using trade-offs. It does not have any bearing on the use of the trade-off tables.

We then call the AnnexFTables method

console_output = AnnexFTables.iGRC_tradeoff_tables(tradeoff_type,
                                                   show_integer_iGRC, show_relative)

and run

for s in console_output:
    print(s)

to show the result

Modified raw iGRC table for trade-off scenario T2
             Max dim   |   2.0 m     6.0 m    16.0 m    40.0 m    80.0 m
Dpop         Max speed |  25 m/s    35 m/s    75 m/s   150 m/s   200 m/s
-----------------------+--------------------------------------------------
         Controlled    |      1         2         3         4         5
        12 ppl/km^2    |      3         4         5         6         7
       125 ppl/km^2    |      4         5         6         7         8
      1250 ppl/km^2    |      5         6         7         8         9
     12500 ppl/km^2    |      6         7         8         9        10
    125000 ppl/km^2    |      7         8         9        10        11
    125000 ppl/km^2    |      7        n/a       n/a       n/a       n/a

1

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

Example 5: Isoparametric critical area plot

This example makes a colored plot of the lethal area for varying impact angles and impact speeds. The size of the aircraft as well as the other parameters are fixed. The idea is to visualize the relation between angle and speed, since this is one of the challenges in understanding the transition from the JARUS model to the iGRC table.

The target for this example is the second column in the iGRC table, where the defined lethal area is 80 m^2. Therefore, the isoparametric curve for 80 m^2 is shown in yellow along with other iso curves in white for comparison.

We start by setting up the critical area class.

CA = CriticalAreaModels()

We instantiate the AnnexFParms class, since we want to use the parameters for the smallest size aircraft in the iGRC table, and we can find them in that class.

AFP = AnnexFParms()

We then setup the aircraft based on the parmaters for the first column in the iGRC table (thus the index [0] on AFP.CA_Parms). Note that this is not specific for fixed-wing, since this information is not currently being used by the critical area computation.

aircraft_type = enums.AircraftType.FIXED_WING
aircraft = AircraftSpecs(aircraft_type, AFP.CA_parms[1].wingspan, AFP.CA_parms[1].mass)

We do not use any fuel.

aircraft.set_fuel_type(enums.FuelType.LION)
aircraft.set_fuel_quantity(0)

And we use the friction coefficient from Annex F.

aircraft.set_friction_coefficient(AFP.friction_coefficient)

We want to plot over the full range of impact speed and impact angles. However, since very shallow impact angles are not handled well by the model, we start at 5 degrees. For speed, we cap it at 40 m/s. Each axis will be 100 steps, which is fine for a relatively smooth plot.

y_speed = np.linspace(0, 50, 100)
x_angle = np.linspace(5, 70, 100)
X_angle, Y_speed = np.meshgrid(x_angle, y_speed)

We then compute the critical area for all combinations of speed and angle. Note that we could have replaced one of the loops with an array input, but since we cannot replace both (since the critical_area method does not support 2D array input), we have chosen to run both dimensions as loops for code clarity.

Z_CA = np.zeros((X_angle.shape[0], Y_speed.shape[0]))
for i, y in enumerate(Y_speed):
    for j, x in enumerate(X_angle):
        Z_CA[i, j] = CA.critical_area(aircraft, y_speed[i], x_angle[j])[0]

The plot is setup with room for a colorbar on the right.

fig = plt.figure()
ax = plt.axes()
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.05)

The matrix is added to the plot.

im = ax.imshow(Z_CA, extent=[x_angle[0], x_angle[-1], y_speed[0], y_speed[-1]],
               aspect='auto', origin='lower')

Conturs are added to the plot.

CS = ax.contour(X_angle, Y_speed, Z_CA, [5, 10, 15, 20, 30, 40, 60, 100, 120], colors='white')
ax.clabel(CS, inline=1, fontsize=10, fmt="%u m^2")

A yellow contour is added for the target critical area for the first column in the iGRC table.

CS2 = ax.contour(X_angle, Y_speed, Z_CA, [AFP.CA_parms[1].critical_area_target], colors='yellow')
ax.clabel(CS2, inline=1, fontsize=10, fmt="%1.1f m^2")

And the colorbar is added along with axes labels and title.

fig.colorbar(im, cax=cax, orientation='vertical')
ax.set_xlabel('Angle [deg]')
ax.set_ylabel('Speed  [m/s]')
ax.set_title('Lethal area [m^2] for {:d} m, {:d} kg, {:2.2f} friction'.format(AFP.CA_parms[1].wingspan,
                                                                          AFP.CA_parms[1].mass,
                                                                          AFP.friction_coefficient))

The contours in the output image show where there is a constant critical area for varying combinations of speed and angle.

There is a significant “drop” visible for the lower speeds. This is coming from lethal kinetic energy limit described in Annex F, i.e., the disregard of the slide part of the critical area where the slide speed is small enough to have a kinetic energy less than the lethal limit. When the kinetic energy at the start of the slide is less than the lethal kinetic energy, only the glide part is lethal, and this does not depend on speed, only on impact angle. This results in the CA remaining constant for changing speeds for sufficiently low speeds. Consequently, the iso curve becomes vertical in this area.

Example 6: Ballistic descent

This example shows how to compute a series of values for a ballistic descent. The same approach is used in Annex F for computing the ballistic values found in Appendix A.

The class BallisticDescent2ndOrderDragApproximation is based on [2], which also describes the detail on how the ballistic descent model works.

We first instantiate the class.

BDM = BallisticDescent2ndOrderDragApproximation()

We set the standard aircraft parameters

aircraft_type = enums.AircraftType.FIXED_WING
width = 2.8
mass = 90

and instantiate the aircraft class.

aircraft = AircraftSpecs(aircraft_type, width, mass)

We need to set the drag coefficient and the frontal area, since they are used in the computation of the ballistic descent. The drag coefficient is typically in the range from 0.7 to around 1 for an aircraft which is not descending “nicely”, i.e., in its usually direction of travel. Here we assume a normal helicopter-like rotorcraft, which may be tumbling during descent. We guestimate the coeffient to be 0.8. A lower coefficient would give a higher impact speed, and thus be more conservative. But 0.8 is not unreasonable for a tumbling rotorcraft. See for instance the wiki page on drag coefficient for more detail on drag for various shapes.

The frontal area is the area covered by the aircraft in the direction of travel during descent. Here we guess that it will be 60 cm by 60 cm, which is reasonable for a 90 kg aircraft.

aircraft.set_ballistic_drag_coefficient(0.8)
aircraft.set_ballistic_frontal_area(0.6 * 0.6)

The ballistic descent class must “be aware” of the aircraft.

BDM.set_aircraft(aircraft)

We also set the initial values for the descent, namely the altitude above the ground and the velocity of the aircraft at the beginning of the descent.

altitude = 100
initial_velocity_x = 28
initial_velocity_y = 0

We can now compute the ballistic descent.

p = BDM.compute_ballistic_distance(altitude, initial_velocity_x, initial_velocity_y)

Note that p is a list with various values about the descent. Note also that

BDM has additional values (attributes) available, relating to the intermediate

computations as decribed in [2]. All available values are printed to the screen here.

print("Distance:     {:1.1f} m    ".format(p[0]))
print("Impact speed: {:1.1f} m/s".format(p[1]))
print("Angle:        {:1.1f} degs".format(p[2] * 180 / np.pi))
print("Time :        {:1.2f} s".format(p[3]))
print("Distances:    {:1.3f}  {:1.3f}  {:1.3f}".format(BDM.distance1, BDM.distance2, BDM.distance3))
print("Times:        {:1.3f}  {:1.3f}  {:1.3f}".format(BDM.time_top, BDM.time_cross, BDM.time_impact))
print("Speeds:       {:1.3f}  {:1.3f}".format(BDM.velocity_x, BDM.velocity_y))

which gives

Distance:     115.7 m
Impact speed: 45.7 m/s
Angle:        61.9 degs
Time :        4.66 s
Distances:    0.000  74.245  41.461
Times:        -0.000  2.854  4.664
Speeds:       21.478  40.288

It is possible to do ballistic descent computations with vector input. We do this by setting the initial velocity to an array.

altitude = 150
initial_velocity_x = np.linspace(0, 80, 100)
initial_velocity_y = np.linspace(-10, 10, 100)

However, we cannot have both horizontal and vertical velocities be arrays at the same time, so we use them individually.

p_vel_x = BDM.compute_ballistic_distance(altitude, initial_velocity_x, -2)
p_vel_y = BDM.compute_ballistic_distance(altitude, 30, initial_velocity_y)

The drag coefficient can also be an array.

drag_coef = np.linspace(0.7, 1, 100)
aircraft.set_ballistic_drag_coefficient(drag_coef)
BDM.set_aircraft(aircraft)
p_drag_coef = BDM.compute_ballistic_distance(altitude, 30, -2)

Now p contains the output from each of the three calculations, and this can be plotted to show the variation of the output. Here are some of the possible combination of input and output parameters.

Finally, this example also shows how to get the ballistic computations found in Annex F Appendix A. This is done by instantiating the AnnexFParms class, since all the ballistic values are computed internally in the class during instantiation.

AFP = AnnexFParms()

Then we can easily access the associated attributes in the class to get the various ballistic values for each category.

for c in range(5):
    print("{:2d} m      {:3d} m/s     {:4d} m    {:3.0f} m/s   {:3.0f} m/s   {:1.0f} deg    {:4.0f} m     "
          "{:4.1f} s   {:6.0f} kJ".format(
        AFP.CA_parms[c].wingspan, AFP.CA_parms[c].cruise_speed, AFP.CA_parms[c].ballistic_descent_altitude,
        AFP.CA_parms[c].terminal_velocity, AFP.CA_parms[c].ballistic_impact_velocity,
        AFP.CA_parms[c].ballistic_impact_angle, AFP.CA_parms[c].ballistic_distance,
        AFP.CA_parms[c].ballistic_descent_time, AFP.CA_parms[c].ballistic_impact_KE / 1000))

The printed values

Ballistic descent computations
------------------------------
Class   Init horiz     From     Terminal   Impact   Impact   Distance   Descent      KE
          speed      altitude   velocity   speed    angle    traveled    time      impact
 1 m       25 m/s       75 m     25 m/s    24 m/s   76 deg      63 m      4.7 s        1 kJ
 3 m       35 m/s      100 m     45 m/s    39 m/s   64 deg     123 m      4.9 s       39 kJ
 8 m       75 m/s      200 m     63 m/s    60 m/s   57 deg     335 m      6.9 s      718 kJ
20 m      150 m/s      500 m    112 m/s   106 m/s   51 deg    1044 m     10.8 s    27952 kJ
40 m      200 m/s     1000 m    120 m/s   121 m/s   59 deg    1690 m     16.0 s    73038 kJ

are the same as found in the ballistic table in Appendix A of Annex F [1].

1

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

2(1,2)

Anders la Cour-Harbo. Ground impact probability distribution for small unmanned aircraft in ballistic descent. In 2020 International Conference on Unmanned Aircraft Systems (ICUAS), 1442–1451. Aalborg University, IEEE, sep 2020. doi:10.1109/ICUAS48674.2020.9213990.

Example 7: JARUS model vs the iGRC values

In this example we compare the iGRC table in Annex F with the actual output of the JARUS model. This is done by plotting the critical area for a variety of wingspan and population densities. This is similar to the iGRC figures that can be computed with the Figures class, which do the same thing, but for the simplified model.

First, we need to set the impact angle. This is not used in the simplified model, where the impact angle has been simplified away. However, for the JARUS model we need to set this, since the critical area depends on it. Here we choose default value from AnnexFParms (35 degrees) since this angle has a specific meaning in Annex F. But any angle will do, but obviously produce a slightly different plot. Also note that the impact speed may vary with the impact angle.

impact_angle = AnnexFParms.scenario_angles[1]

We instantiate the CA model class.

CA = CriticalAreaModels()

We set the sampling density for the two axes.

pop_density_samples = 400
wingspan_samples = 200

We set the arrays for the sampling on the two axes. Note that the wingspan axis is linear, while the population density axis is logarithmic.

wingspan = np.linspace(0, 21, wingspan_samples)
pop_density = np.logspace(-3 + np.log10(5), 6, pop_density_samples)

We instantiate the Annex F parmeters class AnnexFParms.

Initialize the background matrix with zeros.

GRC_matrix = np.zeros((pop_density_samples, wingspan_samples))

We run a for loop over the range of the wingspan.

for j in range(len(wingspan)):

Now, for every value of the wingspan, we want to use the aircraft values associated with the wingspan, i.e., the appropriate column in the iGRC table.

    if wingspan[j] <= AFP.CA_parms[0].wingspan:
        column = 0
    elif wingspan[j] <= AFP.CA_parms[1].wingspan:
        column = 1
    elif wingspan[j] <= AFP.CA_parms[2].wingspan:
        column = 2
    elif wingspan[j] <= AFP.CA_parms[3].wingspan:
        column = 3
    else:
        column = 4

Since we previously set the impact angle to 35 degrees, we will use the cruise speed as impact speed.

    impact_speed = AFP.CA_parms[column].cruise_speed

However, if we had chosen another angle, we would consider using another speed. Specificallty, in Annex F it is assumed that for an impact at 10 degrees, the speed is reduced equivalent to multiplying with a factor 0.65. This glide speed value can also be obtained from the AFP class as AFP.CA_parms[column].glide_speed.

We set the wingspan according to the maximum value from each class.

    AFP.CA_parms[column].aircraft.width = wingspan[j]

We then compute the size of the critical area using the standard values from Annex F. Since there is no deflagration, the overlap is set to zero.

    M = CA.critical_area(AFP.CA_parms[column].aircraft, impact_speed, impact_angle)[0]

Finally, we populate the background matrix for the plot, that is the iGRC values associated with the combination of population density and critical area.

    for i in range(len(pop_density)):
        GRC_matrix[i, j] = AnnexFParms.iGRC(pop_density[i], M)[0]

Note that the [0] at the end selects the raw iGRC value, while a [1] would select the rounded up iGRC value.

This is followed by numerous lines for creating the plot. This includes showing the matrix along with contours, setting the tick locations and labels, and white lines and text for overlaying the iGRC. The result is as follows.

The iso-parametric curves in black show where the iGRC value has exactly integer values in the background plot. So for instance the curve for iGRC 9 is the boundary between values below 9 and above 9.

The percent of area in Denmark which corresponds to the different population density bands is also shown. These values are based on locations of addresses and a grid size of 1 km by 1 km. This information does not come from CasEx, and is just hardcoded into this example.

Example 8: Obstacles

The use of obstacles for reducing the critical area is demonstrated in this example. The nominal critical area for the 3 m category is 140 \(\mathrm{m}^2\), and we want to see how much the critical area is reduced for an obstacle density of 800 houses per \(\mathrm{km}^2\).

The critical area of 140 \(\mathrm{m}^2\) is associated with the second column in the iGRC, which is for the 3 m wing span. Therefore, we set the width of the CA to 3.

CA_width = 3

Consequently, the length of the CA must be 140 divided by the width.

CA_length = 140/CA_width

We set the density to 800 houses per \(\mathrm{km}^2\), and convert to a density measure in obstacles per square meter.

houses_per_square_km = 800
obstacle_density = houses_per_square_km / 1e6

The size of the houses follows a 2D normal distribution for width and length. We set the mean values and standard deviations for width

obstacle_width_mu = 23
obstacle_width_sigma = 6

and for length

obstacle_length_mu = 9
obstacle_length_sigma = 2

We instantiate the Obstacles class. The third input argument is only used for simulation, which we don’t do here. So it is just set to zero.

We want to compute the probability for varying length of the critical area, so we generate an array from 0 to the length of the nominal CA. Since the resulting curve is quite smooth, it is not necessary to use all that many points. We have chosen 10 here, which is usually sufficient.

x_resolution = 10
x = np.linspace(0, CA_length, x_resolution)

The computations of the cumulatative density function there are three integrals that have to be evaluated. We must specify at what resolution this happens. A value of 25 is typically sufficient. Increasing it will not produce notiable different results.

pdf_resolution = 25

We now compute the CDF curve.

p_x, EX, beta_analytical, sanity_check = OS.cdf(x, obstacle_density, obstacle_width_mu, obstacle_width_sigma,
                                            obstacle_length_mu, obstacle_length_sigma, pdf_resolution)

The resulting curve is plotted.

The blue curve shows the probability of the CA being reduced to a given value or less. The curve noticable does not go to 100%, which may seem odd, since obviously the probability of the CA being less then 140 is 100%. This is because the graph in reality should go vertically up from the current ending point around 76% to 100%, but this jump is not implemented in the code. This is caused by the fact that a certain percentage of the crashes do not impact an obstacles, and thus stays at nominal CA size. This percentage in this case would be around 24%.

Example 9: Recreating Annex F figures

THIS EXAMPLE IS NOT YET UPDATED TO SORA 2.5 AND WILL NOT WORK PROPERLY.

This example reproduces the figures in Annex F [1].

The various figures can be separately generated by calling the functions with the appropriate boolean settings as shown in this example. However, other variations of the figures can also be generated. The documentation for each figure is provided in Figures.

The four different figures produced are shown below.

1

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

Bug fixes and updates

Version 1.2.3

• Updated documentation for examples 1 through 8, including adjustments to the examples to fit good documentation.

• Example 9 on figures is still work in progress.

Version 1.2.2

• Updated documentation for examples. Still not complete.

Known issues still to be fixed:

• figures.py still needs more work to fit SORA 2.5.

• Documentation has only partially been updated to fit all the changes following SORA 2.5.

Version 1.2.1

• Removed models for RTI, RCC, NAVCAD, and FAA. These are no longer avaiable for computation. Also removed the associated enum.

• Length has been removed from AircraftSpecs, since it is no longer needed.

• Changed gravity from 9.82 to 9.81 to comply with MKS system.

• Rearrange the order of the examples.

Version 1.2.0

• Updated Casex to correspond to SORA 2.5. Most the following updates are adjustments to SORA 2.5, while most of the bug fixes are resulting from a general update of CasEx.

• Changed numerous values in annex_f_tables.py to correspond to the coming public version of SORA 2.5.

• Added ballistic_descent_table() to annex_f_tables.py (originally in figures.py) and expanded to show the full table.

Smaller changes

• Change upper proper CoF value in aircraft_specs.py from 1.5 to 1.0.

• Moved COR_from_impact_angle() from aircraft_specs.py to annex_f_parms.py and made it static.

• Added default values for person_radius (0.3) and person_height (1.8) to AnnexFParms.

• Change LA_inert to CA_inert and LA_deflagration to CA_deflagration in critical_area() method to reflect critial area instead of lethal area.

• Change slide_distance_lethal to slide_distance_non_lethal in critical_area() method, since it is in fact the non-lethal distance.

• Updated example 1 to comply with changes.

• Deleted example 2.

• Updated example 3 to comply with changes.

• Updated example 4 to include obstacle reduction.

• Updated example 7 to match new population bands.

• Updated JARUS model equations in critical_area() to accommodate concession of 35 degree impact angle.

• Added population bands to AnnexFParms.

• Development status change from beta to production/stable.

Version 1.1.11

• The disc-shaped end of the critical area was counted double, i.e., a full disc was added to both glide and slide. Now only half a disc is added to glide and slide.

• In the web interface for the online tool “Person width” is changed to “Person radius”.

• Updated values for calculation of CoR to fit SORA 2.5.

• The Obstacles simulation currently does not work. A fix is in progress.

• Thanks to Janis Mauch for point out these bugs.

Version 1.1.10

• First deployed version.

Reference

AircraftSpecs

Class to hold parameters on the aircraft used in area computations.

class casex.aircraft_specs.AircraftSpecs(aircraft_type, width, mass, fuel_type=<FuelType.GASOLINE: 1>, fuel_quantity=0)

This class is designed to hold all the parameters on a specific aircraft for which a critical area is to be computed. Many of the these parameters are not used in computations of critical area, but are reserved for future use.

Parameters

• aircraft_type (enums.AircraftType) – Type of aircraft.

• width (float) – [m] Width of aircraft (wingspan, characteristic dimension).

• mass (float) – [kg] Mass of the aircraft.

• fuel_type (:class:'enums.FuelType, optional) – Fuel type, such as fossil fuels or batteries (the default is FuelType.GASOLINE).

• fuel_quantity (float, optional) – [L] The quantity of fuel in liters. For batteries the quantity is also given in L, i.e. the volume of the battery (the default is 0, which means that no deflagration is assumed upon crash).

width

[m] The width of the aircraft is the horizontal size of the aircraft orthogonal to the direction of travel. This value is used to determine the width of the glide and slide areas. Therefore, this value is the wingspan for fixed wing aircraft, the rotor diameter for rotorcraft, and the rotortip to rotortip distance for multirotor aircraft.

Type

float

mass

[kg] Mass of the aircraft in kg. This is the total mass at the time of crash, including fuel.

Type

float

aircraft_type

The type of aircraft.

Type

enums.AircraftType

fuel_type

Fuel type, such as fossil fuels or batteries (the default is FuelType.GASOLINE).

Type

enums.FuelType

fuel_quantity

[L] The quantity of fuel in liters. For batteries the quantity is also given in L, i.e. the volume of the battery (the default is 0, which means that no deflagration is assumed upon crash).

Type

float

friction_coefficient

[-] Coefficient of friction between aircraft and ground. Appropriate values can be found using FrictionCoefficients (the default is 0.6).

Type

float

coefficient_of_restitution

[-] Coefficient of restitution expresses the loss of energy on impact (the default is 0.7).

Type

float

ballistic_frontal_area

[m^2] Frontal area of the aircraft during ballistic descent. This is the area size of the aircraft as projected in the direction of descent.

Type

float

ballistic_drag_coefficient

[-] Drag coefficient of the aircraft during ballistic descent. For future use.

Type

float

glide_drag_coefficient

[-] Drag coefficient of the aircraft during glide descent.

Type

float

max_flight_time

[s] Maximum flight time on the given amount of fuel. For future use.

Type

float

cruise_speed

[m/s] Cruise speed of the aircraft.

Type

float

glide_speed

[m/s] Glide speed of the aircraft when descending in a glide without thrust.

Type

float

glide_ratio

[-] Glide ratio of the aircraft when descending in an optimal glide angle without thrust.

Type

float

parachute_deployment_time

[s] Deployment time for the parachute, measured from the time deployment is signalled to full deployment.

Type

float

parachute_area

[m^2] Area of the parachute generating drag during descent and full deployment.

Type

float

parachute_drag_coef

[-] Drag coefficient.

Type

float

reset_values()

Reset all aircraft parameters.

Returns

Return type

None

set_aircraft_type(aircraft_type)

Set aircraft type.

The type of aircraft such as fixed wing or rotorcraft.

Parameters

aircraft_type (enums.AircraftType) – Type of aircraft.

Returns

Return type

None

set_ballistic_drag_coefficient(ballistic_drag_coefficient)

Set drag coefficient for ballistic descent.

Parameters

ballistic_drag_coefficient (float) – [-] Drag coefficient for ballistic descent.

Returns

Return type

None

set_ballistic_frontal_area(ballistic_frontal_area)

Set frontal area for computation of terminal velocity in ballistic descent.

Parameters

ballistic_frontal_area (float) – [m^2] Frontal area of the aircraft during ballistic descent.

Returns

Return type

None

set_coefficient_of_restitution(coefficient_of_restitution)

Set coefficient of restitution.

Parameters

coefficient_of_restitution (float) – [-] Coefficient of restitution for the ground impact.

Returns

Return type

None

set_cruise_speed(cruise_speed)

Set cruise speed.

Parameters

DOC (MISSING) –

Returns

Return type

None

set_friction_coefficient(friction_coefficient)

Set coefficient of friction between aircraft and ground.

Parameters

friction_coefficient (float) – [-] Coefficient for the glide resistance (the default is 0.6).

Returns

Return type

None

set_fuel_quantity(fuel_quantity)

Set the fuel quantity.

Parameters

fuel_quantity (float) – [L] Quantify of fuel. Set to 0 for no fuel in case no deflagration is desired in the computations.

Returns

Return type

None

set_fuel_type(fuel_type)

Set the type of fuel.

Sets the type of fuel onboard the aircraft. For a list of options, see enums.FuelType.

Parameters

fuel_type (enums.FuelType) – Type of fuel.

Returns

Return type

None

set_glide_drag_coefficient(glide_drag_coefficient)

Set drag coefficient for glide.

Parameters

glide_drag_coefficient (float) – [-] Drag coefficient for glide descent.

Returns

Return type

None

set_glide_speed_ratio(glide_speed, glide_ratio)

Set glide ratio.

Parameters

• glide_speed (float) – [m/s] Glide speed in the direction of travel.

• glide_ratio (float) – [-] The horizontal distance per vertical distance for the aircraft at glide speed.

Returns

Return type

None

set_mass(mass)

Set the aircraft mass.

Parameters

mass (float) – [kg] Mass of aircraft.

Returns

Return type

None

set_max_flight_time(max_flight_time)

Set max flight time.

Parameters

max_flight_time (float) – [s] Maximum flight time.

Returns

Return type

None

set_parachute(parachute_deployment_time, parachute_area, parachute_drag_coef)

Set parachute parameters.

Parameters

• parachute_deployment_time (float) – [s] Time from initiation of parachute deployment to a fully deployed parachute.

• parachute_area (float) – [m^2] Size of the parachute.

• parachute_drag_coef (float) – [-] Drag coefficient for the parachute.

Returns

Return type

None

set_width(width)

Set the aircraft width.

Parameters

width (float) – [m] Width of the aircraft

Returns

Return type

None

terminal_velocity()

Compute terminal velocity for free falling aircraft.

Returns

terminal_velocity – [m/s] The terminal velocity for the aircraft.

Return type

float

width_mass_check()

Out of range check on width and mass.

Performs an out of range check of width and mass. Warnings are issued for parameters out of range.

Returns

Return type

None

AnnexFParms

This class provides support for redoing some of the computations found in Annex F [1].

class casex.annex_f_parms.AnnexFParms

This class contains the following parameters for the 5 size classes in the iGRC table:

Parameters

impact_angle (float) – [deg] The impact angle of the descending aircraft, relative to horizontal.

aircraft

The class containing information about the aircraft.

Type

AircraftSpecs

aircraft_type

The type of aircraft. This parameters is not currently used.

Type

enums.AircraftType

critical_area_target

[m^2] Size of the largest critical area for each size class.

Type

float

ballistic_descent_altitude

[m] Assumed altitude for beginning of ballistic descent.

Type

float

ballistic_descent_time

[s] Computed descent time for ballistic descent.

Type

float

ballistic_distance

[m] Computed horizontal distance traveled during ballistic descent.

Type

float

ballistic_drag_coefficient = 0.8

[-] Drag coefficient used for ballistic descent.

Type

float

ballistic_frontal_area

[m^2] Assumed frontal area used in ballistic computations.

Type

float

ballistic_impact_angle

[deg] Computed impact angle with 0 being horizontal.

Type

float

ballistic_impact_KE

[J] Computed kinetic energy of aircraft just prior to impact.

Type

float

ballistic_impact_velocity

[m/s] Assumed horizontal velocity for beginning of ballistic descent.

Type

float

cruise_speed

[m/s]Maximum cruise speed for each size class.

Type

float

friction_coefficient = 0.65

[-] The friction coefficient is assumed constant at 0.65 throughout Annex F [1].

Type

float

glide_reduce = 0.65

[-] Reduction in glide speed relative to cruise speed.

Type

float

glide_speed

[m/s] The glide speed resulting from multiplying the cruise speed by glide_reduce.

Type

float

impact_angle

[deg] The impact angle of the aircraft when crashing, measure relative to horizontal.

Type

float

KE_critical

[J] Non-lethal energy during slide.

Type

float

mass

[kg] Assumed biggest mass for each size class.

Type

float

population_bands = [0.25 25 250 2500 25000 250000]

[ppl/km2] The population density bands used for the iGRC.

Type

float array

scenario_angles = [10, 35, 62]

[deg] The three impact angles for the three descent scenarios. The 62 degrees is not actually used but recomputed for each ballistic descent.

Type

float array

terminal_velocity

[m/s] Terminal velocity for aircraft.

Type

float

wingspan

Characteristic dimension of the aircraft. See Annex F [1] for more detailed explanation on what that is.

Type

float

static CoR_from_impact_angle(impact_angle, angles=None, CoRs=None)

Compute a coefficient of restitution for a given impact angle.

This method assumes an affine relation between impact angle and CoR. Therefore, two angles and two CoR values are used to determine this relation. The default is as described in Annex F that the CoR is 0.8 at a 10 degree impact and 0.6 at a 90 degree (vertical) impact. This values are used as defaults, but others can be specified.

Parameters

• impact_angle (float) – [deg] The impact angle between 0 and 90.

• angles (float array, optional) – [deg] Array with two different angle of impact values (the default is [10, 90]).

• CoRs (float array, optional) – [-] Array with two COR values corresponding to the two angles (the default is [0.8, 0.6]).

Returns

coefficient of restitution – [-] The coefficient of restitution for the given impact angle.

Return type

float

static applied_obstacle_reduction_factor(width)

Compute the obstacle reduction factor used in the iGRC in Annex F [1].

Parameters

pop_dens (float) – [ppl/km^2] Population density

Returns

obstacle_reduction_factor – The reduction factor for use in iGRC.

Return type

float

static iGRC(pop_dens, CA, TLOS=1e-06, use_conservative_compensation=False)

Compute the finale integer iGRC as described in Annex F [1].

This method computes the integer and the raw iGRC values for a given population density and size of critical area. The TLOS, target level of safety, can also be set, but the default value is \(10^{-6}\) as described in Annex F [1].

Note

This method converts the population density to ppl/m^2 as needed for the equation. This is because the unit for the input is ppl/km^2, since this is typically how the density is known.

Parameters

• pop_dens (float) – [ppl/km^2] Population density

• CA (float) – [m^2] Size of the critical area.

• TLOS (float, optional) – [fatalities per flight hour] Target level of safety (the default is 1e-6). This value is described in more detail in Annex F [1].

• use_conservative_compensation (bool, optional) – if True, the 0.3 reduction in iGRC value is applied.

Returns

• iGRC (integer) – The intrinsic ground risk class (as an integer). This is the raw value rounded up to nearest integer.

• raw iGRC (float) – The raw iGRC before rounding up.

1(1,2,3,4,5,6,7)

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

AnnexFTables

This class allows for recreating some of the tables in Annex F [1].

class casex.annex_f_tables.AnnexFTables

This class contains methods for generating a variety of tables in Annex F [1].

Parameters

none. –

none.

static ballistic_descent_table()

Compute the ballistic descent table in Annex F [1].

Produces console output to show the ballistic descent table.

Parameters

none. –

Returns

console_output – The output to show.

Return type

list of strings

static iGRC_tables(show_with_obstacles=True, show_ballistic=False, show_with_conservative_compensation=True, show_glide_angle=False, show_additional_pop_density=False)

Compute the iGRC tables Annex F [1].

Produces console output to show the iGRC tables in a variety of forms. The default setup will produce the nominal iGRC table.

Parameters

• show_with_obstacles (bool, optional) – If true, show iGRC values when obstacles are present. Default True.

• show_ballistic (bool, optional) – If true, show iGRC values for ballistic descent (typically rotorcraft). This superseeds show_glide_angle. Default False.

• show_with_conservative_compensation (bool, optional) – If true, subtract 0.3 from all values. The concept is explained in Annex F. Default True.

• show_glide_angle (bool, optional) – If true, show for 10 degree impact angle instead of 35 degree. Default False.

• show_additional_pop_density (bool, optional) – If true, show additional population density rows. Default False.

Returns

console_output – The output to show.

Return type

list of strings

static iGRC_tradeoff_tables(tradeoff_type, show_integer_iGRC, show_relative)

Compute the iGRC tradeoff tables Annex F [1].

Produces console output to show the iGRC tradeoff tables.

The tables can be shown with integer iGRC values (rounded up to nearest integer) and with one decimal. The latter makes it easier to determine how close the iGRC value is to the lower integer.

The tables can also be shown relative to the nominal table. Here positive values are iGRC values above the nominal and negative below nominal.

Parameters

• tradeoff_type (integer) – If 0, produces the nominal iGRC table. For n = 1 through 6 produces the corresponding Tn table, as described in Annex F. For any other value, an error text is produced.

• show_integer_iGRC (bool) – If True, show the iGRC values rounded up to nearest integer. Otherwise, show with one decimal. This only applies to actual iGRC values. Relative values are always shown with one decimal.

• show_relative (bool) – If True, the table shows the difference between nominal iGRC values and values for the tradeoff table. The numbers are always shown with one decimal.

Returns

console_output – The output to show.

Return type

list of strings

static scenario_computation_table(scenario)

Compute table of intermediate values for the three descent scenarios in Annex F [1].

Produces console output.

Parameters

scenario (integer) – Number of scenario, must be either 1, 2, or 3.

Returns

console_output – The output to show.

Return type

list of strings

1(1,2,3,4,5,6)

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

BallisticDescent2ndOrderDragApproximation

Class supports computation of a ballistic descent under the influence of gravity and drag.

class casex.ballistic_descent_models.BallisticDescent2ndOrderDragApproximation

This class allows for computation of ballistic descent using a variation of the standard second order drag model. In this implementation it has been modified to support fast calculations. Details on this can be found in [1].

aircraft

Class holding information about the aircraft.

Type

:class:’AircraftSpecs’

distance1

[m] Horizontal distance travelled between start time of descent and time time_top.

Type

float

distance2

[m] Horizontal distance travelled between time time_top and time time_cross.

Type

float

distance3

[m] Horizontal distance travelled between time time_cross and time time_impact.

Type

float

velocity_x

[m/s] Horizontal part of the impact velocity.

Type

float

velocity_y

[m/s] Vertical part of the impact velocity.

Type

float

time_top

[s] Time from start of descent to reaching the top of the descent curve. This is 0 if the descent starts with a descent.

Type

float

time_cross

[s] Time from start of descent to reaching the point, where the vertical velocity exceeds the horizontal velocity.

Type

float

time_impact

[s] Time from start of descent to impact.

Type

float

c

This represents the multiplication of frontal area, air density, and drag coefficient, which often appears together in the computations. The c attribute is just a placeholder for this multiplication.

Type

float

gamma

A placeholder for \(\sqrt{m \cdot g / c}\), where \(m\) is the mass of the aircraft, and \(g\) is the gravitational constant.

Type

float

compute_ballistic_distance(altitude, initial_velocity_x, initial_velocity_y)

Compute the distance, time, angle, and velocity of a ballistic descent impact.

One of the following parameters can be an numpy.array

• altitude

• initial_velocity_x

• initial_velocity_y

• aircraft.ballistic_drag_coefficient

• aircraft.ballistic_frontal_area

and the output will also be an numpy.array for the outputs that depends on the parameters given as an numpy.array. Note that aircraft refers to the variable set with the methods set_aircraft.

The following requirements apply to the inputs

• initial_velocity_x must be positive.

• initial_velocity_x must be larger than initial_velocity_y.

• absolute value of initial_velocity_y must be less than the terminal velocity.

Parameters

• altitude (float) – [m] Altitude of aircraft at time of event.

• initial_velocity_x (float) – [m/s] Horizontal velocity as time of event.

• initial_velocity_y (float) – [m/s] Vertical velocity as time of event.

Returns

• distance_impact (float) – [m] Horizontal distance to impact point relative to event point.

• velocity_impact (float) – [m/s] Impact velocity.

• angle_impact (float) – [deg] Impact angle (relative to horizontal).

• time_impact (float) – [s] Time from event to impact.

Raises

• NegativeHorizontalVelocityError – If the initial horizontal velocity is negative.

• HorizontalSmallerThanVerticalVelocityError – If the initial horizontal velocity is smaller than the initial vertical velocity.

set_aircraft(aircraft)

Set the aircraft.

Parameters

aircraft (AircraftSpecs) – Class holding information about the aircraft.

Returns

Return type

None

1

A. la Cour-Harbo and H. Schioler. How obstacles may reduce the impact of a crashing unmanned aircraft. Preprint, pages, 2021.

Constants

Constant used throughout the package.

casex.constants.GRAVITY

[kg/m/s^2] Gravitational constant: 9.81

casex.constants.AIR_DENSITY

[kg/m^3] Density of air at sea level: 1.225

Conversion

Provide functionality for converting between imperial and metric system.

class casex.conversion.Conversion

Provide functionality for converting between imperial and metric system.

Note that all methods are static, and can be called without instantiating Conversion, i.e.

length = casex.Conversion.ft_to_m(32)

static ft_to_m(length)

Converts feet to meter. Simply a multiplication with 0.3048 ft/m.

Parameters

length (float) – Length in [ft].

Returns

length – Length in [m].

Return type

float

static ftlb_to_J(energy)

Converts foot-pound to joule. Simply a multiplication with 1.355818 J/ft-lbs.

Parameters

energy (float) – Energy in [ft-lbs].

Returns

energy – Kinetic energy in [J].

Return type

float

static kg_to_lbs(mass)

Converts kg to pounds. Simply a division with 0.45359237 lbs/kg.

Parameters

mass (float) – Mass in [lbs].

Returns

mass – Mass in [kg].

Return type

float

static lbs_to_kg(mass)

Converts pounds to kg. Simply a multiplication with 0.45359237 lbs/kg.

Parameters

mass (float) – Mass in [kg].

Returns

mass – Mass in [lbs].

Return type

float

static m_to_ft(length)

Converts meter to feet. Simply a division with 0.3048 ft/m.

Parameters

length (float) – Length in [m].

Returns

length – Length in [ft].

Return type

float

CriticalAreaModels

This class provide methods for computing critical area for a variety of models.

class casex.critical_area_models.CriticalAreaModels(buffer=0.3, height=1.8)

This class is used for computing the critical area using the JARUS model for a crash when the basic parameters for the aircraft is known. The math behind the model is described in Appendix B in Annex F [1].

The main method in this class is critical_area, which computes the size of the critical area given a number of input arguments relating to the aircraft.

There are two attributes in this class for describing a standard person, namely buffer and height.

Parameters

• buffer (float, optional) – [m] Radius of a standard person as seen from above (the default is 0.3 m).

• height (float, optional) – [m] The altitude above the ground at which the aircraft can first impact a person (the default is 1.8 m).

buffer

[m] Radius of a standard person as seen from above (the default is 0.3 m).

Type

float, optional

height

[m] The altitude above the ground at which the aircraft can first impact a person (the default is 1.8 m).

Type

float, optional

static check_glide_angle(glide_angle)

Checks the glide angle.

Issues a warning if the glide angle is out of range, or close to zero (which produces unrealistic results in the CA model). It also flips the glide angle if it is between 90 and 180 degrees.

Parameters

glide_angle (float) – [deg] The glide angle to be checked.

Returns

glide_angle – [deg] The glide angle, which is either the same as the input, or flipped if needed.

Return type

float

critical_area(aircraft, impact_speed, impact_angle, critical_areas_overlap=0, lethal_kinetic_energy=- 1, use_obstacle_reduction=True)

Computes the lethal area as modeled by different models.

This function supports one of the following input parameters to be a vector, which will give a vector of the same size as output:

• impact_speed

• impact_angle

• critical_area_overlap

• aircraft.width

• aircraft.fuel_quantity

This vector is given as numpy.array, and only one of the parameters can be a vector for each call. The return values are then also numpy.array IF the input parameter that is a numpy.array is used in the computation.

Parameters

• aircraft (casex.AircraftSpecs) – Class with information about the aircraft.

• impact_speed (float) – [m/s] Impact speed of aircraft (this is speed along the velocity vector).

• impact_angle (float) – [deg] Impact angle relative to ground (90 is vertical, straight down). Note that when using the JARUS model and the width of the aircraft is <= 1 m, the impact angle used is minimum 35 degrees, regardless of input. See Appendix A.5 in Annex F [1] for details.

• critical_areas_overlap (float, optional) – [0 to 1] Fraction of overlap between lethal area from glide/slide and from explosion/deflagration. Default is 0.

• lethal_kinetic_energy (float, optional) – The lethal kinetic energy threshold in J. If not specified (or set to -1), the standard approach as described in Annex F Appendix A [1] is used. If set to a positive value, this value is used independently of all other parameters. If set to a negative value other than -1, this value is used following the Annex F standard approach, but with absolute of the given value.

• use_obstacle_reduction (bool, optional) – If set to true, the Annex F obstacle reduction is applied as described in the Annex. If set to false, no obstacle reduction is applied. Default is true.

Returns

• critical area (float) – [\(\mathrm{m}^2\)] Size of the critical area for the selected model.

• estimated glide area (float) – [\(\mathrm{m}^2\)] The glide and slide areas are estimated as the relation between the glide and slide distances multiplied by the glide+slide area.

• estimated slide area (float) – [\(\mathrm{m}^2\)] The estimated slide area.

• critical area inert (float) – [\(\mathrm{m}^2\)] The inert part of the critical area.

• deflagration area (float) – [\(\mathrm{m}^2\)] The deflagration area as given by the deflagration model. Always 0 for non-JARUS models.

• glide distance (float) – [\(\mathrm{m}\)] The glide distance.

• slide distance non lethal (float) – [\(\mathrm{m}\)] The slide distance until non-lethal. Always 0 for non-JARUS models.

• velocity min kill (float) – [\(\mathrm{m/s}\)] The speed at which the slide is no longer lethal. Always 0 for non-JARUS and NAWCAD models.

• t safe (float) – [\(\mathrm{s}\)] The time between impact and when slide speed is non-lethal. Always 0 for non-JARUS and NAWCAD models.

Raises

• InvalidAircraftError – If the aircraft is not of type AircraftSpecs.

• OnlyOneVetorInputError – There is more than one vector input, which is not supported.

static glide_angle_from_glide_ratio(glide_ratio)

Compute glide angle from glide ratio.

Parameters

glide_ratio (float) – [-] The ratio between vertical and horizontal speed during glide.

Returns

glide_angle – [deg] The glide angle for the given glide ratio.

Return type

float

glide_distance(glide_angle)

Compute glide distance based on glide angle.

Glide distance is the distance an aircraft will glide through the air for a given glide angel from altitude height until it impacts the ground. Thus, the glide starts at altitude Height and continues until the aircraft impacts the ground.

Parameters

glide_angle (float) – [deg] The angle of the aircraft relative to the ground as is impacts the ground. Must be between 1 and 179 degree. Values above 90 degrees are used as ‘180 - GlideAngle’.

Returns

distance – [m] The glide distance.

Return type

float

static horizontal_speed_from_angle(impact_angle, impact_speed)

Compute horizontal speed component for a given impact angle and impact speed.

Parameters

• impact_angle (float) – [deg] Impact angle of the aircraft.

• impact_speed (float) – [m/s] Impact speed of the aircraft (speed in the direction of travel).

Returns

horizontal_speed – [m/s] The horizontal compotent of the impact speed.

Return type

float

static horizontal_speed_from_ratio(glide_ratio, impact_speed)

Compute horizontal speed from glide ratio.

Parameters

• glide_ratio (float) – [-] The ratio between vertical and horizontal speed during glide.

• impact_speed (float) – [m/s] The impact speed of the aircraft (in the direction of travel).

Returns

horizontal_speed – [m/s] The horizontal compotent of the impact speed.

Return type

float

static slide_distance_friction(velocity, friction_coefficient)

Computes slide distance based on initial velocity and friction.

Sliding distance computed based on the assumption

\[F = -f \cdot w,\]

where \(F\) is the frictional force, \(f\) the frictional coefficient, and \(w\) the body weight. The slide distance is the length of the slide between impact and the body coming to rest.

This is a standard assumption found in most sources that includes friction. See for instance [2].

Parameters

• velocity (float) – [m/s] Horizontal component of the impact velocity.

• friction_coefficient (float) – [-] Friction coefficient, typically between 0.4 and 0.7.

Returns

distance – [m] Distance from impact to rest.

Return type

float

static speed_from_kinetic_energy(KE, mass)

Compute speed from kinetic energy.

Parameters

• KE (float) – [J] Kinetic energy of the aircraft.

• mass (float) – [kg] Mass of the aircraft.

Returns

speed – [m/s] The speed associated with the given kinetic energy and mass.

Return type

float

static vertical_speed_from_angle(impact_angle, impact_speed)

Compute vertical speed from impact angle.

Parameters

• impact_angle (float) – [deg] Impact angle of the aircraft.

• impact_speed (float) – [m/s] Impact speed of the aircraft (speed in the direction of travel).

Returns

vertical_speed – [m/s] The vertical compotent of the impact speed.

Return type

float

1(1,2,3)

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

2

John A. Ball, Michael Knott, and David Burke. Crash Lethality Model - Report no. NAWCADPAX/TR-2012/196. Technical Report, Naval Air Warfare Center Aircraft Division, 2012.

Enums

class casex.enums.AircraftMaterial(value)

Enum of the aircraft materials.

ALUMINUM = 3

CARBONFIBER = 2

GLASSFIBER = 1

RUBBER = 7

STEEL = 4

STYROFOAM = 6

WOOD = 5

class casex.enums.AircraftType(value)

Enum for aircraft type.

FIXED_WING = 2

GENERIC = 1

LIGHTER_THAN_AIR = 5

MULTI_ROTOR = 4

ROTORY_WING = 3

class casex.enums.FuelType(value)

Enum for fuel type.

AVGAS = 5

DIESEL = 2

GASOLINE = 1

JETA1 = 3

LIFE = 7

LION = 8

LIQUID_BUTANE = 10

LIQUID_HYDROGEN = 9

METHANOL = 6

PETROL = 1

class casex.enums.GroundMaterial(value)

Enum of the ground materials.

ASPHALT = 2

CONCRETE = 1

GRASS = 3

SAND = 4

SOIL = 5

class casex.enums.Wrapping(value)

Enum for wrapping type for aircraft heading.

NONE = 2

PI2PI = 1

Custom exceptions

exception casex.exceptions.HorizontalSmallerThanVerticalVelocityError

The horizontal velocity is smaller than the vertical velocity.

exception casex.exceptions.InvalidAircraftError

The aircraft is not of type AircraftSpecs.

exception casex.exceptions.NegativeHorizontalVelocityError

The horizontal velocity is negative.

exception casex.exceptions.OnlyOneVetorInputError

There are more than one vector input, which is not supported.

ExplosionModels

Class for computing the explosion and deflagration models.

class casex.explosion_models.ExplosionModels

This class implements a model for the lethal area for an explosion, the lethal thermal area for explosion as well as deflagration, and the size of a deflagration fireball (which is also considered lethal).

The main sources for the models are [1], [2], [3], and a brief review of explosion and deflagration is given in Annex F [4].

The models are all based on TNT equivalent mass, since this is how the literature does it. This means that for any of the models it is necessary to convert the fuel amount to a given TNT mass which has the same energy density.

static TNT_equivalent_mass(type_of_fuel, fuel_quantity)

Compute the equivalent mass of TNT for a number of fuels.

Compute the amount of TNT in kg equivalent in energy density to a given type of fuel. This value is used for the explosion and deflagration computations.

Available types of fuels are:

Fuel type	Energy density	Mass density
TNT	4.184 MJ/kg	–
Gasoline	46.4 MJ/L	0.75 kg/L
Diesel	45.6 MJ/L	0.83 kg/L
Jet A1	43 MJ/L	0.80 kg/L
AvGas	44.7 MJ/L	0.69 kg/L
Methanol	19 MJ/L	0.79 kg/L
LiFe batteries	1.8 MJ/kg	–
LiOn batteries	0.8 MJ/kg	–
Liquid hydrogen	142 MJ/kg	–
Liquid butane	27.8 MJ/kg	–

Note

The above values are drawn from a variety of sources of unconfirmed reliability. As such, these are the authors’ best guesses at appropriate values.

Parameters

• type_of_fuel (enums.FuelType) – The type of fuel.

• fuel_quantity (float) – [L] The amount of fuel.

Returns

TNT_mass – [kg] TNT equivalent mass.

Return type

float

static fireball_area(TNT_mass)

Compute lethal area for fireball (deflagration).

Compute the size of a fireball for a given amount of propellant. The model is found in [2] page 62.

Note

See note for lethal_area_explosion.

Parameters

TNT_mass (float) – [kg] The equivalent TNT mass for the propellant.

Returns

area – [m^2] Size of area of fireball.

Return type

float

static lethal_area_explosion(TNT_mass, K=7.14)

Compute lethal area for explosion.

Determines the lethal area for an explosion. The model used is

\[A_{explosion} = \pi D^2\]

where \(D\) [m] is the distance from the explosion, where it is no longer lethal. It is given as

\[D = K W^{1/3}\]

where \(K\) [m/kg \(^{1/3}\)] is a scaling factor for the acceptable risk and W [kg] is the TNT equivalent mass [1], [2], [3].

A recommended value for the scaling factor for 3.5 psi overpressure is 18 ft/lb \(^{1/3}\) (also called K18, for unprotected persons), which is equal to

\[K = 7.14~\mathrm{m/kg}^{1/3}\]

in SI units. For details on this value, see [1].

Note

The area given by the used model assumes a near-perfect combustion of the fuel, which typically requires a close to ideal mixing of fuel and oxidizer. As this normally do not happen during an aircraft crash, the model tends to be rather conservative in the estimate of the lethal area. For more detail on this, please consult Annex F [4].

Warning

The model for lethality caused by explosions has not been deeply investigated, and the model is provided as is, with references to it origin. Ultimate responsibility for determining the lethal area at an explosion rests with the user.

Parameters

• TNT_mass (float) – [kg] Equivalent TNT mass to the explosive material. Use TNT_equivalent_mass to determine this value.

• K (float, optional) – [kg/m^(1/3)] Scaling factor for acceptable risk (the default is 7.14).

Returns

area – [m^2] Size of the lethal area, which is in the form of a disc.

Return type

float

static lethal_area_thermal(TNT_mass, p_lethal)

Compute lethal area for thermal radiation (deflagration).

Computes the lethal area for deflagration based on the thermal radiation. Model is taken from [2] page 61.

Note

See note for lethal_area_explosion.

Parameters

• TNT_mass (float) – [kg] The equivalent TNT mass for the propellant.

• p_lethal (float) – [-] A number between 0 and 1 indicating the target lethality.

Returns

area – [m^2] Size of lethal area.

Return type

float

1(1,2,3)

Department of Defense. DoD Ammunition and Explosives Safety Standards: General Explosives Safety Information and Requirements (DOD 6055.09-M, Volume 1). Technical Report, 2012.

2(1,2,3,4)

John A. Ball, Michael Knott, and David Burke. Crash Lethality Model - Report no. NAWCADPAX/TR-2012/196. Technical Report, Naval Air Warfare Center Aircraft Division, 2012.

3(1,2)

Meredith J. Harwick, John Hall, John W. Tatom, and Robert G. Baker. Approved Methods and Algorithms for DoD Risk-Based Explosives Siting. Technical Report, DDESB, jul 2009.

4(1,2)

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

Figures

Recreates some of the figures in Annex F [1].

class casex.figures.Figures

This class provides some static methods that will were used to make some of the figures in Annex F [1]. The methods take a number of boolean variables to control the content of the figures. The specific choice of parameters for each figure is shown in example 8.

static figure_angle_vs_speed(show_matrix=False, show_contours=False, save_fig=False)

Recreates the figure in Annex F relating impact angles and impact speed for the different size classes.

This method also outputs a variety of computations for the ballistic descent as listed in Annex F Appendix A [1].

Parameters

• show_matrix (bool, optional) – Set to False to get plots in Annex F (the default is False).

• show_contours (bool, optional) – Set to True to see the CA matrix (the default is False).

• save_fig (bool, optional) – If True save the figure to a PNG (default is False).

Returns

Return type

None

static figure_iGRC_CA_vs_PopDensity(filename, show_with_obstacles=False, show_reduced_CA_axis=True, show_old_quantization=True, show_iGRC_prefix=True, show_additional_grid=False, show_colorbar=False, show_x_wingspan=True, show_x_velocity=True, show_x_CA=False, show_x_CA_above=False, show_title=True, save_fig=False, return_fig=False, show_descriptors=False)

Recreates the figures showing iGRC values and iGRC table in Annex F Section 1 [1].

Parameters

• filename (str) – Name of the output file (only applicable if save_image is True).

• show_with_obstacles (bool, optional) – If True show the calculated iGRC values with a reduction due to obstalces (see Annex F Appendix A) (default is False).

• show_reduced_CA_axis (bool, optional) – If True use a reduced granularity on the CA axis (default is True).

• show_old_quantization (bool, optional) – If True uses the SORA V2.0 quantization (only applicable if ?show_reduced_CA_axis? is True) (default is True).

• show_iGRC_prefix (bool, optional) – If True show the iGRC numbers as iGRC-X instead of just X (default is True).

• show_additional_grid (bool, optional) – If True show additional grid lines. Makes it a bit cluttered, but assists in reading the table (default is False).

• show_descriptors (bool, optional) – If True and if ?show_old_quantization? is True, then the descriptors of population density from SORA V2.0 is added to the second axis.

• show_colorbar (bool, optional) – If True show colorbar instead of numbers for the background ISO plot (default is False).

• show_x_wingspan (bool, optional) – If True show wingspan on the x axis (default is True).

• show_x_velocity (bool, optional) – If True show velocities on the x axis (default is True).

• show_x_CA (bool, optional) – If True show size of CA on the x axis (default is False).

• show_x_CA_above (bool, optional) – If True and CA is shown, show it above the figure instead of below (default is False).

• show_title (bool, optional) – If True show the title (default is True).

• save_fig (bool, optional) – If True save the figure to a PNG with the given filename (default is False).

• return_fig (bool, optional) – If True return the created matplotlib figure instead of showing it (default is False).

Returns

Return type

None

static figure_obstacle_critical_area_reduction(save_fig=False)

Recreates the figure in Annex F Appendix B for the CA reduction used in the iGRC table..

Parameters

save_fig (bool, optional) – If True save the figure to a PNG (default is False).

Returns

Return type

None

1(1,2,3,4)

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

FrictionCoefficients

Friction coefficients between a variety of aircraft materials and ground types.

class casex.friction_coefficient.FrictionCoefficients

This class provides some help to determine an appropriate friction coefficient for a sliding aircraft. The values provided here are guidance only, and not all combinations of aircraft material and ground type is known in this class.

The coefficient takes the following values.

	Ground type
Aircraft material	Concrete	Asphalt	Grass	Sand	Soil
Glass fiber	0.2	n/a	0.15	n/a	n/a
Carbon fiber	n/a	n/a	n/a	n/a	n/a
Aluminum	0.4	n/a	n/a	n/a	n/a
Steel	0.45	n/a	n/a	0.2	0.4
Wood	0.6	n/a	n/a	n/a	n/a
Styrofoam	n/a	n/a	n/a	n/a	n/a
Rubber	0.7	0.9	0.35	0.5	n/a

Note

The above coefficients are drawn from a variety of sources of unconfirmed reliability. As such, these are the authors’ best guesses at appropriate coefficients. The true coefficient depends on the state of the materials (wet, dry, greased, etc.) and the actual material (e.g. there are numerous types of concrete, soil, etc.).

Note

There are combinations of materials for which there is not coefficient available. This is du to lack of information. If anyone has such information, please contact the authors, and addition coefficient can be included in a update of this package.

Warning

The friction coefficients provided are guidance only! Ultimate responsibility for correct choice and use of friction coefficient rests with the user!

get_coefficient(aircraft_material, ground_material)

Provide friction coefficient for aircraft sliding over ground.

Parameters

• aircraft_material (enums.AircraftMaterial) – Type of the aircraft material.

• ground_material (enums.GroundMaterial) – Type of the ground material.

Returns

friction_coefficient – The friction coefficient between the aircraft and ground material. Returns -1 if the coefficient is not available (see table above). Returns -2 if either material is not recognized.

Return type

float

misc

Miscellaneous classes for supporting the package.

class casex.misc.InitialSpeeds(initial_speed_x_mu, initial_speed_x_sigma, initial_speed_y_mu, initial_speed_y_sigma)

Class for holding initial speeds for ballistic descent.

Parameters

• initial_speed_x_mu (float) – [m/s] The mean value of the normal distribution for the initial horizontal speed.

• initial_speed_x_sigma (float) – The standard deviation of the normal distribution for the initial horizontal speed.

• initial_speed_y_mu (float) – [m/s] The mean value of the normal distribution for the initial vertical speed.

• initial_speed_y_sigma (float) – The standard deviation of the normal distribution for the initial vertical speed.

initial_speed_x

[m/s] The initial horizontal speed.

Type

float

initial_speed_y

[m/s] The initial vertical speed.

Type

float

class casex.misc.NormalDistributionParameters(mu=0.0, sigma=1.0, wrapping_type=<Wrapping.NONE: 2>)

Class for provide support for generating and using normal distributions.

Parameters

• mu (float, optional) – Mean of the normal distribution (the default is 0).

• sigma (float, optional) – Standard deviation of the normal distribution (the default is 1).

• wrapping_type (enums.Wrapping, optional) – The wrapping type for mu. When set to EWrapping.PI2PI, mu is wrapped to the interval -\(\pi\) to \(\pi\) (the default is EWrapping.NONE).

input_set

The domain for the sampling (i.e. the input to the distribution, or the x axis values).

Type

float array

output_set

The value set for the sampling (i.e. the output of the distribution, or the y axis values).

Type

float array

mu

Mean of the normal distribution (the default is 0).

Type

float, optional

sigma

Standard deviation of the normal distribution (the default is 1).

Type

float, optional

wrapping_type

The wrapping type for mu. When set to EWrapping.PI2PI, mu is wrapped to the interval -\(\pi\) to \(\pi\) (the default is EWrapping.NONE).

Type

enums.Wrapping, optional

compute_sampling(times_sigma, num_of_samples)

Computes a sampling of the normal distribution.

The normal distribution can be plotted using output set against input set. This method computes the input set as a linear function and the output set as the normal distribution from the input set. Both sets are parameters in the class.

Parameters

• times_sigma (float) – This value is multiplied onto sigma and the results plus/minus is the interval for the sampling.

• num_of_samples (int) – Number of samples in the sampling.

Returns

Return type

None

Obstacles

Support both computation and simulation of the reduction of critical area.

class casex.obstacles.Obstacles(CA_width, CA_length, num_of_obstacles, trial_area_sidelength)

This class has methods for computing the theoretical reduction in the size of the critical area when there are obstacles in the ground area as well as for simulating this reduction.

The theoretical reduction is based on the work [1].

Examples of how to MISSING DOC

num_of_obstacles

Number of obstacles in the simulation.

Type

int

trials_count

Number of trials in the simulation.

Type

int

CA_width

[m] Width of critical area.

Type

float

CA_length

[m] Length of critical area.

Type

float

intersected_obstacles

List of all obstacles that have caused reduction in a CA.

Type

List of Polygon

closest

The point closest to the beginning of a CA for each reduced CA.

Type

Point

CA_cut_off_coords

Coordinates of CAs where they are cut off as a result of impact with an obstacle.

Type

List of Point

obstacles_rtree

Intermediate variables for increasing computations speed.

Type

STRtree

CA_lengths

List of the length of every CA after reduction

Type

Array of floats

total_obstacle_area

[m^2] The total area of all obstacles not considering any overlap (sum of the area of every obstacle).

Type

float

total_coverage

[m^2] The total area covered by obstacles. This means that overlapping areas are only counted once.

Type

float

trial_area_sidelength

[m] Length of each side of the square trial area.

Type

float

obstacles

A list of all obstacles in the simulation.

Type

List of Polygon

CAs

A list of all nominal CAs in the simulation (before potential reduction).

Type

List of Polygon

CAs_reduced

A list of all CAs after potential reduction.

Type

List of Polygon

num_of_empty_CA

The count of how many CAs have become empty in the simulation.

Type

int

num_of_reduced_CA

The count of how many CAs are reduced in the simulation.

Type

int

class DistributionParameters(distribution_type: ‘DistributionType’, loc: float = 0.0, scale: float = 1.0)

class DistributionType(value)

An enumeration.

static Minkowski_difference_convex_polygons(A, B)

Compute the polygon that is the Minkowski difference of two convex polygons A and B.

Compute the Minkowski difference of two convex polygons. This methods is based on the Minkowski sum for convex polygons, and it is therefore required that both input polygons are convex, otherwise the result is not correct.

The result is returned as a MultiPoint type.

Parameters

• A (Polygon) – The one polygon in the Minkowski difference.

• B (Polygon) – The other polygon in the Minkowski difference.

Returns

C – The result of the computation as a Multipoint (collection of points) and not a polygon.

Return type

MultiPoint

static Minkowski_sum_convex_polygons(A, B)

Compute the polygon that is the Minkowski sum of two convex polygons A and B.

This methods is based on convex hull for computing the Minkowski sum, and it is therefore required that both input polygons are convex, otherwise the result is not correct.

The result is returned as a MultiPoint type.

Parameters

• A (Polygon) – The one polygon in the Minkowski sum.

• B (Polygon) – The other polygon in the Minkowski sum.

Returns

C – The result of the computation as a Multipoint (collection of points) and not a polygon.

Return type

MultiPoint

static Minkowski_sum_convex_polygons_area(w, x, a, b, theta1, theta2)

Compute the area of the Minkowski sum of two rectangles polygons.

This is a fast method for computing the Minkowski sum of two polygons that are both rectangles.

For details on how this is done, see [1].

Parameters

• w (float) – Width of rectangle 1.

• x (float) – Length of rectangle 1.

• a (float) – Width of rectangle 2.

• b (float) – Length of rectangle 2.

• theta1 (float) – [deg] Angle of rectangle 1.

• theta2 (float) – [deg] Angle of rectangle 2.

Returns

area – Area of the Minkowski sum of the two rectangles.

Return type

float

class ObstaclesSizeProperties(width_mu: float, width_sigma: float, length_mu: float, length_sigma: float)

cdf(x, obstacle_size_resolution=10, CA_orientation_resolution=10, obstacle_orientation_resolution=10, probability_threshold=0.0044318484119380075, show_progress=False)

Compute the CDF for the length of the critical area when rectangular obstacles are present.

This is the CDF for the length of the critical area when there are a given obstacle density of rectangular obstacles with dimension given by normal distributions. To draw the full CDF, a typical input for x is an array ranging in value from 0 to the nominal length of the CA. Since this is usually a rather smooth curve, it can be approximated well by relatively few x values (typically 10 or 15).

The resolution should not be lower than

Note that this function relies on a previous call to generate_rectangular_obstacles_normal_distributed().

For a more detailed explanation of the CDF, see [1].

Parameters

• x (float array) – [m] The length of the critical area for which the CDF is computed. This can be a scalar or an array.

• obstacle_size_resolution (int (default is 10)) – The number of points for the discretization of the two integrals for width and length of obstacles. A good starting value is 10. For high resolution, 20 is an approprite choice.

• CA_orientation_resolution (int (default is 10)) – The number of points for the discretization of integral for CA orientation. A good starting value is 10. For high resolution, 20 is an approprite choice.

• obstacle_orientation_resoltion (int (default is 10)) – The number of points for the discretization of the integral over the obstacle orientation. The value depends highly on the chosen distribution. If the orientation type is FIXED, this value is not used, as the integral does not need to be evaluated.

• probability_threshold (float (default is PDF for normal distribution evaluated at 3 sigma, approx 0.0044)) – This value is used to speed up computation of the integrals. For any probability below this threshold, the contribution to the integrals is ignored. To include all values, set this to zero. Adjusting this value will affect the sanity checking, so if it becomes too large (i.e., too much of the integral contributions are ignored), a warning will be issued.

• show_progress (bool (default False)) – Write the progress in percent to the prompt for the multple integral computation.

Returns

• p_x (float array) – The CDF value for the given x. This return parameter has the same type as input x.

• expected_value (float) – [m] The expected value of the length.

• beta (float) – The beta values as computed in [1].

• acc_probability_check (float) – A sanity check on the triple integral. This values should be relatively close to 1, especially for high value of pdf_resolution.

• pdf_width (np.array) – The PDF for the obstacle width as used in the integral.

• pdf_length (np.array) – The PDF for the obstacle length as used in the integral.

• pdf_CA_orientation (np.array) – The PDF for the CA orientation as used in the integral.

compute_CA_lengths()

Make a list of the actual lengths of the CAs.

Parameters

None –

Returns

Return type

None

compute_coverage(show_progress=False)

Determine total obstacle coverage.

Returns

Return type

None

compute_reduced_CAs(show_progress=False)

Compute the reduction for each CA

Any CA that intersects with an obstacles is reduced such as to no longer intersect with any obstacles. This method runs through all CAs and all obstacles and produces a list of CAs that.

Returns

Return type

None

generate_CAs(trials_count)

Generate a number of critical areas for simulation.

A number of critical areas are generated with 2D uniformly distributed location and uniformly distributed orientation between 0 and 360 degrees. They are have the width and length as set at initialization of the class.

Parameters

trials_count (int) – Number of trials to perform.

Returns

Return type

None

generate_rectangular_obstacles_along_curves(width_mu, width_sigma, length_mu, length_sigma, houses_along_street, rows_of_houses, distance_between_two_houses)

Generate a set of obstacles that follows a curve.

Rectangular obstacles are generated so that they follow a vertical curve and located in parts. This is to simulate houses along a road. The size and density of houses can be adjusted.

An area larger than the trial area is covered with obstacles, but only the ones inside the trial area are preserved. The number of obstacles is set to the number of obstacles preserved.

Parameters

• width_mu (float) – [m] The mean of the normal distribution of the width of the obstacles.

• width_sigma (float) – The standard deviation of the normal distribution of the width of the obstacles.

• length_mu (float) – [m] The mean of the normal distribution of the length of the obstacles.

• length_sigma (float) – The standard deviation of the normal distribution of the length of the obstacles.

• houses_along_street (float) – The number of houses along the road (22 is a good starting number).

• rows_of_houses (float) – The number of rows of houses (12 is a good starting number).

• distance_between_two_houses (float) – [m] The distance between two houses that make up the pair of houses that shares a common border, but are on two different streets (20 m is a good starting number).

Returns

Return type

None

generate_rectangular_obstacles_normal_distributed(width_mu, width_sigma, length_mu, length_sigma)

Generate a set of uniformly distributed rectangular obstacles.

This method generates a number of rectangular obstacles which a co-linear with the axes, and with width and length varying according to normal distributions with the mean and standard deviation as given by input parameters. The position of the obstacles are uniformly distributed in 2D in the trial area.

Parameters

• width_mu (float) – The mean of the normal distribution of the width of the obstacles.

• width_sigma (float) – The standard deviation of the normal distribution of the width of the obstacles.

• length_mu (float) – The mean of the normal distribution of the length of the obstacles.

• length_sigma (float) – The standard deviation of the normal distribution of the length of the obstacles.

Returns

Return type

None

static mirror_polygon_in_origin(polygon)

Compute the mirror of a polygon by negative all corner coordinates.

Parameters

polygon (Polygon) – A polygon to be mirrored.

Returns

m – A polygon that is the mirror of the original polygon

Return type

Polygon

missed_obstacle_CA_intersections()

Identifies missed intersections metween obstacle and reduced CA.

In the simulation, it may happen that a critical areas is reduced, and the part of the critical area that has been reduced away intersects another obstacle that the one causing the reduction. If this overlap is such that the nominal critica area is not completely separated in two (i.e., the obstacle has a corner inside the CA without going all the way through), the overlapping area will not contribute correctly to the simulation results.

This is a rare event, and is therefore ignore. To visualize these problematic overlaps, this function provide a list of the affect obstacles and CAs, plus compute the area missed in the simulation.

Returns

• intersection area (MISSING DOC) – MISSING DOC

• problematic obstacles (MISSING DOC) – MISSING DOC

• problematic CAs (MISSING DOC) – MISSING DOC

obstacle_density()

Return the density of obstacles.

static set_limits(ax, x0, xN, y0, yN, step=1)

Set the axes limits for a axis in a plot.

Parameters

• ax (Handle MISSING DOC) – MISSING DOC

• x0 (float) – MISSING DOC

• xN (float) – MISSING DOC

• y0 (float) – MISSING DOC

• yN (float) – MISSING DOC

• step (int, optional) – MISSING DOC

Returns

Return type

None

show_CDF(ax, show_CA_as_size=True, line_label=None, line_color='blue', line_width=3)

Plots the CDF

Plots a histogram of simulation results to the provided axis.

Parameters

• ax (handle) – Handle to the plotting axis.

• show_CA_as_size (bool (default is True)) – If True, show the first axis as size of CA. If False, show the first axis as CA length (i.e., CA size divided by width).

• line_label (text (default is None)) – If None, a label text is generated. If not None, the provided text is used for label text.

• line_color (color (default is 'blue')) – Color of CDF graph (passed directly to plot).

• line_width (float (default 3)) – Width of CDF graph (passed directly to plot).

Returns

Return type

None

show_simulation(ax, problematic_obstacles=None, problematic_CAs=None, show_CAs=True, show_CA_first_point=False, show_CAs_reduced=True, show_obstacles=True, show_obstacles_intersected=True, show_debug_points=False)

Visualize simulation with obstacles and critical areas.

This functions makes it easy to visualize the result of simulations. It will show the simulated square area with obstacles (possible showing which are impacted and which are not), critical areas (both which are “full” and which have been reduced).

It can also show a list of problematic obstacles and critical areas. For more details on what that is and how to detect them, see the help for sanity_check().

Parameters

• ax (axis to plot in) – Handle to the plot axis

• problematic_obstacles (list of polygon (default None)) – Show the provided problematic obstacles (in yellow).

• problematic_CAs (list of polygon (default None)) – Show the provided problematic CAs (in yellow).

• show_CAs (bool (default True)) – Show the nominal critical areas.

• show_CA_first_point (bool (default False)) – Show a point at the starting end of each CA.

• show_CAs_reduced (bool (default True)) – Show the reduced critical areas.

• show_obstacles (bool (default True)) – Show all obstacles (in green).

• show_obstacles_intersected (bool (default True)) – Show the obstacles intersected by critical areas (in orange instead of green).

• show_debug_points (bool (default False)) – Show debug points on the CAs to help verify that intersections and polygonal reduction is corrected. Typically only used for debugging.

Returns

Return type

None

singleton_objects_CDF(x)

CDF for singleton objects.

This implements equation (15) in the paper. MISSING DOCS

Note that this function requires objects to have been specified first.

Parameters

x (float, scalar or np.array) – The length of the CA, typically starting from 0 to the full length (nominel CA size divded by width)

Returns

Return type

CDF with same number of samples as input x (but always as np.array)

1(1,2,3,4)

A. la Cour-Harbo and H. Schioler. How obstacles may reduce the impact of a crashing unmanned aircraft. Preprint, pages, 2021.

Indices and tables

• Index

• Module Index

• Search Page

1

John A. Ball, Michael Knott, and David Burke. Crash Lethality Model - Report no. NAWCADPAX/TR-2012/196. Technical Report, Naval Air Warfare Center Aircraft Division, 2012.

2

Meredith J. Harwick, John Hall, John W. Tatom, and Robert G. Baker. Approved Methods and Algorithms for DoD Risk-Based Explosives Siting. Technical Report, DDESB, jul 2009.

3

JARUS. JARUS guidelines on SORA - Annex F - Theoretical Basis for Ground Risk Classification. Technical Report, 2022. JAR-WG6-QM.01.

4

A. la Cour-Harbo and H. Schioler. How obstacles may reduce the impact of a crashing unmanned aircraft. Preprint, pages, 2021.

5

Anders la Cour-Harbo. Ground impact probability distribution for small unmanned aircraft in ballistic descent. In 2020 International Conference on Unmanned Aircraft Systems (ICUAS), 1442–1451. Aalborg University, IEEE, sep 2020. doi:10.1109/ICUAS48674.2020.9213990.

6

Robert M. Montgomery and James A. Ward. Casualty Areas from Impacting Inert Debris for People in the Open - RTI Report No. RTI/5180/60-31F. Technical Report, Research Triangle Institute, 1995.

7

Department of Defense. DoD Ammunition and Explosives Safety Standards: General Explosives Safety Information and Requirements (DOD 6055.09-M, Volume 1). Technical Report, 2012.

8

Federal Aviation Administration. Flight Safety Analysis Handbook Federal Aviation Administration. Technical Report, Federal Aviation Administration, sep 2011.

9

Range Commanders Council. Range safety criteria for unmanned air vehicles - rationale and methodology supplement. Supplement to document 323-99. Technical Report, Range Commanders Council, 2001.