#### Abstract

*Aims*. The complex dynamics of bodies, originating from the interplanetary matter and passing through Earth’s atmosphere, defines their further position, velocity, and final location on Earth’s surface in the form of meteorites. One of the important factors that affect the movement of a body in the atmosphere is its shape and orientation. Our goal is to model the interaction of real shape meteoroids with Earth’s atmosphere and compare the results with the standard spherical body approach. *Methods*. In the simulation, we use 3D models of fragments of the Košice meteorite with different sizes and shapes. Using a 3D model of fragments, we consider the real shape of the body to define its resistance properties during atmospheric transition more specifically. The simulation is performed using virtual wind tunnel in the MicroCFD (Computational Fluid Dynamics) software to obtain more realistic drag coefficients and using the µ(m)-Trajectory software to model the particle trajectory in the atmosphere including the wind profile. The final outputs from these programs are the drag coefficient as a function of the altitude and the particle orientation. Using these parameters we get the more realistic body trajectory and the impact area coordinates. Comparison of the results for real and spherical model meteorite impact location is discussed. *Results*. Simulation showed significant differences in trajectory and the impact area for the different real body orientations compared to the spherically symmetric body. Also, an important result is a difference in the impact area of the real body with a specific orientation without rotation and the body with considered rotation. The significant difference between the modeled impact of a real shape body and its real place of finding compared to a spherically symmetric body indicates the importance of the method used.

#### 1. Introduction

The components of interplanetary matter are rich sources of information that serves to describe and better understand the processes associated with the origin, evolution, structure, and dynamics of the Solar System bodies. Material from meteorites falls is still rare and its collection is important for subsequent study of physical, mineralogical, and chemical properties. The study of the meteorites is dependent on the instrumental records and following the effective determination of impact area and their recoveries.

