This paper introduces a phase tracking method for planetary radio science research with computational algorithm implemented fo r NVIDIA GPUs. In contrast to the phase-locked loop (PPL) phase counting method used in traditional Doppler data processing, this method fits the tracking data signal into the shape expressed by the Taylor polynomial with optimal phase and amplitude coefficients. The Differential Evolution (DE) algorithm is employed for polynomial fitting. In order to cope with high computational intensity of the proposed phase tracking method, the graphics processing units (GPUs) are employed. As a result, the method estimates the instantaneous phase, frequency, derivative of frequency (line-of-sight acceleration) and the total count phase of different integration scales. This data can be further used in planetary radio science research to analyze the planetary occultation and gravitational fields. The method has been tested on MEX (Mars Express, ESA) and Chang'E 4 relay satellite (China) tracking data. In a real experiment with 400K data block size and $\sim$80,000 DE solver objective function evaluations we were able to acheive the target convergence threshold in 6.5 seconds and do real-time processing on NVIDIA GTX580 and 2$\times$ NVIDIA K80 GPUs, respectively. The precision of integral Doppler (60s) is 2 mrad/s and 4 mrad/s for MEX(3-way) and Chang'E 4 relay satellite(3-way) respectively.

Electron density profiles of the Venus' ionosphere are inverted from the Venus Express (VEX) one-way open-loop radio occultation experiments carried out by Shanghai 25 m antenna from November 2011 to January 2012 at solar maximum conditions and by New Norcia 35 m antenna from August 2006 to June 2008 at solar intermediate conditions. The electron density profile (from 110 km to 400 km) retrieved from the X-band egress observation at Shanghai station, shows a single peak near 147 km with a peak density of about $2 \times 10^4 \rm{cm}^{-3}$ at a solar zenith angle of 94$^{\circ}$. As a comparison, the VEX radio science (VeRa) observations at New Norcia station were also examined, including S-, X-band and dual-frequency data in the ingress mode. The results show that the electron density profiles retrieved from the S-band data are more analogous to the dual-frequency data in the profile shape, compared with the X-band data. Generally, the S-band results slightly underestimate the magnitude of the peak density, while the X-band results overestimate that. The discrepancy in the X-band profile is probably due to the relatively larger unmodeled orbital errors. It is also expected that the ionopause height is sensitive to the solar wind dynamical pressure in high and intermediate solar activities, usually in the range of 200 km - 1000 km on the dayside and much higher on the nightside. Structural variations ("bulges" and fluctuations) can be found in the electron density profiles in intermediate solar activity, which may be caused by the interaction of the solar wind with the ionosphere. Considerable ionizations can be observed in the Venus' nightside ionosphere, which are unexpected for the Martian nightside ionosphere in most cases.