A New Open-source Geomagnetosphere propagation tool (OTSO) and its applications

We present a new open-source tool for magnetospheric computations, that is modelling of cosmic ray propagation in the geomagnetosphere, named the “Oulu - Open-source geomagneToSphere prOpagation tool” (OTSO). A tool of this nature is required in order to interpret experiments and study phenomena within the cosmic ray research ﬁeld. Within this work OTSO is applied to the investigation several ground level enhancement events. Here, we demonstrated several applications of OTSO, namely computation of asymptotic directions of selected cosmic ray stations, eﬀective rigidity cut-oﬀ across the globe at various conditions within the design, general properties, including the magnetospheric models employed. A comparison and validation of OTSO with older widely used tools such as MAGNETOCOSMICS was performed and good agreement was achieved. An application of OTSO for providing the necessary background for analysis of two notable ground level enhancements is demonstrated and the their spectral and angular characteristics are presented.


Introduction
The Earth is under constant bombardment by high-energy particles known as cosmic rays (CRs).Primary CRs can have a solar, galactic or extra-galactic origin and are composed of ≈ 90% protons, 9% helium nuclei, and 1% heavier element nuclei (e.g.Gaisser et al., 2016, and references therein).CRs with a solar origin are produced during and following solar eruptions, such as solar flares and/or coronal mass ejections (CMEs) (e.g.Desai & Giacalone, 2016, and references therein), whilst CRs produced outside of the solar system are believed to come primarily from supernova remnants (e.g.Blasi, 2013, and references therein).CRs were first discovered in the early 20th century by Dr. Victor Hess and since then our knowledge of CRs has been constantly developing through the application of groundbreaking experiments, recent examples include the PAMELA (Payload for Antimatter Exploration and Light-nuclei Astrophysics) detector (Adriani et al., 2017), AMS-02 (Alpha Magnetic Spectrometer) (Aguilar et al., 2021) space-probes, and a plethora of ground-based experiments.
If a CR reaches the Earth's atmosphere it collides with atmospheric constituents.This collision produces numerous secondary particles, which then proceeds to collide or decay further into more secondary particles creating a cascade, the process developing until a threshold energy is reached.This phenomenon is known as an extensive air shower and is widely exploited by ground based detectors as a mechanism for study CRs.Within this work neutron monitors (NMs) are used as an example for CR detection, specifically of solar origin.NMs are fixed to a single location on the Earth making them especially good at studying CRs of solar origin, typically revealing anisotropy (e.g.Moraal & Mc-Cracken, 2012;Bütikofer et al., 2009).When NMs detect an increased flux of CRs with solar origin, it is dubbed a ground level enhancement (GLE), a relatively rare event that occurs only a few times per solar cycle (Shea & Smart, 2012).Introduced during the 1957 -1958 International Geophysical Year, NMs are standardised, nowadays assembled in a global network (Simpson, 2000;Mavromichalaki et al., 2011).
CRs are charged particles and as such experience the Lorentz force when travelling within a magnetic field.Thus, when a CR encounters the Earth's magnetic field there is the potential for the CR to be deflected by the magnetic field or penetrate it.The ability of the Earth's magnetic field to deflect certain CRs is known as magnetic shielding.
Whether a CR is able to penetrate the magnetosphere depends greatly on the CR's energy, the geomagnetic conditions at the time of arrival, the location of the CR's arrival, and its incidence.Only once a CR has penetrated the magnetosphere can it proceed to reach the Earth's atmosphere.An important characteristic of CRs is their rigidity, this value is typically used instead of CR energy as it is independent of the CR charge and species (Cooke et al., 1991).Rigidity quantifies the impact that magnetic fields have on the propagation of the CR, the larger the rigidity value the less deflected the particle is by a magnetic field.The rigidity of a CR is calculated using equation: where P is the rigidity, p is the CR's momentum, c is the speed of light, Z is the atomic number, and e is the elemental charge.Rigidity is important when considering CRs arriving at Earth as it can tell us which CRs are able to penetrate the magnetosphere and which are deflected away as a result of magnetic shielding.The rigidity needed by the CR to penetrate the magnetosphere ranges from 0 GV, at the magnetic poles, to ≈ 17 GV at the magnetic equator, this is due to the increase in magnetic shielding at lower latitudes (Gerontidou et al., 2021).Knowing the rigidity needed by a CR to arrive at different locations on the Earth (known as the cutoff rigidity) is important to analyse both ground-based and inside the geomagnetosphere space-borne experiments.
Determining an exact cutoff rigidity can be a difficult task due to the complex nature of CR propagation in the magnetosphere.It is typical to see a collection of CRs with sequential rigidities having a mixture of trajectories that can and can not penetrate the magnetosphere, referred to allowed and forbidden respectively (for details see Cooke et al., 1991).This region of allowed and forbidden CR trajectories is known as the penumbra.In order to get a useful quantitative value for the cutoff at given points on the Earth's surface an effective cutoff rigidity (R c ) is found, which accounts for the effects of the penumbra.As mentioned previously, the trajectories of CRs can be very complex, increasingly so at lower rigidities, this means that CRs that are detected at e.g.NMs do not necessarily arrive from directly above the station, therefore each station has its own asymptotic direction of acceptance for CRs, that is which part of the sky the detector is actually observing (Rao et al., 1963).
The complex nature of the Earth's magnetic field structure and CR propagation within it makes modelling CR trajectories very computationally intensive.The equations of motion that describe the trajectory of a charged particle in the Earth's magnetosphere currently have no known closed form solution.As such the trajectory of said particle must be determined using numerical integration (e.g.Bütikofer, 2018).It is almost impossible to predict where an arriving CR will encounter the Earth based on its point of arrival at the magnetosphere, therefore, the trajectory is typically computed backwards from just above the Earth's surface, around 20 km in altitude, to the CR's point of entry into the magnetosphere.
During GLEs, SEPs can have energies ranging from 10 MeV/nucleon up to about several GeV/nucleon (Biswas, 2000), relativistic effects should thus be accounted for in the model.To resolve any issues that can arise from this, the computation of the particle's trajectory must be done in small steps to avoid the model breaking, however this exacerbates the computational intensity of the modelling process.Finding a good balance between maintaining accuracy of the simulation and time efficiency is one of the main tasks of creating a magnetosphere computation tool.
The magnetosphere is a complex and dynamic environment, constantly changing in response to external conditions, which makes the accurate modelling very challenging.Empirical observations made by spacecraft have been used historically to create models describing the magnetic field structure (Jordan, 1994).As of present the field is best described as a combination of the inner magnetic field (created by the dynamo process in the Earth's core) and external magnetic fields (created by the various different currents within the magnetosphere).For the internal field, models such as the IGRF (Alken et al., 2021) and dipole models (Nevalainen et al., 2013) can be used and for external fields the Tsyganenko models are typically used (Tsyganenko, 1989(Tsyganenko, , 1995(Tsyganenko, , 1996(Tsyganenko, , 2002a(Tsyganenko, , 2002b)).
A tool that can compute the trajectories of charged particles within an accurate model of the magnetosphere under various conditions is highly valuable within the CR research field.The usefulness of such a tool has lead to the creation of multiple tools in the past (see (Bütikofer, 2018) and references therein), some examples are the developed by Smart and Shea (2001), COR by Gecášek et al. (2022), and MAGNETOCOSMICS by Desorgher (2006), the latter taken as a reference tool in this work, being widely used over the years.We emphasise that MAGNETOCOSMICS was designed within the framework of the Geant4 toolkit (Agostinelli et al., 2003), was released in 2006, and is in practice outdated nowadays.The only way to resolve this issue is to either update MAGNE-TOCOSMICS to be compatible with the newer Geant4 versions or create a new tool.Here we have chosen the latter approach and present the newly developed Oulu -Open-source geomagneToSphere prOpagation tool (OTSO), which can be tailored by the scientific community to meet the corresponding needs.

Formalism of CR Propagation
The trajectories of charged particles are influenced by the magnetic field generated by the Earth's core.This is due to the Lorentz force generated perpendicular to a charged particle moving through a magnetic field.The Lorentz force is described by: where F is the force [N], q is the charge [C], E and B are the electric [V m −1 ] and magnetic fields [T] respectively, and v is the particle's velocity [ms −1 ].
When considering magnetosphere calculations we can neglect E from the equation as its influence is negligible due to the high electrical conductivity of the region [for details see Bütikofer (2018) and the discussion therein].It is important to note that the bulk of CRs are travelling at relativistic speeds and as such the impact this has on the particle's mass must be considered when calculating the acceleration.This is achieved by incorporating the Lorentz factor γ: Combining equation 2 and the Lorentz factor within Newton's second law we can determine the equations of motion for a relativistic particle as a result of the Lorentz force within the magnetosphere in Cartesian coordinates: with m 0 being the rest mass of the particle.Knowing the acceleration of the CR at a specific point in the magnetosphere allows for the trajectory to be determined by preforming numerical integration of the equations of motion, as there is no solution in their enclosed form.There are multiple methods that can be implemented to achieve this, such as the Runge-Kutta, Euler, Boris, and Vay methods.All of these methods have their own benefits and drawbacks when considering computation speed and accuracy [see R.P. and Leigui (2019) and references therein].The most widely used method in CR simulations is the 4 th order Runge-Kutta method discussed by Smart et al. (2000), and as such it has been incorporated into the present work.This method offers a good balance between approximation accuracy and calculation time.Other integration methods can be added to OTSO later if required.
The trajectory of the particle is then determined up until the allowed or forbidden condition for the trajectory is met.Similarly to Smart and Shea (1981), we employ the approach of starting the particle propagation at 20km above the surface with the particle's velocity being directed vertically at the zenith, yet incidences of varying zenith and azimuth are also possible within OTSO.The trajectory is considered allowed if the particle is then able to reach the model magnetopause boundary and forbidden if it returns below the 20km starting altitude.We emphasise that 20km is selected as this is the typical altitude that atmospheric cascades start as a result of CR collisions with atmospheric particulates (Grieder, 2001).While collisions are possible above this point the model assumes the atmosphere is collisionless until 20km for simplicity.In addition, the tool is ending the simulation once the CR has travelled more that 100 Earth radii, to avoid endless simulation of a trapped particle that neither escapes or returns to Earth, and in this instance the trajectory is assumed forbidden.
OTSO can calculate individual trajectories of CRs with any given initial rigidity and start position on the Earth's surface.The user can select a range of rigidity values to test across as well as a rigidity step value, ∆R.OTSO will then repeat the trajectory calculation over all iterations of rigidity within the given range determining allowed and forbidden rigidity values.Due to the penumbra, that is encountered around the region in which rigidities change from allowed to forbidden, it is important to know the upper most accepted rigidity before the first forbidden value, R U , and the last allowed value, R L .These values are recorded during computations and to account for the effect of the penumbra an effective cutoff rigidity, R c , is calculated using a method described in Cooke et al. (1991).In which the sum of the number of allowed rigidities multiplied by ∆R is subtracted from R U , seen in equation 5.

Determining Time-step
One of the most important parts of the computation is determining the time-step (∆t) to use during the numerical integration of the equations of motion.If ∆t is too large then errors will accumulate leading to incorrect trajectories or the simulated particles accelerating to speeds faster than light, breaking the simulation.Vice versa, if ∆t is too small then the computation can take an irrationally long time to complete, making the tool impractical to use.
A convenient way to determine ∆t involves using the CR's properties and position in the magnetosphere.This work uses the method developed within Smart and Shea (1981) in which ∆t is the time taken for the particle to travel 1.0% of its gyration distance, making the assumption the magnetic field is uniform over the step.In order to optimise the computation further, the adaptive time-step method also outlined within Smart and Shea (1981) was utilised.This allows ∆t to grow by a maximum of 10% between Runge-Kutta iterations, only if the previous iteration was completed within an accepted error range, and sets the maximum value of ∆t to be 1.5% of the gyration time.This growth limit prevents any regions of sudden acceleration being skipped by large ∆t values.However, through testing of the new program the 1.5% limit was found to lead to extended computation times and was changed to 15 % with marginal impact on the results of the calculations.This limit can be edited or disabled depending on the accuracy of results required by the user and the desired computation time.The error between Runge-Kutta iterations is determined by checking the β value of the CR before and after the step, where β is the speed of the CR in units of c.As a charged particle in a magnetic field, with no other external force, should have a constant speed we can take the change in β to represent the error between steps.Within the tool this is implemented so that if β has grown by over 0.0001% during a Runge-Kutta step ∆t is assumed too large and the same iteration is repeated using ∆t 2 .This β error check is quite conservative and reduced error values have been able to reproduce similar results in a vastly reduced computation time, especially at higher latitude stations.The value of this error check parameter can be selected at the leisure of the user.

Employed Magnetosphere Models
In order to model the magnetic field OTSO uses an internal and external magnetic field model.There are only two main internal field models included in OTSO, these being the IGRF (Alken et al., 2021) and geodipole field models.The external component of the magnetic field is modelled using the Tsyganenko models: TSY87, TSY89, TSY96, TSY01, and TSY01S (Tsyganenko, 1987(Tsyganenko, , 1989(Tsyganenko, , 1995(Tsyganenko, , 1996(Tsyganenko, , 2002a(Tsyganenko, , 2002b;;N. Tsyganenko et al., 2003).We plan to include other models in the future, based on a convenient parameterisation.This allows for easier comparison with e.g.MAGNETOCOSMICS and/or other similar tools.The latter Tsyganenko models get increasingly complex and computationally intensive, leading to long simulation times.The use of later Tsyganenko models should be considered during periods of intense geomagnetic activity (e.g.periods of Kp index above 6).Generally a combination of the IGRF and TSY89 models is sufficient to provide fast and reliable results (e.g.Kudela & Usoskin, 2004;Nevalainen et al., 2013, and the discussion therein).Unless stated differently, the future calculations using OTSO will be performed using this combination of models.
In order to determine whether a CR has escaped the magnetosphere, this tool constantly checks the CR's position in relation to the model magnetopause chosen for the simulation.If the CR reaches the magnetopause boundary it is then assumed to have escaped.The TSY96, TSY01, and TSY01S models for the external magnetic field contribution use their own model magnetopause described within Tsyganenko (1995Tsyganenko ( , 1996Tsyganenko ( , 2002aTsyganenko ( , 2002b)).However, TSY89 has no such empirical magnetopause model used within it and therefore a "de-facto" boundary must be selected by the user which is then applied to the simulation when using the TSY89 model.
Models of the magnetopause have historically been produced using empirical methods, utilising data from satellite magnetopause crossings to best fit the shape.The models that have been included in the tool currently are a sphere with a 25 Earth radii radius centred around the Earth (for use when not considering any external field models) as well as the Formisano, Sibeck, and Kobel models (Formisano et al., 1979;Sibeck et al., 1991;Kobel, 1992;Flückiger & Kobel, 1990).Due to the differing assumptions made during the creation of these models the magnetopause shape can vary significantly between them, leading to slightly different simulation outcomes.When using TSY89 within this work the Kobel (1992) model has been used.
There are many more more advanced magnetopause models that take into account different variables, such as the solar wind conditions.Some examples of newer models are Lin et al. (2010) and Shue et al. (1998).While these models may provide a more accurate portrayal of the magnetopause they are not included in this tool at present.While these models may provide a more accurate portrayal of the magnetopause they are not included in this tool at present.These models require many more input variables to function, complicating and increasing the computational strain of the simulation.Their inclusion is planned to be accommodated in future versions of the tool if the need arises for the extra accuracy they may provide.

Programming Languages
OTSO has been developed within the framework of both the python and fortran programming languages.A precompiled language, such as fortran, was crucial in the development of this tool as the processing speed offered by compiled languages help complete the computationally intensive CR trajectory simulations within a reasonable time frame.Fortran is an old and dated language, with limited utility when compared to other more modern complied languages, such as C++.However, fortran does benefit from being an older language by being relatively simple as well as having many freely accessible and verified libraries previously written for it, which are already extensively used by the CR community, such as the Tsyganenko models, geopack library, and IRBEM library (https://prbem.github.io/IRBEM/).For these reasons fortran was chosen to utilise these libraries and speed up the development of this tool.
Python was also used as the way to initialise the tool and input the parameters.
Python is a very simple programming language and was picked to allow anyone with a basic understanding of programming to use the tool.The installation of the python language is also simple, being easily achieved through the download of the anaconda software.Anaconda also includes all the needed python modules needed for the tool to run, such as F2PY (which allows python and fortran to transfer information between each other).
The result of both these decisions is that the tool is simple to obtain, use, and edit.

Example Results, Comparison, and Applications
As this tool is being constructed to be a possible alternative to older programs, namely MAGNETOCOSMICS, the analysis of this new tool relies on the comparison of results between the two programs.To achieve this, several cases were selected, specifically related to GLE analysis, and both programs were used to conduct computations for said GLE(s).The first case is a well studied, within not complicated as magnetospheric conditions, derived characteristics GLE event, namely GLE # 70 (Vashenyuk et al., 2006;Bütikofer et al., 2009;Mishev & Usoskin, 2016).NM stations at various latitudes were used to test OTSO by conducting cutoff and asymptotic cones computations.All computations were done using a rigidity step of 1×10 −3 GV, manuscript submitted to JGR: Space Physics increasing the precision of the results and allowing the penumbra to be shown in greater detail, and employing combination of IGRF and TSY 89 models.

Cutoff Rigidity
The results for the vertical cutoff computations can be seen in Table 1, where a general good agreement between the OTSO and MAGNETOCOSMICS is found.The stations with the greatest difference between the two tools were Oulu and Tixie Bay.The slight variation in results can be attributed to the accuracy of the integration methods used within the two tools.

Asymptotic Cones
Once the trajectory of a CR has been simulated the asymptotic longitude and latitude are computed.The CRs with accepted trajectories then have these values plotted in order to construct the asymptotic cone of acceptance.Figure 1 shows the high rigidity value region of the cones created by OTSO and MAGNETOCOSMICS for three of the NM stations considered: Cape Schmidt, Oulu, and Rome respectively, encompassing the case of anti-sunward NM (CAPS), polar sunward (Oulu), and low latitude station (Rome).
One can see that the cones are in good agreement with each other, this is particularly true at the higher end of the rigidities.The cones calculated for the Oulu, and Rome NMs (see Figure 1) are almost identical with minor variations, however Cape Schmidt's cone shows that there can be some deviations in the cone shape between the two tools at lower rigidity values, namely the width of the OTSO cone increased (left panel of Fig- ure 1).This is because the trajectories of lower rigidity CRs are more complex, especially when their rigidity is close to the cutoff value, making simulation of these CRs more difficult.The accuracy of the integration method is important in these circumstances, and is likely the cause of the difference.

Global Cutoff
OTSO can compute the vertical cutoff rigidity on a global scale.Due to the mixed language nature of the tool multi-core processing was implemented using python to con- duct the large number of computations required for this operation in a time efficient manner.The global cutoff map was created by conducting cutoff calculations at regular intervals of 1 • in latitude and longitude.The same computation was done using MAGNE-TOCOSMICS and, in order to compare the two results, the absolute value for the difference between the computed cutoff rigidities at each point on the Earth was found and plotted in Figure 2.There are two clear results that can be inferred from Figure 2. Firstly, in general the difference between the two tools is minor, this is especially evident in the polar and equatorial regions.Secondly, there are anomalous regions on the Earth where the difference between the two tools is noticeable, with the most prominent region being found over the south pacific ocean, Figure 3

Ground level enhancement analysis
OTSO was employed for GLE analysis, namely for computations of the cut-off rigidity and asymptotic directions (ADs) for NMs used as input for the method deriving spectral and angular distribution of SEPs.Here, we used a method based on neutron monitor records analysis [e.g.Shea and Smart (1982); Cramp et al. (1997)], whose details and applications are given elsewhere (Mishev et al., 2018;Mishev, Koldobskiy, Usoskin, et al., 2021;Mishev et al., 2022).The method is an unfolding procedure, that is mod-  (Tsyganenko, 1989) or Tsyganenko 01 (Tsyganenko, 2002).The former combination allowed straightforward computation of ADs and rigidity cut-offs with reasonable precision [e.g.Kudela and Usoskin (2004); Kudela et al. (2008); Nevalainen et al. (2013)], whilst the latter is usually employed in the case when the K p index is greater than 6 [for details see Smart et al. (2000)].
Here we analysed two notable GLEs: GLE # 66 and GLE # 71.The GLE #66 occurred during one of strongest geomagnetic storms when the 3-hour planetary K p index was 9, on October 29, 2003.The event was the second in the sequence of three GLEs, the so-called Halloween events [e.g.Gopalswamy et al. (2005); Liu and Hayashi (2006); Gopalswamy et al. (2012)], recorded by the global neutron monitor network, the count rate increases are given in (http://gle.oulu.fi).In addition to the complicated geomagnetospheric conditions, a strong Forbush decrease was also observed prior to and during this event, which was explicitly considered in our analysis similarly to Mishev, Koldobskiy, Kocharov, and Usoskin (2021).Hence, the complicated geomagnetospheric conditions and accompanying Forbush decrease, make the analysis of this particular GLE specifically challenging.After computing the ADs with OTSO (see the left panel of Fig. 4), and employing the method described above, we derived the spectral and angular characteristics of the GLE producing SEPs.
The best-fit of the derived SEP spectra was obtained by a modified power-law rigidity spectrum [e.g.Vashenyuk et al. (2008); Mishev, Koldobskiy, Kocharov, and Usoskin (2021)], i.e.,: where the flux of particles with rigidity P in [GV] is along the axis of symmetry identified by geographic latitude Ψ, longitude Λ and the power-law exponent is γ with the steepening of δγ, J 0 is the particle flux at 1 GV in [m −2 s −1 sr −1 GV −1 ], for SEPs with rigidity P > 1 GV.Accordingly, for SEPs with P ≤ 1 GV, the rigidity spectrum is: For the angular distribution the best-fit was obtained by Gaussian-like distribution: where α is the pitch angle, σ accounts for the width of the distribution.
Details of derived spectra and pitch angle distribution (PAD) are given in Table 2 and 3 for the application of OTSO using IGRF+TSY 89 and IGRF+TSY 01 respectively.The merit function (Equation 9), that characterised the quality of the fit, that is the residual according to Himmelblau (1972); Dennis and Schnabel (1996) is defined as: Normally D ≤ 5 % for strong events [e.g.see Vashenyuk et al. (2006)], whilst for weak events it is ≈ 10 -15%, in some cases it can even approach 20% [for details see Mishev et al. (2018see Mishev et al. ( , 2022))].Here, the angular distribution of the arriving SEPs was fitted by complicated pitch an-

394
-13-manuscript submitted to JGR: Space Physics gle distribution (PAD) with a shape similar to that considered by Cramp et al. (1997), namely superposition of two Gaussians: where α is the pitch angle, σ 1 and σ 2 are parameters corresponding to the width of the pitch angle distribution, B and α ′ are parameters corresponding to the contribution of the second Gaussian, including the direction nearly opposite to the derived axis of symmetry.The best fit for the spectra was obtained by employing modified power-law, details given in Table 4.The derived characteristics of the SEPs during GLE # 71 are in practice the same as by Mishev, Koldobskiy, Usoskin, et al. (2021), and are in good agreement with the PAMELA direct measurements [for details see Adriani et al. (2015)].

Conclusion
A new open-source tool for conducting magnetospheric computations, called OTSO, has been developed at the request of the wider CR research community.The primary aim of which is to provide a user friendly alternative to older tools that fulfil the same purpose, such as MAGNETOCOSMICS.
OTSO has a good agreement with other magnetospheric computation tools, with the variations in results being likely due to differences in the integration methods used.
New models, integration methods, and optimisations will be incorporated into OTSO over time by the community upon scientific goals requests.Some of the additions to OTSO The creation of OTSO bolsters the CR research field's arsenal of tools that can be used to study CRs in the Earth's magnetosphere, providing the basis for detailed analysis of various CR experiments including GLEs and space weather service(s).
Within this work OTSO was successfully used for the analysis of two GLEs, namely the event that occurred during one of the strongest geomagnetic storms, GLE #66 on October 29, 2003, and the widely studied complex event, used for verification of NM data analysis using PAMELA measurements, GLE #71 that occurred on May 17, 2012.OTSO was able to obtain a good agreement with prior studies and in-situ space-borne measurements for these two events, proving that it is capable of being used to study complex events, such as those with high anisotropy like GLE #71, as well as events during intense and complicated magnetospheric conditions, such as GLE #66.Hence, it has been demonstrated that OTSO can be used as a reliable tool for geomagnetospheric computations under various conditions and circumstances, providing the necessary basis for strong SEP analysis.
As such OTSO represents, a community requested new generation tool, with the possibility for constant improvement, providing reliable geomagnetospheric computations related to CR research.
occurred on 13 December 2006 during the declining phase of solar cycle 23 as the result of a X3.4/4 B solar flare, with an associated GLE being detected around 03:00 UTC.The magnetic field distortion is described using an IOPT value within the Tsyganenko (1989) model used for these tests.For GLE # 70 the IOPT was set to 5, this corresponds to a planetary Kp index value of 4-, 4, 4+ at the time of the GLE. 12

Figure 1 .
Figure 1.Computed asymptotic cones for three selected NMs during GLE # 70 using both OTSO and MAGNETOCOSMICS, as denoted in the legend.The NMs shown are: Cape Schmidt (left), Oulu (middle), and Rome (right).
looks into this anomalous region in more detail.Within Figure3OTSO shows a gradual decrease in rigidity values with a significant penumbra present across the south pacific anomaly.In contrast MAGNETOCOS-MICS' plot is much more sporadic with sudden changes in R u and R l with a small penumbra is some regions, leading to the difference in R c seen in Figure2within this region.

Figure 2 .
Figure 2. Absolute difference in calculated effective vertical cutoff rigidities for the entire Earth during GLE # 70 between OTSO and MAGNETOCOSMICS.

Figure 3 .
Figure 3. Cross sections of the cutoff rigidity values over the region of largest difference between MAGNETOCOSMICS (left) and OTSO (right), taken at a longitude of −140 • .
elling the global NM network response and optimisation of the model response over experimental data, which involves computation of the ADs and cut-off rigidity of NM stations used for the data analysis; assuming a convenient initial guess for the optimisation [e.g.Cramp et al. (1995); Mishev et al. (2017)]; selection of model parameters and the optimisation itself (Mishev & Usoskin, 2016).The method was recently verified by direct space-borne measurements [for details see Koldobskiy et al. (2019); Mishev, Koldobskiy, Kocharov, and Usoskin (2021); Koldobskiy et al. (2021)].The ADs were computed employing the aforementioned two field combination: internal, namely IGRF geomagnetic model (Alken et al., 2021), and external model, Tsyganenko 89 The derived GLE spectral and PAD employing different combinations of magnetospheric models, namely IGFR+TSY 89 and IGRF+TSY 01 are comparable, despite several differences in asymptotic directions, specifically for MCMD and TERA NMs, for details see the right panel of Fig.4.Virtually the same spectra and PAD can be explained by the complexity of the unfolding procedure [see the discussion inHimmelblau (1972);Mishev, Koldobskiy, Usoskin, et al. (2021)].In general, the employment of IGRF+TSY 01 resulted on slightly harder spectra, wider PAD and reduced D. Note, that the asymptotic directions of SOPO, the NM with the greatest count rate increase, which is the station with maximal weight during the optimisation, are in practice the same.OTSO was also used for the analysis of GLE # 71, which occurred on May 17, 2012.The event was observed as a weak enhancement of the count rates at several NMs with greater signals recorded by APTY, OULU, and SOPO/SOPB NMs, while the other stations registered marginal count rate increases.This implied large anisotropy of the SEPs, confirmed by the following analysis [e.g.Mishev et al. (2014);Kocharov et al. (2018)].

Figure 4 .
Figure 4. Left panel: Asymptotic directions (IGRF+TSY 89) of selected NM stations during GLE #66 at 21:00 UT.The small circle depicts the derived apparent source position, and the cross the interplanetary magnetic field (IMF) direction obtained by the Advanced Composition Explorer (ACE) satellite.The lines of equal pitch angles relative to the derived anisotropy axis are plotted for 30 • , 60 • , and 90 • for sunward directions, and 120 • , 150 • for anti-Sun direction.Right panel: Comparison of computed asymptotic directions of selected NM stations during the GLE #66 employing TSY 89 (solid lines) and TSY 01 models (dashed lines).
will allow it to more accurately recreate older tools.The open-source and community driven element of this new tool will allow it to evolve into a robust magnetospheric computation tool that can facilitate the many needs of the CR research community, including space weather service(s), latitude surveys etc... [e.g.Mavromichalaki et al. (2018);Nuntiyakul et al. (2020)].OTSO has been designed to be as user friendly as possible, for both those wishing to edit the program and those with little programming knowledge.The main tool being accessed via python opens the tool up to computer novices and the inclusion of libraries such as IRBEM provides a strong foundation for OTSO's further development.As such the new tool provides a good starting point for a community driven magnetospheric computation tool.

Table 1 .
Data for the calculated effective vertical cutoff rigidity for selected NMs using both OTSO and MAGNETOCOSMICS.

Table 2 .
Derived spectral and angular characteristics during GLE # 66 on October 29, 2003 fitted with modified power-law rigidity spectrum employing AD computed with IGRF and TSY

Table 4 .
Derived spectral and angular characteristics during GLE # 71 on May 17, 2012 fitted with modified power-law rigidity spectrum employing AD computed with IGRF and TSY 89 models.