1. INTRODUCTION
Driven by the increasing interest in and demand for both military and civil applications of satellite-based navigation systems, as well as the nation's quest for independence, China began developing its BeiDou Navigation Satellite System (BDS) in 1990 (Teunissen and Montenbruck, Reference Teunissen and Montenbruck2017). With 16 satellites in orbit, an independent service providing relative and standalone positioning with centimetre-level precision throughout the Asia-Pacific region has already been provided (Gu et al., Reference Gu, Shi, Lou, Feng and Ge2013). More recently, the development of the BDS reached its final stage, the construction of the global service system BDS-3, as scheduled (Yang et al., Reference Yang, Xu, Li and Yang2018).
The deployment of the BDS-3 constellation began in 2015, with the launch of five BDS-3 in-orbit validation satellites, comprising two inclined geosynchronous orbit (IGSO) satellites, I1-S and I2-S, and three medium Earth orbit (MEO) satellites, M1-S, M2-3 and M3-S (Yang et al., Reference Yang, Xu, Li and Yang2018). Based on this demonstration, a wide range of valuable studies have been published, providing a comprehensive analysis of signal performance (Zhang et al., Reference Zhang, Li, Lu, Wu and Pan2017; Yang et al., Reference Yang, Xu, Li and Yang2018), satellite precise orbit determination (POD) (Xie et al., Reference Xie, Geng, Zhao, Liu and Wang2017; Hu et al., Reference Hu, Wang, Wang and Hernández Moraleda2018; Ye et al., Reference Ye, Yuan and Ou2018) and inter-satellite measurement evaluation (Ren et al., Reference Ren, Yang, Zhu and Xu2017). Notably, satellite-induced elevation-dependent pseudo-range multipath error, one of the main factors limiting the performance of the previous generation of BDS satellites, is no longer present in the signals of the new BDS-3 satellites. After the positive assessment of the design of BDS-3 based on the demonstration system, the first two operational new-generation satellites were launched (Space Vehicle numbers 19 and 20) on 5 November 2017, marking the beginning of China's construction of a global BDS system (Yang et al., Reference Yang, Xu, Li and Yang2018). In the last 8 months, the construction of BDS-3 has increased in pace, with 11 more satellites launched.
Of the various benefits of the signals of the new generation of satellites, the contribution of BDS-3 to precise point positioning (PPP) is of special interest, as PPP has been proven to efficiently provide centimetre-level positioning with standalone stations near and above the Earth (Zumberge et al., Reference Zumberge, Heflin, Jefferson, Watkins and Webb1997; Kouba and Héroux Reference Kouba and Héroux2001). PPP has also been widely used in a number of geodesy and geodynamics studies investigating, for instance, tectonic plate motion (Heki et al., Reference Heki, Miyazaki, Takahashi, Kasahara, Kimata, Miura, Vasilenko, Ivashchenko and An1999), positioning and navigation (Shi et al., Reference Shi, Zheng, Lou, Gu, Zhang, Dai, Li, Guo and Gong2017), early warning systems for earthquakes and tsunamis (Shi et al., Reference Shi, Lou, Zhang, Zhao, Geng, Wang, Fang and Liu2010) and Global Navigation Satellite System (GNSS)-based space weather (Zheng et al., Reference Zheng, Lou, Gu, Gong and Shi2017).
With the modernisation of the Global Positioning System (GPS) and Globalnaya Navigazionnaya Sputnikovaya Sistema (GLONASS), as well as the progress of the newly developed BDS and Galileo, multi-constellation PPP has received increasing attention from those seeking to improve positioning performance, and especially to accelerate PPP convergence. Cai and Gao (Reference Cai and Gao2007) presented a combined GPS/GLONASS PPP solution soon after the recovery of GLONASS. However, their introduction of GLONASS observations failed to improve PPP performance. Recently, Shi et al. (Reference Shi, Yi, Song, Lou, Yao and Zhang2013) argued that as the GLONASS signal is division multiple access modulated, the receiver inter-frequency bias (IFB) should be carefully calibrated in GLONASS-only PPP. The results suggested an improvement of almost 50% during the initialisation period. BDS-2 was designed as a regional navigation system with five geostationary orbit (GEO) and five IGSO satellites. Therefore, although it can provide centimetre-level PPP results independently, its convergence time is usually much longer than that of GPS or GLONASS, so its potential to provide multi-constellation PPP is limited (Ge et al., Reference Ge, Zhang, Jia, Song and Wickert2012). It should be noted that the ionosphere-free (IF) combination observation was typically enabled in the abovementioned studies.
To access the full capabilities of multi-frequency signals, an undifferenced and uncombined observation model has been promoted and demonstrated as an efficient method of simultaneously estimating ionosphere delay and multi-frequency phase and code biases (Schönemann et al., Reference Schönemann, Becker and Springer2011; Gu et al., Reference Gu, Shi, Lou, Feng and Ge2013; Gu et al., Reference Gu, Shi, Lou and Liu2015b; Zhang et al., Reference Zhang, Teunissen, Yuan, Zhang and Li2018). Based on the deterministic plus stochastic ionosphere model for GNSS (DESIGN), Lou et al. (Reference Lou, Zheng, Gu, Wang, Guo and Feng2015) presented a comprehensive analysis of multi-constellation PPP for both single- and dual-frequency undifferenced and uncombined observations. Zhou et al. (Reference Zhou, Dong, Ge, Li, Wickert and Schuh2018) compared PPP performance in terms of convergence and accuracy with different estimation schemes of GLONASS code IFBs, and argued that the relationship between GLONASS code IFBs and frequency number may not be strictly linear or quadratic.
As a global navigation system, BDS-3 consists of 2 IGSO and 10 MEO satellites, as listed in Table 3 (later), and is thus expected to provide better PPP performance and more rapid convergence than BDS-2. To reveal the full potential of the BDS in positioning and navigation, this paper presents an initial evaluation of the PPP performance of mixed BDS-2 and BDS-3 satellites, using an undifferenced and uncombined model.
The remainder of this paper is organised as follows. First, a general PPP model with undifferenced and uncombined observations from both BDS-2 and BDS-3 satellites is proposed. Second, based on this model, PPP performance is analysed in detail in terms of accuracy and convergence. The final section concludes the paper.
2. DATA PROCESSING MODEL
As suggested by Gu et al. (Reference Gu, Lou, Shi and Liu2015a) and Zhao et al. (Reference Zhao, Wang, Gu, Zheng, Shi, Ge and Schuh2018), a basic model of the BDS with undifferenced and uncombined observations as individual measurements can be written as follows:
$$\left.\begin{array}{@{}c@{}} P_{r,f}^s =\rho_r^s +t_r -t^s+\alpha_r^s \cdot T_z +\beta_f \cdot I_{r,Z}^s +b_{r,f}^s +\varepsilon_P \\ \Phi_{r,f}^s =\rho_r^s +t_r -t^s+\alpha_r^s \cdot T_z +\beta_f \cdot I_{r,Z}^s -N_{r,f}^s +\varepsilon_{\Phi} \end{array}\right\}$$
where 
$P_{r, f}^{s}$ and 
$\Phi_{r, f}^{s}$ are the pseudo-range and carrier phase from receiver r to satellite s on frequency f in length units, respectively, ρ is the geometric distance with antenna phase centre corrections and phase windup corrections applied, t r and t s are the receiver and satellite clock error, respectively, T z is the zenith tropospheric delay (ZTD), which can be converted to slant delays with the mapping function α, 
$I_{r, Z}^{s}$ denotes the zenith total electron content with the frequency and mapping function M s-dependent factor, i.e., βf = (40 · 3/f 2)(1/M s), N is the float ambiguity by definition in cycle units with the corresponding wavelength λf and 
$b_{r, f}^{s}$ and 
$b_{f}^{s}$ denote the frequency-dependent code bias for receiver and satellite, respectively.
It should be emphasised that the superscript s here may indicate either BDS-2 or BDS-3 satellites, i.e., 
$s\in \left(\begin{smallmatrix} s_{\text{BDS-2}} &s_{\text{BDS-3}} \end{smallmatrix}\right)$. As signal frequency differs between the two types of satellite, the frequency-dependent code bias terms are given special attention. The code bias for BDS-2 is generally expressed as follows:
$$b_{r,f}^s =b_{r,f} -b_f^s +b_f^s (E)$$ The satellite elevation-dependent variation 
$b_{f}^{s} ( E) $ can be removed in advance with a third-order polynomial model, as suggested by Lou et al. (Reference Lou, Gong, Gu, Zheng and Feng2016):
$$b_f^s (E)=a_1 \cdot E+a_2 \cdot E^2+a_3 \cdot E^3$$
where E denotes the elevation in radius and 
$a_{i} \left(i\in \left(\begin{smallmatrix} 1 & 2 & 3 \\ \end{smallmatrix}\right)\right)$ denotes the frequency-dependent coefficients of the GEO, IGSO and MEO satellites for BDS-2, as listed in Table 1 (Lou et al., Reference Lou, Gong, Gu, Zheng and Feng2016).
Table 1. Coefficients of GEO, IGSO and MEO in code bias variation model for different signal frequencies.
We ignore the variation in BDS-3 code bias because it is reduced to a negligible level, as mentioned previously. The precise satellite orbit and clock products obtained by Wuhan University (WHU) (Xie et al., Reference Xie, Geng, Zhao, Liu and Wang2017) are applied to Equation (1) for BDS-3 satellites. The precise clock products are estimated with a B1/B3 IF combination in which the satellite clock and code bias are lumped together, i.e.,
$$t^s:=t^s+\frac{f_1^2 \cdot b_1^s}{f_1^2 -f_3^2}-\frac{f_3^2 \cdot b_3^s}{f_1^2 -f_3^2}$$Therefore, we have the pseudo-range observation
$$\left.\begin{array}{@{}l@{}} P_{r,1}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_1 \cdot I_{r,Z}^s+b_{r,1} +\dfrac{f_3^2 \cdot b_{13}^s }{f_1^2 -f_3^2 }+\varepsilon_P \\ P_{r,2}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_1 \cdot I_{r,Z}^s+b_{r,2} +\dfrac{f_3^2 \cdot b_{13}^s }{f_1^2 -f_3^2 }+b_{12}^s +\varepsilon_P \\ P_{r,2}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_1 \cdot I_{r,Z}^s+b_{r,2} +\dfrac{f_3^2 \cdot b_{13}^s }{f_1^2 -f_3^2 }+b_{12}^s +\varepsilon_P \end{array}\right\}$$
where 
$b_{12}^{s}=b_{1}^{s}-b_{2}^{s}$ denotes the differential code bias (DCB) of frequency 1 and frequency 2, and 
$b_{13}^{s}=b_{1}^{s} -b_{3}^{s}$ denotes the DCB of frequency 1 and frequency 3. Although the DCB of BDS-2 (C01–C14) can be retrieved from the International GNSS Service (IGS), no official product for BDS-3 (C19–C22 and C27–C30) satellite DCB correction is available. Thus, we use the relevant values obtained by WHU, as listed in Table 2.
Table 2. Correction of satellite DCB for B1B2 and B1B3 (unit: ns).
Similarly, we define b r,12 = b r,1 − b r,2 and b r,13 = b r,1 − b r,3, and denote t r: = t r + b r,1. Next, Equation (5) is simplified the following equation by correcting the satellite DCB for B1B2 and B1B3:
$$\left.\begin{array}{@{}c@{}} P_{r,1}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_1 \cdot I_{r,Z}^s+\varepsilon_P \\ P_{r,2}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_1 \cdot I_{r,Z}^s+b_{r,12} +\varepsilon_P \\ P_{r,3}^s =\rho_r^s +t_r +\alpha_r^s \cdot T_z +\beta_3 \cdot I_{r,Z}^s+b_{r,13} +\varepsilon_P \end{array}\right\}$$Importantly, as not all of the BDS precise satellite ephemeris data are available in WHU's BDS-3 product, the GBM product from the German Research Centre for Geosciences (GFZ) is used for BDS-2 precise orbit and clock correction. In this case, the time reference differs between the GBM and WHU satellite clock solutions, so two receiver clock parameters for BDS-2 and BDS-3 satellites, respectively, are estimated simultaneously in Equation (6).
In addition, for the ionosphere delay 
$I_{r}^{s}$, we adopt the DESIGN method promoted by Zhao et al. (Reference Zhao, Wang, Gu, Zheng, Shi, Ge and Schuh2018):
$$I_{r,Z}^s =a_0 +a_1 dL+a_2 dL^2+a_3 dB+a_4 dB^2+r_r^s$$
where a 0 is the average value of ionospheric delay over the station; a 1, a 2, a 3 and a 4 are the coefficients of the two second-order polynomials along the east–west and south–north direction, respectively; 
$r_{r}^{s}$ is the residual ionospheric effect for each satellite; and dL and dB are the differences in longitude and latitude, respectively, between the ionospheric pierce point and the approximate location of the station. Furthermore, the daily variation of the deterministic part, i.e., a i, is expressed as a Fourier series with fixed frequency, i.e., 
$f_{j} \in \left(\begin{smallmatrix} 1/24 & 1/12 & 1/8 \end{smallmatrix}\right)$, and the amplitude, i.e., x 0, y j, z j, is expressed as a daily constant updated every 24 h:
$$a_i (t)=x_0 +\sum_{j=1}^3 ( y_j \cdot \sin ( 2\pi tf_j) +z_j \cdot \cos ( 2\pi tf_j) ) $$The stochastic part is estimated as a random walk with the variogram summarised in Equations (9)–(12):
$$\gamma (h)=\begin{cases} \displaystyle c_s \cdot \left(\left(\frac{3}{2}\right)\left(\frac{h}{a_s}\right)- \left(\frac{1}{2}\right)\left(\frac{h}{a_s}\right)^3\right),& 0\le h< a_s \\ c_s,& h\ge a_s \end{cases}$$where the maximum correlation distance a s = 9, 000 s, h represents the distance (time interval in our case) between two samples and c s varies with geomagnetic latitude B as a Gaussian function:
$$c_s=c_{s,\min} +c_{s,\max} \cdot e^{-(|B|-15)^2/128}$$The variation in ionospheric activity is modelled by c s, min and c s, max with the Epstein function as
$$\left.\begin{array}{@{}c@{}} c_{s,\min}=0.6+(7.5-0.6)\cdot (1/(1+e^x))\\ c_{s,\max}=6.0+(75.0-6.0)\cdot (1/(1+e^x)) \end{array}\right\}$$with x dependent on the sunspot number s n:
$$x(s_n)=\frac{s_n -100}{20}$$ In addition, owing to the progress of ionosphere delay modelling, an a priori ionosphere delay correction 
$I_{r, {\rm corr}}^{s}$ is usually available from, for instance, global ionosphere maps (GIMs) or regional ionosphere products, which can be used as pseudo-observations for each line of sight:
$$I_{r,{\rm corr}}^s =a_0 +a_1 dL+a_2 dL^2+a_3 dB+a_4 dB^2+r_r^s +\varepsilon_I$$where ɛI is the corresponding noise.
3. EXPERIMENTAL ANALYSIS
To demonstrate the performance of the new generation of BDS satellites, both BDS-2 and BDS-3 observations are collected in the PPP experiment and analysed in terms of positioning accuracy and convergence. In this section, we begin with an overview of the BDS-3 system, then introduce the data and strategy, and finally present the results.
3.1. System overview
Currently, the global BDS system consists of 2 IGSO and 10 MEO new-generation satellites. Table 3 presents the status of BDS-3 as of 25 May 2018. IGSO01-S, IGSO02-S, MEO01-S and MEO02-S are experimental satellites.
Table 3. Status of BDS-3 satellites (as of 25 May 2018)
The corresponding track of sub-satellite points is drawn in Figure 1. The corresponding BDS satellite visibility is presented in Figure 2. As illustrated by Figure 1 and 2, more than 10 BDS satellites are available for the region, from approximately 55°E to 150°E and 40°S to 45°N, with a mask angle of 5°. The four stations involved in the experiment, i.e. the International GNSS Monitoring and Assessment System (IGMAS) stations CANB, DWIN, KNDY and PETH, are also shown in Figure 1. These stations are selected because they have both BDS-2 and BDS-3 tracking capability and are located in the Asia-Pacific region, where the BDS can provide positioning independently.
Figure 1. Ground track of sub-satellite points for BDS-2 (in blue) and BDS-3 (in red) as of 25 May 2018. The four stations involved in the experiment, i.e., CANB, DWIN, KNDY and PETH, are plotted with yellow stars
Figure 2. Satellite visibility of BDS-2 and BDS-3 as of 25 May 2018 (http://www.beidou.gov.cn/xt/jcpg/201805/t20180530_14900.html).
3.2. Experiment details
The processing strategy is realised based on the Fusing in GNSS (FUSING) software developed at WHU for validation. At present, the FUSING software is capable of multi-GNSS, i.e., GPS, GLONASS, BDS and Galileo, positioning; POD of GNSS; high-frequency satellite clock estimation; ionosphere and troposphere modelling; and multi-frequency precise positioning (Shi et al., Reference Shi, Guo, Gu, Yang, Gong, Deng, Ge and Schuh2018; Zhao et al., Reference Zhao, Wang, Gu, Zheng, Shi, Ge and Schuh2018; Yang et al., Reference Yang, Gu, Gong, Song, Lou and Liu2019). Measurements are collected from the 4 stations shown in Figure 1, which have both BDS-2 and BDS-3 tracking capability, for the period DOY 133–DOY 168, 2018, with an interval of 30 s. It should be noted that the data obtained from these stations are unavailable for some days. For more details of the PPP strategy, please refer to Table 4.
Table 4. Details of the PPP strategy.
3.3. Analysis of results
Compared with the reference coordinates based on the weekly static PPP solution, Table 5 lists the overall root mean square (RMS) of the PPP at each station. Note that unless otherwise stated, the statistics presented in the experiment are based on the full 24 h period, i.e., samples both before and after convergence are included.
Table 5. RMS of PPP for stations averaged over the period DOY 133 to DOY 168, 2018.
As demonstrated by Table 5, the PPP at the stations DWIN, KNDY and PETH is more accurate than that at CANB, regardless of whether BDS-2 alone or BDS-2 and BDS-3 are involved. This is reasonable, as DWIN, KNDY and PETH are located in the ‘core’ service region of the current BDS system, as suggested by Figures 1 and 2. Comparing the BDS-2-only solution with the BDS-2 and BDS-3 solution reveals that positioning accuracy is significantly improved by the addition of the new BDS-3 satellites, with improvements of 28%, 21% and 5% for the up (U), north (N) and east (E) components, respectively. In addition, the three-dimensional RMS values for different days are plotted in Figure 3. Compared with BDS-2 alone, BDS-2 and BDS-3 decrease the RMS of PPP from 0·56 m to 0·39 m, 0·31 m to 0·25 m, 0·31 m to 0·23 m and 0·36 m to 0·27 m for CANB, DWIN, KNDY and PETH, respectively.
Figure 3. Three-dimensional PPP RMS of BDS-2 (black bars) and BDS-2 and BDS-3 (red bars) at sites CANB, DWIN, KNDY and PETH from DOY 133 to DOY 168, 2018.
In addition to accuracy, convergence time is an important factor limiting the applications of PPP. To evaluate the initialisation performance, the 68% quantile convergence series of the experimental period for PPP are plotted as black lines for BDS-2 alone and red lines for BDS-2 and BDS-3 in Figure 4. Suppose that the set of absolute values of positioning difference is 
$X_{t}=\left(\begin{matrix} x_{1} & x_{1} &\cdots & x_{n} \end{matrix} \right)$ for epoch t; then the 68% quantile q t for this epoch is defined as
Figure 4. Upper panel: 68% quantile horizontal convergence series of BDS-2 (black line) and BDS-2 and BDS-3 (red line) PPP for DOY 133 to DOY 168, 2018 over a 6-hour pass. Bottom panel: 68% quantile vertical convergence series of BDS-2 (black line) and BDS-2 and BDS-3 (red line) PPP for DOY 133 to DOY 168, 2018 over a 6-hour pass.
$$P(x\le q_t)=68\% (x\in X_t),$$where P denotes the probability function. Obviously, q t can be derived by sorting the set X t incrementally. As shown in Figure 4, with more MEO satellites the initialisation of PPP significantly speeds up. With the threshold of 0·15 m as the indicator, the time needed for convergence is about 2·3 h in both the horizontal and vertical directions for BDS-2 PPP. For BDS-2 and BDS-3 PPP, the time is shortened to about 1·5 h and less than 1 h, respectively. As shown in the figure, the improvement in the vertical direction is more significant than that in the horizontal direction, and a slight improvement is found even after convergence in the vertical direction. This is reasonable, because as suggested by Zheng et al. (Reference Zheng, Lou, Gu, Gong and Shi2017), the poor geometry of the BDS-2 satellite constellation results in a high correlation between the ZTD and the vertical component, limiting the corresponding accuracy. As the precision in the horizontal direction remains roughly the same, as suggested by the upper panel in Figure 4, the improvement in RMS indicated in Table 5 can be attributed mainly to the convergence period for the horizontal direction.
To provide more details of the PPP series, the UNE (Up/North/East) errors of PPP in a single day are further presented as an example, and Figures 5–8 compare the results for CANB, DWIN, KNDY and PETH, respectively. The corresponding number of satellites involved in the solution is also plotted. With about eight more new-generation (BDS-3) MEO satellites, about two additional satellites can be tracked for each station and each epoch.
Figure 5. Positioning difference series at CANB for BDS-2 PPP (upper panel) and BDS-2 and BDS-3 PPP (bottom panel) for DOY 142, 2018. The corresponding number of satellites is plotted as a grey line.
Figure 6. Positioning difference series at DWIN for BDS-2 PPP (upper panel) and BDS-2 and BDS-3 PPP (bottom panel) for DOY 142, 2018. The corresponding number of satellites is plotted as a grey line.
Figure 7. Positioning difference series at KNDY for BDS-2 PPP (upper panel) and BDS-2 and BDS-3 PPP (bottom panel) for DOY 142, 2018. The corresponding number of satellites is plotted as a grey line.
Figure 8. Positioning difference series at PETH for BDS-2 PPP (upper panel) and BDS-2 and BDS-3 PPP (bottom panel) for DOY 144, 2018. The corresponding number of satellites is plotted as a grey line.
The above analysis presents the performance of kinematic PPP with BDS-2 and BDS-3 satellites. In addition, Table 6 lists the repeatability of daily static PPP in terms of standard deviation (STD) from DOY 133 to DOY 168, 2018. Obviously, the contribution of additional BDS-3 satellites to positioning repeatability is rather limited. In addition, the STD is around 24 mm for the vertical direction and 10 mm for north and east regardless of whether BDS-2 alone or BDS-2 and BDS-3 are involved in the daily static PPP solution.
Table 6. Repeatability of daily static PPP from DOY 133 to DOY 168, 2018.
4. CONCLUSION
This paper has presented an investigation of PPP with measurements from both BDS-2 and BDS-3 satellites. To simultaneously process observations of different signals, i.e., B1/B1C, B2/B2a and B3C, an undifferenced and uncombined observation model with ionosphere delay constrained by DESIGN has been used as the basic model. Using the FUSING software developed at WHU, measurements collected at four stations in the Asia-Pacific region with both BDS-2 and BDS-3 tracking capability have been compared across the experimental period of DOY 133 to DOY 168, 2018.
The results suggest that compared with the reference coordinates based on a weekly static PPP solution, the new model significantly improves positioning precision. For BDS-2 PPP, the average RMS values are 0·29 m, 0·17 m and 0·17 m for the up, north and east directions, respectively; when both BDS-2 and BDS-3 are used, these values decrease to 0·21 m, 0·13 m and 0·16 m. Analysis of the convergence time also reveals that the initialisation of PPP is significantly faster for BDS-2 and BDS-3 than for BDS-2 alone. With the threshold of 0·15 m as an indicator, the time needed for convergence is about 2·3 h in both the horizontal and vertical directions for BDS-2 PPP. This time is shortened to about 1.5 h and less than 1 h, respectively, for the BDS-2 and BDS-3 solution.
It should be emphasised that this paper is only an initial evaluation of the PPP performance provided by BDS-3. As China intends to establish its global satellite navigation system by 2020, BDS-3 satellites will be launched very frequently this year, inevitably improving PPP performance.
ACKNOWLEDGEMENTS
This study was supported by the State Grid Corporation Science and Technology Project (Grant No. 5211XT180047) and the State Key Research and Development Program (2017YFB0503401). The authors thank the anonymous reviewers for their valuable comments and the IGS, the iGMAS and WHU for providing data.}