The interaction of a meteoroid with Earth’s atmosphere is associated with several phenomena related to the 4 regimes of the atmospheric motion: preheating, ablation, dark flight, and impact. An important phenomenon accompanying preheating and ablation is ionization. The dark phase of flight and impact are characterized by nonemission of radiation [1]. Based on this is a description of dynamic, geometric, and physical properties of meteoroid, during the interaction with the atmosphere, determined from instrumental observations. Detailed physical procedures and models are given in Ceplecha [2] which includes a description of a meteorite fall originating from the observation of a fireball. Despite the fact that we have tens of thousands of cataloged meteorites, there are only 38 meteorites with known heliocentric orbit, by October 2021 (https://www.meteoriteorbits.info/, Meteoritical Bulletin Database). Since 1959, there have been significant cases such as Prˇibram, the first meteorite with the known heliocentric orbit [3, 4], Almahata Sitta, observed before entering to the atmosphere [5], Križevci, the first meteorite fall determined from video observations with the most accurately known orbit [6], Chelyabinsk, the largest meteorite with the known heliocentric orbit [7], Neuschwanstein, almost identical orbit with the meteoroid Přibram [4], and others.

In this study, we used a detailed documented and analyzed case of the Košice meteorite fall [8] on February 28, 2010, in the area west of the town Košice in eastern Slovakia. Its trajectory, fragmentation model, orbit definition, and impact area of fragments are given in Borovička et al. [8]. A detailed analysis of the impact area is given in Tóth et al. [9]. Subsequently, the material and magnetic properties (density, porosity, magnetic susceptibility) of the Košice meteorite fragments were defined [10]. Spectroscopic analysis to determine the presence of iron-containing compounds was discussed in Sitek et al. [11]. Statistical analysis of the mass distribution for fragments can be found in Gritsevich et al. [12] and Betzler and Borges [13].

The aim of our study is to obtain more realistic shape of trajectories and impact locations of individual fragments, using dynamic interaction of the lowest part of the planetary atmosphere (the dark phase of flight) with 3D models of fragments of the Košice meteorite. The first analyses of this issue were discussed in our previous work [14, 15]. The study was done by using two software tools. To determine the resistance properties of a real body model we use the virtual wind tunnel MicroCFD. To define the shape of the trajectory during the dark phase of flight and the individual places of impact of the fragments we use the *µ*(m)-Trajectory program, our own software.

In the paper, we analyze the resistance properties, the drag coefficient Γ of the three real 3D scans of meteorites of different sizes and shapes (for masses of different orders, ∼2000 g, ∼200 g, and ∼20 g). We compare the dynamic properties of the used real shapes of meteorite fragments with the standard approach from Ceplecha [2]; in which a symmetrical body shape with corresponding resistance properties is used. At the same time, we compare the places of impact obtained by the Borovička [8] model for the Košice case with the impact area of bodies with real shapes obtained by our simulations.

#### 2. The Dark Phase of Flight

Defining the shape of meteoroid trajectory and impact area is based on the dark phase of flight modeling. Modeling of the lowest flight phase or dark phase of flight is based on a mathematical model by Ceplecha [2]. During this part of trajectory, the body and its fragments cease to emit radiation and move on ballistic trajectories under the influence of gravity and interaction with Earth’s atmosphere. The body’s velocity during the dark phase of flight is in the order of hundreds to tens of ms^{−1} and is significantly lower than its entry velocity. The wind field in individual level of altitudes has a significant effect on the body motion during this phase. The result of this effect is a change of the original direction of flight. The typical impact areas have dimensions of several km^{2}. In addition, the dark flight phase is affected by the body shape [1].

#### 3. The Drag Force and Determination of the Drag Coefficient Γ

The study of fluid flow around the body was first investigated by Sir Isaac Newton in “Principia Mathematica” [16], where he expected a direct dependence of the fluid resistance and velocity of a moving body. Empirical studies by Stokes and Navier expanded this issue and created a set of laws of motion in a real fluid [17]. One of these laws defines the term drag force. The drag force is a resistance of the liquid or gas acting opposite to the relative motion of the body with respect to a surrounding fluid. It is well described by Newton’s law (in the case of supersonic velocities by Modified Newton’s law, Figure 1) and it depends on the fluid flow velocity or the relative velocity of the body [18]. The drag force is proportional to the velocity of the body for a laminar flow and the squared velocity for a turbulent flow [19, 20]. Laminar or turbulent fluid flow around a body is defined by a dimensionless variable, the Reynolds number Re, which correlates the ratio of inertial and viscous forces acting in the flow [21]. During the dark phase of flight, there is a significant change in the velocity of the meteoroid and it is necessary to take into account a broad interval of the Reynolds number. Based on experiments that analyzed the supersonic flow of fluid or air around a smooth sphere, the Reynolds number is in the interval 10^{5}–10^{6} which represents strong turbulent flow and low viscosity [22]. Based on Newton’s theory, the drag force for turbulent flow is defined as where *F*_{d} is a drag force that is identical to the relative direction flow of the fluid, *ρ* is a density of the fluid present, *u* describes the relative velocity of the body with respect to a surrounding fluid, *S* is effective cross section, and Γ is the drag coefficient. The drag coefficient is a dimensionless quantity that is used to quantify the resistance of an object in a fluid and it is always associated with a particular surface area. The value of the drag coefficient depends on the relative velocity, orientation, and size of the body, density of the fluid present, and its viscosity [21]. We can rewrite equation (1) in the form

Equation (2) expresses the calculation of the value of the drag coefficient using the drag force. In the case of a meteorite flight through the atmosphere, the resistance properties are described by the value of the drag coefficient Γ as a function of the velocity in Mach (M) (equation (2)).

The standard model is based on certain approximations where the symmetrical shape of the body with a tabulated function of the drag coefficient (Table 1) and without rotation is considered [2]. However, most models use a constant drag coefficient Γ of 0.7, when describing meteorite falls with relatively good results (e.g., [8]). It should be remembered that the experimental data indicate the drag coefficient as a strong function of the Mach number, especially in the low supersonic and subsonic regimes [23, 24]. The search for a useful fit for this experimental data is discussed in Carter et al. [25] (see Figure 1). Experiments that define the dependence of the drag coefficient with the Reynolds number for a smooth sphere [22] are fitted in the article of Morrison [26].

In our study, in the first step we focus on the use of a spherically symmetric body and constant value of the drag coefficient Γ of 0.7, and then using the steps given in Section 5, we obtain real values of the drag coefficient for real shapes of meteorites.

#### 4. Used Software

An important factor in describing of the meteoroid dynamics during the dark phase of flight is the body shape and its rotation and orientation. The result of the specific shape of the body for a given altitude and orientation in the flight’s direction is a specific value of the drag coefficient Γ. Its value or the value of the drag force can be obtained by simulating real physical conditions during the body flight through the atmosphere. For this purpose is used a virtual wind tunnel and in this paper we use the commercial program *MicroCFD* (Computational Fluid Dynamics, https://www.microcfd.com/). The functionality of the program, the physical background, practical testing, and its use are described in detail in the articles of Rohde [27–30]. The calculated values of the drag coefficient Γ as a function of altitude enter into calculation of the dark phase of flight model.

The virtual wind tunnel *MicroCFD* (Figure 2) simulates the interaction of flowing air (with realistic viscosity) with 3D models of bodies. In the interaction, it is assumed a high Reynolds number and turbulence. The input files are 3D stereolithographic models of type.stl (standard CAD format). The surface of the body in this format is made up of individual triangulating shapes. A 3D model of a real body is created by scanning of meteorite fragment over a full 360° circle with the specific angular step using tomography equipment [31]. The input parameters to the MicroCFD program are the velocity of air flow in Mach units, the fluid pressure in hPa, the temperature in K, the gas constant that describes the type of gas used for the simulation, and its average specific heat. The values of the input parameters are directly related to atmospheric conditions at selected altitudes.

The simulations in the virtual wind tunnel are taking place in lower atmosphere under 30 km. As the output of the simulation are the values of the drag force acting on the body resolved into three orthogonal components *X*, *Y*, and *Z*. The resulting force is the sum of aerodynamic forces acting on individual triangulating areas of the 3D body model, where their mutual orientation is considered [32]. We also obtain moment of forces acting in three orthogonal components *X*, *Y*, and *Z* and effective cross section of the body in two perspectives (*XY*, *YZ*).

The second software used in this work is the *μ*(m)-Trajectory program (Figure 3). It is a complex tool with a user interface, created in Python by the first author of the paper. The program computes the motion of a body in the atmosphere, in three orthogonal components based on the aerodynamic model of the dark phase of flight [2]. The program also works with the changing values of the drag coefficient as a function of altitude and thus enables use of the wind tunnel results.

While the standard aerodynamic model is based on the change of the body motion as a function of altitude, for our needs we use equations rewritten as a function of time. The primary part of the model is the equations of motion defining the change of the three orthogonal components of velocity over time. The equations of motion include gravitational and wind effects on the body, considering the presence of the Coriolis force. For each integration step, the length of the body’s trajectory, the change of its geographical coordinates, and the change of its azimuth and zenith distance are calculated [2]. Our tests showed that an integration step size of 0.1 seconds is sufficient for the accuracy of the calculation. The simulation of the body motion, based on aerodynamic equations, uses the numerical method of Runge-Kutta of 4th order.

The program uses data from the MSIS-E-90 atmosphere model [33] and from University of Wyoming server (http://weather.uwyo.edu/) downloads wind altitude profiles nearest to the meteorite fall. Initial variables define location in the atmosphere (latitude, longitude, and altitude), azimuth, zenith angle, and physical conditions of the meteoroid in the beginning point of the dark phase of flight (initial velocity, density, and radius of the meteorite body). At the same time, the integration step, the body rotation velocity (if considered), and value of the drag coefficient are defined (constant or variable value). The output of the program provides data such as body position and its velocity for each integration step, the geographic coordinates of the impact area, and a graphical representation of trajectory and impact directly in Google Earth application.

#### 5. Procedure

The simulation goes in two iterations steps, which we will call first (I.) and second (II.):(1)We enter to the first iteration step (I.) with a symmetrical body in the form of a sphere. From the known volume and density of the fragment [10], the radius *r* for the spherical shape of the body with a mass equivalent to the real meteorite entered to the *μ*(m)-Trajectory program is determined.(2)In the *μ*(m)-Trajectory program, parameters of the dark phase of flight for the spherical model are calculated using the input values (wind field, atmosphere, position of the body, etc.), the velocity of the body in the beginning point of the dark phase of flight, and the drag coefficient Γ for the spherical body.(3)Due to unusually long calculation time in MicroCFD software (in the order of hours), three to four altitudes are selected from the dark phase model based on point 2. For each altitude, based on the model, we know the velocity of the body and we can determine the density of the atmosphere. For these selected parameters, new values of the drag coefficient Γ are determined (Section 3) using the MicroCFD wind tunnel, where 3D models of real fragment are used instead of the spherical body. This method provides more accurate values of the drag coefficient Γ for different orientations of the body during a flight.(4)The new values of the drag coefficient Γ, for each orientation of the body, are plotted as the function of altitude, Γ(h) (e.g., Figure 4), fitted by the least squares method.(5)We enter to the second iteration step (II.), where the new parameters of the dark phase with more precise velocity at specific altitude are calculated again using the *μ*(m)-Trajectory program. In this step, into the calculation enter the values of the drag coefficient obtained from the function Γ(h) (point 4). Optionally, we can define the period of the body rotation.(6)The simulation process is repeated once again in points 3 to 5, where to the point 3 are entered the new velocities of the body in selected altitudes.

This procedure ensures the improving of the value of the drag coefficient Γ, obtaining the function Γ(h) for all defined orientations of the real 3D model of the body and impact area for the individual orientations of the fragments on the trajectory.

In the MicroCFD we create simulations of the drag coefficient in six orientations for each fragment. The first three orientations are in the direction of the *X*, *Y*, and *Z* axes of the body, and the remaining are in the opposite directions, which will be referred to as *IX*, *IY*, and *IZ*. For conditions in the specific altitude of the atmosphere, we simulate the value of the drag coefficient for four altitudes and in different orientations. So for each meteorite fragment we performed 24 simulations (6 orientations × 4 altitudes).

At the output of the *µ*(m)-Trajectory we get an impact area composed of the individual impacts relating to the orientations of the body *X*, *Y*, *Z*, *IX*, *IY*, *IZ*. In addition, we obtain two impact locations for the simulated rotation of the body (links with the period of the body angular velocity, point 5) in gradual (REG) and random (RAN) setting of individual orientations in one-second interval. In this procedure, for each meteorite, we obtain at the output one impact point for the spherical body from the first iteration step (I), six impact points for individual orientations of the real body model, and two impact points for REG and RAN rotation of the real body model from the second iteration step (II).

#### 6. Input Parameters and Used 3D Models

In our simulation we use three 3D models of the fragments of Košice meteorite representing three orders of mass: catalog no. 53 (23.2 g), no. 21 (208 g), and one uncatalogued fragment of 2370 g [9] (Figure 5).

The calculation of the dark phase of flight begins at the altitude in which the Košice meteoroid ceased emitting radiation (beginning point of the dark phase of flight, 17.4 km). We do not consider changing the shape of the body or its mass after entering the dark phase of flight. However, there is a continuing ablation (a very short time) after crossing the terminal point of the fireball visible by particular observing system [34]. In our simulation, we consider the beginning point of the dark phase of flight to be fixed and we do not work with additional ablation. In the future, it is necessary to define this beginning point in the altitude where the ablation process is definitively terminated. The velocity of the largest fragment in the beginning point has been estimated at 4.5 km/s [8]. For the meteorite models no. 21 and no. 53 we set the velocity in the beginning point to 2.5 km/s. This value is based on the fact that we cannot associate the origin of smaller fragments with a specific fragmentation, so based on Borovička et al. [8] we assume their origin at higher altitudes of the Košice fireball trajectory. Naturally, their velocity decelerated to smaller value than the velocity of the largest fragment in the beginning point of the dark phase of flight. Also, we note that the value 2.5 km/s comes from the meteoroid fragmentation model as the velocity limit for meteoroid ablation [8]. However, the presumption of the value of the initial velocity can be set arbitrary in our case, because we consider only the relative shifts of the impact points of the real body individual orientations compared to the spherical body and do not search for an absolute solution.

For the largest fragment model, we start from the original geographical position of the beginning point (48.746N, 21.083E) [8]. For the models no. 21 and no. 53, we slightly modify the real geographical location of the beginning point (at a given altitude) in order to get the closest impact point to the real place of finding the meteorite. Then we compare the relative change of dark flight simulation for spherical and real shape fragments. To determine the function of the drag coefficient Γ(h) four altitudes were chosen: 15 km, 13 km, 9 km, and 5 km.

#### 7. Simulation Results

After performing the simulation using 3D models and passing through both iterations steps, we obtain the impact area for our meteorite fragments. For each orientation of the fragment in combination with the selected altitudes in which simulation was ongoing, we can create a function of changing the drag coefficient Γ(h) (for six different orientations). The shape of the function is obtained by fitting using the polynomial of the 2nd or 3rd degree. The choice of the fitting function in the form of a polynomial is linked with the shape of the predicted function of the drag coefficient depending on the velocity in Mach units (Figure 1). In the case of the meteorites used, the initial velocity of 4.5 km/s or 2.5 km/s represents approximately a value of 13.5 Mach or 6.5 Mach (for the altitude 17.4 km). However, in the range of altitudes of 15–13 km, the velocity of the bodies decreases to a level of approximately 1 Mach or below, which represents the left side of the function on Figure 1 and a sudden decrease in the value of the drag coefficient Γ. This is in qualitative agreement with our data depending on each individual meteorite shape and orientation. The overall physical behavior corresponds to theoretical expectations. The fitting function is based on the minimum value of root mean square deviation. Now we focus on the results of 3D models of the larger fragments in the following paragraphs.

##### 7.1. Relevance of Fitting Function Shapes

Before moving on to simulating individual fragments, let us look at the degree of reliability of the results, using a smaller number of altitudes. To verify the objectivity of the curve shapes of individual orientations, we performed a test simulation of the body at different altitudes than the default ones (5, 9, 13, and 15 km) used. We focused on fragment no. 21 (208 g), specifically on the *X* orientation. Simulations were made at eight selected altitudes: 5, 7, 9, 11, 13, 14, 14.5, and 15 km. The resulting curves of the drag coefficient Γ(h) with default altitudes and new altitudes, using a polynomial of the 3rd degree (equation (3)), can be found in Figure 4.

The obtained curves in Figure 4 are slightly different. The largest relative difference between the curves of the drag coefficient Γ(h) is 0,1. This difference may be affected by the shape of the trajectory when simulating the dark phase of flight and thus change the final shape of the impact area. From the acquired results, we consider that it would be more suitable to define at least twice the number of altitudes, but based on the time-consuming simulation in the wind tunnel, we worked with default four altitudes.

##### 7.2. Košice Meteorite, Uncatalogued Largest Fragment (2370 g)

The simulation of the dark phase of flight of the largest Košice fragment was performed without complications. The resulting values of the drag coefficient Γ for 3D model of the largest fragment are given in Table 2.

For all orientation, the polynomial fitting of the 3rd degree was used (in general form, equation (3)). The obtained shape of the function is used for calculation of the values of the drag coefficient Γ in each step of the integration of the dark phase of flight (Section 5, points 4 and 5).

The variable *x* represents the altitude, in which the body is located at a given time. The obtained functions of the drag coefficient Γ(h) for individual orientations are plotted in Figure 6.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The shape of the drag coefficient Γ(h) curve for each orientation of the fragment has a similar character. Decreasing shape of the curve can be explained by a significant change of the fragment velocity from 1.4 km/s (15 km) to 0.3 km/s (13 km). Subsequently, at an altitude of approximately 7 km, the value of the drag coefficient increases. This increase can be explained by the transition of the fragment into a denser atmosphere. The velocity on the ground level was from 46 m/s to 53 m/s depending on the orientation.

The localisation of individual impacts for the spherical model, the model computed by Borovička et al. [8], and the body model with real shape for different orientations and rotation in prograde and retrograde mode are in Figure 7. The impact distance of the spherical model compared to the real place of the recovered fragment impact is 640 m. The impact distances of individual orientations of the body with real shape, gradual (REG), and random (RAN) rotation ranges from 340 m (orientation IY) to 490 m (orientation RAN) compared to the real place of impact. Naturally, rotation somehow averaged the atmosphere interaction with specific orientation shapes and cross sections. That is way, the RAN and REG simulations ended in the middle of 6 oriented flights. The orientation *Z* has a significantly different geographical position of the impact. In this case, we assume that the fragment was oriented to the side with the largest effective cross section. The given result points to an important influence of the orientation and the real shape of the body on the change of the trajectory and position of the impact place. It is also necessary to mention that the individual impacts of the real body keep the flight direction (with a small shift to the south).

##### 7.3. Košice Meteorite, No. 21 (208 g)

Simulation of the dark phase of flight of the Košice fragment no. 21 was successful in all aspects. The resulting values of the drag coefficient Γ for 3D model of this fragment are given in Table 3 and plotted in Figure 8.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

For all orientations, the polynomial fit of the 3rd degree was again used (equation (3)). Similar to the largest fragment of the Košice meteorite, the obtained shape of the function Γ(*h*) is used to calculate the values of the drag coefficient in each step of the integration of the dark flight (see Section 5, points 4 and 5).

The shape of the function for the drag coefficient Γ(h) in each orientation of the body follows a similar pattern. A significant difference between the curves can be seen in case of the *Z* and *IZ* orientation. This change in the shape of the curve and the consequence change in the values of the drag coefficient are related to a significant change of the cross section for these two orientations. We did not expect a lower value of the drag coefficient for the *IZ* orientation at an altitude of about 15 km compared to other orientations, but we assume that the conical shape for the given orientation plays a significant role. As a result there is a significantly different position of impact for *IZ* orientation (Figure 9). The value of the drag coefficient Γ from an altitude of 14 km is less than 0.5 (except for *Z* and *IZ* orientation). We see the connection in a smaller cross section compared to the study of the largest fragment. Considerable deceleration of the fragment from 2.5 km/s (altitude 17,4 km) to 0.1 km/s (altitude 15 km) resulted in decrease of the drag coefficient. The gradual transition to the denser part of the atmosphere leads to subsequent increase of the drag coefficient, which is observed approximately from an altitude of 11 km. The velocity of the fragment on the ground level was from 36 m/s to 46 m/s depending on the orientation.

The distance of the modeled impacts compared with the real position of the recovered fragment in impact area (Figure 9) is variable depending on the orientation of the fragment. The final values of individual orientations and gradual (REG) and random (RAND) rotation ranges from 440 m (orientation RAN) to 720 m (orientation *IZ*) compared to spherical body of equivalent effective cross section. The absolute differences depend on input conditions of simulations. The spread of impact area is within 1620 m with minimal side spread (few tens of meters), but the majority of impact position are within 450 meters. In Figure 9 we see maintaining of the line of the predicted direction of flight. The impacts of individual orientations in one line are generally related to the greater effect of the wind field on the less massive fragment.

##### 7.4. Košice Meteorite, No. 53 (23.2 g)

Simulation of the dark phase of flight of the Košice fragment no. 53 was also successful. For this fragment, which is close to the spherical shape, we performed 48 simulations. The resulting values of the drag coefficient Γ for 3D model of Košice fragment no. 53 are given in Table 4.

In the case of the Košice fragment no. 53, in terms of minimum value of the root mean square deviation, the polynomial fitting of the 2rd degree was used (in general form, equation 4).

Similar to the previous fragments of the Košice meteorite, the variable *x* in equation (4) represents the altitude of the fragment. The final drag coefficient values as a function of the altitude for individual orientation of the fragment are plotted in Figure 10.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The shape of the function for drag coefficient Γ(h) is different compared with curve shape for fragment no. 21 or the largest fragment. The drag coefficient curve has an increasing character and the value of Γ is < 0.5 in each orientation. This result is due to more efficient deceleration of the smallest fragment and based on this, our simulations define the area behind the sudden decrease of the function in Figure 1. At the same time, the nearly spherical shape of the fragment leads to a lower drag coefficient value and less effective deceleration in the lower atmosphere. At an altitude of 17.4 km, the body enters into the dark phase of flight with a speed of 2.5 km/s. This value is identical to the fragment no. 21. However, the velocity 0.3 km/s was reached at an altitude of 13.6 km (for gradual orientation). The velocity on the ground level was from 31 m/s to 42 m/s, which is only approximately 5 m/s less than fragment no. 21.

The impact area for the individual orientations of the fragment is shown in Figure 11. Impact distances in comparison with a real position of the fragment impact (or spherical body of equivalent effective cross section) ranges from 560 m (orientation *X*) to 1210 m (orientation *Z*). In some cases, the distance between impact points for individual orientations is only a few tens of meters. This is caused by the fact that the real shape of the meteorite is close to the sphere. Similarly like in the case of the fragment no. 21, the relative spread of the individual impacts is within 680 m along the flight trajectory with minimal side spread (few tens of meters).

#### 8. Conclusion

The aim of the study was to process, evaluate, and analyze the obtained results from simulation of the dark flight for three real shape fragments of the Košice meteorite. The simulation was fully successful for meteorite fragments, where using the MicroCFD virtual wind tunnel we have made the description of the drag coefficient Γ in 6 orientations and in different altitudes. The results of the dark phase of flight, using the *µ*(m)-Trajectory program, brought significant shifts of the impact area compared to the spherical shapes of the fragments, which proved the necessity of such precise simulation.

However, we must keep in mind the inaccuracies entering into the simulations of the dark phase of flight. On the one hand, from the end point of the light phase of flight, physical and location parameters enter into the simulation with standard deviations: geographical position of the body, velocity, zenith angle, azimuth, etc. On the other hand, the atmospheric conditions such as air density, wind field, and the associated drag coefficient Γ also have an impact on the overall measurement inaccuracy. These parameters come from real measurements using meteorological balloons, often tens to hundreds of kilometers away from the event, or from the meteorological models that extrapolate real data. Figure 4 also shows the different shape of the drag coefficient curves Γ(h) and as we showed it is sensitive on velocity, real shapes of meteoroids, and the choice of the number of altitudes where the simulation is performed and thus can impact results.

In general, we pointed out that used methodology in the form of simulation through several iteration steps has its importance. However, it is necessary to gradually improve this methodology and find a solution for its more effective processing. Also, it is important to process simulations for a larger number of fragments of different types, shapes, and masses. We aim to create a dataset that will contribute to the effective definition of impact areas for a wide range of meteorite falls.

In the future, to obtain the more realistic model of the dark phase of flight and the effective definition of impact area, we plan to create more functions and modifications in the program. One possibility is to generate a statistically large number of particle clones in combination with using the probability density function of random variables, kernel density estimation (KDE) [35], with the purpose of estimating the subarea of impact with the statistically highest concentration of particles and, at the same time, the use of random variations in wind speed and direction for each particle at a level ±10–20%, which are based on the uncertainties of meteorological models [34]. Finally, in the simulation, except for the use of real shapes of fragment, it is appropriate to test several empirical or theoretical models of the drag coefficient [26, 36].

#### Data Availability

The authors can provide background data that support the results of their study upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the Slovak Research and Development Agency under Grant nos. APVV-16-0148 and APVV-18-0103, by the Slovak Grant Agency for Science, Grant no. VEGA 1/0596/18, and by the Grants for PhD students and young researchers at Comenius University in Bratislava, Grant no. UK/174/2019.