1. Introduction
For the last few decades, geographic information systems (GIS) techniques have developed and been used in archaeology (Conolly & Lake 2006). A good example of the use of GIS in spatial analyses is viewshed analysis (Tandy 1967), which is often applied in archaeological studies to improve understanding of ancient landscapes (Kay & Sly 2001), such as intervisibility between long-gone settlements (Wheatley 1995) and visibility of cultural objects (Ogburn 2006). Key terrain, Observation and fields of fire, Concealment and cover, Obstacles and Avenues of approach (KOCOA) is an analytical method employed by war commanders and leaders (USMC 1988) to assess the battlefield features of terrain and maximise gained advantages (USMC 2018). In conflict archaeology, researchers apply the KOCOA method to identify and compare terrain features with accounts of historical events (Neumann 2021). KOCOA is often applied to the battlefields of the First and Second World Wars (Spennemann 2020) as well as less recent sites (Brown et al. 2017; Brown 2021). Sometimes, KOCOA analysis aids in identifying features such as fortifications and weapon positions that would be difficult to detect by other methods (Holas 2022; Sivilich & Sivilich 2015). Currently, fields of fire are still identified with known weapon ranges, but this approach becomes insufficient in terrain with high elevation differences. Additionally, areas within range can be obstructed by mountainous terrain or anthropogenic obstacles, so although the results of viewshed analyses can identify visible areas, fields of fire may cover larger areas, especially for weapons with high elevation angles and high trajectories to which invisible areas are accessible.
Therefore, in this paper, we suggest using a ‘throwshed’ analysis with the KOCOA method to identify fields of fire and cover in conflict archaeology. Compared to viewshed, throwshed does not refer to visible areas of the terrain but to areas that are accessible by a projectile thrown or shot from any given place in accordance with the physical laws of external ballistics. Similarly to viewshed tools that process digital elevation models (DEMs) with a few parameters to create raster results, we introduce a tool for GIS users that is freely available for archaeological purposes. The tool not only augments the KOCOA methodology but also enriches GIS techniques with a new spatial analysis using DEMs (Mokarram & Hojati 2017; Rocha et al. 2023). The creation of such spatial analyses would not be possible without the DEMs that have arisen in recent decades as a result of progress in new approaches to topographical mapping and improvements to their resolution, accuracy, and accessibility both globally (Hawker et al. 2022; NASA 2015) and regionally (Leitmannová et al. 2022; Swisstopo 2022). Favourable performance of these DEMs is often proved by their comparison to elevations with higher accuracy (at least three times better) (Marsh, Harder & Pomeroy 2023; Maune 2007) or to hydrological data (Archer et al. 2018; Boulton & Stokes 2018).
Throwshed is a relatively new analysis within GIS. Previous works have used ‘fireshed’ analysis (Lacey 2003). However, fireshed addresses the issue of fields of fire and cover with a different method that involves significant approximations, as acknowledged by Lacey (2003) and others (Athanson 2010). A study conducted by Silliman and Batt (2015), although employing exterior ballistics, approached the issue in a distributional manner. An approach using 3D domes is also of interest (Spennemann 2020). The current paper presents an advanced methodology that was partially inspired and developed anew. The aim of this work is to demonstrate throwshed’s contribution to the KOCOA method due to its innovative approach to solving the issues of the standard approach. The current study also aims to demonstrate the functionalities of the tool with diverse simulated results with modified settings and secure its credibility by validating simulated results from current, historical, and archaeological spatial data.
2. Methods
The tool was constructed with several inspirations from existing concepts of similar spatial analyses and physical models. The analysis can be carried out via a Python script with associated modules defining the model that comprise the whole tool.
The overall method is described by a flowchart in Figure 1 and as follows. Input parameters are used to construct the ballistic trajectories sequentially. The shooting angle rises in increments over a selected range with an optimally set step. The resultant trajectories together form a set in one vertical plane that broadly characterises the projectile’s reach. This is shown in Figure 2. The new method finds a trajectory with exactly the right shooting angle to intersect a specific cell in a raster DEM. Due to the numerical character of the trajectory computation, this cannot be done analytically; instead, the trajectory that intersects each cell of DEM has to be acquired iteratively by searching for a specific shooting angle, constructing the trajectory, and then testing whether it intersects the target cell. Because it would be much more intensive to iterate through the whole shooting angle range, which can be set from –90° to +90°, the trajectory set helps to divide the whole surrounding space into chunks. After finding the chunk that corresponds with the target cell, two delineating trajectories and their known shooting angles are used to find the cell-intersecting trajectory. This iterative process is shown in Figure 2. Once the intersecting trajectory has been found, the model assesses whether it experiences another terrain crossing ahead of the target cell’s intersection. If not, the target cell is considered accessible. Otherwise, it is obstructed. Another situation arises when the cell is out of reach and thus inaccessible. All the scenarios are presented graphically in Figure 3. After all the cells of the DEM have been addressed, a throwshed raster file is created for browsing. These factors are described in greater detail in the following subsections.
Figure 1
Flowchart of throwshed general method.
Figure 2
Graphical description of throwshed general method. The procedures shown are the construction of the trajectory in detail (A), the generation of the trajectory set (B), the selection of the target chunk (C), and the iterative search for the target cell’s intersecting trajectory (D).
Figure 3
Three scenarios of the intersection: target cell intersected and accessible (A), target cell intersected but inaccessible (B), and target cell out of reach (C). Note that the first and last trajectories of the set are also depicted.
2.1 Ballistic trajectory solution
A ballistic trajectory can be constructed by several methods. Analytical methods encompass the parabolic solution (Esparza 1984), Siacci’s model (Siacci 1880) and others (Chudinov, Eltyshev & Barykin 2018; Turkyilmazoglu 2016; Weinacht, Cooper & Newell 2005). Although these may not be computationally demanding, they carry several drawbacks. These include ignoring impactful air resistance when the parabolic model considers only gravitational force and limiting input parameters to small ranges, such as elevation angle limited to a maximum of 15° in Siacci’s solution (McCoy 2012). In contrast, numerical methods, although laborious due to the number of computations, are capable of trajectory formation within a fraction of a second with modern processors. In addition, numerical solutions tend to create more authentic trajectories and, therefore, are applied in the physical model.
The numerical procedure divides the whole trajectory into small elements. These are constructed gradually from the shooting point to the end of the trajectory with the differential formulae of Newton’s laws of motion, initial conditions and other parameters, some of which change across the trajectory. Computations result in a list of 3D coordinates of points with a resolution derived by the incrementation approach. Formulae (1) to (3) introduce Newton’s equations of motion for each axis, where the change of independent variable from time to distance is applied (McCoy 2012):
where ρ corresponds to air density in kg/m3, A is an area of perpendicular cross-section of the projectile in m2, CD is the nondimensional aerodynamic drag coefficient, m is the mass in kg, is overall velocity in m/s, Vx, Vy, and Vz are its components in m/s, Wx, Wy, and Wz are wind velocity components in m/s and g is gravitational acceleration in m/s2.
Distance integration is more effective when compared to time integration. In the latter, the first elements of the trajectory may acquire an enormous length, leading to a high level of negligence of force impact, and, at the end of the trajectory, a redundant number of tiny elements can be formed as the projectile rapidly loses velocity. Hence, the distance incrementation is adopted and adapted to slant distances. This is because purely horizontal differentiation of the trajectory can become inaccurate at high shooting angles. Both horizontal and slant distance steps can be applied and set to the desired size; otherwise, these steps are set to the resolution of the DEM by default.
Two numerical solutions are considered in the model. The Euler method is elemental, and many other solutions are inspired by it. Another method is Heun’s solution, a modification of the Euler method that is also called the ‘second-order’ method. The difference between these two lies in the mechanism that addresses the derivatives. Whereas Euler’s method recognises only the derivative at the beginning of the trajectory element, Heun’s method considers derivatives from the beginning and the end of the element. The Euler method creates a tangent with a slope, which is used with the element’s size to find subsequent points on the trajectory. The problem is that due to the concavity and convexity of the real solution, a straight tangent will create a point slightly off the real trajectory, and the more elements in the trajectory or the larger the element size, the greater the deviation of the trajectory points. Nonetheless, this method is applied within the physical model because it suffices for shorter-range weapons. In contrast, Heun’s method computes tangents on both element ends easily by computing one more step ahead with the Euler method. This provides a corrector to the values from the first, predictor, element. Hence, Heun’s method is also identified as the predictor–corrector method and is equivalent to the second-order Runge-Kutta method. Higher-order Runge-Kutta methods iterate through this averaging mechanism (McCoy 2012) and produce even more accurate results. However, these are not used in the model due to the high number of computations required.
The model presented in the method uses a 2-degree of freedom (DOF) point mass trajectory solution because it considers the position of a projectile only within the vertical plane composed of X and Y axes. Due to the approaches used to search for trajectories intersecting specific target DEM cells, the Z axis, thus equation (3), as well as the Coriolis and wind effects, are omitted. Nevertheless, gravity and air resistance forces shape the trajectory in the model. Fortunately, these are found to be the most impactful ones, even at short distances, unlike the Coriolis effect, whose influence increases only at very long distances over kilometres. To close this, up to 6-DOF solutions exist. These also incorporate three rotations of the spin-stabilised projectile around the XYZ axes and create very authentic trajectories. Unfortunately, they require input data (Elsaadany & Wen-jun 2014) for all the physical forces and moments involved and thus are unavailable to users. Thus, the solution presented here compromises accuracy for a limited quantity of inputs.
2.2 Envelope of trajectories
Fixing all input parameters to static constants and permitting only the elevation, or shooting angle, to vary will produce a family of trajectories (Butikov 2015) that create the envelope. This envelope is also identified as a ‘parabola of safety’ for parabolic trajectories because it separates an outer ‘safe’ zone from an inner ‘unsafe’ zone (Pascoal, Castro & Rosa 2011). The envelope can be derived analytically (Baće et al. 2002; Butikov 2003), but incorporating the complex trajectories from the current physical model renders the equations inaccurate. The solution can be acquired by creating a large number of trajectories with different elevation inputs that result in a dense trajectory set from which marginal fragments of individual trajectories contribute to the construction of the whole envelope. Such a curve, if drawn from the shooting point, can then be utilised to determine the area of the DEM that is potentially accessible. Everything outside the envelope is considered inaccessible. This is in accordance with the varying range for any altitude difference between the shooting and target point. Figure 4 shows the computed envelope with a theoretical envelope, which would be created with an infinitesimal step of shooting angle and, therefore, with an infinite number of trajectories without any area omitted. The model computes such an envelope where the omitted area is minimal.
Figure 4
Theoretical (green) and computed (red) envelopes of trajectories with omitted (grey) and accessible (yellow) areas depicted on the trajectory set (black).
2.3 Intersecting the target cells
Although discerning what is within the envelope of the trajectory set and what is outside, in the safety zone, can be satisfactory for basic use cases and relatively flat terrain, natural or anthropogenic obstacles may render areas in the envelope inaccessible due to shadowing effects in defilade positions (Bellamy 1990). To assess the accessibility of a target cell of the DEM, the trajectory that intersects this cell must be constructed and compared to the terrain profile. Thus, a specific elevation angle is sought. This model applies an iterative method of gradual trajectory nearing from surrounding arbitrary trajectories towards the middle of the cell. Algorithms of previous approaches (McCoy 2012) computed vertical differences between the target point and the surrounding trajectories above and below. These differences can aid in adjusting the elevation angle within the range of angles of neighbouring trajectories. This method is effective for flat and direct trajectories with low elevation angle values. However, the ratio of the vertical differences within tall trajectories of high shooting angles would be heavily corrupted by the shape of the trajectories and would complicate the computation process. Hence, the method of searching normals was developed. Perpendicular distances from the closest points of surrounding trajectories towards the target point create an efficacious ratio, and although normals are algorithmically more difficult to compute than simple vertical differences, they reduce the number of iterations in the nearing process. Additionally, this approach is applicable for trajectories of all elevation angles from –90° to 90°.
The drawback, again due to the numerical nature of the model, is that a trajectory intersecting the very central point of the cell is seldom found, because finding such a perfect trajectory requires extensive computation. Thus, trajectories that are as close as half of the cell size from its central point are considered intersecting. This value is considered an indirect indicator of the accuracy of this and a few more algorithms within the overall model.
Assessing the trajectory intersecting the exact centre of the cell would entail calculating the accessibility of the whole cell area. Instead, only the half of the intersected cell is determined, whether in front of the central point or behind it. This information is valuable for cells located on the edges of the obstacle shadows, because they are only partially accessible. Among well-known rasterisation methods (Biagi & Negretti 2001), the hierarchy method is applicable if the trajectories intersecting both sides of the cell are allocated. A statistical sample of two trajectories is not ideal, but this is employed to cope with the situation and numerical solution in compromise with the duration of computations. If at least one side of the cell is proven to be accessible, the whole cell is deemed accessible. The disadvantage is that the cell might be assessed as inaccessible even though, in reality, a tiny fraction of it is reachable. This can occur if the trajectory intersects the inaccessible part of the cell’s half. In summary, it is difficult to obtain statistically sufficient outcomes and to subdue the procedure’s numerical nature.
2.4 Trajectory assessment
To assess whether the target cell is accessible by a projectile, its trajectory is compared to the terrain profile from the shooting point to the target cell. If the trajectory intersects terrain sooner than the target cell, this cell is considered inaccessible; otherwise, it is considered to be accessible. In a viewshed analysis, the relation of line of sight (LOS) to terrain is defined by comparing the slope of LOS with the slopes to all the intermediate cells in the terrain profile (Petrasova et al. 2015). If any of the intermediate cell-related slopes exceeds the slope of the LOS, the target cell is considered invisible. Due to the character of the trajectories, in the throwshed model the heights of the ballistic trajectory points are compared to the heights of respective points on the DEM. If any intermediate terrain height exceeds the trajectory point’s height, the target cell is considered inaccessible.
Two trajectories intersect each cell located within the envelope of trajectories. A direct trajectory has a low elevation angle and is more prone to intersect obstacles, whereas an indirect trajectory has a high elevation angle, as is often applied in artillery. Therefore, two bands of result raster are applied: one for direct firing and the other for indirect. Cells outside the envelope can be hit zero times, and cells right on the envelope can be hit only once (Donnelly 1992). Thus, any cell within the DEM can be hit 0–2 times.
2.5 Input parameters
Several initial conditions are required by the tool to be initiated. Some of them are mandatory, whereas others are required only for specific features that the throwshed tool offers. Parameters are described in Table 1.
Table 1
Input parameters for throwshed tool.
| MANDATORY | |
|---|---|
| Initial height h (m) | Height above the terrain from which the projectile starts its trajectory |
| Elevation θ (°) | Also known as shooting angle; it is set as a range in interval <–90, +90> |
| Gravitational acceleration g (m/s–2) | Causes projectile’s attraction vertically down, usually around 9.80 |
| Initial velocity v0 (m/s) | Or muzzle velocity; the speed at which the projectile leaves the weapon |
| Temperature T (°C) | Important to recalculate the temperature at different relative heights |
| Air density ρ (kg/m3) | One of the factors causing air resistance, a standard value of 1.225 |
| Mass m (kg) | Mass of the projectile |
| Area A (m2) (d (m)) | The area of the projectile perpendicular cross-section; can also be set as the diameter of the rotationally symmetric body |
| Drag coefficient CD | Nondimensional aerodynamic drag coefficient unique for projectile |
| Azimuth α (°) | Set as a range in interval <0, 360> |
| OPTIONAL | |
| Shooter’s eyes height hE (m) | Applicable for viewshed analysis |
| Target height hT (m) | |
| Wall height hW (m) | Applicable for burning artificial wall obstacles into DEM, negative wall height value is considered as a depth of artificial ditch |
| Wall width wW (m) | |
| Oscillation frequency f (Hz) | Frequency of oscillating/rotating projectile (e.g. wobbling arrow) |
| Oscillation distance dO (m) | The distance at which the oscillating effect ceases to be present |
| Peak drag coefficient | Drag coefficient of oscillating projectile at its max |
| Peak area AP (m2) | Area of oscillating projectile at its max |
2.6 Aerodynamic drag coefficient
When using the tool for a specific situation and projectile, most of the mandatory input parameters are relatively simple to obtain and apply. Conversely, acquiring drag coefficient may be difficult because this value is unique for projectiles and needs to be implemented as a function for transonic and supersonic projectiles of velocities close to and above Mach 1. It cannot be neglected because it plays an important role in determining the deceleration due to aerodynamic drag force. The magnitude of deceleration aAD can be described by equation (4).
The coefficient depends on the air’s physical characteristics, such as its density, viscosity, and speed relative to the projectile, as well as the shape and size of the projectile and even its roughness (Shin et al. 2022).
Older techniques of obtaining the coefficient values involved experimental shooting of the projectiles at spark photography ranges, such as the proving ground of the Ballistic Research Laboratory (BRL) (McCoy 2012). Measured data were transformed into drag-to-Mach tables and graphs providing coefficient values for velocities in Mach numbers. BRL reports from such experiments exist and are available online for many projectiles (McCoy 1990a; McCoy 1990b). Even reports of older experiments are accessible (Schmidt 1954). These utilised coefficients in accordance with the older nomenclature (McShane, Kelley & Reno 1953). Transformation to modern CD from older KD is provided by (McCoy 2012) in the following equation (5).
Another method of obtaining the coefficient experimentally would be to use a wind tunnel.
Technological progress allowed the rise of programmes such as Ansys Fluent and its computational fluid dynamics (CFD) mechanism, which enables users to calculate drag and other coefficients owing to numerical airflow simulations around modelled bodies. This task can prove difficult without knowledge of aerodynamics and the behaviour of various kinds of flows. However, many studies have computed the values of the coefficients for specific projectiles, and these are often available for readers. Sometimes, authors even compare the simulated data with the experimental ones, which lends more credibility to their data (Merda & Magier 2017). Moreover, some studies (Ko et al. 2020) provide coefficients calculated by a semi-empirical method, PRODAS, developed by Arrow Tech Associates, and compare them with values obtained by the classical CFD method.
Note that the model currently does not differentiate between total and zero-yaw drag coefficients. The total value of the coefficient is obtained by equation (6) (McCoy 2012).
Here, corresponds to zero-yaw drag coefficient, is a yaw-drag coefficient, and δ2 is an element calculated by pitch and yaw angles. These correspond to the orientation of the projectile’s axis of symmetry towards the direction of flight. Although projectiles without spin are less affected by this approximation, it becomes problematic for spinning projectiles within high elevation trajectories, such as modern artillery shells, where the direction of flight changes drastically while maintaining the projectile’s absolute orientation stable due to the gyroscopic effect. This angled position of the projectile increases the second component of equation (6). Although many studies provide values for total coefficients and respect this issue, these data can depend on specific shooting angles. Finally, besides providing the drag coefficient in the manner of constant or list values by the user, drag tables of a few projectiles are built into the throwshed tool.
2.7 Atmospheric models
When considering indirect artillery trajectories with apexes at much greater heights than the shooting points, atmospheric properties can change rapidly, impacting the magnitude of the air resistance effect. Temperature and air density gradients influence the speed of sound, which inversely affects the projectile’s velocity in Mach numbers. As previously described, the drag coefficient depends strongly on the Mach number, which differs by altitude even though the projectile’s velocity remains the same. Several models of standard atmosphere address this issue for various altitudes. The throwshed physical model employs the older US Army Standard Metro and the newer International Civil Aviation Organization (ICAO) standard atmospheres that are described and applied by McCoy (2012). Because errors in models performance proliferate above the altitude of 12 km, users must approach high-altitude projectile trajectories with caution.
2.8 Digital terrain and surface model
Until now, reference has only been made to DEM. DEM is a term that incorporates both the digital terrain model (DTM) and the digital surface model (DSM). The former addresses only natural terrain, but the latter also includes natural and artificial features above the terrain, such as vegetation and buildings. These may present a barrier to projectile trajectories, so the DSM must be used carefully. A DEM of any projected, metric coordinate system is applicable with any resolution.
Additionally, the throwshed tool can use line vector layers included into the input DEM at specific height or depth and width. This function may help to recreate long-gone anthropogenic features such as walls or ditches. Nonetheless, their presence influenced the fields of fire aspects of the KOCOA analysis to some extent.
2.9 Probable throwshed
Due to the numerical character of the physical model, it is very difficult to determine the accuracy of the result. Addressing uncertainty propagation through all the elements of computed trajectory is costly but can be overcome by introducing the Monte Carlo method to simulate error-affected input parameters to create a ‘probable throwshed’ by summing the outcomes of repetitive calculations. This is employed analogously to the ‘probable viewshed’ (Fisher 1993), but additionally to the simulation of DEM errors, input parameter errors are simulated as well. The result is then expressed as a probability value between 0 and 1 instead of a Boolean value. The Monte Carlo method allows the generation of pseudorandom values, and the incorporation of the central limit theorem accomplishes a normal distribution of the values (Brandt 2014). Using the mean and standard deviation values defined for the input parameters and the permuted congruential generator (PCG64) algorithm (O’Neill 2014), the analysis tool addresses the issue of uncertainty of DEM and geometric and physical parameters.
2.10 Other features
Other functionalities of the tool may help provide more authentic and practical results, some of which are explained here. One of them addresses the oscillation or rotation of the projectiles. This is accomplished by considering the cross-sectional area and drag coefficient of a subsonic projectile at regular and maximum positions. Therefore, these are extreme values. This does not include supersonic spinning projectiles, which follow other laws in 6-DOF trajectory models. Instead, it refers to the rotation of asymmetrical projectiles or the well-known wobble effect of bow arrows related to the archer’s paradox. An arrow released from a bow performs a complex bend that repeats in a decaying manner across a definite distance. Although the transition from regular position to maximal bend is addressed by a goniometric function sine, the transition from maximal to zero bend at a given distance after the wobble effect fades is approached linearly. Besides the slant distance over which the oscillation or rotation effect is cancelled, the frequency is also required. The magnitude of an arrow’s bend depends on its stiffness, dimensions, and energy delivered by the bowstring (Kooi & Sparenberg 1997). This can be computed by mathematical models (Kooi 1998) or derived experimentally (Ando et al. 2014; Park 2011) and then used to model the projectile and obtain the values required with airflow simulation programmes. Although this approach deals with a phenomenon that affects the drag force, it must be treated cautiously due to its limitations and approximative character. Other features incorporate the application of the viewshed analysis that is used to clip the throwshed result directly or inversely, and analogously to viewshed analysis, the option of cumulative throwshed is also possible for multiple shooting points.
2.11 Main premises for analysis application and validation
Considering a few presumptions and properties can help to understand the limitations of the method and its applications. The throwshed model can compute trajectories for any spin and fin projectile whose physical parameters are known in addition to other parameters characterising the 3D position and atmospheric properties. This includes subsonic to supersonic projectiles that do not experience any additional propulsion after being shot. From ballistics theory, only the projectile’s exterior component is incorporated into the solution. Thus, no interior nor terminal ballistics, therefore, any maximum terminal-effect range is modelled by the tool. Instead, a minimum-to-maximum theoretical range is calculated around the shooting place regardless of the energy carried and the potential damage that could be caused. The lack of terminal ballistics cannot be overlooked in specific cases, but this can be mitigated by setting the proper elevation angle range by the user. Areas shaded by obstructions are interpreted within these ranges as holes. The throwshed tool can be used for any area and landscape whose terrain or surface is provided in a DEM raster. The two-band output raster inherits the resolution and coordinate system of the DEM, and cell values in each of the bands reaching 0–n, where n is the number of weapons used. Archaeological studies often focus on one age or event within the historical era. However, this analysis is not bound to a specific epoch and can be used for any event or battlefield for which broad knowledge of the weapons and projectiles involved is available.
Validation sections 4.2 and 4.3 require the muzzle velocities of the cannons, otherwise unknown, to initiate the analysis. Robins’s (1742) book, which was a huge contribution to the ballistics of guns, established the principles of interior smooth bore ballistics, which thanks to Collins’s (2015) work, can be formulated into equation (7).
Here, v0 stands for velocity in m/s, and R is the ratio of hot gas pressure to atmospheric pressure in the cannon’s bore after igniting the gunpowder. Robins (1742) calculated this constant to 1000, but later, Hutton (1805) examined and recalculated this constant to a more precise 1600, which is used in the current paper. Patm stands for atmospheric pressure, which is set to average sea-level 101325 Pa, p is the mass of a gunpowder charge used to fire the cannonball, m is the mass of the cannonball, η is the density of the gunpowder, and L is the full length of the cannon bore. The diameter of the bore, d, can be utilised to compute c, which stands for the length of the gunpowder charge in m. Formula (8) explains this simple calculation.
Drag coefficients in validation sections 4.2 and 4.3 were interpolated from work done by Miller and Bailey (1979), which collects data from 200 years of drag and other sphere aerodynamic characteristics measurements, mainly by Bashforth (1870). As in most historical scientific experiments, these were conducted with spheres limited in Reynolds numbers and reaching diameters far smaller than the diameters of the cannonballs in question. Bashforth’s measurements with equally spaced distances from the muzzle with an electrical chronographic system were considered a precursor of modern aeroballistic ranges (Miller & Bailey 1979) and an advancement over older measurements, which used less precise ballistic pendula (Didion 1857; Hutton 1812).
3. Results
The main aim of this section is to demonstrate the functionalities of the throwshed tool and interpret the results computed with it. Two different projectiles are considered: an idealised spherical stone projectile released by a sling and a 0.308 Sierra 168 gr bullet shot by a rifle. The temperature, air density and gravitational acceleration values were chosen as standard for the study area described below. Parameters describing the bullet shot are conventional, and the physical parameters of the sling projectile were set to deliver long ranges to showcase the advantages and functionalities of the tool. However, note that the values and shape of the sling stone come from recorded intervals and descriptions; for example, the initial velocity of the projectile reaches 30 to 90 m/s (Denny 2025; Harrison 2006) depending on the throwing technique, the shape of the projectile, and up to 450 g and 5.0 to 6.3 cm in mass and diameter, respectively (Dohrenwend 2002; Korfmann 1973). Lastly, the drag coefficient for such a sphere is empirically well known. Table 2 contains all input parameters, and Table 3 shows drag-to-Mach data for the second projectile. A DTM in a mountainous region of Slovakia was selected to perform the analysis. Figure 5 shows a 1200 × 1200 m location for the sling (A) and a 10 × 10 km location for the rifle (B). It also shows the terrain conditions with the shooting points superimposed. The DTM was acquired from the website of the Geodesy, Cartography and Cadastre Authority of the Slovak Republic (ÚGKK SR 2019) and was resampled from a 1-meter resolution to a 10-meter resolution for the rifle.
Table 2
Values of input parameters of sling and rifle projectiles for throwshed demonstration.
| PARAMETER | SLING PROJECTILE | RIFLE PROJECTILE |
|---|---|---|
| Initial height (m) | 1.700 | 1.700 |
| Elevation (°) | <–90, +90> | <–90, +90> |
| Gravitational acceleration (m/s–2) | 9.810 | 9.810 |
| Initial velocity (m/s) | 70.000 | 808.000 |
| Temperature (°C) | 15.0 | 15.0 |
| Air density (kg/m3) | 1.058 | 1.058 |
| Mass (kg) | 0.2500 | 0.0109 |
| Diameter (m) | 0.0500 | 0.0078 |
| Drag coefficient | 0.47 | Drag-to-Mach data in Table 3 |
| Azimuth (°) | <0, 360> | <0, 360> |
Table 3
Drag-to-Mach coefficient values of 0.308 Sierra 168 gr bullet.
| Mach | 0.000 | 0.800 | 0.850 | 0.900 | 0.950 | 1.000 | 1.050 | 1.100 |
| CD | 0.140 | 0.140 | 0.142 | 0.160 | 0.240 | 0.430 | 0.449 | 0.447 |
| Mach | 1.200 | 1.400 | 1.600 | 1.800 | 2.000 | 2.200 | 2.500 | |
| CD | 0.434 | 0.410 | 0.385 | 0.365 | 0.350 | 0.339 | 0.320 |
Figure 5
DTM with shooting points, location for sling throwshed (A), and location for rifle throwshed (B).
The first of the basic functionalities is simple throwshed, which only ascertains the relation of the DTM cells to the trajectory family envelope. The trajectory is not compared to the terrain. Figure 6 presents the sling (A) and rifle (B) results. Note that only the first band of the raster result is displayed, representing only the cells accessible to direct shooting. Due to the maximum span of the elevation angle range and the method applied, the second band would be identical. Areas within the unsafe zone of just one shooting point are much smaller for the sling projectile due to its much smaller initial velocity. The maximum theoretical range of slingshot varies between 300 and 430 m, which seems exaggerated but is still in agreement with both the experience of modern sling hobbyists and historical and ethnographic accounts that report and estimate the ranges of professional sling shooters up to 500 or even 700 m (Brown & Craig 2009; Harrison 2006). The rifle range varies from 4770 to 5170 m. The range of the weapons on a horizontal plane was superimposed to depict the differences with real range, which varies due to the terrain. Circular range also suggests how the fields-of-fire aspect is usually addressed in the KOCOA method. The smaller circles indicate the range of projectiles shot at an elevation of 0°, and larger circles indicate the range of projectiles shot at 45° angles. Differences between the throwshed range and that of a horizontal plane at 45° are up to 130 m for the slingshot and up to 500 m for the rifle.
Figure 6
Simple 1-point throwshed for sling (A) and rifle (B). Each smaller circle represents the shooting range of a weapon with a shooting angle of 0° in terrain shaped like a horizontal plane. The larger circles represent each weapon’s range with the shooting angle set to 45°.
The second result differs in the number of shooting places incorporated, which was increased to three. Again, simple throwshed was computed both in binary mode in Figure 7 and cumulative mode in Figure 8. Now, the unsafe zone is much larger. Although the binary result represents true or false values, the cumulative throwshed depicts the overlaps of the results of each shooting point and, therefore, the number of weapons that can hit the area.
Figure 7
Simple 3-point binary throwshed for sling (A) and rifle (B).
Figure 8
Simple 3-point cumulative throwshed for sling (A) and rifle (B).
The following simple throwshed was calculated for one shooter with the visibility analysis, which clipped the result to the visible area. Figure 9 displays the clipped analysis result for both projectiles. In both cases, the area of the unsafe zone shrank. This is due to the position of the shooting point and the complex terrain.
Figure 9
Simple 1-point throwshed for sling (A) and rifle (B) clipped by viewshed.
The probable throwshed was also computed with the standard deviations of input parameters presented in Table 4. This result was computed as a simple throwshed with 50 repetitions for one shooting place. Figure 10 depicts the result. Cell values now express the 0–1 probability of each cell being hit. Range differences from lowest to highest probable areas vary from 100 m to 200 m for a sling (A) and between 900 m and 1200 m for a rifle (B). These differences are not negligible when compared to the ranges of the ideal throwshed in Figure 6.
Table 4
Standard deviations of input parameters and DTM for probable throwshed.
| STANDARD DEVIATION OF PARAMETER | SLING PROJECTILE | RIFLE PROJECTILE |
|---|---|---|
| Initial height (m) | 0.100 | 0.100 |
| Gravitational acceleration (m/s–2) | 0.010 | 0.010 |
| Initial velocity (m/s) | 2.500 | 5.000 |
| Temperature (°C) | 5.0 | 5.0 |
| Air density (kg/m3) | 0.050 | 0.050 |
| Mass (kg) | 0.0100 | 0.0005 |
| Diameter (m) | 0.0050 | 0.0001 |
| Drag coefficient | 0.05 | 0.00 |
| DTM elevation (m) | 0.160 | 0.160 |
Figure 10
Simple 1-point probable throwshed for sling (A) and rifle (B).
Fields of fire that are shaped by obstructions on the unsafe zone are addressed in the subsequent classic throwshed. In this mode, the cells in the envelope are assessed by comparing their intersecting trajectories to the terrain profile. Figures 11 and 12 contain classic throwshed results for both projectiles and both bands. The gaps in band 1 of the raster are much larger than in band 2, with only small gaps present.
Figure 11
Classic and direct 1-point throwshed for sling (A) and rifle (B).
Figure 12
Classic and indirect 1-point throwshed for sling (A) and rifle (B).
An analysis was conducted for multiple shooters in the cumulative mode of the classic throwshed. Figure 13 shows the result from band 1 for sling (A) and rifle (B). The number of overlaps can be derived from the image.
Figure 13
Classic and direct 3-point cumulative throwshed for sling (A) and rifle (B).
Another analysis was conducted with walls providing obstructions. The height of the walls was set to 6 m and their width was set to 3 m. These walls were set in the DTM to create the DSM used for the analysis. Figure 14 displays the result with walls represented as lines. The walls created more inaccessible areas for the sling (A). For the rifle (B), this effect is almost negligible due to the scale of the area affected. The azimuth range was changed to 45–90°.
Figure 14
Classic and direct 1-point throwshed using walls for sling (A) and rifle (B) with adjusted azimuth interval.
Finally, a throwshed analysis was conducted for one shooting point in classic mode for both projectiles, in which the elevation angle range was adjusted to 0–60°. Figures 15 and 16 depict how much this affected the results for both bands. The inaccessible areas surrounding the shooting point appear in both direct and indirect cases.
Figure 15
Classic and direct 1-point throwshed for sling (A) and rifle (B) with adjusted shooting angle range.
Figure 16
Classic and indirect 1-point throwshed for sling (A) and rifle (B) with adjusted shooting angle range.
4. Validation
Several methods can be used to validate the throwshed tool results. These include comparing the simulated data with experimental ranges from firing tables, ballistic trajectories simulated in other models, and spatial data for projectiles excavated from archaeological studies. Matching the results of the throwshed tool with these data would indicate the tool’s performance and reliability. Three validation methods were selected.
4.1 Match with the data modelled in the existing literature
Three literature sources computing the range of projectiles were consulted. McCoy (2012) gives several examples of computation of ballistic trajectories. One of the examples calculates the range of a 105 mm high-explosive M1 shell. The result agrees precisely with that obtained by throwshed tool and, in addition, with the characteristics of this shell claimed by the manufacturers for commercial purposes.
Silliman and Batt (2015) determined the range of horizontally shot musket ball numerically to almost 149 m, which varies by 8 m from the range of almost 157 m obtained in the throwshed analysis. Ranges were computed for the same projectile by Miller (2010). For the same input parameters except for the drag coefficient considered as a function, the range is evaluated at 146 m, which varies by 3 m from the range of almost 149 m computed by the throwshed tool. These variations may be caused by different integration methods using time instead of distance steps or by tiny deviations in the drag coefficient used.
4.2 Match with experimental data from firing tables
This validation inspects the match between historical experimental and simulated data from shooting ranges. The experimental data represent distance ranges of the cannons used in the military ordnance of Imperial and Royal Austria at the beginning of the 19th century. These data were chosen for validation because the following subsection is a case study that uses similar technical parameters. Experimental data are derived from firing tables in the Handbook for Austrian Artillery Officers (Smola 1839). Simulated ranges are calculated by the throwshed tool, for which input parameters were acquired from the Handbook, which elaborates on the technical parameters of the cannons, munition, and gunpowder. Each of these parameters warrants individual discussion.
One can easily derive shooting ranges by finding the appropriate elevation angles corresponding to the offsets set by the elevation screw of the cannon, which are given alongside the ranges in the Handbook. Although the clarity of these ranges is indisputable, the surrounding area is only assumed to be a horizontal plane, because the Handbook states that the experiments took place in 1826 at Pest, Hungary, without any site terrain characteristics. The Handbook also specifies the deviation of the range for various types of cannons and elevations of 0–5° as a ratio of the average range. These may be helpful when determining the tolerance allowed for differences between simulated and experimental ranges. Table 5 contains the deviations and the minimal and maximal offsets with corresponding elevation angles.
Table 5
Deviations of experimental ranges and limiting parameters of the cannons.
| CANNON TYPE | CALIBRE (lb) | DISTANCE DEVIATION OF RANGE (%) | MIN–MAX ELEVATION OF THE REAR RING (zoll) | MIN–MAX ELEVATION ANGLE (°) |
|---|---|---|---|---|
| Field cannons | 6 | 11 | 0–10 | 0.6000–11.1000 |
| 12 | 8 | 0–18 | 0.6000–15.4167 | |
| 18 | 6 | 0–20 | 0.6000–15.0000 | |
| Battery cannons | 12 | 7 | 0–30 | 0.5667–16.2333 |
| 18 | 7 | 0–32 | 0.5833–15.8167 | |
| 24 | 8 | 0–34 | 0.6000–15.9500 | |
| Defence cannons | 6 | 9 | 0–20 | 0.4833–13.2333 |
| 12 | 10 | 0–22 | 0.4500–12.0667 | |
| 18 | 7 | 0–22 | 0.5000–11.5500 |
The initial height is roughly 1 m, and the air density is set at a standard 1.225 kg/m3. The diameter and mass of the cannonballs depend on the calibre and sometimes on the cannon used, and four cannonball types are addressed: 6, 12, 18, and 24 lb. These values were derived from the Handbook, converted into metric units, and added to Table 6. Temperature is assumed to be 15 °C, and gravitational acceleration is assumed to be 9.81 m/s2. The elevation angle is set to multiple exclusive values as only a single trajectory is modelled to simulate the distance range. The elevation value depends on the offset of the elevation screw on the rear cannon ring, which can vary 0–8 inches. Some of the values had to be interpolated as did the corresponding experimental ranges. Elevations in inches and degrees and the corresponding ranges are presented in Table 7. As section 2.11 indicates, the drag coefficients were derived from the graphs in Miller and Bailey (1979) as drag-to-Mach tables for all diameters because some of the cannonballs reach up to 2 Mach velocities.
Table 6
Technical parameters of the cannons and cannonballs from the Handbook (Smola 1839).
| CANNON TYPE | CALIBRE (lb) | BORE DIAMETER (m) | BORE DEPTH (m) | CANNONBALL DIAMETER (m) | CANNONBALL MASS (kg) | GUNPOWDER CHARGE (kg) | INITIAL VELOCITY (m/s) |
|---|---|---|---|---|---|---|---|
| Field cannons | 6 | 0.0946 | 1.3568 | 0.0905 | 2.7475 | 0.8400 | 499.68 |
| 12 | 0.1183 | 1.7096 | 0.1139 | 5.4995 | 1.4000 | 471.81 | |
| 18 | 0.1348 | 1.9732 | 0.1305 | 8.2558 | 2.5200 | 497.61 | |
| Battery cannons | 12 | 0.1183 | 2.7424 | 0.1139 | 5.4995 | 2.2400 | 597.10 |
| 18 | 0.1348 | 3.0087 | 0.1305 | 8.2558 | 3.0801 | 575.42 | |
| 24 | 0.1480 | 3.1681 | 0.1436 | 11.0208 | 3.9201 | 561.56 | |
| Defence cannons | 6 | 0.0946 | 2.2572 | 0.0905 | 2.7475 | 1.1200 | 603.35 |
| 12 | 0.1183 | 2.7068 | 0.1139 | 5.4995 | 2.2400 | 595.56 | |
| 18 | 0.1348 | 2.8376 | 0.1305 | 8.2558 | 3.0801 | 568.84 |
Table 7
Simulated and experimental shooting ranges of the cannons according to elevation angles.
| PARAMETERS P: ELEVATION ANGLE EA (°), EXPERIMENTAL RANGE ER (m), SIMULATED RANGE SR (m) | ELEVATION OF REAR RING (ZOLL) | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| CANNON TYPE | CAL. (lb) | P | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| Field cannons | 6 | EA | 0.6000 | 1.7000 | 2.7833 | ||||||
| ER | 379 | 683 | 910 | ||||||||
| SR | 428 | 763 | 1012 | ||||||||
| 12 | EA | 0.6000 | 1.4667 | 2.3167 | 3.2000 | 4.0500 | |||||
| ER | 455 | 721 | 910 | 1138 | 1290 | ||||||
| SR | 423 | 715 | 944 | 1144 | 1312 | ||||||
| 18 | EA | 0.6000 | 1.3667 | 2.1167 | 2.8000 | 3.5333 | |||||
| ER | 455 | 759 | 986 | 1176 | 1366 | ||||||
| SR | 465 | 759 | 989 | 1167 | 1334 | ||||||
| Battery cannons | 12 | EA | 0.5667 | 1.1000 | 1.6500 | 2.2000 | 2.7333 | 3.2833 | 3.8333 | ||
| ER | 379 | 759 | 910 | 1100 | 1252 | 1396 | 1517 | ||||
| SR | 547 | 782 | 974 | 1135 | 1275 | 1272 | 1514 | ||||
| 18 | EA | 0.5833 | 1.0667 | 1.5667 | 2.0667 | 2.5667 | 3.0667 | 3.5500 | 4.0500 | ||
| ER | 379 | 683 | 936 | 1037 | 1214 | 1366 | 1487 | 1631 | |||
| SR | 551 | 774 | 964 | 1124 | 1266 | 1393 | 1504 | 1611 | |||
| 24 | EA | 0.6000 | 1.0667 | 1.5500 | 2.0167 | 2.4833 | 2.9667 | 3.4333 | 3.9000 | 4.3667 | |
| ER | 531 | 759 | 910 | 1062 | 1214 | 1366 | 1467 | 1669 | 1770 | ||
| SR | 556 | 778 | 969 | 1127 | 1268 | 1399 | 1514 | 1621 | 1721 | ||
| Defence cannons | 6 | EA | 0.4833 | 1.1500 | 1.8000 | 2.4667 | 3.1167 | ||||
| ER | 379 | 645 | 860 | 1062 | 1214 | ||||||
| SR | 478 | 738 | 931 | 1096 | 1235 | ||||||
| 12 | EA | 0.4500 | 1.0000 | 1.5333 | 2.0833 | 2.6167 | 3.1667 | 3.7167 | |||
| ER | 379 | 645 | 835 | 986 | 1163 | 1305 | 1517 | ||||
| SR | 486 | 740 | 934 | 1101 | 1241 | 1370 | 1488 | ||||
| 18 | EA | 0.5000 | 1.0333 | 1.5500 | 2.0667 | 2.5833 | 3.1000 | 3.6167 | 4.1333 | ||
| ER | 379 | 683 | 910 | 1100 | 1252 | 1404 | 1593 | 1644 | |||
| SR | 500 | 750 | 946 | 1112 | 1257 | 1387 | 1506 | 1614 | |||
Lastly, muzzle velocity is an essential parameter that differs between cannons of differing calibres. Muzzle velocity is not provided in the Handbook. It is unknown whether this parameter has ever been computed for the cannons in question. However, it may be modelled as section 2.11 suggests. At the time (1826), other nations typically used gunpowder charges approximately a third of the cannonball mass (McConnell 1988), and Imperial and Royal Austria was no exception. Detailed charges can be found in the Handbook (Smola 1839). Table 6 shows these masses converted to kg. The density of the gunpowder is set to 940 kg/m3 according to the Handbook, which states the densities for various types of gunpowder, even cannon-firing gunpowder. Commonly, the length of the bore is provided as a multiplication of the diameter of the bore. However, this is not correct because it refers to the cannon barrel length. The Handbook also provides dimensions, even the length of the inner bore, which is important. Table 6 presents these with the bore diameters d converted into m and computed muzzle velocities. These values are incorporated into the physical model to simulate the shooting range of the cannons.
Simulated ranges are presented in Table 7, where they can be compared to the experimental ranges. It is clear for both types of shooting ranges that the higher the elevation angle is, the larger is the range. Most of the differences could be considered negligible because the corresponding simulated ranges would fit with the dispersion of the real ranges. Some, such as 1272 m with a 12 lb battery cannon with 5-inch (zoll) elevation and 934 and 1101 m ranges with 12 lb defence cannons with 2- and 3-inch elevations, respectively, vary slightly over the 7% and 10% deviations of their experimental pairs. This may be caused by the rounding to hundreds of steps of the experimental ranges in the Handbook (Smola 1839), especially with the battery and defence cannons, where the rounding of the ranges turns from 100 to 200 steps from around 1 inch of elevation. Five of nine simulated ranges with an elevation of the rear ring equal to 0 even deviate by approximately 25% from the paired real ranges. This may be caused due to the initial height, determined primarily by the carriage on which the cannon barrel was mounted. These carriages were unique for each type of cannon. Nevertheless, terrain imperfections in the surrounding testing area could also contribute to this situation. Even slight variations in height, either at the beginning or end of the trajectory, could cause large horizontal shifts for the results, both simulated and experimental, as the trajectories are flatter with lower elevation angles. Diminishing differences with increasing shooting angles could test this premise. Despite these factors, the similarities in ranges may be considered significant.
4.3 Match with archaeological data
This section validates the performance of the throwshed physical model with archaeological data from remains discovered during a 2023 excavation from a real historical event: the 1794–1795 siege of Luxembourg (Waersegers 2023). The French First Republic laid siege to the Fortress of Luxembourg with many artillery batteries, while troops of the Habsburg monarchy defended the Fortress (Kreins 2003). Figure 17(A) shows a map of the Fortress from that time with four forts: Fort Neipperg from the detached works of the Thionville Front, Fort Wallis from the exterior works of the same front, and Forts Rubamprez and Rumigny from the detached works of the Trier Front. Figure 17(B) also shows the remaining ditch of the battery, dug by the French to build an earthen breastwork during the siege. It was then backfilled with the Austrian cannonballs, considered here as the archaeological data. If the simulated throwshed areas overlap with the ditch, the validation may be considered successful.
Figure 17
1794 map of the Fortress of Luxembourg (Lefort 1905) with four forts of interest (A). DTM with the ditch and forts (B).
Several parameters are required to simulate the fields of fire. Four forts, Rumigny, Rubamprez, Neipperg, and Wallis, were selected for their appropriate position within the Fortress facing the target ditch (Waersegers et al. 2023, Waersegers et al. 2025). Their position was determined either by the map by Lefort (1905) after second-degree polynomial transformation within QGIS (Neipperg, Wallis) or by a location available on the internet (Rumigny, Rubamprez). The altitudes of these forts were derived from the report (Ulveling 1869). Only the altitude of Fort Rubamprez was estimated from the difference between the original and current terrain altitudes of neighbouring Fort Rumigny. The coordinates of the shooting places are given in Table 8. The altitudes were also used to modify the original DTM of the surrounding area into DSM. The DTM of the Grand Duchy of Luxembourg, produced by the Administration de la Navigation Aérienne du Luxembourg in 2018 with a resolution of 0.5 m, was resampled to 1 m resolution. The DTM is also shown in Figure 17(B).
Table 8
Coordinates of shooting places in LUREF and NGL data and calibres used within forts (Waersegers).
| FORT | E (m) | N (m) | H (m) | CALIBRE (lb) |
|---|---|---|---|---|
| Rumigny | 78376.00 | 75163.00 | 307.64 | 6 |
| Rubamprez | 78186.00 | 75108.00 | 297.37 | 6, 12, 18 |
| Wallis | 77536.00 | 74403.00 | 305.00 | 6, 12, 18 |
| Neipperg | 78040.00 | 74418.00 | 312.12 | 6, 12 |
The types of the cannons are unknown, but due to their standardisation, with reference to the Handbook (Smola 1839), the character of the event, and archaeological finds, 6, 12, and 18 lb calibres of field, battery, and defence cannons of types predating 1823 were chosen. Information from the Handbook (Smola 1839) might risk inaccuracy considering the date of the siege and changes in cannon type production over almost 30 years. Identifying similar technical parameters for the cannons from older literature (Demian 1812; Unterberger 1807) reduces the risk. Almost 200 cannonballs were recovered during the archaeological research at the location of the ditch. Cannonballs of 6, 12, and 18 lb were recognised, and with the permission of the researchers (Waersegers & Theis, unpublished work), the cannonballs’ average mass and diameter are noted in Table 9 with the cannon parameters. These were used to compute the muzzle velocities. Specific calibres used within forts are also known from Waersegers’s unpublished work. These are presented in Table 8. The remaining parameters are initial height of 1 m, gravitational acceleration of 9.81 m/s2, temperature of 15 °C, air density of 1.225 kg/m3, drag coefficients derived from graphs by Miller and Bailey (1979), and elevation ranges also shown in Table 9. Due to the situation and demanding computational time, the azimuth range was diminished to a range of 110–160°.
Table 9
Technical parameters of the cannons and recovered cannonballs. Mean technical parameters of the cannonballs were provided by Waersegers and Theis (unpublished work).
| CANNON TYPE | CALIBRE (lb) | BORE DIAMETER (m) | BORE DEPTH (m) | CANNONBALL DIAMETER (m) | CANNON-BALL MASS (kg) | GUNPOWDER CHARGE (kg) | INITIAL VELOCITY (m/s) | MIN–MAX ELEVATION ANGLE (°) |
|---|---|---|---|---|---|---|---|---|
| Field cannons | 6 | 0.0946 | 1.3568 | 0.0900 | 2.6426 | 0.8400 | 509.50 | 0.6000–11.1000 |
| 12 | 0.1192 | 1.7096 | 0.1125 | 5.3056 | 1.4000 | 481.78 | 0.6000–15.4167 | |
| 18 | 0.1364 | 1.9732 | 0.1305 | 8.1608 | 2.5200 | 503.00 | 0.6000–15.0000 | |
| Battery cannons | 12 | 0.1192 | 2.7424 | 0.1125 | 5.3056 | 2.2400 | 609.72 | 0.5667–16.2333 |
| 18 | 0.1348 | 3.0087 | 0.1305 | 8.1608 | 3.0801 | 581.41 | 0.5833–15.8167 | |
| Defence cannons | 6 | 0.0946 | 2.2572 | 0.0900 | 2.6426 | 1.1200 | 615.21 | 0.4833–13.2333 |
| 12 | 0.1192 | 2.7068 | 0.1125 | 5.3056 | 2.2400 | 608.16 | 0.4500–12.0667 | |
| 18 | 0.1364 | 2.8376 | 0.1305 | 8.1608 | 3.0801 | 574.82 | 0.5000–11.5500 |
Figures 18, 19, 20, 21 present the simple throwshed; Figure 18 shows the outcome for a 6 lb field (A) and defence cannons (B) situated on all forts. Figure 19 shows the results for field cannons shooting 12 lb cannonballs from Forts Rubamprez, Neipperg, and Wallis (A) and 18 lb cannonballs from Forts Rubamprez and Wallis (B). Figure 20 depicts the throwshed computed for battery cannons firing 12 lb (A) and 18 lb cannonballs (B) from the forts, and Figure 21 presents the simulated fields of fire for defence cannons firing 12 lb (A) and 18 lb cannonballs (B). The throwshed covers the French ditch and its surroundings in all scenarios. This demonstrates a clear overlap of simulated fields of fire with the real data. Additionally, the throwshed confirms the accessibility of the ditch to all types of cannons from all suggested forts.
Figure 18
Cumulative throwshed for the 6-pounder field (A) and defence cannons (B).
Figure 19
Cumulative throwshed for 12-pounder (A) and 18-pounder field cannons (B).
Figure 20
Cumulative throwshed for 12-pounder (A) and 18-pounder battery cannons (B).
Figure 21
Cumulative throwshed for 12-pounder (A) and 18-pounder defence cannons (B).
5. Discussion
The Results section presented two distinctions between shooting on a horizontal plane and shooting in complex terrain. The first distinction is when shooting at a surrounding downhill-inclined area enables a much larger shooting range than the range of the weapon on the horizontal plane, the usual way of addressing the fields of fire in KOCOA analysis, and vice versa for the uphill-inclined area. Hence, simply superimposing the horizontal plane range around the shooting point is not always sufficient, especially with complex terrain. The second distinction is in fields of fire that can be strongly affected by natural or artificial obstacles. Large areas that are inaccessible can influence the tactics of the attacking side. In contrast, the areas which might be invisible to the defending party could be accessible by its projectiles bypassing the obstacles in indirect and, less likely, direct trajectories. This is where the efficacy of the basic range and viewshed clipping combination decreases.
Differences between the bands of the classic throwshed result are present; direct trajectories in band 1 tend to be flatter and, hence, more prone to the cover effect of obstacles than the indirect trajectories in band 2. That is also a reason why the gaps in band 1 will always be larger than those in band 2: because the indirect trajectories access invisible areas more efficiently due to their larger slope.
Multiple shooting points are also analysed because many archaeological studies address numerous places where the weaponry was situated. Here, the results are interpreted in two ways. The binary approach addresses the question of whether the target area is accessible to any of the shooters. Therefore, the cell values are 0 or 1, representing true and false. In addition, the cumulative approach provides information on how many shooters can hit the area in question. The higher the cell value, the higher the risk of damage to an adversary entering the unsafe zone. Archaeological excavations and subsequent analyses of battlefields can indicate the existence of weathered or destroyed walls or ditches, which also alter the fields of fire, as seen in the results.
Outcome clipped by viewshed can be summarised as ‘I shoot at what I see’. This function integrates the visibility analysis into the throwshed computational workflow; therefore, no need arises to compute a viewshed analysis over the DEM separately and later merge the result with the throwshed raster. An unclipped classic throwshed covers more areas than the clipped result. These areas are invisible and yet have been shown to be accessible to projectiles. This reveals another consequence of the terrain’s complexity. To visualise this contrast, an inverse clip by a viewshed was included within the functionalities to delineate invisible yet accessible areas, as opposed to the results shown in this paper.
The probable throwshed analysis played a significant role, because the range differences were not negligible. Not only does it address the uncertainty of questionable projectile parameters, but it can also be used for projectiles whose parameters are represented by mean values with standard deviations derived statistically from numerous samples of same-type projectiles with varying characteristics. This is particularly relevant for pre-industrial weapons, such as slings and bows. If the lacking or diverse physical parameters of sling stones and arrows are expressed in a form applicable to the probable throwshed computation, the uncertainty and diversity of the parameters can still be interpreted probabilistically. The most challenging parameter to obtain is the drag coefficient, because this is not straightforward to measure. This issue is addressed in the section on the aerodynamic drag coefficient. In practice, an approximation by adopting the values reported for comparable projectile parameters in existing literature or an estimation through an airflow simulation around a modelled mean projectile can be made to determine the value of the coefficient. The initial velocity of the projectile is also difficult to derive as it depends not only on the projectile but also on the weapon. Again, approximations can be made from existing literature and by modelling from knowledge of the internal ballistics, similar to the computation of the initial velocity for the cannons in section 2.11. A final option is to fine-tune a single missing parameter by fitting the throwshed model to archaeological data where a shooting point and a sufficient sample of projectiles similar to the mean projectile have been recovered.
Inaccessible areas appear closer to the shooting point when the elevation angle range is not at its full span from –90 to +90°. The higher the minimal elevation angle, the larger this inaccessible area is for direct trajectories. The lower the maximal elevation angle, the larger the inaccessible area is for indirect trajectories, and if it is lowered sufficiently, covered areas can disappear completely. This behaviour depends strongly on the conditions of the terrain.
The validations of throwshed with a few types of weapons and projectiles, especially cannons, show a match with simulated data from recent literature, experimental data in firing tables from older literature, and archaeological data describing a real conflict event for which all shooting places analysed achieved access to the target. Section 4.3 may serve as a foundation for following use cases. In the event of preliminary archival research on the conflict period, collecting and implementing as many consistent parameters as possible into throwshed (proven or unproven exchanges of fire, positions, types of artillery, projectiles, powder, practical and maximum ranges, fixed or free orientation of artillery, topography, terrain accessibility, natural or built cover) already makes it possible to confirm or refute the ballistic potential of the attacker or defender. From a predictive perspective, by systematizing and cross-referencing variables between attackers and defenders, probability calculations would then allow for a clearer delimitation of new areas of archaeological potential within the landscape, as seen in the case of an artillery battery besieging a city. Throwshed could serve as a specific tool complementary to LiDAR or geophysical prospecting techniques to more precisely guide research zones or areas of high archaeological potential. More empirically, the discovery of projectiles or artillery pieces on a site without any known historical or archaeological context would also allow for the approximate orientation of archaeological monitoring zones to determine the location of the opposing camp. This latter scenario is obviously highly dependent on other initial guiding factors to further the process, such as the identification of the combatant side or a chronological clue.
The accuracy of the raster is influenced in several ways. The properties of the physical model do not match the performance of the advanced 6-DOF models, which also exhibit some degree of inaccuracy. The quantity of variables and the spatial extent of the trajectory impose high demands on the estimation of all the conditions and will always cause uncertainty. The propagation of input parameter errors is laborious due to the model’s numerical character, though the procedure of the probable throwshed analysis endeavours to mitigate this issue. Moreover, a few measures that seek to reduce the duration of calculations degrade the real image. This encompasses the construction of an envelope by a definite number of trajectories, searching for the trajectory with the most distant reach and intersecting the cell near its centre. Subjecting all these techniques to the proposed half-cell stopping criterion leads to some degree of omission. Finally, the accuracy of the physical model and of the procedures of the method are related to 3D space, but the outcome of the tool is in 2D. It is difficult to derive specific values that describe the overall accuracy. Obscurity of accuracy required the validation of the analysis.
Lastly, several summaries and use cases for conflict archaeology can be outlined to demonstrate where throwshed analysis may outperform standard approaches depending solely on linear distances. One key factor is that, in complex terrain with features such as hills, valleys, terraces, and cliffs, theoretical range can vary and inaccessible areas can occur. Similar effects may also arise in flat landscapes due to anthropogenic obstacles such as palisades, walls, ramparts, and towers. When the location of such obstructions is known but their dimensions are uncertain, the distribution of recovered projectiles, together with fine-tuning of throwshed parameters, can assist in approximating the height and width of obstacles. Artificial barriers at or near the shooting position can also influence the range by increasing the shooter’s elevation despite the terrain being flat. Furthermore, projectile finds considered anomalous in excavations or field survey contexts may be explained through throwshed modelling. The KOCOA methodology used in conflict archaeology could benefit from throwshed outputs because these can help define the advantages and drawbacks of both sides of a conflict, for example by identifying the effective firing positions, constraining the azimuth range that should cover the avenues of approach, and delineating safe areas for an attacker where obstructed blind spots are present. Such information can resolve otherwise inexplicable decisions and enrich understanding of the tactics used in the battlefield. When combined with archaeological spatial data, throwshed analysis may also extend and delimit the broader area of archaeological interest in projectile remains and potential impact targets. Because projectiles were often collected for reuse after a conflict, throwshed modelling can also indicate areas that may be underrepresented in the archaeological record and therefore appear ‘negative’ during prospection. Overall, the proposed methodology advances the narrative from a simple question of range, ‘how far?’, to a spatial question, ‘where?’, by mapping plausible firing and throwing impact zones as a spatial shed across a 3D surface: hence throwshed.
6. Conclusion
The method for throwshed analysis presented above locates areas potentially accessible to projectiles thrown or shot from any given place where the terrain conditions are available as a DEM as well as the physical parameters of the projectile and the atmosphere. The graphical results demonstrate how the terrain’s features and relative elevations were recognised and rendered as distinct idealised accessible and obstructed areas and varying ranges. These constitute the fields of fire and cover aspects of the KOCOA methodology used in conflict archaeology. Additionally, it may be concluded that throwshed can be used in several ways within the archaeological field: before a potential or future operation, to guide a research project or delimit new zones of archaeological vigilance and potential in the context of planned or preventive archaeology; or after an operation, once archaeological remains are contextualized, for the purpose of factual clarification and refining historical sources.
Moreover, validation affirmed the simulation of the external ballistics. Validation with archaeological data also demonstrated that the technical parameters of weapons and projectiles, although demanding, may be researched or computed to initiate the simulation. Omitting various physical effects and introducing methodological resolutions that solve the problems distort the overall outcome to varying extents. These resolutions were deployed to optimise the program, make the tool available for users with ordinary computational power, and ease the acquisition of suitable input data.
In a terrain that closely resembles a horizontal plane, the simple superimposition of a weapon range perimeter and generation of a viewshed may suffice within the KOCOA analysis to address the fields of fire and cover elements. The aim of this paper is to introduce an analysis with evidence in graphical form demonstrating how the range of a weapon may change with high elevation differences on a battlefield. The evidence also suggests that more complex terrain leads to large areas that are invisible to viewshed analysis but were still accessible to the weapons used. Moreover, should the firing tables for the weapon be unavailable, no shooting range can be determined. Throwshed tool addresses these issues and conveys binary, cumulative and probabilistic results. Functionalities demonstrated by multiple setting options could be of help in various use cases, as discussed above. These facts and the validation of the physical model by other simulated, historical, and archaeological data support the conclusion that throwshed analysis can advance the KOCOA methodology.
Certainly, there are limitations, such as no multiple points and no drag-to-Mach tables available for supersonic projectiles in the probable throwshed and the significant issue with the oscillation effect over short distances. With new procedures, these deficiencies are expected to be reduced or eliminated. One anticipated development is the deployment of the tool as a GIS plugin with its own graphical interface in the QGIS program. Additionally, although an extension to the 3-DOF model with the Coriolis effect and the influence of the wind has been created, handling 3D trajectories has yet to be implemented in the whole model. A method for computing terminal ballistics is also anticipated, because the maximum terminal-effect range of projectiles varies from the maximum theoretical range. Optimisation of the code will always remain the primary task. Lastly, updating to the 6-DOF model will remain an option.
Data Accessibility Statement
The throwshed tool scripted in Python language is openly available at a public GitHub repository: https://www.github.com/Tadeo98/throwshed2. Throwshed is licensed under CC BY-NC-SA 4.0.
Archaeological data for cannonball diameters and masses are not currently accessible because they appear in another archaeological article yet to be published.
Author Contributions
Conceptualisation, methodology, programming, validation, research, and writing: Červík, T.
Archaeological data supply and help with research: Waersergers, Y.
Supervision, review, and editing: Zischg, A., Lieskovský, T.