MNRAS 000, 1–11 (2020) Preprint 9 June 2021 Compiled using MNRAS LATEX style file v3.0 Accurate 3D fireball trajectory and orbit calculation using the 3D-FireTOC automatic Python code Eloy Peña-Asensio1,2★ , Josep Maria Trigo-Rodríguez1,3 , Maria Gritsevich4,5,6 and Albert Rimola2 1 Institut de Ciències de l’Espai (ICE, CSIC), Campus UAB, C/ de Can Magrans s/n, 08193 Cerdanyola del Vallès, Catalonia, Spain 2 Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain 3 Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Catalonia, Spain 4 Finnish Geospatial Research Institute (FGI), Geodeetinrinne 2, FI-02430 Masala, Finland 5 Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2a, P.O. Box 64, FI-00014 Helsinki, Finland 6 Institute of Physics and Technology, Ural Federal University, Mira str. 19. 620002 Ekaterinburg, Russia arXiv:2103.13758v2 [astro-ph.EP] 8 Jun 2021 Delivered 2020 November 26. ABSTRACT The disruption of asteroids and comets produces cm-sized meteoroids that end up impacting the Earth’s atmosphere and producing bright fireballs that might have associated shock waves or, in geometrically-favorable occasions excavate craters that put them into unexpected hazardous scenarios. The astrometric reduction of meteors and fireballs to infer their atmospheric trajectories and heliocentric orbits involves a complex and tedious process that generally requires many manual tasks. To streamline the process, we present a software package called SPMN 3D Fireball Trajectory and Orbit Calculator (3D-FireTOC), an automatic Python code for detection, trajectory reconstruction of meteors, and heliocentric orbit computation from video recordings. The automatic 3D-FireTOC package comprises of a user interface and a graphic engine that generates a realistic 3D representation model, which allows users to easily check the geometric consistency of the results and facilitates scientific content production for dissemination. The software automatically detects meteors from digital systems, completes the astrometric measurements, performs photometry, computes the meteor atmospheric trajectory, calculates the velocity curve, and obtains the radiant and the heliocentric orbit, all in all quantifying the error measurements in each step. The software applies corrections such as light aberration, refraction, zenith attraction, diurnal aberration and atmospheric extinction. It also characterizes the atmospheric flight and consequently determines fireball fates by using the 𝛼 − 𝛽 criterion that analyses the ability of a fireball to penetrate deep into the atmosphere and produce meteorites. We demonstrate the performance of the software by analyzing two bright fireballs recorded by the Spanish Fireball and Meteorite Network (SPMN). Key words: methods: analytical – methods: data analysis – Earth – meteorites, meteors, meteoroids – planets and satellites: atmospheres 1 INTRODUCTION better understanding of the impact hazard associated with large mete- oroids (Jenniskens 1998; Brown et al. 2002b; Trigo-Rodriguez et al. Meteor networks provide valuable scientific information about mm- 2007; Gritsevich et al. 2012; Trigo-Rodriguez & Williams 2017; San- to m-sized meteoroids crossing the Earth’s orbit thanks to the con- som et al. 2019; Moreno-Ibáñez et al. 2020; Moilanen et al. 2021). tinuous monitoring of the night sky (Ceplecha 1987). First meteor The recovery and the analyses of new meteorites and the study of their networks were based on classic photography, but after the first appli- dynamic association with comets, asteroids or planetary bodies give cation of CCD and video techniques to meteor observations (Trigo- new clues on the physical processes delivering space rocks to Earth Rodriguez et al. 2005; Trigo-Rodríguez et al. 2006b; Madiedo & (Whipple & Jacchia 1957; Trigo-Rodriguez et al. 2007; Jenniskens Trigo-Rodríguez 2007) great progress has been made. Nowadays, just & Vaubaillon 2008; Trigo-Rodríguez et al. 2009; Trigo-Rodríguez a few decades after this digital revolution, CCD and video cameras et al. 2015; Dmitriev et al. 2015). Moreover, the characterization of produce enough meteor recordings to provide an accurate depiction the atmospheric flight and the study of mechanical properties of the of bright fireballs. meteoroids contribute to impact hazard assessment (Trigo-Rodríguez Meteor detection provides information about the origin of mete- et al. 2006b; Tanbakouei et al. 2019). The analysis of cm- to m-sized oroids and about the continuous decay of asteroids and comets, their meteoroids ablating in the Earth’s atmosphere gives valuable clues main parent bodies (Murad & Williams 2002). By studying the helio- on the delivery of volatiles to Earth by using meteor spectroscopy centric orbits of meteoroids, identifying meteorite-dropping events, (Trigo-Rodriguez et al. 2003; Trigo-Rodríguez et al. 2004; Trigo- and developing the skills to reconstruct their strewn fields, one gains a Rodríguez 2019), but also is relevant to test their ability to penetrate deep into the atmosphere and quantify the consequences of small ★ E-mail: eloy.pena@uab.cat © 2020 The Authors 2 E. Peña-Asensio et al. asteroids coming from similar sources for planetary defense (Brown sula. The main steps of the analyses are set out below: 1) Meteor et al. 2002a; Boslough & Crawford 2008; Silber et al. 2018). trace detection, 2) Star identification and photometry, 3) Pixel to Meteor monitoring differs from most other types of astronomical real-world transformation, 4) Atmospheric trajectory reconstruction, observations since these luminous events cannot be predicted either 5) Parametrisation of the atmospheric flight, and 6) Heliocentric orbit in time or space (Vinkovic & Gritsevich 2020). For this reason, it computation. is important to monitor the sky with full-time and maximum spa- Finally, as examples of reduction procedures two meteoric events tial coverage. That is the foremost goal of multiple stations systems, are discussed: the Taurid Fireball: SPMN251019B and the sporadic often referred to as a meteor network. Some detection networks are superbolide: SPMN160819. The first study case involves typical tuned to very bright meteors, called fireballs when they exceed the recordings from video monitoring stations, while the second case brightness of Venus or superbolides when they are brighter than the involves one recording station, an occasional picture of the trail and magnitude −16 (Trigo-Rodríguez et al. 2015). Over the years, me- data extracted from a US government sensor recording. Using these teor and fireball detection networks have been built in many parts data the software was able to reconstruct the meteoroid trajectories, of the world, for instance, the Harvard Meteor Project (Jacchia & computing their masses, luminosities and, from the computation of Whipple 1956), the European Fireball Network (Ceplecha 1957), the their radiant and initial velocities, obtaining the heliocentric orbits continental scale Desert Fireball Network (DFN) (Bland 2004), the they had prior to the impact. SPanish Meteor Network (SPMN) (Trigo-Rodriguez et al. 2005), the Southern Ontario All-Sky Meteor Network (Weryk et al. 2007), the Finnish Fireball Network (FFN) (Gritsevich et al. 2014), the French 2.1 Fireball trace detection Fireball Recovery and InterPlanetary Observation Network and Me- A key step in developing an automatic astrometry is ensuring the soft- teorite Network (FRIPON) (Colas et al. 2015), the Italian network ware’s ability to detect the meteors appearing into the field of view for meteors observations and trajectory studies (PRISMA) (Gardiol of the video detection systems. For this purpose, we used the open- et al. 2016) and the Global Fireball Observatory (Devillepoix et al. source CV2 OpenCV library (Bradski & Kaehler 2000). Computer 2020). vision techniques are applied to obtain the corresponding pixel co- Since 1999 the SPMN has been continuously monitoring the sky ordinates with the meteor moving from frame to frame. Only frames over the Iberian Peninsula by setting about 34 stations distributed below a fixed mean pixel value will be processed in the most common throughout Spain (Trigo-Rodriguez et al. 2005; Trigo-Rodríguez cases, but saturated frames can appear during the recording of the et al. 2006b). All the data used in this work were obtained by the brightest bolide flares. These sudden increases in the meteor mag- SPMN and the data were processed from the database created and nitude often saturate the images, so they cannot be properly treated. operated by the Meteorite, Minor Bodies and Planetary Sciences The first step in processing a video frame uses the typical method group at the Institute of Space Sciences (CSIC-IEEC). The network of Gaussian blur smoothing of the grayscale image to reduce noise stations consist of two operational systems: 1) All-sky CCD cam- (Wells 1986). This enhances image structure and reduces details by eras (180◦ ) with fish-eye lenses and detectors of 4096𝑥4096 pixels convolutions. Each frame is compared to a reference frame (the frame (Trigo-Rodriguez et al. 2005), and 2) wide-field video systems (90◦ to prior to the beginning of the event), so we can extract the pixel value 120◦ ) working at 25 frames per second (transformed into 50 frames difference above a threshold using the Absdiff function from the CV2 per second by deinterlacing) (Madiedo & Trigo-Rodríguez 2007). package. To improve the detection, some morphological operations For the first system, the entire sky can be recorded without interrup- are applied such as erosion (removing isolated pixels) and dilation tion and reaching stellar magnitude between +6 and +10, depending (expanding the pixel size). The next step is to contour the differenced on the zenith angle and the night sky background luminosity. In the pixels, as done by Suzuki & Abe (1985), to outline the halo of bright case of the second instrumentation, the typical configuration uses 3 pixels. These contours in subsequent frames can be used to determine cameras per station covering 120𝑥90 square degrees up to a limiting if the feature has moved. In each frame, the contour should be the magnitude of +4. meteoroid or its trail. If the feature is determined to be the meteoroid, the centroid of the contour can be used as the location of the me- teoroid (in pixel coordinates). Figure 1 shows a selection of frames from the processed event SPMN300319B, where fireball detections 2 AUTOMATIC METEOR DETECTION AND SOFTWARE are shown in chronological order from left to right, and their subse- ANALYTICAL PROCEDURES quent processing steps top to bottom (see also Table 1 for observer The astrometric reduction of fireballs involved hitherto a complex and information). The first shown frame precedes the appearance of the tedious processes that generally required several manual tasks. The meteor. The following frames demonstrate: detection of the meteor, future of fireball analysis is oriented towards complete automation a false positive due to glare, a rejected frame when brightness of the of the process, as demonstrated by the latest efforts (Colas et al. meteor has saturated the image, another detection of the meteor and 2020). To streamline the study of meteors, we developed a software the detection of the meteor trail. called 3D Fireball Trajectory and Orbit Calculator (3D-FireTOC), Due to the changing nature of fireball recordings, three different an automatic Python code for detection, trajectory reconstruction of methods to avoid false positives are implemented: 1) Discriminating meteors and heliocentric orbit computation from CCD recordings. by contour area size excluding excessively small and large contours, Thanks to the rendering engine VTK integrated into the MayaVI 2) Using the first detected points, a Kalman filter predicts an ex- package and using the NASA visible Earth catalogue for rendering the pected area for the next bolide position restricting the contours search surface (http://visibleearth.nasa.gov/), a realistic 3D model generator (Welch et al. 1997; Sansom et al. 2015) (function implemented in was implemented. Furthermore, we developed a friendly graphical CV2). However, if the first detections are not correct, the filter pro- user interface based on the toolkit Qt. duces wrong predictions, and 3) In parallel, all detected points are This software was developed in the framework of the analyses of saved (including those discarded with the Kalman filter). If the above Spanish Meteor and Fireball Network (SPMN) automatic recordings method does not give a result consistent with a nearly continuous of fireballs from optical stations distributed across the Iberian Penin- straight trajectory, after the detection process the clustering algo- MNRAS 000, 1–11 (2020) 3D-FireTOC 3 Figure 1. Frames from SPMN300319B fireball, an intermediate step in the processing and the result. Event recorded from Observatorio Astronómico Ramón Maria Aller (OARMA). rithm DBSCAN is applied to rule out incorrect points (Ester et al. ties. Because the atmosphere is thin compared to Earth’s radius, it 1996). The cluster associated with the meteor path will present a very may be treated as a plane parallel to the surface. This simplification low point dispersion, unlike obstacles or glare. See Figure 2. allows to easily apply Snell’s law and obtain a relation between the It is worth noting that SPMN control computers are synchronised true zenith distance 𝑧 and the apparent zenith distance 𝜁: either using GPS controlled systems or using a known software called NetTime that guarantees a minimum time accuracy of 0.1 𝑠, often sin 𝜁 cos (𝑧 − 𝜁) + cos 𝜁 sin (𝑧 − 𝜁) = 𝑛 sin 𝜁 . (1) slightly better 0.01 𝑠. The velocity is determined from the typical Assuming 𝑧 − 𝜁 is negligible and dividing by cos 𝜁: 1/25 𝑠 frame rate video frequency that could be improved to 1/50 𝑠 when deinterlacing the imagery. 𝑧 − 𝜁 = (𝑛 − 1) tan 𝜁 . (2) Although the refractive index at sea level may change with pres- sure and temperature, we assume that on average 𝑛 ≈ 1.0003. The 2.2 Star identification and photometry apparent position of the star can be expressed then as: To convert pixel coordinates into equatorial coordinates, it is first   sin 𝐴 tan 𝜁 necessary to identify reference stars with known declination and right 𝛼 0 = 𝐿𝑆𝑇 − arctan , ascension. To highlight the stars and reduce noise and spontaneous cos 𝜙 − sin 𝜙 cos 𝐴 tan 𝜁 (3) fireball flashes, all frames without detection are overlapped and a 𝛿 0 = arcsin (sin 𝜙 cos 𝜁 + cos 𝜙 sin 𝜁 cos 𝐴), logarithmic correction is applied to improve the process. For the where 𝐴 is true azimuth, 𝜙 is the observer’s latitude, 𝐿𝑆𝑇 is the local automated identification of the star coordinate on the image, the sidereal time, 𝛼 0 is the apparent right ascension, 𝛿 0 is the apparent Oriented FAST and Rotated BRIEF (ORB) descriptor is used (Rublee declination (Tatum 2019). et al. 2011). Once again, DBSCAN clustering algorithm is used: since The aberration of light is a phenomenon that occurs due to the stars appear in the sky far from each other in a random distribution, vector difference between the velocity of the Earth and the starlight’s the data labelled as noise by the algorithm will be the one of interest. velocity. This effect displaces the star towards the Earth’s apex and Figure 3 shows some of the most relevant steps of this process applied may be corrected by using Lorentz transformations: to the fireball event SPMN300319B. It is necessary to apply different corrections to the reference stars cos 𝜒 + 𝑣𝑐 that will be used as a comparison to properly estimate the fireball cos 𝜒 0 =  , 1 + 𝑣𝑐 cos 𝜒 magnitude. The star magnitude found in the catalogue is not the (4) sin 𝜒 same magnitude as observed from Earth. This is due to the different sin 𝜒 0 =  , physical phenomena produced by the atmosphere and the Earth’s 𝛾 1 + 𝑣𝑐 cos 𝜒 motion, which must be corrected to obtain proper results. √︁ where 𝛾 is the Lorentz factor 1/ 1 − (𝑣/𝑐) 2 , 𝜒 is the true apical The more air mass a star’s light passes through to the observer, distance, 𝜒 0 is the apparent apical distance, 𝑣 is the Earth’s speed and the more its brightness is reduced due to absorption and scattering 𝑐 is the speed of light. processes. We correct this atmospheric extinction using the table Assuming 𝑣/𝑐 << 1 and applying trigonometric operations it made by Green (1992), based on theoretical values for different at- follows that: mospheric conditions taking into account the observer altitude and  the zenith angle of the star. − 𝑣𝑐 sin 𝜓 csc 𝜒 The atmosphere of the Earth exhibits a non-uniform optical dis- 𝛼0 = + 𝛼, sin 𝜒 sin 𝜔 (5) tribution as a function of altitude. This means that starlights are 𝛿𝜒 refracted as their velocity changes from layers with different densi- 𝛿0 = (− cos 𝜓 sin 𝜒 + sin 𝜓 cos 𝜔 cos 𝜒) + 𝛿, cos 𝛿 MNRAS 000, 1–11 (2020) 4 E. Peña-Asensio et al. Figure 2. Clustering algorithm and statistical calculations for discarding false positives and automatically selecting the points corresponding to the meteor trail. From left to right: All detected points, clusters found and noise, and selected fireball cluster track. Applied to the SPMN300319B event recorded from OARMA. Figure 3. Sequence of the process of obtaining the coordinates of the stars in the image. Depicted temporarily from left to right. It shows the first frame of the video, the overlapping of all valid frames without detection, the application of the ORB algorithm after a logarithmic correction, the classification with the clustering algorithm, and the result. Applied to the SPMN300319B event recorded from OARMA. where 𝛼 is the true right ascension, 𝛿 is the true declination and 𝜔 is puted to convert the plate coordinates (𝑥, 𝑦) into equatorial ones the angle between the Earth’s apex, the star, and the north polar. (𝛼, 𝛿). This yields the apparent trajectory of the fireball from each Each star moves in the field of view with a relative speed depend- station. However, since the optical system introduces a lens distortion ing on its declination. This causes the stars near the poles to have and possible misalignments, there is a displacement of the stellar im- smaller angular velocities and, hence, to activate the pixels longer. ages from the plate centre. So, pixels cannot be converted directly to Consequently, as suggested by Rendtel (1993) it is appropriate to cor- equatorial coordinates. Therefore, the transformation needs interme- rect the magnitude of the stars to a reference declination as follows: diate steps. It is necessary to transform measured plate coordinates or pixel coordinates into standard coordinates (𝜉, 𝜂), that is, stere- 1 ographic projected coordinates or true coordinates, since: 1) The 𝑚(𝛿 𝑠 ) = 𝑚(𝛿0◦ ) − 2.5𝑝 log , (6) stereographic projection and the camera sensor are not necessarily cos 𝛿 𝑠 aligned, 2) The photographic objective provides a distorted image where 𝑝 is the Schwarzschild exponent (typically between 0.7 and of the celestial sphere, as it is the result of the sphere projection on 0.8), 𝛿 𝑠 is the star declination, 𝛿0◦ is the reference declination and 𝑚 the focal plane, and 3) The wide-field and fish-eye lenses typically is the apparent magnitude. produce pincushion or barrel distortion. Once the magnitudes of the reference stars are corrected, we per- After the first empirical proposal of an absolute astrometric model form an aperture photometry by counting the pixel intensity of both by Ceplecha (1987), some refinements and improvements to the pa- the stars and the fireball. In this way, a logarithmic regression can rameter estimation were suggested (Borovička 1992; Borovicka et al. be made to obtain the magnitude of the fireball. To standardize the 1995). However, these models present a high non-linearity, so they luminosity, we corrected its magnitude as if it had been observed at are hardly reversible and the convergence of estimation algorithms the zenith and calculated its absolute magnitude, i.e., its magnitude is not easily achieved. For this reason, new parametrisation based at 100 km distance: on polynomial representation was proposed (Bannister et al. 2013; ℎ Barghini et al. 2019; Jeanne et al. 2019). Following these latest works 𝑀 = 𝑚 − 5 log , (7) 100 and since SPMN’s stations are equipped with very diverse lenses, we where ℎ is the height and 𝑀 is the absolute magnitude. model, as a first approximation, the distortion due to the lens with a quadratic expression as suggested by Hawkes (1993), which can be expanded to higher orders if the number of reference stars allows it: 2.3 Pixel to real-world transformation A key process concerns the transformation from the digital chip sys- 𝜉 − 𝑥 = 𝑎𝑥 2 + ℎ𝑥𝑦 + 𝑏𝑦 2 + 𝑔𝑥 + 𝑓 𝑦 + 𝑐, (8) tem coordinates into the equatorial coordinates that we will carry out 𝜂 − 𝑦 = 𝑎 0 𝑥 2 + ℎ 0 𝑥𝑦 + 𝑏 0 𝑦 2 + 𝑔 0 𝑥 + 𝑓 0 𝑦 + 𝑐 0 , by comparing reference stars in the detecting field of view (FOV). Once the apparent positions of the reference stars are known, to- where 𝑎, ℎ, 𝑏, 𝑔, 𝑓 , 𝑎 0 , ℎ 0 , 𝑏 0 , 𝑔 0 , 𝑓 0 are parameter to fit and the gether with their pixel positions, a transformation matrix is com- plate constants. Finally, the transformation of standard coordinates MNRAS 000, 1–11 (2020) 3D-FireTOC 5 Figure 4. Graphical representation of the real meteor trajectory calculation by intersecting the planes and obtaining the radiant by projecting backwards Figure 5. A schematic diagram for radiant error computation. until the collision with the celestial sphere. The vertical projection is shown as well. where the negative root of 𝑡 𝑅 is chosen since it is the closest into equatorial coordinates (Steyaert 1990) is computed: point to (𝑋1 , 𝑌1 , 𝑍1 ). The cartesian coordinates of the radiant are   (𝑋 𝑅 , 𝑌𝑅 , 𝑍 𝑅 ). 𝜉 𝛼 = 𝐴 + arctan , The presence of the Earth’s gravity not only disturbs the velocity 𝜂 · sin 𝐷 − cos 𝐷 of the meteoroid but also modifies its velocity vector, having con- ! (9) 𝜂 · cos 𝐷 + sin 𝐷 sequences in the determination of its radiant in the sky (Dmitriev 𝛿 = arctan √︁ . et al. 2015). The method proposed by Andreev (1990) corrects this 𝜉 2 + (𝜂 · sin 𝐷 − cos 𝐷) 2 shift of the radiant towards the zenith, the so-called zenith attraction. being 𝛼 the right ascension, 𝛿 the declination and ( 𝐴, 𝐷) the Furthermore, the diurnal aberration has to be taken into account. unknown optical axis. Since the Earth rotates around its axis, the position of the radiant Since there is no analytical method to find the position of the moves away. The diurnal aberration is caused by the velocity of the optical axis, the simplex algorithm is used to find the solution for observation point on the rotating surface of the Earth. Therefore, it the system 8 and 9 that minimizes the mean squared error (Motzkin depends not only on the moment at which the observation is made, 1952). but also on the latitude and longitude of the observer as the Earth’s rotation around its axis moves the position of the radiant as well. It is corrected using the approximation suggested by Bellot-Rubio 2.4 Atmospheric trajectory reconstruction (1992). Numerous methods have been proposed for meteoroid triangula- Likewise, by performing geometric operations, the height of the tion (Ceplecha 1987; Borovicka 1990; Gural 2012), some of them meteoroid ℎ, the distance to each station, the length travelled and the very recent (Jansen-Sturgeon et al. 2020). We follow the Method of angle between the fireball trajectory and the local horizon 𝛾 can be Planes proposed by Ceplecha (1987), the average plane containing deduced. the apparent trajectory and the geographic coordinate from each ob- The calculation of errors consists of assuming the worst scenario servation point is obtained and then the stereoscopic intersection of from the simplex method uncertainties, that is to say, that each point the apparent trajectories calculated. See Figure 4. of the apparent trajectory arranges in the way that most deviates from We remark that the radiant is computed, as usual in related bibliog- the radiant. This will occur when the points are aligned crosswise raphy, by doing a backward projection of the atmospheric trajectory along the path, as shown in Figure 5: on the right of that figure, from each station until their encounter with the celestial sphere at the four possibilities of deviation assuming the worst case for each 𝑅∞ . Once the first meteor point is known (𝑋1 , 𝑌1 , 𝑍1 ) and using any right ascension and declination are depicted; on the left, the two other point contained on the trajectory (𝑋2 , 𝑌2 , 𝑍2 ), the radiant can largest possible deviations for each apparent trajectory are shown, be computed by a line-sphere intersection: which delimits the radiant error. The standard deviation assumed in the apparent trajectory comes from how our pixel to real-world 𝑎 𝑅 = (𝑋2 − 𝑋1 ) 2 + (𝑌2 − 𝑌1 ) 2 + (𝑍2 − 𝑍1 ) 2 , transformation matches the reference stars.  𝑏 𝑅 = 2 · (𝑋2 − 𝑋1 ) · 𝑋1 + (𝑌2 − 𝑌1 ) · 𝑌1 + (𝑍2 − 𝑍1 ) · 𝑍1 , (10) In a similar way, using the clone trajectories for the worst-case scenario and following the plane intersection method, we obtain the 𝑐𝑅 = 𝑋12 + 𝑌12 + 𝑍12 − 2 𝑅∞ average observed velocity in the first 10 % of the luminous trajectory being 𝑎 𝑅 , 𝑏 𝑅 and 𝑐 𝑅 the parameters of the resulting equation from for each of them, as suggested by Whipple & Jacchia (1957). Veloci- substituting the equation of the line into the sphere. ties are calculated using the closest and/or most reliable observation. This yields to the parametric line: Then, we perform a linear regression to estimate the pre-atmospheric √︃ velocity at an instant prior to the first detection, specifically, the time −𝑏 𝑅 − 𝑏 2𝑅 − 4 · 𝑎 𝑅 · 𝑐 𝑅 interval corresponding to one frame. In this way, we obtain the me- 𝑡𝑅 = , teoroid velocity at atmospheric impact and its associated maximum 2𝑎 𝑅 (11) error. Using this velocity, we apply the aforementioned diurnal aber- 𝑋 𝑅 = 𝑋1 + 𝑡 𝑅 · (𝑋2 − 𝑋1 ), ration and zenith attraction corrections propagating the error by de- 𝑌𝑅 = 𝑌1 + 𝑡 𝑅 · (𝑌2 − 𝑌1 ), riving the equations involved and taking into account the astrometric 𝑍 𝑅 = 𝑍1 + 𝑡 𝑅 · (𝑍2 − 𝑍1 ), errors. MNRAS 000, 1–11 (2020) 6 E. Peña-Asensio et al. 2.5 Parametrisation of the atmospheric flight 2009, 2012; Turchak & Gritsevich 2014; Sansom et al. 2019; Moreno- Ibáñez et al. 2020). In this regard, we implemented latest modifica- One of the most complex parts of the meteor reduction analysis is tion of the method proposed by Sansom et al. (2019) for determining to develop a mathematical model that properly describes the atmo- fireball fates using 𝛼 − 𝛽 criterion. Physically meaningful parametri- spheric flight, allowing us to extract physical information. Following sation of the luminous flight allows the pre-atmospheric and final the classical dynamic third-order time dependent system for charac- mass (corresponding to the terminal height) to be computed: terizing meteor deceleration and assuming that the body does not suffer any kind of fragmentation (Hoppe 1937) developed the well- !3 known Single Body Theory (SBT). We adopted this SBT approach 1 𝑐 𝑑 𝐴𝑒 𝜌 0 ℎ 0 𝑀𝑒 = , which treats the heat exchange and drag coefficients as constants 2 𝛼𝜌 2/3 sin 𝛾 𝑚 (16) along the luminous flight.   2 !! 𝛽 𝑉 By introducing convenient dimensionless quantities, the trajectory 𝑀 𝑓 = 𝑀𝑒 exp − 1− , 1−𝜇 𝑉𝑒 equations acquire the form of (Gritsevich 2009): where 𝜌 𝑚 meteoroid bulk density and 𝐴𝑒 is the pre-atmospheric 𝑑𝑣 1 𝜌0 ℎ0 𝑆 𝑒 𝜌𝑣𝑠 shape factor (usually ranges between 1.21, for an ideal sphere, and 𝑚 = 𝑐 , 𝑑𝑦 2 𝑑 𝑀𝑒 sin 𝛾 1.8) (Trigo-Rodríguez et al. 2015; Meier et al. 2017; Gritsevich et al. (12) 𝑑𝑚 1 𝜌0 ℎ0 𝑆 𝑒 𝑉𝑒2 𝜌𝑣 2 𝑠 2017; Lyytinen & Gritsevich 2016). = 𝑐ℎ , From the initial mass and approximating the shape of the mete- 𝑑𝑦 2 𝑀𝑒 𝐻 ∗ sin 𝛾 oroid to a sphere, the initial size can be estimated. This value can be the scale factor ℎ0 = 7.16 𝑘𝑚, sublimation heat 𝐻 ∗ , atmospheric contrasted with the calculation of the diameter from the radiated en- density near the sea level 𝜌0 = 1.29 · 103 𝑔/𝑐𝑚 3 , dimensionless ergy. Assuming that the kinetic energy value is the registered impact mass 𝑚 = 𝑀/𝑀𝑒 , velocity 𝑣 = 𝑉/𝑉𝑒 , air density term 𝜌 = 𝜌 𝑎 /𝜌0 energy 𝑇𝐸 , the equivalent meteoroid diameter 𝐷 is computed as: and cross-sectional area 𝑠 = 𝑆/𝑆 𝑒 . The subscript "𝑒" indicates the √︄ parameters at the entry to the atmosphere. 3𝑇𝐸 To find an analytical solution, it is assumed the isothermal atmo- 𝐷=23 , (17) 2𝜋𝜌𝑣 2 spheric model 𝜌 = 𝑒 −𝑦 and according to Levin (1956) the body mass and its middle section are connected by introducing the shape change where 𝜌 is the meteoroid bulk density and 𝑣 the velocity of the coefficient 𝑠 = 𝑚 𝜇 . The dimensionless parameter 𝜇 is treated as a meteor. constant and can be inferred in each case by studying the meteor light It is worth noting that since observed velocities have marked in- curve (Gritsevich & Koschny 2011; Bouquet et al. 2014; Drolshagen accuracies, data must be pre-processed before it can be used to fit et al. 2020). The first integrals for the system (Eq. 12) was proposed the parameters and compute velocities and deceleration. Experience by Stulov et al. (1995) with the initial conditions 𝑚 = 1, 𝑣 = 1, says that an optimal way to approximate these velocities is to adjust 𝑦 = ∞: the distances with a least-squares to the following equation:   𝛽 𝑚 = 𝑒𝑥 𝑝 − (1 − 𝑣 2 ) , 𝐿 = 𝑎 + 𝑏𝑡 + 𝑐𝑒 𝑘𝑡 , (18) 1−𝜇 (13) Δ 𝑦 = ln 𝛼 + 𝛽 − ln , where 𝐿 is the path length, 𝑡 is the time and 𝑎, 𝑏, 𝑐 and 𝑘 are 2 variables to be determined in the curve fitting (Whipple & Jacchia where: 1957; McCrosky & Posen 1968). Once the adjustment is made, by Δ = 𝐸𝑖(𝛽) − 𝐸𝑖(𝛽𝑣 2 ), deriving the previous expression velocities and decelerations are obtained in a trivial way. ∫ 𝑥 𝑡 𝑒 𝑑𝑡 (14) 𝐸𝑖(𝑥) = 𝑑𝑥, −∞ 𝑡 showing that the trajectory depends on two dimensionless parame- ters: 2.6 Heliocentric orbit computation 1 𝜌0 ℎ0 𝑆 𝑒 The last step to know the origin of the meteoroid in the Solar Sys- 𝛼= 𝑐 , 2 𝑑 𝑀𝑒 sin 𝛾 tem is to reconstruct its heliocentric orbit. Once the radiant has been (15) obtained and corrected and the atmospheric flight velocity curve has 𝑐 ℎ 𝑉𝑒2 𝛽 = (1 − 𝜇) . been computed, the orbital elements that define the meteor’s orbit 2𝑐 𝑑 𝐻 ∗ can be calculated. Following the steps of Ceplecha (1987) and Jen- where, in this section only, 𝛼 symbol refers to the ballistic coefficient niskens & de Lignie (1987), first the coordinates of the geocentric and 𝛽 to the mass loss parameter. radiant (𝛼𝐺 , 𝛿𝐺 ) are transformed into ecliptical longitude and lati- The parameter 𝛼 characterizes the aerobraking efficiency since tude (𝐿 𝐺 , 𝐵𝐺 ). Thus, the heliocentric ecliptic system of rectangular it is proportional to the mass of a trajectory-aligned atmospheric coordinates can be defined as: column of cross-section divided by the body mass. The parameter 𝛽 is proportional to the fraction of the kinetic energy supplied to a 𝑋 = 𝑟 · cos 𝐿 · cos 𝐵, unit mass of the body as heat divided by the effective destruction 𝑌 = 𝑟 · sin 𝐿 · cos 𝐵, (19) enthalpy. 𝑍 = 𝑟 · sin 𝐵, These parameters bring great simplicity to the characterization of the atmospheric flight and can also be used to estimate how likely where 𝑟 is the distance to the Sun. a fireball produces meteorites (Gritsevich 2008a,b; Gritsevich et al. Then the ecliptical longitude of the Earth’s Apex 𝐿 𝐴𝑃 and the MNRAS 000, 1–11 (2020) 3D-FireTOC 7 Earth’s velocity 𝑉 𝐴𝑃 are extracted from 𝐽𝑃𝐿 𝐻𝑜𝑟𝑖𝑧𝑜𝑛𝑠 ephemerides. intermediate part, which was only filmed from the Eivissa station. The heliocentric velocities of the meteoroid can be expressed as: Despite this, from the astrometric measurements of the video frames and the integration of the data we achieved the trajectory reconstruc- 𝑣 𝐻𝑥 = −𝑣 𝐺 · cos 𝐿 𝐺 · cos 𝐵𝐺 + 𝑉 𝐴𝑃 cos 𝐿 𝐴𝑃 , tion. The fireball light was first detected at a height of 79.0 ± 0.1 𝑘𝑚 and the end occurred at 58.3 ± 0.1 𝑘𝑚 having a trajectory angle of 𝑣 𝐻𝑦 = −𝑣 𝐺 · sin 𝐿 𝐺 · cos 𝐵𝐺 + 𝑉 𝐴𝑃 sin 𝐿 𝐴𝑃 , (20) 𝛾 = 28.7 ◦ , which indicates the very remote possibility of being a 𝑣 𝐻𝑧 = −𝑣 𝐺 · sin 𝐵𝐺 . meteorite-dropper since its terminal height was too high (Moreno- Ibáñez et al. 2015). Following the photometry procedure described The specific angular momentum ℎ¯ = (ℎ 𝑥 , ℎ 𝑦 , ℎ 𝑥 ) and the ascend- in Section 2.2, we obtained a magnitude of −13.5 ± 0.5, as bright as ing node vector 𝑛ˆ are needed to determine the parameters, which can the full moon. be computed as: The pre-atmospheric velocity was retrieved from the velocity mea- ℎ¯ = 𝑟¯ × 𝑣, sured at the earliest part of the fireball trajectory by doing a regres- (21) sion and extrapolating with a backpropagation. It was estimated to be 𝑛¯ = 𝑧ˆ × ℎ¯ = (−ℎ 𝑦 , ℎ 𝑥 , 0), 28.0 ± 0.2 𝑘𝑚/𝑠. Assuming a shape change coefficient of 𝜇 = 2/3, where 𝑥,ˆ 𝑦ˆ and 𝑧ˆ are the unit axes of the heliocentric coordinate a shape factor of 𝐴𝑒 = 1.3, a drag coefficient of 𝑐 𝑑 = 1.3 and a system. relatively low density of 𝜌 = 1.6 𝑔/𝑐𝑚 3 (Harmon & Nolan 2005), Finally, each of the orbital elements can be obtained: the initial and final mass were computed using the method detailed in ℎ𝑧 Section 2.5. Using the D-criterion of Southworth & Hawkins (1963) cos(𝑖) = , we obtained a value of 𝐷 𝑆𝐻 = 0.35. This indicates that the orbit of ℎ ℎ𝑦 the event SPMN251019B is suggestive of being dynamically associ- cos(Ω) = √︃ , ated with the established Southern Taurid shower (Jenniskens et al. ℎ2𝑥 + ℎ2𝑦 2016). −ℎ 𝑦 𝑒 𝑥 + ℎ 𝑥 𝑒 𝑦 The calculated radiant and the velocity are shown in Table 2, cos(𝜔) = √︃ , together with the computed orbital parameters and the main fireball 𝑒 ℎ2𝑥 + ℎ2𝑦 (22) parameters. Figure 6 shows the summed frames of the recordings 1  𝐺 𝑀   and the graphic representation of the apparent trajectories over the 𝑒¯ = 𝑣2 − 𝑟¯ − (𝑟¯ · 𝑣¯ ) 𝑣¯ , celestial sphere and the atmospheric flight in real scale. In addition, 𝐺𝑀 𝑟 Figure 7 depicts the orbit of the SPMN160819B progenitor associated 1 𝑎= , to comet 2P/Encke. 2 − 𝑣2 𝑟 𝐺𝑀 Bright fireballs recorded in October are often belonging to one 𝑒¯ · 𝑟¯ of the Taurid streams (Northern or Southern branches). The Taurids cos 𝑣 0 = , 𝑒𝑟 exemplify impact hazard associated with large meteoroids due to the where 𝑖 is the inclination, Ω the longitude of the ascending node, 𝜔 frequency and size of the bodies reaching the Earth’s atmosphere. the argument of perihelion, 𝑒 the eccentricity, 𝑎 the semimajor axis The entire Taurid complex consists of Near Earth Objects (NEOs), and 𝑣 0 the true anomaly (Dubiago 1961). plus several meteoroid streams. The complex itself is considered a potential source of risk related to possible impacts by cosmic objects. In fact, it was proposed that the Tunguska event was produced by an asteroid-size body associated with the Taurid complex (Sekanina 3 STUDY CASES 1998). Several studies have demonstrated the dynamic association The software was successfully applied to study different events as between the Taurid complex and the disruption of a much larger test cases. We tested earlier the computer vision system with the 2P/Encke progenitor comet (see e.g. (Kresak 1978)). SPMN300319B case as it presented notable complications such as obstacles and frame saturation. We chose two other events to exem- plify a complete reduction (see Table 1). 3.2 Sporadic superbolide: SPMN160819 The SPMN251019B fireball is a typical reduction case thanks to the favorable astrometry made based on the recordings from the three On August 16, 2019, a very bright superbolide catalogued as stations, which we propose as belonging to the Taurids complex. The SPMN160819 event occurred (see Table 1). It was an event of con- other studied case is the superbolide SPMN160819 that demonstrates siderable importance due to its magnitude that, unfortunately, was the ability to combine satellite data and video recordings. These only partially recorded from the Eivissa station of the SPMN network events are listed in Table 1 with their corresponding observation (Peña-Asensio et al. 2020b). However, thanks to citizen collaboration, data. we had access to two more records: an image from Costa Brava and a video from Sardinia, which were used in the superbolide analysis. Since casual records of extremely rare events have limited resolu- 3.1 Taurid Fireball: SPMN251019B tion, we had to use the peak brightness coordinates measured by the The first example is the bolide SPMN251019 that occurred on Oc- Center for Near Earth Object Studies (CNEOS) at NASA to perform tober 25th, 2019 at 04:36:46 UTC (Peña-Asensio et al. 2020a). The a correct reduction of this event. event was videotaped by three SPMN monitoring stations: Astro- From the recording from Eivissa, in which the Moon appears at a nomical observatory at Puig des Molins (Eivissa), Montseny Astro- similar altitude, the superbolide was more luminous than the Moon. nomical Observatory (Barcelona) and Folgueroles (Barcelona). The It was estimated to exhibit an absolute magnitude of −16.5 ± 0.5. station coordinates are listed in Table 1. The superbolide from Eivissa was so distant that the first detected One of the complications of this case is that the two closest sta- light was at a height of 67 ± 3 𝑘𝑚 and ended at 23 ± 2 𝑘𝑚. The result tions recorded the beginning and the end of the fireball, but not the for the pre-atmospheric velocity was 15 ± 1 𝑘𝑚/𝑠 and the terminal MNRAS 000, 1–11 (2020) 8 E. Peña-Asensio et al. Table 1. Table with the different events recorded by the SPMN network. The SPMN300318 event is from single station since it was used to illustrate the operation of the code. ∗ The observation point does not belong to the SPMN network. Name Stations Longitude Latitude Altitude Date Start Time (UTC) End Time (UTC) SPMN300319B OARMA 08◦ 33’19”W 42◦ 52’33”N 236 m 2019/03/30 19h46m30.4s 19h46m34.4s SPMN251019B Eivissa 01◦ 25’45”E 38◦ 54’21”N 45 m 2019/10/25 04h36m48.4s 04h36m50.4s Folgueroles 02◦ 19’33”E 41◦ 56’31”N 580 m 04h36m49.976s 04h36m50.657s Montseny 02◦ 32’01”E 41◦ 43’47”N 194 m 04h36m46.279s 04h36m48.310s SPMN160819 Eivissa 01◦ 25’45”E 38◦ 54’21”N 45 m 2019/08/16 20h36m01.3s 20h36m05.6s Costa Brava∗ 03◦ 04’10”E 41◦ 49’03”N 2m 20h36m04s 20h36m04s Sardinia∗ 08◦ 31’43”E 39◦ 54’37”N 30 m 20h36m01s 20h36m06s Radiant Data Observed Geocentric Heliocentric 𝛼 (◦ ) 42.7±0.2 40.5±0.2 346.7±0.4 𝛿 (◦ ) 11.3±0.1 9.5±0.2 -4.2±0.3 V (km/s) 28.0±0.2 26.0±0.2 36.5±0.2 Orbital Parameter a (AU) e q 𝜔 (◦ ) 𝜔 (◦ ) i (◦ ) 1.97±0.07 0.792±0.007 0.410±0.006 109.2±0.9 31.199±10−4 6.0±0.4 SPMN251019B Mag ℎ𝑖 (𝑘𝑚) ℎ 𝑓 (𝑘𝑚) 𝑉𝑖 (𝑘𝑚/𝑠) 𝑉 𝑓 (𝑘𝑚/𝑠) 𝑀𝑖 (𝑔) 𝑀 𝑓 (𝑔) -13.5±0.5 80.0±0.1 58.3±0.1 28.0±0.2 17.59±0.2 43.1 0.003 Table 2. Top: SPMN251019B observed, geocentric and heliocentric radiant and velocities. Middle: SPMN251019B calculated orbital parameters. Bot- tom: SPMN251019B computed atmospheric trajectory, velocity and mass. Figure 7. An orbital projection to the ecliptic plane of different solar system bodies is shown, indicating the Vernal equinox to the right. The Sun is shown in yellow, Earth’s orbit in blue, Mars’s orbit in red, the Main Belt Asteroids in gray, Jupiter in green, the orbit of the SPMN251019B meteoroid in pink, comet 2P/Encke in dashed pink and the orbit of the SPMN160819 meteoroid in orange colour. velocity 11±1 𝑘𝑚/𝑠. The pre-atmospheric velocity is consistent with that recorded by CNEOS (14.9 𝑘𝑚/𝑠). Figure 8 shows the stacked frames of the recordings and the graphic representation of the apparent trajectories over the celestial sphere and the atmospheric flight in real scale. Figure 7 shows the orbit of the SPMN160819 progenitor. The slope between the trajectory and the local horizon is one of the key parameters that define the fate of the meteoroid as a consequence of the ablation. In this case, the trajectory slope was estimated to be 49 ◦ . After performing the fitting of the normalised velocity and the normalised height in order to parametrize the atmospheric flight (see Figure 9) and assuming a mean value of ordinary chondrite’s density of 2.7 𝑔/𝑐𝑚 3 (Consolmagno & Britt 1998; Blum et al. 2006), a shape change coefficient of 𝜇 = 2/3, a shape factor of 𝐴𝑒 = 1.3 and a drag coefficient of 𝑐 𝑑 = 1.3 the masses are calculated from Eq. 16. The initial mass of the meteoroid was estimated to be 2100 𝑘𝑔 Figure 6. Top: SPMN251019B apparent trajectory recorded and reduced corresponding to the initial size of 1.2 𝑚 and the terminal mass com- from Eivissa (orange), Folgueroles (red) and Montseny (green). Bottom: puted is 190 𝑘𝑔. Introducing the radiated energy peak recorded by SPMN251019B atmospheric trajectory with vertical projection (white). CNEOS, 𝑇𝐸 = 0.089 𝑘𝑡, on Eq. 17 gives a diameter of 1.3 𝑚, which is in good agreement with our results. This emitted energy could be compared with the Villalbeto de la Peña superbolide videotaped on MNRAS 000, 1–11 (2020) 3D-FireTOC 9 Radiant Data Observed Geocentric Heliocentric 𝛼 (◦ ) 228.2±1.5 204.0±1.4 226.1±0.4 𝛿 (◦ ) 68.0±0.2 67.6±0.4 16.9±0.4 V (km/s) 15±1 10±1.5 31±0.6 Orbital Parameter a (AU) e q 𝜔 (◦ ) 𝜔 (◦ ) i (◦ ) 1.15±0.06 0.17±0.04 0.953±0.008 126±9 143.43±10−4 17±3 SPMN160819 Mag ℎ𝑖 (𝑘𝑚) ℎ 𝑓 (𝑘𝑚) 𝑉𝑖 (𝑘𝑚/𝑠) 𝑉 𝑓 (𝑘𝑚/𝑠) 𝑀𝑖 (𝑘𝑔) 𝑀 𝑓 (𝑘𝑔) -16.5±0.5 67±3 23±3 15.1±1 11±1 2100 190 Table 3. Top: SPMN160819 observed, geocentric and heliocentric radiant and velocities. Middle: SPMN160819 calculated orbital parameters. Bottom: SPMN160819 computed atmospheric trajectory, velocity and mass. Figure 9. Top: Plot of observational data with velocity normalized to the entry velocity and height normalized to the atmospheric scale height for the SPMN160819 event. Bottom: SPMN160819 flight parametrization. The bounding line for a 50 𝑔 meteorite is shown in black for the case where there is no spin (𝜇 = 0) and in grey where spin allows uniform ablation over the entire surface (𝜇 = 2/3). Parameters 𝛼 and 𝛽 come from Eq. 15. and its ability to penetrate the atmosphere and pose a potential haz- ard. Thanks to the application of these new techniques, the analysis Figure 8. Top: SPMN160819 apparent trajectory recorded and reduced from of the atmospheric deceleration of cm to m-sized bodies penetrating SPMN Eivissa station (Orange), Sardinia (Green) and Costa Brava (pur- ple). Bottom: SPMN160819 atmospheric trajectory with vertical projection the atmosphere at hypervelocity is facilitated. In summary, the main (white). The red dot corresponds to the point of highest radiated energy as conclusions of this work are: registered by CNEOS. • An automatized processing system to perform detection and astrometric reduction of meteor video recordings is presented. The Jan. 4th, 2004 that produced a blast with a kinetic energy of about software uses state-of-the-art vision techniques, image processing 0.09 𝑘𝑡 (Llorca et al. 2005; Trigo-Rodríguez et al. 2006a). The 𝛼 − 𝛽 and motion detection methods to achieve fast astrometry, and a reli- criterion shows that this event was likely to produce meteorites, as it able calculation of astrometric errors. is depicted in Figure 9. The results are shown in Table 3. • New reduction techniques allows for avoiding false positives associated with bright flares experienced during the ablation process, as the Kalman filter is implemented to predict the motion of the object in the image. In addition, to discard incorrect points, a post- 4 CONCLUSIONS processing treatment was developed using clustering algorithms. A software tool for the detection and reconstruction of meteor trajec- • A corner algorithm is applied to automatically identify refer- tories was developed. The entire reduction procedure is automatic, ence stars. Subsequent treatment is performed to avoid possible false this way increasing our capacity to quantify, almost in real time, the positives due to presence of other objects, such as trees or buildings. meteoroid properties, the fireball trajectory, the heliocentric orbit, Also, corrections of the atmospheric extinction and refraction, as well MNRAS 000, 1–11 (2020) 10 E. Peña-Asensio et al. as the light aberration due to the Earth’s motion, was implemented Borovicka J., 1990, Bulletin of the Astronomical Institutes of Czechoslovakia, to improve the photometry. 41, 391 • We implemented a model to approximate the distortion of the Borovička J., 1992, PAICz, 79 lenses produced by wide-field and all-sky cameras using quadratic Borovicka J., Spurny P., Keclikova J., 1995, Astronomy and Astrophysics expressions and the simplex algorithm. The software characterizes Supplement Series, 112, 173 Boslough M., Crawford D., 2008, International Journal of Impact Engineer- the meteor flight and computes the pre-atmospheric mass. It identifies ing, 35, 1441 if a bolide is a meteorite-dropper using the 𝛼 − 𝛽 criterion. Bouquet A., Baratoux D., Vaubaillon J., Gritsevich M. I., Mimoun D., Mousis • A realistic atmospheric trajectory model in 3D was developed. O., Bouley S., 2014, Planetary and Space Science, 103, 238 Astrometric errors are propagated to infer the uncertainty in the Bradski G., Kaehler A., 2000, Dr. Dobb’s Journal of Software Tools, 25, 120 determination of the heights, radiant and inferred velocity along the Brown P. G., Revelle D. O., Tagliaferri E., Hildebrand A. R., 2002a, Mete- luminous path of the fireballs. oritics & Planetary Science, 37, 661 • The performance of the software is demostrated by computing Brown P., Spalding R., ReVelle D. O., Tagliaferri E., Worden S., 2002b, the heliocentric orbits of the two study cases. We found that the Nature, 420, 294 SPMN251019B fireball data are in good agreement with the values Ceplecha Z., 1957, Bulletin of the Astronomical Institutes of Czechoslovakia, associated with the Southern Taurid meteoroid stream. 8, 51 Ceplecha Z., 1987, Bulletin of the Astronomical Institutes of Czechoslovakia, • In reference to the superbolide SPMN160819, the obtained data 38, 222 show that it was produced by a m-sized sporadic meteoroid that, after Colas F., et al., 2015, in Proceedings International Meteor Conference. pp disruption in the atmosphere, might have produced meteorites. 34–38 Colas F., Zanda B., Bouley S., Jeanne S., Malgoyre A., Birlan M., et al., 2020, Submitted to A&A ACKNOWLEDGEMENTS Consolmagno S., Britt D. T., 1998, Meteoritics & Planetary Science, 33, 1231 Devillepoix H., et al., 2020, Planetary and Space Science, 191, 105036 JMT-R, EPA and AR acknowledge financial support from the Span- Dmitriev V., Lupovka V., Gritsevich M., 2015, Planetary and Space Science, ish Ministry (PGC2018-097374-B-I00 funded by MCI-AEI-FEDER, 117, 223 PI: JMT-R; CTQ2017-89132-P,PI: AR). MG acknowledges support Drolshagen E., et al., 2020, Luminous efficiency based on FRIPON meteors from the Academy of Finland project no. 325806, and the Russian (arXiv:2011.06805) Foundation for Basic Research, project nos. 18-08-00074 and 19- Dubiago A., 1961, New York 05-00028. AR is indebted to the “Ramón y Cajal” program and Ester M., Kriegel H.-P., Sander J., Xu X., et al., 1996, in Kdd. AAAI Press, DIUE (project 2017SGR1323). This project has received funding pp 226–231 Gardiol D., Cellino A., Di Martino M., 2016. In Proceedings of the Interna- from the European Research Council (ERC) under the European tional Meteor Conference, Egmond., pp 76–79 Union’s Horizon 2020 research and innovation programme (grant Green D. W. E., 1992, International Comet Quarterly, 14, 55 agreement No. 865657) for the project “Quantum Chemistry on In- Gritsevich M., 2008a, Solar System Research, 42, 372 terstellar Grains” (QUANTUMGRAIN). We thank Prof. José A. Do- Gritsevich M., 2008b, in Doklady Physics. p. 588 cobo and Dr. Pedro P. Campo for the video obtained to exemplify the Gritsevich M., 2009, Advances in Space Research, 44, 323 software (Fig. 1-3) recorded from Observatorio Astronómico Ramón Gritsevich M., Koschny D., 2011, Icarus, 212, 877 Maria Aller (OARMA), Universidad de Santiago de Compostela. The Gritsevich M. I., Stulov V. P., Turchak L. I., 2009, Doklady Physics, 54, 499 casual video recording from Sardinia was kindly provided by Claudio Gritsevich M., Stulov V., Turchak L., 2012, Cosmic Research, 50, 56 Porcu (Fig. 8). We also thank Quico Terradelles for provinding the Gritsevich M., et al., 2014, in Rault J. L., Roggemans P., eds, Proceedings persistent train of the SPMN160819 superbolide recorded from Platja of the International Meteor Conference, Giron, France, 18-21 September 2014. pp 162–169 d’Aro, Costa Brava (Fig. 8). We also thank preliminary feedback re- Gritsevich M., et al., 2017, in , Astrophysics and Space Science Proceedings. ceived from that superbolide by Peter Brown (UWO). The authors Springer International Publishing, pp 153–183, doi:10.1007/978-3-319- thank Dr. Eleanor K. Sansom for the detailed and constructive review 46179-3_8 that helped us to improve this paper. Gural P. S., 2012, Meteoritics & Planetary Science, 47, 1405 Harmon J., Nolan M., 2005, Icarus, 176, 175 Hawkes R. L., 1993, in Meteoroids and their Parent Bodies. p. 227 AVAILABILITY OF DATA Hoppe J., 1937, Astronomische Nachrichten, 262, 169 Jacchia L. G., Whipple F. L., 1956, Vistas in Astronomy, 2, 982 The data underlying this article will be shared on reasonable request Jansen-Sturgeon T., Sansom E. K., Devillepoix H. A., Bland P. A., Towner to the corresponding author. M. C., Howie R. M., Hartig B. A., 2020, The Astronomical Journal, 160, 190 Jeanne S., et al., 2019, Astronomy & Astrophysics, 627, A78 Jenniskens P., 1998, Earth, Planets and Space, 50, 555 REFERENCES Jenniskens P., Vaubaillon J., 2008, AJ., 136 Andreev G., 1990, in Proceedings of the International Meteor Conference, Jenniskens P., de Lignie M., 1987, J. DMS, 9 9th IMC, Violau, Germany, 1990. pp 25–27 Jenniskens P., et al., 2016, Icarus, 266, 331 Bannister S. M., Boucheron L. E., Voelz D. G., 2013, Publications of the Kresak L., 1978, Bulletin of the Astronomical Institutes of Czechoslovakia, Astronomical Society of the Pacific, 125, 1108 29, 129 Barghini D., Gardiol D., Carbognani A., Mancuso S., 2019, Astronomy & Levin B., 1956, Akad. Nauk SSSR, Moscow Astrophysics, 626, A105 Llorca J., et al., 2005, Meteoritics and Planetary Science, 40, 795 Bellot-Rubio L., 1992, Introducción a la Teoría Física de los Meteoros. Inter- Lyytinen E., Gritsevich M., 2016, Planetary and Space Science, 120, 35 science Publishers Madiedo J. M., Trigo-Rodríguez J. M., 2007, Earth, Moon, and Planets, 102, Bland P. A., 2004, Astronomy & Geophysics, 45, 5 133 Blum J., Schräpler R., Davidsson B. J., Trigo-Rodriguez J. M., 2006, The McCrosky R. E., Posen A., 1968, SAO Special Report, 273 Astrophysical Journal, 652, 1768 Meier M. M. M., Welten K. C., Riebe M. E. I., Caffee M. W., Gritsevich MNRAS 000, 1–11 (2020) 3D-FireTOC 11 M., Maden C., Busemann H., 2017, Meteoritics & Planetary Science, 52, APPENDIX A: CALCULATIONS 1561 Moilanen J., Gritsevich M., Lyytinen E., 2021, MNRAS, in revision Moreno-Ibáñez M., Gritsevich M., Trigo-Rodríguez J. M., 2015, Icarus, 250, This paper has been typeset from a TEX/LATEX file prepared by the author. 544 Moreno-Ibáñez M., Trigo-Rodríguez J. M., Gritsevich M., Silber E. A., 2020, Not. Roy. Astron. Soc., pp 316–324 Motzkin T., 1952, Project... scoop Murad E., Williams I. P., 2002, Meteors in the Earth’s Atmosphere Peña-Asensio E., et al., 2020a, in LPI Contribution, 39th Lunar and Planetary Science Conference- No., 1391. Peña-Asensio E., Trigo-Rodriguez Mas-Sanz E., Ribas J., 2020b. EuroPlanet Science Congress 2020. EPSC2020-459 Rendtel J., 1993, Handbook for visual meteor observers (Belgium: Interna- tional Meteor Organization) Rublee E., Rabaud V., Konolige K., Bradski G., 2011, in 2011 International Conference on Computer Vision. pp 2564–2571 Sansom E. K., Bland P., Paxman J., Towner M., 2015, Meteoritics & Planetary Science, 50, 1423 Sansom E. K., et al., 2019, The Astrophysical Journal, 885, 115 Sekanina Z., 1998, Planetary and Space Science, 46, 191 Silber E. A., Boslough M., Hocking W. K., Gritsevich M., Whitaker R. W., 2018, Advances in Space Research, 62, 489 Southworth R., Hawkins G., 1963, Smithsonian Contributions to Astro- physics, 7, 261 Steyaert C., 1990, Photographic Astrometry. International Meteor Organiza- tion Stulov V. P., et al., 1995, Moscow: Nauka, p. 236 Suzuki S., Abe K., 1985, Computer Vision, Graphics, and Image Processing, 30, 32 Tanbakouei S., Trigo-Rodríguez J. M., Sort J., Michel P., Blum J., Nakamura T., Williams I., 2019, A&A, 629, A119 Tatum J. B., 2019, Celestial Mechanics. Chapter 11: Photographic Astrome- try. Trigo-Rodriguez J. M., Williams I. P., 2017, Springer, pp 11–32 Trigo-Rodriguez J. M., Llorca J., Borovička J., Fabregat J., 2003, Meteoritics & Planetary Science, 38, 1283 Trigo-Rodríguez J. M., Llorca J., Fabregat J., 2004, Monthly Notices of the Royal Astronomical Society, 348, 802 Trigo-Rodriguez J., et al., 2005, Earth, Moon, and Planets, 95, 553 Trigo-Rodríguez J. M., Borovička J., Spurnỳ P., Ortiz J. L., Docobo J. A., Castro-Tirado A. J., Llorca J., 2006a, Meteoritics & Planetary Science, 41, 505 Trigo-Rodríguez J. M., Llorca J., Castro-Tirado A. J., Ortiz J. L., Docobo J. A., Fabregat J., 2006b, Astronomy & Geophysics, 47, 6.26 Trigo-Rodriguez J. M., Madiedo J. M., Llorca J., Gural P. S., Pujols P., Tezel T., 2007, MNRAS, 380, 126 Trigo-Rodríguez J. M., Madiedo J. M., Williams I. P., Castro-Tirado A. J., Llorca J., Vítek S., Jelínek M., 2009, Monthly Notices of the Royal Astronomical Society, 394, 569 Trigo-Rodríguez J. M., et al., 2015, Monthly Notices of the Royal Astronom- ical Society, 449, 2119 Trigo-Rodríguez J. M., 2019, in 2053-2563, Hypersonic Meteoroid Entry Physics. IOP Publishing, pp 4–1 to 4–23 Turchak L. I., Gritsevich M. I., 2014, Journal of Theoretical and Applied Mechanics, 44, 15 Vinkovic D., Gritsevich M., 2020, Journal of the Geographical Institute Jovan Cvijic, SASA, 70, 45 Welch G., Bishop G., et al., 1997, Chapel Hill, NC, USA Wells W. M., 1986, IEEE Transactions on Pattern Analysis and Machine Intelligence, pp 234–239 Weryk R. J., Brown P. G., Domokos A., Edwards W. N., Krzeminski Z., Nudds S. H., Welch D. L., 2007, Earth, Moon, and Planets, 102, 241 Whipple F. L., Jacchia L. G., 1957, Smithsonian Contributions to Astro- physics, 1, 183 MNRAS 000, 1–11 (2020) 12 E. Peña-Asensio et al. Table A1. SPMN251019B data reduction for the station Eivissa, Folgueroles and Montseny. SAO number, plate coordinates in pixels (𝑥, 𝑦), standard coordinates ( 𝜉 , 𝜂), right ascension and declination and their respective errors are shown. Station Ref SAO x (px) y (px) 𝜉 𝜂 RA (◦ ) DE (◦ ) err. RA (%) err. DE (%) Eivissa 1a 30631 457.2 497.85 0.09 0.54 268.302 57.01 0.009 0.009 2a 17074 538.1 390.8 0.31 0.73 245.92 61.57 0.018 0.053 3a 17365 493.6 390.2 0.18 0.78 257.125 65.769 0.067 0.06 4a 18222 409.2 404.0 -0.09 0.81 287.974 67.721 0.023 0.04 5a 8220 511.9 305.0 .35 1.13 230.092 71.859 0.032 0.032 6a 8102 504.7 282.2 0.37 1.27 222.576 74.174 0.017 0.023 7a 19019 333.85 432.5 -0.34 0.76 310.523 61.908 0.03 0.094 8a 34137 257.6 411.3 -0.65 0.9 332.736 58.253 0.049 0.01 9a 19302 311.65 414.0 -0.43 0.85 319.53 62.642 0.116 0.061 10a 20268 276.5 345.85 -0.62 1.23 342.449 66.23 0.009 0.012 Folgueroles 1b 60198 393.05 45.75 0.57 -0.1 113.477 31.89 0.034 0.076 2b 79666 405.2 91.35 0.56 -0.18 116.042 28.03 0.037 0.078 3b 95895 614.1 48.65 1.16 -0.35 99.092 15.758 0.004 <0.001 4b 115456 575.5 210.75 0.9 -0.64 111.678 8.299 0.029 0.024 5b 115756 574.55 258.85 0.85 -0.75 114.641 5.236 0.024 0.026 6b 61414 210.5 185.1 0.13 -0.14 140.163 34.395 0.012 0.063 7b 81064 201.95 302.25 0.02 -0.31 148.141 26.014 0.038 <0.001 8b 81004 226.85 313.9 0.04 -0.35 146.409 23.782 0.066 0.103 9b 98967 241.15 442.5 -0.06 -0.6 152.053 11.981 0.012 0.037 Montseny 1c 94027 386.1 234.85 0.15 -0.24 68.982 16.52 0.014 0.062 2c 77168 363.4 420.8 -0.06 -0.01 81.566 28.612 0.078 0.453 3c 39955 261.9 466.0 0.02 0.23 75.488 43.823 0.086 0.285 4c 40186 262.5 496.7 -0.02 0.28 79.165 45.996 0.02 0.136 5c 58636 363.8 510.75 -0.18 0.12 89.914 37.213 0.06 0.455 6c 40750 308.25 547.6 -0.16 0.28 89.864 44.945 0.021 0.101 7c 40756 299.9 553.75 -0.15 0.3 89.965 45.934 0.042 0.306 8c 39053 154.8 400.15 0.28 0.36 55.739 47.789 0.028 0.037 9c 38787 124.0 397.95 0.33 0.42 51.091 49.863 0.038 0.09 10c 23789 87.95 410.15 0.37 0.52 46.212 53.509 0.033 0.052 Table A2. SPMN160819 data reduction for the station Eivissa. SAO number, plate coordinates in pixels (𝑥, 𝑦), standard coordinates ( 𝜉 , 𝜂), right ascension and declination and their respective errors are shown. Ref SAO x (px) y (px) 𝜉 𝜂 RA (◦ ) DE (◦ ) err. RA (%) err. DE (%) 1 91781 339.0 424.8 -0.54 -0.32 3.309 15.184 0.416 0.541 2 108378 400.05 349.75 -0.26 -0.39 346.19 15.205 0.395 0.327 3 54471 175.55 413.15 -0.68 0.15 17.433 35.621 0.186 0.091 4 54058 228.2 406.25 -0.61 0.01 9.832 30.861 0.194 0.122 5 73765 266.35 379.4 -0.5 -0.06 2.097 29.09 0.464 0.101 6 90981 326.7 299.65 -0.21 -0.13 345.944 28.083 0.279 0.142 7 90816 357.5 293.5 -0.16 -0.2 342.501 24.602 0.18 0.04 8 90734 329.2 262.85 -0.11 -0.09 340.751 30.221 0.026 0.363 9 90238 399.05 229.1 0.06 -0.23 331.753 25.345 0.618 1.342 10 127029 492.05 240.0 0.15 -0.49 326.046 9.875 0.179 0.06 11 22268 62.45 272.6 -0.46 0.62 21.454 60.235 0.652 0.283 12 11482 79.75 252.75 -0.39 0.59 14.177 60.717 0.182 0.08 13 21609 107.7 267.8 -0.39 0.5 10.127 56.537 0.153 0.258 14 21133 116.65 230.6 -0.29 0.51 2.295 59.15 0.281 0.365 15 125122 557.9 48.2 0.75 -0.45 297.696 8.868 0.133 0.474 16 105223 547.6 36.1 0.77 -0.42 296.565 10.613 0.138 0.227 17 105500 482.4 38.6 0.67 -0.26 299.689 19.492 0.347 0.776 18 20268 118.45 142.55 -0.08 0.58 342.42 66.2 0.342 0.072 19 34137 168.25 129.25 0.01 0.46 332.714 58.201 0.162 1.578 20 19302 166.85 81.7 0.12 0.5 319.645 62.586 0.186 0.166 21 19019 178.55 53.35 0.2 0.49 311.322 61.839 0.291 0.582 MNRAS 000, 1–11 (2020)