Layout 6 1 ANNALS OF GEOPHYSICS, 64, 4, GM440, 2021; doi:10.4401/ag-8507 O P E N ACC E S S Impact of approximate implementation manner in Klobuchar model Aghyas Aljuneidi*,1 and Hala Tawfek Hasan2 (1) Higher Institute for Applied Sciences and Technology (HIAST), Damascus, Syria (2) Higher Institute for Earthquake Studies and Research (HIESR), University of Damascus, Syria Article history: received June 5, 2020; accepted May 11, 2021 Abstract This paper focuses on the approximations that John A. Klobuchar made in mid 70s in his famous algorithm of ionospheric correction model for single frequency GPS receiver. At that time Klobuchar used a system of fixed geomagnetic north pole coordinates which are not accurate nowadays according to the International Geomagnetic Reference Field and to the World Magnetic Model because the geomagnetic poles move slowly. In addition, Klobuchar had to do other trigonometry simplifications in his implementation to avoid sophisticated computations. In order to evaluate this approximate implementation in a single frequency GPS receiver, ionospheric time and range delay are estimated on the entire day of January 1st 2010, using a different implementation in MATLAB. The required GPS data is obtained from recorded RINEX files at UDMC near DAMASCUS, SYRIA. In this comparative study, we reformulated the standard equations of Klobuchar model and examined the influence of its approximations on the ionospheric range delay and found a non- negligible bias in order of ten centimeters, whereas the influence of the movement of the geomagnetic poles was in order of few centimeters. Keywords: GPS; Geomagnetic North Pole; Total Electron Content; Ionospheric Correction; Klobuchar Model. 1. Introduction Global Navigation Satellite Systems (GNSSs) began after the launching of Sputnik, the first artificial satellite, by the Soviet Union in 1957, which was very motivating for the US scientists. Tracking of Sputnik by using the Doppler shift inspired the idea to track a receiver’s position on Earth by acquiring a certain modulated radio frequency emitted by satellites. Starting by the first satellite navigation program, named TRANSIT, and because of its numerous weaknesses, the US in mid 70s released its first Global Positioning System in the world, named GPS, which still operate today and is subject to continuous modernizations. Four years after the launch of GPS in 1978 the Soviet Union’s counterpart called GLONASS also started launching first satellites. Other systems like the European Galileo or the Chinese BEIDOU/COMPASS were started decades later. The fundamental principle of satellite positioning is based on the measurements of the signal transit time between the satellite and the receiver, the GPS system uses the L band frequencies 1575.42 MHz, 1227.6 MHz and recently 1176.45 MHz for the transmission of data. The observation set for a single channel GPS receiver are the pseudorange, the Doppler shift, the C/N0, and recently the carrier phase measurements [Parkinson, 1996]. The pseudorange is the time difference, when the data was sent by the satellite and when it was received, multiplied by the speed of light. The satellites also transmit data in their navigation message, which includes information like the health status of the satellite, rough clock and orbit correction terms and eight parameters for ionospheric error correction. Table 1 presents the most important GPS errors for single frequency receivers. Table 1. Statistical ranging error budget (1𝜎) adopted from [Kaplan and Hegarty, 2006]. After the end of Selective Availability (SA) in May 1st 2000, the ionosphere became the highest source of ranging error for GPS C/A code receivers, as a result, reducing this error is fundamental [Kaplan and Hegarty, 2006]. Different methods can be adopted to minimize the ionosphere effect in single receiver: • use of multi-frequency technique; • use some kind of augmentation system; • use of ionospheric model. Due to the dispersive nature of the ionospheric refractivity, the ionospheric error can be corrected by the ionosphere-free combination with simultaneous measurements at two or more frequencies, [Wang, 2016], this method is the most effective, but it cannot be used in a single frequency receiver. Instead, single frequency receivers can acquire the ionospheric correction as differential corrections from Differential GPS (DGPS) or from Satellite Based Augmentation System (SBAS), by transmitting corrections to the GPS receivers either via satellite or terrestrial radio. Augmentation systems are not available around the whole globe and involve the use of complex architectures and expensive receivers, while the use of broadcast ionospheric models, such as the Klobuchar Model, is cheaper and easy to implement in a single frequency GNSS receiver. 2. Problem statement The estimation and removal of ionospheric delay errors in the real-time applications of the single frequency users continued to remain a potential challenge to be addressed. The Klobuchar model was the first ionospheric correction model implemented in near real time. Over few decades it has become the standard ionospheric correction used by virtually all single frequency GPS receivers in near Earth for mid.-latitude. It is still used by a very large class of commercial low-cost single frequency receivers around the world till now. In spite of its numerous advantages, like not requiring additional parameters outside those contained in the Aghyas Aljuneidi et al. 2 Segment Source Error Source GPS1σ Error Space Ephemeris data Satellite clock 2.1 m 2.1 m User Ionospheric delay Tropospheric delay Multipath Receiver noise 4.0 m 0.7 m 1.4 m 0.5 m Sum User equivalent error RMS 5.3 m GPS navigation messages, fastness and ease of implementation, it is considered an outdated model from an accuracy point of view relative to nowadays classifications. However this model can reduce the RMS position error about 50% to 60% [Klobuchar, 1975]. When it was implemented and because of computational load concerns, John Klobuchar, the inventor, used some approximate equations. Finally, John Klobuchar fixed the geodetic coordinate of the geomagnetic north pole in his algorithm according to some fitting in that time as 69W longitude and 78.3N latitude, which is not accurate nowadays because the poles wander with time. Based on the International Geomagnetic Reference Field document (IGRF-13) in 2020, the geomagnetic north pole’s 72.7W longitude and 80.7N latitude are slightly different from their counterparts in the World Magnetic Model (WMM2020) located at 72.68W longitude and 80.65N latitude for the same year. In reality Klobuchar adopted the geodetic coordinates of the geomagnetic north pole as 69W Longitude and 78.3N Latitude, [Klobuchar, 1975; Davies, 1965], for more details see page 20. 2.1 Objectives of study The main objective of this paper is to analytically study some defects in the implementation of the Klobuchar model and to do a new implementation with different considerations, then compare the results, in two cases: a) standard implementation of Klobuchar equations without any trigonometry approximations; b) set the exact coordinate of geomagnetic north pole. 2.2 Importance of study Because the GPS satellite message has a place for only eight coefficients to describe the worldwide behavior of the Earth’s ionosphere and no global coverage of Augmentation Systems, it is virtually impossible for most of the single frequency GPS users around the world to improve the performance of their receivers concerning ionospheric corrections, which is the largest error in GPS. 3. Ionosphere and geomagnetism The ionosphere is defined as that part of the upper atmosphere where the density of free electrons and ions is high enough to affect the propagation of electromagnetic radio frequency signals [Kaplan and Hegarty, 2006]. The signals are delayed instead, and this delay is part of the travel time, which represents an error into the computed range. The error is negative for the carrier phases (phase advance; that is, the measured range is smaller than the real range), and positive for the pseudorange (group delay; that is, the measured range is longer than the real range). The phase advance and group delay are equal in value and opposite in sign. For a good first order approximation, the delays that introduced in radio signals by the ionosphere (T����), are correlated with the Total Electron Content (TEC) which is defined as the total number of electrons presents in a cross section unit area of 1 𝑚² along the propagation path of the satellite signal to the receiver and inversely proportional to the square of the carrier phase frequency (f ), consequently: (1) Thus, all electron density values along the propagation path of the signal are integrated in order to obtain STEC, or slant TEC, which different from Vertical TEC or VTEC, can be mathematically represented as this integral: 𝑆𝑇𝐸𝐶 = ∫�𝑁�(𝑠). 𝑑𝑠 (2) T����(second) = 40.3 ���� ��² 3 Impact of implementation manner in Klobuchar model Aghyas Aljuneidi et al. 4 Where, 𝑁�(𝑠) is the electron content per volume unit and 𝑝 is the propagation path between the source and the receiver. STEC is measured in TEC unit’s i.e. TECU, which equals to 1016 electrons/𝑚². Accordingly, estimating and modelling TEC is a necessary step to correct the ionospheric effects on electromagnetic waves propagating from satellites to Earth, On the other hand, Total Electron Content (TEC) measurements provide a rich source of information about the Earth’s ionosphere. thus, having GPS-derived TEC observations in comparison with another expensive method like satellite altimetry, incoherent scatter radar… etc., GPS became a promising remote sensing sensor to model the Ionosphere and to do many studies about related geophysical phenomenon [Shuanggen et al., 2014]. Worth to mention, that TEC is highly correlated with the critical frequency of the F2 layer (fo F2 or NmF2). The NmF2 is the maximum electron density at the F2 layer of the ionosphere, approximately 300 - 500 km altitude, which traditionally was measured by Ionosonde [Kouris et al., 2008]. The Solar activity follows a regular periodic variation with eleven years periodicity, commonly known as solar cycle. The solar cycle is the periodic recurrence of sunspots at Sun’s surface, therefore TEC and NmF2 are affected by solar activity (as measured by solar radio flux f10.7) as a long-term ionospheric variations, on the other, during the solar cycle there occurs periodic change in the solar radiation and the ejection of solar material which induce geomagnetic storms, characterized by the fluctuations of the Earth magnetic field under the influence of solar wind and interplanetary magnetic field, which is a short-term ionospheric variations (including TEC). Particularly, TEC changes with time in the forms of diurnal, seasonal, and long-term yearly variations, Spatially, it changes with geographic longitude and latitude of the GPS receiver (user) location, since the ionospheric properties are based mainly on the interactions between the solar radiation and the Earth’s geomagnetic field [Davies, 1965], as well as TEC is related by the geomagnetic coordinates of the point of ionosphere where its estimation is done. The Klobuchar model assumes the Earth’s magnetic field to be an Earth-centered dipole, so the relationship between the geomagnetic coordinates (𝜙� the latitude and 𝜆� the longitude) and the corresponding geographic coordinates (𝜙₁, 𝜆₁) at a certain point P is given by (3) where 𝜙� and 𝜆� are the geographical latitude and longitude of the north geomagnetic pole, [Davies, 1965]. 3.1 Geomagnetic poles A certain shift or displacement exists between the geomagnetic poles and the Earth’s poles in the two directions of the south and the north, and moves in certain acceleration to the west. Its change will affect the sector of navigation, communication and military, and other aspects [e.g., Guo et al., 2011]. The Earth’s magnetic field is in fact a composite of several magnetic fields generated by variety of sources. These fields are overlaid on each other and interact with each other through inductive processes. The most important of these geomagnetic sources are: a) the Earth’s outer core; b) magnetized rocks in Earth’s cover; c) fields generated outside the Earth by electric currents which are flowing in the ionosphere and magnetosphere; d) electric currents flowing in the Earth’s crust (usually induced by varying external magnetic fields); e) Ocean current effects. All of these parameters vary with time of measures, ranging from milliseconds to millions of years (i.e. magnetic reversals phenomena) [Prö�lss, 2004]. The World Magnetic Model (WMM) and the International Geomagnetic Reference Field (IGRF) are estimated from the most recent data and they are of a comparable quality. The difference between WMM and IGRF are within expected model of inaccuracy. The WMM is a predictive only �sin𝜙� = sin𝜙� sin𝜙� + cos𝜙� cos𝜙� cos(𝜆� ‒ 𝜆�) sin𝜆� = cos 𝜙� sin(𝜆�+𝜆�) cos 𝜙� model and the WMM2020 is valid for the current epoch (2020-2025). The IGRF is retrospectively updated and the latest update IGRF-13 is valid for years 1900-2025. Models need to be revised at least every five years because of the varying nature of the magnetic field with both time and location, this change in the geomagnetic field called “the geomagnetic secular variation”, [Chulliat et al., 2020; Thébault et al., 2015], and each epoch model comprises a set of spherical harmonic coefficients (called Gauss coefficients) in a series expansion of the global magnetic field. The geomagnetic poles (dipole poles) are the intersections of the Earth’s surface and the axis of a bar magnet theoretically placed at the center the Earth by which we approximate the geomagnetic field and can be computed from the first three coefficients from a main field model, On the other hand, the magnetic poles are the points at which magnetic needles become vertical and can be computed from all the Gauss coefficients using an iterative method. In Table 2, we show predicted geodetic locations of the geomagnetic north pole by “International Geomagnetic Reference Field” (IGRF-13) including prediction. 4. Ionospheric modelling There have been numerous researches on ionospheric modeling, whether global, regional and local, or near real time and post processing models, these researches include: • Empirical models based on ionospheric on extensive world-wide data measurements. Usually, the data are collected over an extended period of time and then fitted with simple analytical functions. 5 Impact of implementation manner in Klobuchar model Table 2. Locations of geomagnetic north pole in geodetic coordinate based on IGRF-13 including prediction adopted from [Kyoto University, World Data Center for Geomagnetism, Kyoto]. Year IGRF-13 Latitude Longitude 1965 78.6N 69.9W 1970 78.7N 70.2W 1975 78.8N 70.5W 1980 78.9N 70.8W 1985 79.0N 70.9W 1990 79.2N 71.1W 1995 79.4N 71.4W 2000 79.6N 71.6W 2005 79.8N 71.8W 2010 80.1N 72.2W 2015 80.4N 72.6W 2016 80.4N 72.6W 2017 80.5N 72.6W 2018 80.5N 72.7W 2019 80.6N 72.7W 2020 80.7N 72.7W 2021 80.7N 72.7W 2022 80.7N 72.7W 2023 80.8N 72.7W 2024 80.8N 72.6W 2025 80.9N 72.6W • Global Ionospheric Map (GIM) and regional ionospheric maps are Numerical maps, which provide values of ionospheric parameter in a grid. • Analytical models are based on orthogonal function fits to the output that obtained from numerical models. • Physical models are developed based on typically solution of the continuity equation or the momentum and energy equations for the electrons and ions. 4.1 Empirical ionospheric models There are four types of empirical model: • Bent model; • International Reference Ionospheric (IRI) model; • NeQuick model; • Klobuchar model. 4.1.1 Bent model Is an empirical, global model which was created in 1973 using soundings and satellite measurements taken from 1962 to 1969, a period which covered both a solar minimum and a solar maximum. Bent model was developed for ground-to-satellite communications by using hundreds of parameters, this global model provides the TEC in the altitude range from 150 – 2000 km as a function of position of the observer, time, solar flux and the sunspot number … etc., Bent predicts that the model is 75–80% accurate [Bent et al., 1973]. 4.1.2 IRI model The first version of International Reference Ionospheric (IRI) model was formulated by International Union of Radio Science (URSI) and Committee on Space Research (COSPAR) in 1978. With the development of techniques and the accumulation of observation data, IRI model has been updated frequently. Data sources of IRI model include Ionosonde, radars and several satellites and rockets. [Bilitza et al., 2017]. In this model, Ionospheric electron density profile can be derived easily after specifying date, time, location, sunspot number and some other parameters, then the TEC and ionospheric delay can be calculated, worth mentioning, IRI models suffer from lack of knowledge of the upper part of the ionosphere, since it is limited at 2000 km height. 4.1.3 NeQuick model NeQuick is a three-dimensional and time dependent ionospheric electron density model, taking advantage of the increasing amount of available data, there are many NeQuick versions aiming to update better representations of the ionosphere for the whole globe. [Nava et al., 2008]. The NeQuick1 has been used by the European Space Agency (ESA) European Geostationary Navigation Overlay Service (EGNOS) project and has been chosen to single-frequency positioning applications in Galileo, it has also been adopted by the International Telecommunication Union, Radio communication Sector (ITU-R) as an appropriate method for (TEC) modeling. 5. Klobuchar model basics This is the standard correction used by almost all GPS receivers, which operate near Earth, with fewer parameters. The Klobuchar model is a simple ionospheric model that gives the vertical ionospheric delay at a given time and Aghyas Aljuneidi et al. 6 location for the GPS single-frequency users in near real-time. It was developed in 1975 by John A. Klobuchar at the Air Force Geophysics Laboratory, U.S. Due to quite limited resources of computer microprocessor systems used in navigation receivers at that time, this system had to be extremely undemanding for such resources. Klobuchar model is an ionospheric correction algorithm (ICA). Over many years of operation, this model demonstrated to be efficient. Since the TEC converted into correction to pseudo-range can vary from a few meters to tens of meters depending on solar activity, time, user location, satellite’s azimuth and elevation relative to user location and many others parameters [Klobuchar, 1987]. 5.1 Single layer ionospheric assumption Ionospheric Klobuchar model was designed based on the Bent model. It is defined as a single layer ionospheric model (SLM-Single Layer Model), because the TEC is assumed as concentrated in an infinitesimal layer placed at an average altitude of 350 km from the Earth’s surface (i.e. maximum electron density at the F2 peak). In a SLM the VTEC (vertical value of TEC) is calculated at the Ionospheric Pierce Point (IPP) which is a geographic point, obtained by the intersection between the thin layer, (SLM) and the line of sight between the satellite at ~22500 km height and the receiver (user) on earth surface (Figure 1). Where: • 𝜙�, 𝜆�: geodetic location of user. • 𝑍: The zenith angle at the user location. • 𝛽: The zenith angle at the IPP. • 𝜙�, 𝜆�: geodetic location of the Ionospheric Pierce Point (IPP), it is worth to mention that the geomagnetic latitude of the same point is 𝜙�. • 𝐸: is equal to �/₂ ‒ 𝑍, the elevation of the GPS satellite relative to user location, at the current moment,, local time (seconds) noticing that the Azimuth well be 𝐴𝑧. • RE: radius of the Earth, 6378137.0m according to WGS84, worth to mention. • hS: altitude of SLM, equal to 350 km. • 𝜓pp: Earth-centered angle, equal to (𝑍 ‒ 𝛽), then 𝜓pp = 𝑍 ‒ 𝛽 = �/₂ ‒ 𝐸 ‒ 𝛽. (4) with a simple trigonometry 7 Impact of implementation manner in Klobuchar model Figure 1. Single layer ionospheric model adopted from [Bolaji et al., 2016]. (5) (6) 5.2 Ionosphere-mapping function All slant TEC measurements are converted to an equivalent vertical value of TEC, i.e. VTEC by using the thin shell elevation mapping function because Mapping to vertical is necessary if the radial electron density structure is not modeled [Mannucci et al., 1993], thus for any arbitrary line of sight, STEC (the slant angle) of the satellite must be taken into account: (7) As illustrated in Figure 2. The GPS Klobuchar model uses an elevation mapping function, i.e. obliquity factor to convert vertical to slant delay at user level. (8) which leads to: (9) the simplified elevation mapping function [Klobuchar, 1975] proposed by Klobuchar in function of elevation E (semicircles) is: 𝐹 = 1 + 16. (0.53 + 𝐸)³ (10) 𝑠𝑖𝑛 𝛽 = � 𝑠𝑖𝑛 𝑍� = � 𝑐𝑜𝑠 𝐸��� �� + ���� �� + �� 𝜓pp = 𝑍 ‒ 𝛽 = �/₂ ‒ 𝐸 ‒ 𝑠𝑖𝑛⁻¹ � 𝑐𝑜𝑠 𝐸��� �� + �� 1 ����𝑆𝑇𝐸𝐶 = 𝑉𝑇𝐸𝐶 𝐹 = = 1 �������� ���� 𝑇���� = 40.3 * 𝐹 * ���� ��² Aghyas Aljuneidi et al. 8 Figure 2. STEC and VTEC in Single layer ionospheric model adopted from [Janssen, 2012]. Bear in mind that STEC is equal to TEC of equation (1). 9 Impact of implementation manner in Klobuchar model While according to [Mannucci et al., 1993], the elevation mapping function, with a simple trigonometry, became: (11) 5.3 Radius of the Earth Klobuchar articles [Klobuchar, 1975; Klobuchar,1987] do not state explicitly the length of the Earth radius used in the approximated formulas. However, spherical geometry was applied in the mid 70s. We could not use WGS84 formula directly without justification, and the value of radius the earth could be inferred from the numerical equations provided by Klobuchar himself. So, by analyzing [Klobuchar, 1975] page 10 and comparing these two equations: Where • 𝑍 is the zenith angle at the IPP; • 𝑒𝑙 is the elevation angle, • ℎ is the altitude of SLM layer, which equal to 350 km, and • 𝑟₀ is the length of earth radius in Klobuchar Implementation. By solving the equation = 0.948 with taking in account the approximated value 0.948, we found 𝑟₀ = 6380.76 km which is greater than the equatorial radius = 6378.16 km and the polar radius = 6356.16 km and the main radius = 6371.03 [Allen, 1976; Flock, 1983], where all these previous values when applied to will give the same value, i.e. 0.948, thus, as a decision, we will update as a medium and accepted value of the main radius of the Earth, 𝑅� = 6378137.0 m according to WGS84 in all formulas of this article regardless of Latitude dependence (considering it as the mean radius of the spherical earth as Klobuchar has done in its model of mid 70s). 5.4 Model layout Klobuchar model assumes that there is a constant delay of 5 ns during night time equal to 9.2 TECU i.e. 1.5 m range error at L1 frequency, and a cosine function (in the very basic variant, expanded in a Taylor series in order to save the computing resources) in daytime that is centered at the 14th hour of local time. Therefore, the maximum semi period for the cosine function is 20 hours, because the value at the end of the day must be the same at the beginning of the following day (i.e. 5 ns according to the model, Figure 3). Accordingly, the time delays 𝑇���� (in seconds) at observer’s location [Klobuchar, 1975]: (12) where • DC (seconds) is the night-time constant offset term (set to 5 ns corresponding to 1.5 meters ionospheric delay on L1). • A is the amplitude term (seconds). • 𝜙 is the phase of the maximum vertical ionospheric delay which is assumed to be at 14:00hours local mean solar time (50400 seconds). • P is the period (seconds). • t is local time in IPP (seconds), where local time is equal to Time Zone (related to the longitude) + UT (the 𝑍 = 𝑠𝑖𝑛⁻¹ ( cos (𝑒𝑙)) and 𝑍 = 𝑠𝑖𝑛⁻¹(0.948cos (𝑒𝑙))�₀ �₀ � � 𝐹 = 1��1�� cos 𝐸�²��� �� + �� 𝑇���� = 𝐹�𝐷𝐶 + 𝐴𝑐𝑜𝑠( )��²�⁽�⁻�⁾� � Universal Time), as there are several versions of UT, [IERS, 2018], the most commonly used being Coordinated Universal Time (UTC), GPS TIME was zero at 0h 6-Jan-1980 and since it is not perturbed by leap second, GPS is now ahead of UTC by 18 seconds, so according to Klobuchar: thus Hence 𝑡(𝑠𝑒𝑐𝑜𝑛𝑑) = 12 * 3600 * 𝜆��(semi circles) + 3600 * 𝑈𝑇(ℎ𝑜𝑢𝑟). (13) and • 𝑭 is the Elevation mapping function. The A and P terms each have four time and solar-activity dependent coefficients to represent the behavior of the global ionosphere. To analytically represent the amplitude and period terms, a cubic polynomial in geomagnetic latitude is used to represent ionospheric delays which was the best available to model the whole world TEC at the time of development. [Klobuchar, 1975]. Because the Bent model did not cover the low latitudes appropriately, Klobuchar did not adopt higher than the cubic polynomial, as might be required in the equatorial anomaly region, thus, it is correct to say that the Klobuchar model is less accurate in the low-latitude regions that in the middle latitudes. The model is also less accurate in the highly variable polar regions. 5.5 Broadcast ionospheric model The Klobuchar model is the official ionospheric broadcast model (IBM) for single frequency GPS receivers in real time, it is derived from the empirical knowledge of the long-term ionosphere behavior. It is designed for single frequency users in which the required model parameters are broadcast via the navigation message from the satellites, the GPS Ground Control Segment updates daily these coefficients according to season and solar activity level. 𝑡(ℎ𝑜𝑢𝑟) = + 𝑈𝑇(ℎ𝑜𝑢𝑟)���₍���₎ ¹⁵ 𝑡(ℎ𝑜𝑢𝑟) = 12 + 𝑈𝑇(ℎ𝑜𝑢𝑟) = 12𝜆��(semi circles) + 𝑈𝑇(ℎ𝑜𝑢𝑟).���₍���₎ ¹⁸⁰ Aghyas Aljuneidi et al. 10 Figure 3. Ideal time delay half cosine function in day time adopted from [Sanz et al., 2013], simplified from [Klobuchar, 1987]. 11 Impact of implementation manner in Klobuchar model A few broadcast parameters and coefficient are decoded from UTC and IONO blocks from page 18 sub frame 4 in ALMANAC file, these eight coefficients are determined by a competent processing center in the USA and transmitted to GPS satellites, after which they are translated to users as part of GPS navigation messages, these parameters are: 𝛼₁, 𝛼₂, 𝛼₃, 𝛼₄, and 𝛽₁, 𝛽₂, 𝛽₃, 𝛽₄. As a rule, the coefficients are updated incomprehensibly once in a few days depending on the ionosphere variability, this makes the model intrinsically insensitive to short-term ionospheric variations (e.g. geomagnetic storms, interplanetary magnetic field effects, etc.). Many efforts had been done in the aim to improve the consistency of Klobuchar model taking into consideration that only eight coefficients are introduced to describe the global variation of the ionosphere due to the communication limitation in GPS navigation message. accordingly, many studies which operate in post processing time attempted to re_estimate these eight parameters [Yuan et al., 2008], or to extend the structure of Klobuchar model (Klobuchar-like model) [Orus et al., 2002]. In this study we will take in consideration a set of Improved Klobuchar-Style Ionospheric, [IGS, 2010], in purpose for comparison with the set of coefficients downloaded from ALMANAC in the same date. While the Chinese System, BeiDou, ionospheric model (BIM) is based on the Klobuchar model, it has eight parameters although it has some modifications, the difference between them is that BIM is situated in terms of geographic coordinates while Klobuchar is situated in terms of geomagnetic coordinates, another difference is BIM uses the ionosphere height as 375 km, compared with Klobuchar model which uses 350 km. [BeiDou Navigation Satellite System Signal In Space Interface Control Document, 2016]. The NeQuick-G model is a specific implementation of the NeQuick model for its use in the Galileo system for single-frequency ionospheric correction, [European GNSS (Galileo) open service, 2016], the model coefficients are broadcast in the navigation message since many years, three world-wide coefficients are to be broadcast to the user, these values are computed at the Galileo system level, using a set of world-wide distributed monitoring stations to evaluate slant TEC needed to determine the coefficients on a current day for use on the following day. The proposed GLONASS model, [Ivanov, 2017], is a semi-empirical model of vertical distribution of the concentration of ionospheric free electrons, the control input parameters are the solar activity index 𝐹₁₀.₇ (the intensity of RF radiation of the Sun at the wavelength of 10.7 cm), the geomagnetic activity index 𝐴�, and the correction parameter of electrons concentration at the peak of iono phere region 𝐹₂. The transmission of the first and the third parameters to users in the GLONASS navigation messages started in 2017-2018. Although the index 𝐴� is available on public Internet sites. Merited to mention, the GEMTEC model, [Gorbachev, 2015] which is another proposed broadcast ionosphere model overcomes the Klobuchar model and requires only a single external input parameter 𝐹₁₀.₇. 5.6 The model accuracy The high variability in the ionosphere will prevent any model from being more than 75% accurate without the knowledge of the specific conditions [Klobuchar, 1975]. Klobuchar model assumes an ideal smooth comportment of the ionosphere, therefore any significant variations from day to day will not be modeled properly. The model can reduce the RMS position error due to uncompensated ionospheric delay by about 50% to 60% depending on the solar activity or the region. Therefore, its accuracy is often unsatisfactory even for the absolute positioning, especially in the case of the high solar activity [Stepniak, 2014]. Regardless of our concern in this study about the weakness of implementation of Klobuchar model, the model’s accuracy is basically limited due to the following three reasons: • First, the thin layer model currently used in GPS has deficiencies resulting from conversion of slant TEC to effective vertical TEC. The deficiencies come from inappropriate attribution of the thin shell height. This conversion introduces a few errors in the middle latitude where electron density is small. [Norsuzila et al., 2008]. • Second, a table of predefined values is used to update the eight coefficients in the GPS system. These values are predefined for a particular condition in the solar cycle. They are not set to be updated every day. • Third, the diurnal variation of the ionosphere is calculated using the eight parameters in the broadcast ephemeris, and setting the delay values at night as a constant, regardless of the latitude, makes the model fail to model the temporal changes, particularly at night. • Fourth, the use of a mapping function is clearly an approximation, and can lead to errors of several TECU when applied to regions of large horizontal electron density gradients [Klobuchar et al., 1993]. • Sixth reason, seems to be an additive approximation, the value of the Earth radius was fixed independently of latitude because Klobuchar used the spherical model for the Earth. 6. Klobuchar model algorithm The Klobuchar model is a compromise between computational complexity and accuracy, its limited number of inputs makes it an efficient and easy tool of many researches. The TEC must be found at the Ionospheric Pierce Point (IPP), since the TEC is best ordered in geomagnetic, rather than geodetic coordinates, the assumption is that the TEC is only a function of geomagnetic latitude and local time. Any two points on the Earth having the same geomagnetic latitude are assumed to have the same TEC at the same local time. Thus, a conversion to geomagnetic latitude is required, [Klobuchar, 1975]. 6.1 Implementation In this study we considered the coordinate of geomagnetic pole as an input because it is a varying with time. 6.1.1 Inputs • 𝜙�₀, 𝜆�₀: geodetic location of geomagnetic north pole. • 𝜙�, 𝜆𝑢: geodetic location of user(rad). • 𝐸, 𝐴𝑧: elevation and Azimuth of GPS satellite relative to user location(rad), worth to mention that these values could be calculated from the , and the coordinates of Satellite ,,. • 𝑡: local time in the IPP (seconds). • RE, hS: radius of the Earth, 6378137.0 m according to WGS84 and the altitude of SLM, equal to 350 km. • 𝛼₁(second), 𝛼₂(second/semicircle), 𝛼₃�( )²�, 𝛼₄�( )³� and 𝛽₁(second), 𝛽₂�( )�, 𝛽₃�( )²�, 𝛽₄�( )²�: the eight parameters broadcast in the navigation message. 6.1.2 Outputs • 𝑇���� the slant time delay (seconds), thus we can obtain both the range delay 𝑑����(𝑚) and correspond TEC value. • 𝐷���� the range delay (meters). • The correspond TEC value. 6.1.3 Implementation steps In the days of Klobuchar, the calculation complexity was an obstacle for implementation of Klobuchar model equations, such several approximations are made in the geometric calculations to further reduce operational user computer storage and running time requirements, hence we cited in detail, with their units, both the approximate equations i.e. indexed by ₋���� and the standard equations without any approximation with no index i.e. which face no problem nowadays looking at the huge progress in the digital era. Thus: Aghyas Aljuneidi et al. 12 13 Impact of implementation manner in Klobuchar model Step1: the Earth-centered angle (14) Step2: geodetic latitude of IPP (15) as a condition imposed by Klobuchar. If 𝜙�� is less than ‒0.416 rad let 𝜙�� = ‒0.416 (semicircles), and if 𝜙�� is more than 0.416 let 𝜙�� = 0.416 (semicircles). Take into account that 0.416 = 75/180 which is identical for Klobuchar assumption that its algorithm is suitable for only mid. Latitude not greater than absolute (75) deg. [Klobuchar, 1975], noteworthy for a specific satellite, the geodetic latitude of IPP varies proportionally with the geodetic latitude of user. Step3: geodetic longitude of IPP (16) Step4: geomagnetic latitude of IPP 𝜙� = 𝑠𝑖𝑛⁻¹�𝑠𝑖𝑛 𝜙�� 𝑠𝑖𝑛 𝜙�₀ + 𝑐𝑜𝑠𝜙�� 𝑐𝑜𝑠𝜙�₀ 𝑐𝑜𝑠(𝜆�� ‒ 𝜆�₀)� (17) With 𝜙�₀ and 𝜆�₀ are the latitude and longitude of geomagnetic north pole respectively, and 𝜙�_���� = 𝜙��_���� + 0.064𝑐𝑜𝑠 �𝜋(𝜆��_���� ‒ 1.617)� (18) Notice that constants 0.064 and 1.617 are rounded to geomagnetic north pole coordinates in the original Klobuchar algorithm 𝜙� = 78.3 and 𝜆�₀ = 291 , where: 𝑐𝑜𝑠 (𝜙�₀)/𝜋 = 0.0645, and 𝜆�₀ /𝜋 = 1.6167, thus we have: (19) From a mathematical viewpoint on the Eq. (17), (19) and based on the IGRF-13 secular variation, the slowly drifting value of geomagnetic north pole coordinates towards north and east, knowing that 𝜙�₀= 80.9N for the year 2025 from Table 2. we conclude that the geomagnetic latitude of IPP will basically rely slowly with time on geodetic latitude of IPP, thus, for a specific satellite, also on the geodetic latitude of the user, hence in future, geomagnetic latitude of IPP will be greater on high latitude, based on this conclusion and according to sign of 𝛼₁,𝛼2,𝛼3,𝛼4 and 𝛽₁,𝛽2,𝛽3,𝛽4 the slant time delay on high latitude will be highly affected by this deduction. Step5: Amplitude of half cosine (20) With condition, if A < 0, let A = 0. 𝜓pp (𝑟𝑎𝑑) = �/₂ ‒ 𝐸 ‒ 𝑠𝑖𝑛⁻¹ � cos(𝐸)� 𝜓pp���� (semi circles) = � � ‒ 0.022 �� �� + ��� ₀.₀₁₃₇ ₀.₁₁+ �/� � 𝜙pp (𝑟𝑎𝑑) = 𝑠𝑖𝑛⁻¹ �𝑠𝑖𝑛(𝜙u)𝑐𝑜𝑠(𝜓pp) + 𝑐𝑜𝑠(𝜙u)𝑠𝑖𝑛(𝜓pp)𝑐𝑜𝑠(𝐴𝑧)� 𝜙pp₋���� (semi circles) = + 𝜓pp₋���� 𝑐𝑜𝑠(𝐴𝑧)�� � � 𝜆pp (𝑟𝑎𝑑) = 𝜆u + 𝑠𝑖𝑛⁻¹ �𝑠𝑖𝑛(𝜓pp) 𝑠𝑖𝑛(𝐴𝑧)/𝑐𝑜𝑠(𝜙u)� 𝜆pp₋���� (semi circles) = + 𝜓pp₋���� 𝑠𝑖𝑛(𝐴𝑧)/𝑐𝑜𝑠(𝜙u)�� � 𝜙�_���� = 𝜙��_���� + � � 𝑐𝑜𝑠 (𝜋𝜆������ ‒ 𝜆�₀) 𝑐𝑜𝑠 (𝜙�₀) 𝜋 � 𝐴(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝛼₁ + 𝛼₂ + 𝛼₃ � �² + 𝛼₄ � �³ 𝐴���� (𝑠𝑒𝑐𝑜𝑛𝑑) = 𝛼₁ + 𝛼₂𝜙�_���� + 𝛼₃𝜙�_����² + 𝛼₄𝜙�_����³ �� ��� � �� � � ¹⁸⁰ � ¹⁸⁰ Step6: Period of half cosine (21) With condition if P < 72000 second, means less than 20h, let P = 72000. Step7: the mapping function (22) Step8: local time at IPP (23) With condition if t >86400 seconds, means more than 24h, let t = t-86400, and if t <0 let t = 86400. Step9: Compute the phase Klobuchar, in his algorithm, considered that the peak of Ionospheric delay occur at 14h local time, where 14h = 50400 second. (24) Step10: the slant time delay The slant time delay is converted from vertical time delay when multiplying by the mapping function F. Thus, if the phase X <𝜋/2 then it is the day time: (25) otherwise, it is the night time: (26) with 𝐷𝐶 =⁻⁹(second), the nighttime delay constant. Step11: the delay range error (27) where 𝑐 = 299792458 𝑚⁄𝑠 𝑃(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝛽₁ + 𝛽₂ + 𝛽₃ � �² + 𝛽₄ � �³ 𝑃���� (𝑠𝑒𝑐𝑜𝑛𝑑) = 𝛽₁ + 𝛽₂𝜙�_���� + 𝛽₃𝜙�_����² + 𝛽₄𝜙�_����³� �� ��� � �� � 𝐹���� = 1 + 16(0.53 + �/� )³� ��1�� 𝑐𝑜𝑠 𝐸�²� 1𝐹 = �� �� + �� 𝑡(𝑠𝑒𝑐𝑜𝑛𝑑) = 12 * 3600 * � � + 3600 * 𝑈𝑇(ℎ𝑜𝑢𝑟) 𝑡����(𝑠𝑒𝑐𝑜𝑛𝑑) = 12 * 3600 * �𝜆������� + 3600 * 𝑈𝑇(ℎ𝑜𝑢𝑟)� ��� � 𝑋(𝑟𝑎𝑑) = 2𝜋 � � 𝑋����(𝑟𝑎𝑑) = 2𝜋 � �� �₋₅₀₄₀₀ � �����₋₅₀₄₀₀ ����� 𝑇����(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝐹�𝐷𝐶 + 𝐴𝑐𝑜𝑠(𝑋)� 𝑇����_����(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝐹�����𝐷𝐶 + 𝐴����𝑐𝑜𝑠(𝑋����)�� 𝑇����(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝐹.𝐷𝐶 𝑇��������(𝑠𝑒𝑐𝑜𝑛𝑑) = 𝐹����.𝐷𝐶� 𝑑����(𝑚) = 𝑐.𝑇���� 𝑑��������(𝑚) = 𝑐.𝑇��������� Aghyas Aljuneidi et al. 14 15 Impact of implementation manner in Klobuchar model Step12: the TEC in the IPP (28) 6.2 Data acquisition The data acquisition is done on January 2010 from 0h to 24h (UTC) at (UDMC: Latitude/Longitude: 33°30′36′′N / 36°17′06′′E, with Leica GRX1200 Receiver and Ashtech Antenna) near DAMASCUS, SYRIA. The eight ionospheric coefficients downloaded from ALMANAC are: 𝛼₁,₂,₃,₄ = [7.4506e-9, -1.4901e-8, -5.9605e-8, 1.1921e-7]; 𝛽₁,₂,₃,₄ = [9.0112e4, -6.5536e4, -1.3107e5, 4.5875e5]. In purpose of comparison, we downloaded a set of improved Klobuchar-Style Ionospheric Coefficients form CODE’S GLOBAL IONOSPHERE MAPS FOR DAY 001, 2010, as follows [IGS, 2010]: 𝛼��������₁,₂,₃,₄ = [8.8647e-9, -1.4985e-8, -2.5446e-7, -2.9687e-7]; 𝛽��������₁,₂,₃,₄ = [9.1360e4, -6.8880e4, -1.235e6, -5.2002e6]. 6.2.1 Satellite visibility The satellite visibility is related to many parameters, one of these parameters is the receiver quality, worth to mention that this study is oriented to such low cost or commercial single frequency GPS receiver, with no coverage for augmentation systems, this kind of GPS receivers is very widespread for civilian users around the third world. � 𝑆𝑇𝐸𝐶(𝑇𝐸𝐶𝑈) = 𝑆𝑇𝐸𝐶����(𝑇𝐸𝐶𝑈) = 𝑓²𝑑���� ⁴⁰.³𝑓²𝑑����_���� ⁴⁰.³ Figure 4. Visibility of Satellites for a complete day January 1st 2010 at UDMC. Knowing that the constellation includes nominally a minimum of 24 operational satellites and a certain number of backup satellites are on orbit, making available up to 31 satellites for positioning, in Figure 4 GPS satellites orbit twice in a sidereal day in inertial space i.e. a GPS satellite orbits once every 11 hours and 58 minutes, with the Earth rotating once every 24 hours, such a GPS satellite pass up over a point above the Earth in average once a day, in Figure 4. The visibility is very good with eight satellites at least in each moment while all satellites pass twice above UDMC except for the Satellites 12,13,18,20,21,22,23,24,29,30 and 32 which pass once only and the Satellites 1 and 25 never do. 6.2.2 All satellites range error In aim to evaluate the effect of approximation adopted by [Klobuchar, 1975], two implementations were written: • Approx. Implem.: trigonometry approximate implementation with approximate geomagnetic north pole coordinates adopted by Klobuchar, i.e. 𝜙�₀ = 78.3N and 𝜆�₀ = 291E using equation (18) to calculate the geomagnetic latitude of IPP. • Stand. Implem.: standard implementation with geomagnetic north pole coordinates adopted by Klobuchar, i.e. 𝜙�₀ = 78.3N and 𝜆�₀ = 291E using equation (17) to calculate the geomagnetic latitude of IPP. While in order to quantify the effect of geomagnetic north poles coordinates shift, two another standard implementation was adopted: • Stand. Implem_IGRF.: standard implementation with geomagnetic north pole coordinates adopted by the International Geomagnetic Reference Field for 2010, 𝜙�₀ = 80.1N deg, 𝜆�₀ = 287.8E deg. (from Table 2.). • Stand. Implem_WMM.: standard implementation with geomagnetic north pole coordinates adopted by the World Magnetic Model [Maus et al., 2010] for the year 2010: 𝜙�₀ = 80.02N deg, 𝜆�₀ = 287.79E deg. Aghyas Aljuneidi et al. 16 Figure 5. Ionospheric Delay time (left axis in blue) & Ionospheric Delay range (right axis in green) for all Satellites at UDMC on 01/01/2010 with Standard Implementation of Klobuchar Model with Geomagnetic North Pole coordinates at 78.3N and 291E. 17 Impact of implementation manner in Klobuchar model Clearly to get these figures (Figures 5 and 6) and for a purpose of clarification we had to average the vertical TEC values obtained from each of the satellites, while practically for navigation application, the ionospheric correction is applied to the individual pseudorange measurement. The available RINEX files for this Station near DAMASCUS, Syria is dated on 01/01/2010, i.e., in mid. winter season which explain why the values of the Ionospheric delay range and time are moderate because there no solar activity in that region for that time. From Figures 5 and 6, the Ionospheric delay range and time reach its peak at ~11.58h UT i.e. 14h local time because UDMC is at longitude which means 2.42H in advance, and its minimum values for nighttime. Let’s check the Ionospheric Delay range when using a set of eight improved ionospheric coefficients from [IGS, 2010] (Figure 7). Worth mentioning the remarkable difference between the two ranges, delay range with improved Ionospheric coefficients does not act as the Ideal range delay half cosine function proposed by [Klobuchar, 1975] and there is an approach between day values and night value which is logical because these measurements are done in mid-winter on 01/01/2010. From Figure 8 We noticed that there is a small difference in order of centimeters in the values of the Ionospheric delay range difference between standard implementation and the approximation. While for the coordinates adopted by IGRF (Table 2.) for the year 2010: 𝜙�₀ = 80.1N deg, 𝜆�₀ = 287.1E deg. While for the coordinates adopted by World Magnetic Model [Maus et al., 2010] for the year 2010: ‒𝜙�₀ = 80.02N deg, 𝜆�₀ = 287.79E deg. Worth to note, Figures 9 and 10. That there is also a small difference in order of centimeters in the values of the Ionospheric delay range between standard implementation when we adopt the geomagnetic north pole coordinates whether, IGRF for 2010 or WWM for the 2010 or the counterpart’s values used by Klobuchar in 1975 and the behavior in both case IGRF and WMM is very close to each other. Merited to mention, the correction for the geomagnetic pole has a geophysical meaning, that might not translate directly into an operational application, because these algorithms are not known in their current implementation. Figure 6. Fitting curve of 10th degree of Ionospheric Delay range for all Satellites at UDMC on 01/01/2010 with Standard Implementation of Klobuchar Model with Geomagnetic North Pole coordinates at 78.3N and 291E. Aghyas Aljuneidi et al. 18 Figure 7. Ionospheric Delay range (left axis in blue) with Ionospheric coefficients downloaded from ALMANC and with those improved Ionospheric coefficients (right axis in green) downloaded from [IGS, 2010] for all Satellites at UDMC on 01/01/2010 with Standard Implementation of Klobuchar Model with Geomagnetic North Pole coordinates at 78.3N and 291E. Figure 8. Difference of Ionospheric Delay range for all Satellites at UDMC on 01/01/2010 between Standard and Approximation Implementation of Klobuchar Model with Geomagnetic North Pole coordinates at 78.3N and 291E. 19 Impact of implementation manner in Klobuchar model Figure 9. Difference of Ionospheric range delay in Standard Implementation of Klobuchar Model between Klobuchar Adopted values and IGRF for 2010 of Geomagnetic north pole coordinates. Figure 10. Difference of Ionospheric range delay in Standard Implementation of Klobuchar Model between Klobuchar Adopted values and WMM for 2010 of Geomagnetic north pole coordinates. 6.2.3 Single point range error The inputs were: Prn=18. 𝜙𝑝₀ = 78.3N deg, 𝜆𝑝₀ =291E deg. 𝜙� = 33.5102N deg, 𝜆� =36.2852E deg. 𝐸 = 45.032 deg. 𝐴𝑧 = 43.187 deg. 𝑈𝑇 = 11:21:30. 𝑅� = 6378137.0m. ℎ� = 350 km. Table 3. Standard and approximate implementation of Klobuchar Model. While for the coordinates adopted by IGRF for the year 2010: ‒𝜙𝑝₀ = 80.1N deg, 𝜆𝑝₀ =287.8E deg, and for the coordinates adopted by [Maus et al., 2010] for the year 2010: ‒𝜙𝑝₀ = 80.02N deg, 𝜆𝑝₀ =287.79E deg. Table 4. Standard implementation of Klobuchar Model using IGRF and WMM north pole coordinates for 2010. In Table 3 and 4 the results in single point are compatible with previous Figures, i.e. there is a difference in order of centimeters between Standard implementation and the other implementations. Aghyas Aljuneidi et al. 20 Output Unit Stand. Implem IGRF. Stand. Implem1. WWM 𝜓𝑝𝑝 deg 2.9049 2.9049𝜙𝑝𝑝 deg 35.6044 35.6044𝜆𝑝𝑝 deg 38.7304 38.7304𝜙𝑚 deg 31.5726 31.5348𝐹 1.3470 1.3470𝑇𝑖𝑜𝑛𝑜 sec 11.646e-09 11.653e-09 Range delay m 3.4913 3.4934 VTEC tecu 21.5 21.58 Output Unit Stand. Implem. Approx. Implem. 𝜓𝑝𝑝 deg 2.9049 2.8866𝜙𝑝𝑝 deg 35.6044 35.6149𝜆𝑝𝑝 deg 38.7304 38.7153𝜙𝑚 deg 31.3240 32.1210𝐹 1.347 1.3506𝑇𝑖𝑜𝑛𝑜 sec 11.692e-09 11.575e-09 Range delay m 3.5051 3.4700 VTEC tecu 21.58 21.27 21 Impact of implementation manner in Klobuchar model 6.3 Analytical study In this section, we’ll formulate analytically the impact of imperfect implementation of the coordinate of the geomagnetic north pole in the standard equations of Klobuchar model. Let 𝐼(𝜙�, 𝜆�) the vertical time delay, by using Taylor expansion at the coordinate of geomagnetic north pole (𝜙�, 𝜆�) [Guo et al., 2011], we get: (29) Where 𝐴₁ is constant. Define (𝜙�₁, 𝜆�₁) = (𝜙�₀ + 𝛥𝜙�, 𝜆�₀ + 𝛥𝜆�) thus first order Taylor series expansion of the function can be obtained as follows: (30) Such for the second term of right side of (30), let: (31) Thus (31) becomes (32) Let’s start by the first term of (32), we have: (33) Let 𝑓 = 𝑠𝑖𝑛 (𝜙�(𝜙�₀, 𝜆�₀)), thus (34) if 𝐼������₀ = 𝜙�(𝜙�₀, 𝜆�₀), such we got: (35) = �𝐴(𝜙�₀, 𝜆�₀) · 𝑐𝑜𝑠 � = 𝑐𝑜𝑠 · 𝐴(𝜙�₀, 𝜆�₀) + 𝐴(𝜙�₀, 𝜆�₀) · �𝑐𝑜𝑠 �. 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝜕𝛪(𝜙�₀, 𝜆�₀) 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ Let 𝐼�����₀ = 𝐴(𝜙�₀, 𝜆�₀), and 𝐼�������₀ = 𝑐𝑜𝑠 𝜕 𝜕𝜙�₀𝜕 𝜕𝜙�₀ 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝐼�����₀ · 𝑐𝑜𝑠� � + 𝐼�������₀ 𝐴(𝜙�₀, 𝜆�₀).𝑇 𝑃(𝜙�₀, 𝜆�₀)𝜕𝛪(𝜙�₀, 𝜆�₀) 𝜕𝜙�₀ 𝐼�����₀= 𝐴(𝜙�₀, 𝜆�₀) = �𝛼₁ + 𝛼₂ + 𝛼₃ � �² + 𝛼₄ � �³� = 𝜙� (𝜙�₀, 𝜆�₀) + 𝜙� (𝜙�₀, 𝜆�₀) 𝜙� (𝜙�₀, 𝜆�₀) + 𝜙²� (𝜙�₀, 𝜆�₀) 𝜙� (𝜙�₀, 𝜆�₀). 𝜙�(𝜙�₀, 𝜆�₀) 𝜋 𝜙�(𝜙�₀, 𝜆�₀) 𝜋 𝜙�(𝜙�₀, 𝜆�₀) 𝜋𝜕 𝜕𝜙�₀𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀𝜕 𝜕𝜙�₀ 2𝛼₃ 𝜋²𝛼₂ 𝜋3𝛼₄ 𝜋³ 𝜙�(𝜙�₀, 𝜆�₀) = (𝑠𝑖𝑛⁻¹ �𝑠𝑖𝑛𝜙𝐼PP 𝑠𝑖𝑛𝜙�₀ + 𝑐𝑜𝑠𝜙𝐼PP 𝑐𝑜𝑠𝜙�₀𝑐𝑜𝑠�𝜆𝐼PP ‒𝜆�₀��) = · �𝑠𝑖𝑛𝜙𝐼PP 𝑠𝑖𝑛𝜙�₀ + 𝑐𝑜𝑠𝜙𝐼PP 𝑐𝑜𝑠𝜙�₀𝑐𝑜𝑠�𝜆𝐼PP ‒𝜆�₀�� = · �𝑠𝑖𝑛𝜙𝐼PP 𝑐𝑜𝑠𝜙�₀ ‒ 𝑐𝑜𝑠𝜙𝐼PP 𝑠𝑖𝑛𝜙�₀ 𝑐𝑜𝑠�𝜆𝐼PP ‒𝜆�₀��. 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ 1 √1‒𝑓²1 √1‒𝑓² 𝐼�����₀= 𝐼������₀+ 𝜙� (𝜙�₀, 𝜆�₀)𝐼������₀ + 𝜙²� (𝜙�₀, 𝜆�₀)𝐼������₀.2𝛼₃ 𝜋² 3𝛼₄ 𝜋³𝛼₂ 𝜋 𝐼(𝜙�₁, 𝜆�₁) = 𝐼(𝜙�₀ , 𝜆�₀) + 𝛥𝜙� + 𝛥𝜆�.𝜕𝛪(𝜙�₀, 𝜆�₀) 𝜕𝜙�₀ 𝜕𝛪(𝜙�₀, 𝜆�₀) 𝜕𝜙�₀ 𝐼(𝜙�, 𝜆�) = 𝐴₁ + 𝐴(𝜙�, 𝜆�) 𝑐𝑜𝑠� �𝑇 𝑃(𝜙�, 𝜆�) 𝜕 𝜕𝜙�₀ While for the second term of (32), we have: But Where thus we get: (36) The same altitude for the third term of the right side in (30), it can be written as: (37) where (38) and (39) and (40) Finally, the equation (30) became (41) as follows: (41) Finally, we formulate the deviation in vertical time delay 𝛥𝐼(𝜙�₁, 𝜆�₁) as a function of the shift of geomagnetic north pole coordinates 𝛥𝜙� and 𝛥𝜆�. 𝐼�������₀= 𝑇 � �²· 𝑠𝑖𝑛 � � ·1 𝑃(𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀) � 𝐼������₀+ 𝜙� (𝜙�₀, 𝜆�₀)𝐼������₀+ 𝜙²� (𝜙�₀, 𝜆�₀)𝐼������₀.�𝛽₂ 𝜋 2𝛽₂ 𝜋² 3𝛽₄ 𝜋³ = 𝐼�����₀ · 𝑐𝑜𝑠 + 𝐼�������₀ · 𝐴(𝜙�₀, 𝜆�₀)𝜕𝐼(𝜙�₀, 𝜆�₀) 𝜕𝜆�₀ 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝐼�����₀= 𝐼������₀+ 𝜙� (𝜙�₀, 𝜆�₀)𝐼������₀+ 𝜙²� (𝜙�₀, 𝜆�₀)𝐼������₀.2𝛼₃ 𝜋² 3𝛼₄ 𝜋³𝛼₂ 𝜋 𝐼������₀ = · (𝑐𝑜𝑠𝜙𝐼PP 𝑐𝑜𝑠𝜙�₀𝑠𝑖𝑛(𝜆𝐼PP ‒𝜆�₀ ))1 √1‒𝑓² 𝐼�������₀= 𝑇 � �²· 𝑠𝑖𝑛 � � ·1 𝑃(𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀) � 𝐼������₀+ 𝜙� (𝜙�₀, 𝜆�₀)𝐼������₀+ 𝜙²� (𝜙�₀, 𝜆�₀)𝐼������₀�𝛽₂ 𝜋 2𝛽₃ 𝜋² 3𝛽₄ 𝜋³ 𝐼�������₀= 𝑐𝑜𝑠� � = 𝑐𝑜𝑠 = 𝑠𝑖𝑛� �. � � 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀)𝑇 �₁��₂𝜙�(𝜙�₀, 𝜆�₀)��₃𝜙²� (𝜙�₀, 𝜆�₀)��₄𝜙³� (𝜙�₀, 𝜆�₀)𝑇 �₁��₂𝜙�(𝜙�₀, 𝜆�₀)��₃𝜙²� (𝜙�₀, 𝜆�₀)��₄𝜙³� (𝜙�₀, 𝜆�₀)𝜕 𝜕𝜙�₀ � � = ‒ 𝑇 � �² · ��₁��₂𝜙�(𝜙�₀, 𝜆�₀)��₃𝜙²� (𝜙�₀, 𝜆�₀)��₄𝜙³� (𝜙�₀, 𝜆�₀)�. 𝑇 �₁��₂𝜙�(𝜙�₀, 𝜆�₀)��₃𝜙²� (𝜙�₀, 𝜆�₀)��₄𝜙³� (𝜙�₀, 𝜆�₀)1 �₁��₂𝜙�(𝜙�₀, 𝜆�₀)��₃𝜙²� (𝜙�₀, 𝜆�₀)��₄𝜙³� (𝜙�₀, 𝜆�₀) 𝜕 𝜕𝜙�₀ 𝜕 𝜕𝜙�₀ �𝛽₁ + 𝛽₂ + 𝛽₃ � �² + 𝛽₄ � �³� = 𝐼������₀+ 𝜙� (𝜙�₀, 𝜆�₀)𝐼������₀+ 𝜙²� (𝜙�₀, 𝜆�₀)𝐼������₀. 𝜕 𝜕𝜙�₀ 𝜙�(𝜙�₀, 𝜆�₀) � 𝜙�(𝜙�₀, 𝜆�₀) � 𝜙�(𝜙�₀, 𝜆�₀) �𝛽₂ 𝜋 2𝛽₃ 𝜋² 3𝛽₄ 𝜋³ 𝐼(𝜙�₁, 𝜆�₁) ‒ 𝐼(𝜙�₀, 𝜆�₀) = 𝛥𝐼(𝜙�, 𝜆�) = � � � �.𝐼�����₀· 𝑐𝑜𝑠 𝐼�����₀· 𝑐𝑜𝑠 𝐼�������₀ .𝐴 (𝜙�₀, 𝜆�₀) 𝐼�������₀ .𝐴 (𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝑇 𝑃(𝜙�₀, 𝜆�₀) 𝛥𝜙� 𝛥𝜆� Aghyas Aljuneidi et al. 22 23 Impact of implementation manner in Klobuchar model Bear in mind that one can use (41) to estimate the geomagnetic north pole coordinates when having an accurate ionospheric time delay for example with multi-frequency technique as mentioned in [Guo et al., 2011]. Figure 11. Ionospheric Delay range in Standard Implementation of Klobuchar Model with Geomagnetic North Pole coordinates at 78.3N and 291E and IGRF-2010 with values: 𝛥𝜙� = 1.8° and 𝛥𝜆� = 3.2°. Figure 12. Ionospheric Delay range in Standard Implementation of Klobuchar Model with Geomagnetic North Pole coordinates of IGRF-2010 and WMM-2010 with values: 𝛥𝜙� = 0.08° and 𝛥𝜆� = 0.01°. To obtain the slant time delay, 𝛥𝑇𝑖𝑜𝑛 ,we have: 𝛥𝑇𝑖𝑜𝑛(𝜙�₁, 𝜆�₁) = 𝐹.𝛥𝐼(𝜙�₁, 𝜆�₁) (42) and slant range delay became: 𝛥𝐷𝑖𝑜𝑛(𝜙�₁, 𝜆�₁) = 𝑐.𝐹.𝛥𝐼(𝜙�₁, 𝜆�₁) (43) Where 𝐹 is the Elevation Mapping Function. In the same way, these results are compatible with the previous ones (Figures 11 and 12). 7. Future work By force of circumstances, this article suffers from the lack of comparison between ionospheric experimental data and the various implementations of Klobuchar model. This would require a statistical study on a larger number of days and possibly using stations at various geomagnetic latitudes, comparing with experimental TEC would allow to understand if the proposed corrections are participating to improve the model predictions, or if they are below its intrinsic noise level. 8. Results and discussion Numerical applications were done in aim to estimate the impact of approximate implementation of Klobuchar model, and analytical studies were performed to evaluate the impact of the non-exact values of geomagnetic north pole coordinates, we did many comparisons in these two aspects: • the approximate implementation where the used geomagnetic north pole coordinates are those of the original Klobuchar model, i.e., 𝜆�₀ = 291E and 𝜙�₀ = 78.3N. • the uncertainty in the exact geomagnetic north pole coordinates in Klobuchar model, where the used geomagnetic north pole coordinates are those of IGRF and WWM models for the year 2010. For the first aspect, from Table 3 and Figure 5 and 6, we can conclude that the difference between the standard and approximate implementation of Klobuchar model is in general in the level of centimeters and the most value affected by this implementation is 𝜙� the geomagnetic latitude of IPP because it is directly calculated as a function of geomagnetic north pole coordinates where the biggest approximate was done. Generally, an order of few centimeters’ impact is critical for Precise Point Positioning (PPP) systems, which is not the case for the receivers covered by this study, even though we supposed that nowadays only civil users continue to use a single frequency GPS receiver with no augmentation systems coverage, although this type of low cost receiver plays an essential role in some countries of the third world for application like civil aviation and generally the transport domain, bear in mind that the other constellations like European Galileo or the Chinese BEIDOU offer a little improvements over the performance of GPS concerning the ionospheric delay corrections. For the second aspect, according to Figures 8, 9 and 10, we found immeasurably value when comparing IGRF with WMM implementations, while there is a significant difference between them and the standard implementation that uses geomagnetic north pole coordinates adopted by Klobuchar in his days. 9. Conclusions The approximation of the Klobuchar model is justified due to the deficiencies in digital technique in the 70s, as well as the global knowledge in geomagnetism. In this brief study, a level of centimeters difference was estimated to be the impact of such approximation, which affects the precise positioning systems. Obviously it isn’t the case of the class of devices concerned by this paper. Aghyas Aljuneidi et al. 24 25 Impact of implementation manner in Klobuchar model Regardless of the type and end user of GPS receivers which relay only on Klobuchar Model to correct the ionospheric range error, this study indicates and affirms the usefulness of Klobuchar model as an efficient method for correcting the ionospheric effects till nowadays for a large group of users around the world. The weak point of this formulation is using very few parameters (eight) that included in the broadcast ephemeris for the ionospheric corrections, which is not enough to describe the worldwide behavior of the Earth’s ionosphere. Acknowledgements. The authors express their gratitude to the Syrian National Earthquake Center for the opportunity of using their available RINEX files in this paper. References Allen, C.W. (1976). Astrophysical Quantities, third edition 1973, reprinted with correction 1976, University of London, Athlone Press. BeiDou Navigation Satellite System Signal In Space Interface Control Document (2016). Open Service Signal (Version 2.1) China Satellite Navigation Office, November 2016. Bent, R.B., S.K. Llewellyn (1973). Documentation and description of the Bent ionospheric model. Space and Missile Organisation, Los Angeles, CA, USA. Bilitza, D., D. Altadill, V. Truhlik, V. Shubin, I. Galkin, B. Reinisch, and X. Huang (2017). International Reference Ionosphere 2016: from ionospheric climate to real-time weather predictions, Space Weather, 15, 418–429. Bolaji, O.S., P.A. Izang, O.R. Oladosu, F. Koya, R.F. Fayose, A.R Rabio (2016). Ionospheric Time-delay over Akure Using Global Positioning System Observations, Acta Geophysica, 63, 3, 884-899, doi:10.1515/acgeo-2015- 0018. Chulliat, A., W. Brown, P. Alken, C. Beggan, M. Nair, G. Cox, A. Woods, S. Macmillan, B. Meyer and M. Paniccia (2020). The US/UK World Magnetic Model for 2020-2025: Technical Report, National Centers for Environmental Information, NOAA. doi:10.25923/ytk1-yx35. Davies, K. (1965). Ionospheric Radio Propagation, Dover Publications Inc., Manufactured in the United States of America.European GNSS (Galileo) OPEN SERVICE (2016). Ionospheric correction algorithm for Galileo single frequency users, Issue 1,2, September, 2016. Flock, Warren L. (1983). Propagation Effects on Satellite Systems at Frequencies Below 10 GHz: A Handbook for Satellite Systems, NASA reference publication 1108. Gorbachev, O. A., V. T. Zalutskii, V. B. Ivanov, D. V. Khazanov, A.A. Kholmogorov (2015). Estimating the Quality of GEMTEC Total Electron Content Model in Autonomous GNSS Positioning, Gyroscopy and Navigation, 6, 3, 241–245. Guo, Y., P. Cheng (2011). Inversion of Geomagnetic North-Pole Using GPS Signal and Klobuchar Ionospheric Algorithm, in 2011 Second International Conference on Mechanic Automation and Control Engineering, doi:10.1109/MACE.2011.5988770. IERS, International Earth Rotation and Reference Systems Service (2019). Paris, Bulletin C 57. IGS, International GNSS Service, CODE - Analysis Center. (2010). Klobuchar-Style Ionospheric Coefficients, Index of /CODE/2010,CGIM0010.10N, https://www.aiub.unibe.ch/research/code___analysis_center/klobuchar_style_ionospheric_coefficients/index_en g.html Ivanov, V.B., D.A. Zatolokin, O.A Gorbachev (2017). Comparing Models of Total Electron Content in the Ionosphere for GLONASS, Gyroscopy and Navigation, 8, 4295–299. Janssen, V. (2012). Likely Impact of the Approaching Solar Maximum on GNSS Surveys: Be Alert but Not Alarmed, Proceedings of the 17th Association of Public Authority Surveyors Conference (APAS2012) Wollongong, New South Wales, Australia. Kaplan, E. D. and C. J. Hegart (2006). Fundamentals of Satellite Navigation. Understanding GPS: Principles and Applications 2 ed., Ch. 2, Artech House. Klobuchar, J. (1975). A first-order worldwide ionospheric time delay algorithm. Air Force Cambridge Research Laboratories, Hanscom, AFB, MA, AFCRL-TR-75-0502, AD A018862. https://doi.org/10.1109/MACE.2011.5988770 https://www.aiub.unibe.ch/research/code___analysis_center/klobuchar_style_ionospheric_coefficients/index_eng.html https://www.aiub.unibe.ch/research/code___analysis_center/klobuchar_style_ionospheric_coefficients/index_eng.html Klobuchar, J. (1987). Ionospheric Time-Delay Algorithms for Single-Frequency GPS Users. IEEE Trans. Aerospace Electronic Sys. AES-23, 3, 325-331. Klobuchar, J.A., S. Basu, P. Doherty (1993). Potential limitations in making absolute ionospheric measurements using dual frequency waves from GPS satellites, in Proceedings of the Seventh International Ionospheric Effects Symposium, J. Goodman, cd., Alexandria, VA, May 1993. Kouris, S., K. Polimeris, L. Cander, and L. Ciraolo (2008). Solar and latitude dependence of TEC and slab thicknes, J. Atmosph. Solar-Terr. Phys., 70, 10, 1351 – 1365, doi:10.1016/j.jastp.2008.03.009 Mannucci A.J., B.D. Wilson, C.D. Edwards (1993). A new method for monitoring the Earth’s ionospheric total electron content using the GPS global network. In: Proc. of ION GPS-93. Institute of Navigation, 1323-1332. Maus, S., S. Macmillan, S. Maclean, B. Hamilton, M. Nair, A. Thomson, C. Rollins (2010). The US/UK World Magnetic Model for 2010-2015, NOAA Technical Report NESDIS/NGDC. Nava, B., P. Coïsson, and S.M. Radicella (2008). A new version of the NeQuick ionosphere electron density model, J. Atmosph. Solar-Terr. Phys., 70, 15) 1856 – 1862. Norsuzila, Y., A. Mardina, M. Ismail (2008). Determination of GPS Total Electron Content using Single Layer Model (SLM) Ionospheric Mapping Function, IJCSNS Int. J. Computer Sci. Net. Security, 8, 9. Orus, R., M. Hernanfez Pajares, J.M. Juan, J. Sanz, M. Garcia-Fernandez (2002). Performance of different TEC models to provide GPS ionospheric corrections, J. Atmosph. Solar-Terr. Phys., 64, 18, 2055–2062. Parkinson, B.W. and J.J. Spilker, (1996). Global Positioning System: Theory and Applications, I, Ch. 11, American institute of Aeronautics and Astronautics, Inc. Pro�lss, G.W. (2004). Physics of the Earth’s Space Environment: An Introduction, Springer, New York, doi:10.1007/978- 3-642-97123-5. Sanz, J., Juan, J. and M. Hernandez-Pajares (2013). GNSS Data Processing, I: Fundamentals and Algorithms. ESA Communications, ESTEC. TM-23/1., Noordwijk, the Netherlands. Shuanggen J., C. Estel, X. Feiqin (2014). GNSS Remote Sensing, Theory, Methods and applications, Vol19, Remote Sensing and Digital Image Processing, Springer, doi:10.1007/978-94-007-7482-7. Stepniak,K., Wielgosz, P., Paziewski, J. (2014) Accuracy analysis of the Klobuchar ionosphere model transmitted by the GPS system, The 9th International Conference “Environmental Engineering”, 22–23 May 2014, Vilnius, Lithuania. Thébault, E., C. Finlay, P. Alken, C. Beggan, E. Canet, A. Chulliat, B. Langlais, V. Lesur, F. Lowes, C. Manoj, M.Rother and R. Schachtschneider (2015). Evaluation of candidate geomagnetic field models for IGRF-12, Earth, Planets and Space, 67, 112. Wang, N., Yuan, Y., Li, Z., Huo, X. (2016). Improvement of Klobuchar model for GNSS, Time-Delay Algorithm, single-frequency ionospheric delay corrections, Adv, Space Res., 57, 7, doi:10.1016/j.asr.2016.01.010. Yuan Y., X., Huo J. Ou, K., Zhang, Y. Chai, D. Wen and R. Grenfell R. (2008). Refining the klobuchar ionospheric coefficients based on GPS observations, IEEE Trans. Aerosp. Electron. Syst., 44, 4, 1498_1510. *CO R R E S PO N D I N G A U T H O R: Aghyas ALJUNEIDI, Higher Institute for Applied Sciences and Technology (HIAST), Damascus, Syria, e-mail: aghyas.aljuneidi@hiast.edu.sy © 2021 the Author(s). All rights reserved. Open Access. This article is licensed under a Creative Commons Attribution 3.0 International Aghyas Aljuneidi et al. 26 http://dx.doi.org/10.1016/j.asr.2016.01.010