Investigating the variation of the Sun’s visual shape, atmospheric refraction and Einstein’s special relativity considered

By experimental measurements and theoretical analyses, this paper investigates the variation of the Sun’s visual shape and figures out the reasons for the variation of its shape. First, the method of image processing, the method of moments and the least-square method are combined to perform experimental measurements and calculations, and the features of the Sun’s visual shape are extracted from the photos of the Sun. Second, theoretical analyses are conducted based on atmospheric refraction and the Einstein’s special relativity theory. The relationship model is established between the zenith and azimuth angles of the Sun, the velocity of the Sun relative to the Earth, and the observation time and position; the refraction index of the atmosphere is expressed as a function of altitude and wavelength of light; an iterative algorithm is constructed to trace rays of light in the atmosphere; a set of formulas is derived to determine the contraction ratio and contraction direction of the Sun’s visual shape. Finally, the theoretical and experimental results are compared; their relative errors are less than 0.3%, which verifies the theoretical analyses. Both theoretically and experimentally, this research proves that the Sun’s visual shape is an ellipse; its shape variation mainly results from atmospheric refraction effects; and the length contraction effect of the Einstein’s special relativity also contributes a little, except at the time of sunrise and sunset.


Introduction
When you enjoy sunrise and sunset, do you pay attention to the variation of the Sun's visual shape? Comparatively speaking, the Sun's image at sunrise or sunset is flatter than that at noon; however, although the Sun's image at noon is the roundest on a day, it still appears to be an ellipse and not a perfect circle. Very few research papers systemically explain the theoretical reasons for the variation of the Sun's visual shape, and only a few very short articles on Wikipedia and in popular science magazines simply outline the natural phenomenon based on atmospheric refraction.
The effects of atmospheric refraction have been studied widely in the fields of optical communications and weather forecast. For the free-space optical communication systems, by applying astronomical refraction formulas and considering meteorological conditions on ground, Karin and Florian (2004) have simulated the refraction angles of solar rays in the atmosphere, and calculated the vertical deviation of solar beams at various wavelengths. For the earth-to-satellite laser communications, Xiang (2008) analyzed the effect of atmospheric chromatic dispersion on the pointing error of the uplink laser beam. For the sake of improving accuracy of satellite laser ranging, Yuan et al. (2011) used ray tracing to compute the various optical paths caused by atmosphere refraction, and gave a regional distribution of the optical path difference in China. Jiang et al. (2013) investigated the relationship between the optical paths and the satellite zenith angles, and proposed an atmospheric refraction compensation scheme for the different satellite zenith angles. In order to find out about the effect of atmospheric refraction on the propagation of weather radar beams, Wang et al. (2018) explored the vertical variation of the refraction index in the first kilometer of the atmosphere with regional climate and topography. Balal and Pinhasi (2019) evaluated the effect of atmospheric refraction and absorption on the propagation of millimeter and sub-Communicated by: H. Babaie millimeter wavelengths from land to satellite. Chaim and Hall (2000) made use of ray tracing to determine atmospheric refraction, and combined the digital terrain model to calculate the visual sunrise and sunset times at some cities in Israel. Kambezidis and Papanikolaou (1990) researched the relationship between solar position and atmospheric refraction. Kambezidis and Tsangrassoulis (1993) proposed a new correction of right ascension according to solar position. Kambezidis (1997) established a set of appropriate spherical trigonometric formulas to estimate sunrise and sunset hours by considering flat and complex terrains and atmospheric refraction. Instead of spherical trigonometry, Sproul (2007) used vector analysis to research the position relationship between the Sun and the Earth.
This paper investigates the variation of the Sun's visual shape on a day, and figures out the reasons for the shape variation by massive experimental measurements and systemic theoretical analyses. First, hundreds of photos of the Sun were taken from sunrise to sunset; an image processing was performed to extract features of the Sun's visual shape from the photos. The method of moments and the least-square method were combined to fit the periphery of the Sun in the processed images; and a set of formulas was derived to calculate the feature parameters of the elliptic Sun. Error analyses showed that the relative measurements accuracy was about 0.023%, the standard deviation of the fitting curve of the ellipse was only 4 pixels. The experimental results proved that the Sun's visual shape can be accurately approximated by an ellipse. Second, in order to investigate the reasons for the variation of the Sun's visual shape, the atmospheric refraction effects were first researched to simulate the influence on the Sun's visual shape at different observation times and positions. Based on Kepler's laws, the relative position and velocity between the Sun and the Earth were analyzed, and a set of formulas was derived to calculate the zenith and azimuth angles of the Sun for different observation positions and times; according to the differential equation of light propagation in an inhomogeneous medium, an iterative algorithm was developed to trace the rays of light in the atmosphere. Meanwhile, because the refraction index of the atmosphere mainly depends on air density, it is expressed as a function of altitude and wavelength of light. Under this expression, it is proved theoretically that the trajectory of sunlight is a planar curve. Nevertheless, as a result of atmospheric refraction, the Sun's visual shape contracts only in the zenith direction, resulting in an elliptic Sun in the observer's eyes. In addition, the length contraction phenomenon in the Einstein's special relativity theory due to the relative movement between the Sun and the observer was also investigated to simulate its effects on the Sun's visual shape at different observation times and positions. Because the relative velocity between the Sun and the observer varies with time, especially the direction of the motion, the effect of the length contraction phenomenon on the Sun's visual shape also varies with time; therefore, the analysis included the contraction ratio and the contraction direction of the Sun's visual shape, which describe a ratio of the minor to the major axis and the direction of semi-minor axis of the Sun's visual shape. A relationship between the length contraction effect, observation positions and times was established. Although the effect of the Einstein's special relativity theory also makes the Sun's visual shape contract in one direction and results in an elliptic Sun, this differs from the atmospheric refraction effect, because the contraction direction caused by the Einstein's special relativity theory varies with time. Comprehensively considering the effects of atmospheric refraction and the Einstein's special relativity theory on the Sun's visual shape, although different, we can prove that the Sun's visual shape is still an ellipse. Therefore, a formula calculating the shape parameter of the ellipse was derived. Finally, the observation position was chosen to be Dalian, and the observation time mid-December 2018. Then, the variation of the Sun's visual shape was simulated, including two cases: one with atmospheric refraction effects only, and another with both atmospheric refraction and the Einstein's special relativity theory. By comparison, it can be found that the relative errors of the simulation results and the experimental data are less than 0.3%. These results show that the reason for the variation of the Sun's visual shape is mainly atmospheric refraction; and at times except of sunrise and sunset, the length contraction effect of the Einstein's special relativity theory also contributes, but little.

Experimental measurements and calculations
In order to investigate the variation of the Sun's visual shape, a high-resolution digital camera was used to take more than 300 photos of the Sun from sunrise to sunset. The time chosen was before the winter solstice in the year 2018; the position was Dalian, a city in northeast China. In Fig. 1, three photos are given, the projection plane is perpendicular to the direction from the camera to the Sun, and the horizontal direction of the photo is about parallel to the sea level. Here, the Sun's visual shape is assumed an ellipse. In the later Sections, the assumption will be verified by experimental results and theoretical analyses. A shape parameter describing the ellipse is defined as k = b/a, where a stands for the dimension of its semi-major axis, b the dimension of its semi-minor axis. When the sun rises or sets, the vertical direction of the photo or the direction of b is about perpendicular to the sea level. A bigger k means that the Sun's visual shape is rounder. Seeing from Fig. 1, at sunrise, noon and sunset, the major and minor axes of the Sun are basically along the horizontal and vertical directions of the images.
For extracting the shape parameter k from the Sun's photos, image processing, and the methods of moments and leastsquares were employed. First, the Sun's image expressed by RGB was converted to a binary image expressed by two luminance values of 0 and 1; then the periphery of the Sun was extracted from the binary image by removing the interior pixels; next, considering the Sun's binary image as an ellipse, the method of moments was used to evaluate its feature parameters, such as the semi-major axis, the semi-minor axis, the coordinates of the center, and the angle between the major axis of the ellipse and the horizontal direction of the image, which were defined as the initial values of the iterative algorithm for fitting an ellipse to the periphery of the Sun. Finally, the leastsquares method or the iterative algorithm was applied to fit an ellipse to the periphery of the Sun and compute the feature parameters of the Sun's elliptical image.
Processing the Sun's image taken by a camera A photo of the Sun taken by a camera is a true color image expressed in RGB. In order to extract the shape parameter k from the Sun's photo, the Sun's image should be processed.
First, the true color image was converted to a grayscale image by removing the hue and saturation information, while retaining luminance. The mapping relationship between luminance and RGB can be described as where Y denotes the luminance value at every pixel, ranging from 0 to 1, Y = 1 denotes white, Y = 0 stands for black. Then, in order to reduce noise in the grayscale image, a 3-by-3 neighborhood median filtering was performed. Fig. 2a is a true color image taken by a camera, and Fig. 2b is a processed image after median filtering. Next, according to the contrast of Fig. 2b, a luminance threshold e Y was estimated, and the luminance values of all pixels with Y ≥ e Y were replaced with the value 1, while the luminance values of the other pixels were set to 0. In this way, Fig. 2b is converted to a binary image expressed by two luminance values of 0 and 1. Besides, in order to eliminate isolated points and burrs on the periphery, the morphological operations were applied to the binary image, including image erosion and dilation. During the image erosions and dilations, a 3-by-3 matrix SE was used as a structuring element object After performing 3 times erosions, 3 times dilations were implemented. Fig. 2c is the binary image after performing erosions and dilations. The last step was to remove the interior pixels and retain the periphery of the Sun. For a pixel in Fig.  2c, if the luminance values of all its 4-connected neighbors were 1, its luminance value would be set to 0, otherwise, the pixel would be retained as a point on the periphery of the Sun, and its luminance value is 1. Fig. 2d illustrates the periphery of the Sun.
Extracting features from the Sun's binary image by the method of moments Fig. 2c is a binary image, where the Sun can be regarded as an ellipse. The feature parameters of the ellipse include the dimension a of the semi-major axis, the dimension b of the semi-minor axis, the center coordinate C (x c , y c ), and the angle γ between the major axis of the ellipse and the horizontal direction of the image, as shown in Fig. 3.
The binary image is expressed by an M × N matrix consisting of 0 and 1; 1 and 0 denote white and black, respectively. It is assumed that the matrix element is IM i, j , i = 1, 2, ⋯, M; j = 1, 2, ⋯, N, using the method of moments, the Sun's area A and the center coordinate C (x c , y c ) can be calculated by The angle γ between the major axis of the ellipse and the horizontal direction of the image can also be obtained by where, I xx , I xy and I yy are the secondary moments relative to C (x c , y c ), which can be computed by the following equations Applying Eqs. (4) and (5), the principal moments of inertia, I 1 and I 2 , are obtained If the Sun in the binary image is considered as an ellipse with two semi-axes a and b, its area A and the principal moments of inertia I 1 and I 2 can be expressed as Therefore, a and k can be solved by Applying the method of moments, the feature parameters of the Sun's binary image are finally extracted, including a, b, C (x c , y c ), γ and k.
Taking the Sun's image in Fig. 2c for example, the feature parameters are extracted. The result data are given in Table 1, where the length unit is pixel. Comparing the binary image of the Sun with the ellipse obtained by Eq. (7), the area of nonoverlapped domain is calculated, equal to 37,418; and its error relative to the total area is about 0.38%.
Fitting an ellipse to the Sun's periphery by the least-square method Fig. 2d illustrates the periphery of the Sun's visual shape. It is assumed that P i (x i , y i ), i = 1, 2, ⋯, q is a point on the periphery. In order to fit an ellipse to the periphery points, two coordinate systems are established, one is the global coordinate system O − XY, the other is the local coordinate system C − X e Y e , as shown in Fig. 3. In O − XY, X and Y are along the horizontal and vertical directions of the image, respectively. In C − X e Y e , the original point is set as the center of the ellipse; X e and Y e are along the semi-major axis and semi-minor axis of the ellipse, respectively.
In the local coordinate system C-X e Y e , P i (x i , y i ) on the Sun's periphery can be expressed in the form of a vector P ! i , satisfying The equation of fitting a reference ellipse can be written as where, r ! θ ð Þ is the polar radius of an edge point on the fitting reference ellipse; θ is the corresponding curve parameter; i ! e and j ! e stand for the unit vectors along X e and Y e ; (x e , y e ) is the coordinate of the point in C-X e Y e . Its unit tangent vector T ! θ ð Þ and unit normal vector N ! θ ð Þ can be expressed as Setting , T x , T y , N x and N y satisfy Δx c , Δy c , Δa, Δb] T , and the sign "T" means the transpose of a vector or matrix. The least-square method is applied to adjust the five feature parameters. The objective is to make ∑ q i¼1 ε 2 i minimum. Based on Eq. (18), the normal equation of the least-square method can be described as By Solving Eq. (19), the parameter increments Δx c , Δy c , Δγ, Δa and Δb can be figured out.
According to the principle of statistics, by applying Eq. (18), the variance SE 2 of the residuals of the fitting ellipse can be calculated Subsequently, the variances S 2 a and S 2 b are obtained where, S 2 a is the variance of a; S 2 b is the variance of b; V 4, 4 and V 5, 5 stand for the 4th and 5th diagonal elements of By employing the least-square method (LSM) to fit an ellipse to the periphery of the Sun's image, the optimal result parameters are obtained, including the semi-major axis a new , the semi-minor axis b new , the angle γ new and the center coordinate x new c ; y new Taking the Sun's periphery in Fig. 2d for example, the least-square method is used to fit an ellipse to the periphery. The initial values are the data listed in Table 1, which are the feature parameters calculated by the method of moments. The final results are shown in Fig. 4 and Table 2.
Seeing from Table 2, the variances S 2 a and S 2 b are very little. If the confidence level is 99.73%, for the semi-major axis a, its measurement accuracy is ±3S a = ± 0.1971, and the relative accuracy is 0.011%; for the semi-minor axis b, its measurement accuracy is ±3S b = ± 0.2016, and the relative accuracy is 0.012%. For the shape parameter k in Table 2, its relative measurement accuracy can be estimated as 0.023%, the sum of the relative accuracies of a and b. The standard deviation of the fitted ellipse to the periphery of the Sun is only about 4 pixels, while the average diameter of the Sun is about 1765 pixels. The results show that the Sun's image can be described as an ellipse.

Experimental results
The method of moments and the least-square method are combined to calculate the shape parameters k of the Sun's images. The images have been taken from sunrise to sunset in mid-December at Dalian. In the Sun's binary image, the Sun is considered as an ellipse. The method of moments is used to evaluate the feature parameters of the ellipse at first, and the feature parameters are set as the initial values of the iterative algorithm for fitting an ellipse to the periphery of the Sun; then, the least-square method or the iterative algorithm is performed to fit an ellipse to the periphery of the Sun and compute the feature parameters of the ellipse.
The resulting data are listed in Table 3, including the shape parameter k, the semi-major axis a, the semi-minor axis b, and the corresponding errors, where the maximum and minimum values of k, and their corresponding timings, as well as the maximum error of k are marked in bold. Obviously, the Sun's visual shape always appears to be an ellipse instead of a perfect circle. The Sun's image is the roundest at noon, while at 7:11 Beijing time, the Sun's image is a flat ellipse with the minimum shape parameter of 88.13%. It should be noted that the sunset photos were taken at Jinshitan Beach. Because the mountains obstruct some sunset view, the last time in Table 3 is 16:20 Beijing time.
To easily observe the variation of the shape parameter, a variation curve is plotted based on the data of Table 3, as shown in Fig. 5, where the horizontal axis is Beijing time, and the vertical axis is k. Seeing from Fig. 5, it is observed that k varies fast and abruptly around sunrise and sunset, whereas it is almost flat from 8:00 to 15:00 Beijing time.

Theoretical analyses and simulation
When the Sun is just rising in the morning, the Sun's image is flattened in the vertical direction, and its shape parameter is minimal; as time goes on, the Sun's visual shape becomes rounder and rounder, until noon, when the shape parameter reaches its maximum; however, the roundest image of the Sun is still an ellipse, but not a circle. In this Section, based on atmospheric refraction effects and the length contraction phenomenon in the Einstein's special relativity theory, the theoretical analysis and simulation are conducted to investigate the reasons for the variation of the Sun's visual shape.

Atmospheric refraction effects
By considering the atmospheric refraction effects only, the Sun's visual shape is affected by the refractive index of the atmosphere, the Sun's position in the sky and observation time; therefore, a model to estimate the refraction index at an altitude was established, and a set of formulas was derived to determine the zenith and azimuth angles of the Sun for different observation times and positions. Finally, an iterative algorithm was constructed to trace rays of light in the atmosphere.

Model of the atmospheric altitude and the refraction index
The Sun's diameter is 1.392 × 10 6 km, and the average Sun-Earth distance is about 1.496 × 10 8 km. It is assumed that the rays from the Sun travel in straight lines. When an observer on the Earth sees the Sun, his angle of view can be calculated by the Sun's diameter over the average Sun-Earth distance, equal to 32.004 ′ or 0.5334 o . However, the Earth is surrounded by a thick atmosphere, and the rays of light must pass through the atmosphere before reaching the Earth. By increasing the altitude, the air density decreases, resulting in a decrease in the refractive index of the atmosphere. As a result, the rays of light are bent in the atmosphere, which means the angle of view varies. The variation in the angle of view mainly occurs in the vertical direction, while in the horizontal direction, the angle of view almost remains unaffected. Therefore, the Sun in our eyes is like an ellipse, but not a circle. The shape of the ellipse varies with the observation time and position. The altitude influences the refractive index by the varying temperature and pressure. When the altitude is above 50 km, the refractive index approaches 1, so the paper mainly focuses on the propagation of light rays in the atmosphere below this altitude.
The relationship between air temperature T (in°K) and altitude H (in km) can be formulated as (Karin and Florian 2004;Jiang et al. 2013;Wang et al. 2018) where, T 0 is the air temperature at ground level.   (23), a temperature variation curve T(H) is given in Fig. 6, where the horizontal axis is the altitude H, the vertical axis is the temperature T.
According to hydrostatics, the air pressure P and the altitude H satisfy the condition where, ρ is the air density in g/m 3 ; g is the gravitational acceleration, equal to 9.8 m/s 2 . In the atmosphere, the air satisfies the equation of state where, R = 8.3144 J · Mol −1 · K −1 is the gas constant; μ = 28.966 g · Mol −1 is the molar mass of air. Substituting Eq.
(25) into Eq. (24), the relationship between the air pressure P and altitude H is obtained The relationship between air density ρ and altitude H can be expressed as Substituting Eq. (23) into Eqs. (26) and (27) where, P 0 and ρ 0 stand for the air pressure and air density at ground level. P 0 , ρ 0 and T 0 satisfy Eq.(25); P, T and ρ are functions of H. Setting T 0 = 269.15 o K, P 0 = 1.01325 × 10 5 Pa, ρ 0 can be calculated from Eq. (25), that is ρ 0 ¼ P 0 μ RT 0 ¼ 1:3115 kg=m 3 . Based on Eqs. (23), (28) and (29), the variations of atmospheric pressure and air density as a function of altitude are obtained, as shown in Fig. 7.
Neglecting influences of vapor in the air on the refractive index of the atmosphere, based on Eq. (25), the refractive index n can be calculated by (Karin and Florian 2004;Jiang et al. 2013;Xiang 2008) Where, λ (in μm) is the wavelength of light. Based on Eq. (30), the refractive index n(H) is plotted in Fig.  8, where the red line stands for the red light from the Sun at sunrise, at the wavelength of λ = 0.685 μm; the black line stands for the white light from the Sun at noon, at the wavelength of λ = 0.550 μm. Seeing Fig. 8, it can be found that with an increase in the altitude, the refractive index n decreases. When H > 30 km, n approaches 1, i.e. the refractive index in a vacuum. In order to investigate the influence of wavelength of light on the refractive index n, the difference of n between the red and white light is calculated; the results are shown in Fig. 9. Obviously, the influence of wavelength on n mainly occurs in the atmosphere below 30 km, but the influence is very little and can be neglected.

Model of the theoretical zenith angle of the sun, the observation time and position
The Sun's photos were taken in mid-December 2018 at Dalian. Besides the refractive index of the atmosphere, the observation time and position also affect the Sun's visual shape. Therefore, a model of the Sun's zenith angle, the observation time and position should be established. Here, neglecting atmospheric refraction effects, a set of formulas was derived to calculate the theoretical zenith angle of the Sun according to the observation position and time. Figure 10 is a sectional view of the Earth surrounded by a thick atmosphere. Dalian is located at 121.62 o E longitude and 38.92 o N latitude. It is assumed that a line connects Dalian (point A) and the Earth's center (point O), another one connects the Sun and the Earth's center (point O). The angle φ DL between the two lines is the theoretical zenith angle of the Sun.
The Earth rotates around its axis, and simultaneously moves around the Sun on the ecliptic plane. The Earth's axis is not perpendicular to the ecliptic plane since there is the declination angle between the equatorial plane and the ecliptic plane. Therefore, the theoretical zenith angle φ DL of the Sun varies with the Earth's rotation and its position on the ecliptic plane.
For convenience, it is assumed that the Sun moves around the Earth in an elliptical orbit, and the Earth is located at a focal point of the elliptical orbit. The parameters of the elliptical orbit are: the semi-major axis aŜ ¼ 1:4960 Â 10 8 km, the semi-minor axis bŜ ¼ 1:4958 Â 10 8 km, the focal length cŜ ¼ 25 Â 10 5 km, the eccentricity eŜ ¼ 0:0167. The Sun coordinate system O − X s Y s is shown in Fig. 11, where the original point is set as the focus where the Earth is located; X s points to the perihelion; Y s is perpendicular to X s . In terms of X s the elliptical orbit can be described as where, r(θ) is the polar radius; θ is the corresponding polar angle; PŜ ¼ bŜ À Á 2 =cŜ.
According to Kepler's second law, r(θ) satisfies where, C 0 is an undetermined constant; t stands for the time counting from the perihelion.
The area of the elliptical orbit is π Á aŜ Á bŜ, which is defined as Because the period of revolution of the Earth is T 0 = 365.2422 days, and T 0 = t(2π), substituting Eq. (34) into Eq. (33), C 0 can be obtained Substituting Eqs. (31) and (35) into Eq. (33), t(θ) can be rewritten as Letting tan ð φ2 Þ ¼ replacing the variable θ in Eq. (36) yields For a given time t 0 , in order to determine the Sun's polar angle θ, φ̂should be first solved by using Eq. (38); then θ is obtained by substituting φ̂into Eq. (37).
Eq. (38) is a transcendental equation. In order to solve for φ, an iterative algorithm is constructed where, φ̂n is the iteration variable in the n th iteration, with an initial iteration value of φ̂0 ¼ 2πÁt 0 T 0 . The error formula of the iterative algorithm can be expressed as where,φ * is the solution to Eq. (38); L is Lipschitz constant, the maximum value of the derivative ofφ n , satisfying From Eqs. (40) and (41), it is seen that the iterative algorithm converges very fast. During the iteration, when setting the accuracy as 10 −6 ( o ), the number of iterations is not more than 3, that is n ≤ 3.
In Fig. 11, the Earth coordinate system O − X E Y E Z E is also established, where the original point is the Earth's center, which coincides with the original point of the Sun coordinate system O − X s Y s ; Z E is the Earth's axis; Y E is along the line of the intersection of the equatorial plane and the ecliptic plane,  and coincides with Y s ; X E is in the equatorial plane, the angle δ is the declination angle equal to 23.433 o .
According to the Gregorian calendar, when the Sun is located at an intersection of the equatorial plane and the ecliptic plane, or when the Sun is on the Y E , the time in 2018 is 00:15:24 on March 21 (Beijing time); in this paper, this time is taken as reference time. Beijing time is defined as the time at 120 o E longitude. At 00:15:24 on March 21, ξ 0 BJ is used to denote the angle between the projection of the 120 o E longitude on the plane O − X E Y E and the X E axis. At the reference time, the Sun's polar angle θ = 90 o , and there are 11.7433 h from 00:15:24 to 12:00:00 Beijing time; then, ξ 0 BJ can be calculated by where, ω E ¼ 15:0411 o =hour, is the average angular velocity of the Earth around its axis; ω S ¼ 0:0411 o =hour, is the average angular velocity of the Earth around the Sun. Dalian is located at 121.62 o E longitude, with a longitude difference between Dalian and Beijing of 1.62 o . ξ DL is used to denote the angle between the X E axis and the projection of the 121.62 o E longitude on the O − X E Y E , as shown in Fig. 11. At the reference time, using ξ 0 DL to describe ξ DL , and ξ 0 DL can be calculated by Taking 00:15:24 on March 21, 2018 (Beijing time) as reference time, after Δt hours, ξ DL can be obtained In the Earth coordinate system O − X E Y E Z E , η DL is used to stand for the angle between the radius direction from the original point (point O) to Dalian (point A) and the coordinate plane O − X E Y E , as shown in Fig. 11. Because Dalian is located at 38.92 o N latitude, η DL = 38.92 o . OA ! is the unit vector along the radius direction OA, and can be expressed by where, i E ! , j E ! and k E ! are the unit vectors along X E , Y E and Z E , respectively.  37) and (38), setting θ = 90 o , the moving time of the Sun from the perihelion to the reference time can be obtained, that is Δt 0 = 89.367789 days. In the Sun coordinate system O − X s Y s , applying the iterative algorithm of Eq. (39), setting t o = Δt 0 + Δt, the polar angle θ of the Sun after moving Δt hours can be obtained from Eqs. (37) and (38). At the moment t o , it is assumed that OS ! is the unit vector along the radius vector r(θ) of the Sun. In the Earth coordinate sys- The angle between OS ! and OA ! is the theoretical zenith angle φ DL of the Sun relative to Dalian, as shown in Fig. 10. Based on Eqs. (45) and (46), φ DL can be figured out by Projecting OS ! on the equatorial plane O − X E Y E , ξ Sun is used to denote the angle between the projection vector of OS ! and the X E , which describes the direction of the Sun, satisfying 12:00:00 Beijing time is not the moment when the Sun shines over Dalian. The time difference Δt can be calculated by In Fig. 11 In the Dalian coordinate system A − X DL Y DL Z DL , the components of OS ! can be written as Thus, relative to Dalian, the Sun's azimuth angle ψ DL can be figured out by Defining clockwise rotation as positive, ψ DL describes the angle between the Sun's direction OS ! and due north. It is assumed that atmospheric refraction effects are neglected, the observation position is Dalian, and the observation time is Beijing time one day in mid-December 2018. In Fig. 12, the variations of the theoretical zenith angle φ DL and the azimuth angle ψ DL are given, which are calculated from Eqs. (47) and (53).
Besides, the time when the Sun shines over Dalian is calculated by Eq. (49). The simulation results are given in Fig. 13, where the horizontal axis is the date in the whole year of 2018, and the vertical axis is the time when the Sun shines over Dalian. From Fig. 13, it is seen that the time ranges from 11:36 to 12:08.

Iterative algorithm for tracing rays of light in the atmosphere
Because the refractive index of the atmosphere varies with altitude, the rays of light are bent in the atmosphere surrounding the Earth. In the inhomogeneous medium, the propagation of light satisfies (Born and Wolf 2007) where, R ! s ð Þ stands for the trace of light; s is the length of arc; n is the refractive index of the atmosphere; ∇ ! n is the gradient of the refractive index.
According to Eq. (30), the refractive index of the atmosphere can be expressed as n(R), where R is the distance of a point in the atmosphere to the Earth's center. Because the contour surface of n(R) is spherical, its gradient direction is parallel to the radial direction of the Earth, that is R ! Â ∇ ! n ¼ 0. Thus, Eq. (54) can be simplified as   to the plane. Therefore, Eq. (57) proves that the trace of light in the atmosphere is a plane curve, and the plane is determined by the Earth's center, the incidence point and incidence direction of light, as shown in Fig. 14. Usingθ to denote the angle between the propagation direction of light and the radial direction of the Earth, taking the modulus of the two sides of Eq.
In order to trace the rays of light in the atmosphere, the atmosphere below the altitude of 50 km is equally divided into 500 layers. It is assumed that the refractive index of each layer is constant. For convenience, the subscript"i" is used to express the ith layer of the atmosphere, i = 1, 2⋯500. R i stands for the distance of the ith layer of the atmosphere from the Earth's center; n i is the refractive index of the ith layer; θ i is the angle between the radial direction of the Earth and the propagation direction of light in the ith layer, as shown in Fig. 14. Eq. (58) is rewritten as According to Eq. (59), an iterative algorithm for tracing the rays of light can be constructed. In Fig. 14, a coordinate system O − XY was established, where the original point was the Earth's center; the horizontal axis X pointed to the Sun; the Y axis was perpendicular to the X axis in the propagation plane of the rays. Based on the reversibility of the optical path, it is assumed that the rays of light travel from the Earth to the Sun.
can be traced and calculated by where, β i + 1 is the incidence angle of light from the ith to the (i + 1)th layer. Eq. (60) is a one-step explicit scheme with firstorder accuracy, and can be used to simulate the propagation path of light in the atmosphere. During the simulation, the observation position coordinate P ! 0 x 0 ; y 0 ð Þand the angleθ 0  Fig. 15, Dalian is set as the origin, the horizontal and vertical directions are the same as the X and Y axes in Fig. 10, and the time is Beijing time. During the simulation τ Sun = 0, which means when the light reaches the atmosphere at the altitude of 50 km, its propagation direction is parallel to the X axis.
Defining δφ ¼ e φ DL −φ DL , φ DL is the theoretical zenith angle of the Sun, and e φ DL is the observed zenith angle of the Sun. In Fig. 16, the difference between φ DL and e φ DL is given, where the horizontal axis is Beijing time of a day, the vertical axis is δφ. From Fig. 16, it could be found that δφ is bigger at sunrise and sunset than at noon, a result caused by atmospheric refraction.
In the Dalian coordinate system A − X DL Y DL Z DL , the direction of the Sun can be described by the observed zenith angle e φ DL and the azimuth angle ψ DL Based on Eq. (61), the direction of the Sun in A − X DL Y DL Z DL can be calculated. Fig. 17 is the projection of the Sun's path on the Y DL Z DL plane at X DL = − 1km.

The length contraction effects in the Einstein's special relativity theory
Here, atmospheric refraction effects are neglected, and the length contraction phenomenon in the Einstein's special relativity theory is considered to investigate its effect on the Sun's visual shape. The length contraction phenomenon in the Sun's image can be described by the contraction ratio and the Because the Sun is far away from the Earth, and the Earth rotates around its axis and simultaneously around the Sun, the relative velocity between the Earth and the Sun is very big. According to the length contraction phenomenon in the Einstein's theory of special relativity, when an observer on the Earth sees the Sun, the size of the Sun's visual shape in its moving direction will become shorter, while the size perpendicular to the moving direction remains unchanged. Assuming that the Sun is a sphere, the Sun in the eyes of the observer will be an ellipsoid, and its ratio of the minor to the major axis can be expressed as where, c = 3 × 10 8 m/s, the velocity of light in vacuum; v is the relative velocity between the Sun and the observer. Taking into account that the average angular velocity of the Earth around its axis is ω E ¼ 15:0411 o =hour, the average angular velocity of the Earth around the Sun is ω S ¼ 0:0411 o =hour, the average distance from the Sun to the Earth is 1.496 × 10 8 km, the relative velocity between the Sun and the observer can approximate to v≈1:496 Â 10 8 ω E −ω S ≈1:088 Â 10 7 m=s ð63Þ Substituting Eq. (63) into Eq. (62) yields k SRT is the shape parameter of the Sun's visual shape, describing the effect of the Einstein's special relativity theory on the visual shape of the Sun. Eq. (64) proves that even if neglecting atmospheric refraction effects, the Sun's shape in our eyes is an ellipse but not a perfect circle, and the ratio of the length contraction or the ratio of the dimensional deviation between the minor and major axes to the major axis is about 0.07%. However, seeing from Fig. 17, the Sun's direction relative to the observer varies with time; as a result, the direction of the major axis in the ellipse image of the Sun also varies.
Applying Eqs. (31), (32) and (35), the transient angular velocity ω S (θ) of the Sun around the Earth can be expressed as The Sun rotates around the Earth in the elliptical plane, as shown in Fig. 11; its angular velocity vector ω ! S θ ð Þ can be written as The Earth rotates around its axis k E ! ; its angular velocity vector ω ! E can be written as Meanwhile, applying Eqs. (31) and (46), the relative radial where, V ! r θ ð Þ is caused by the variable distance between the Sun and the Earth.
Based on Eqs. (66), (67) and (68), the transient relative velocity V ! SE θ ð Þ between the Sun and the Earth can be calculated by Substituting Eq. (69) into Eq. (62) yields Eq. (70) describes the way that the Einstein's special relativity theory affects the visual shape of the Sun when the relative position varies between the Sun and the Earth.
When a high-resolution digital camera is used to take photos of the Sun, the camera remains horizontal and points to the Sun, which means the Sun's image is taken in the di-

rection OS
! . In the Dalian coordinate system A − X DL Y DL Z DL , the horizontal direction of the camera can be described by the Sun's azimuth angle ψ DL , the unit vectors AN ! DL and AW ! DL along two coordinate axes, that is where, H ! camera is the horizontal direction of the camera or image. The vertical direction V ! camera of the camera or image can be determined by In the Sun's image, the direction of the semi-minor axis caused by the Einstein's special relativity theory is along the projection of the Sun's velocity V ! SE θ ð Þ relative to the Earth on the image. Defining ψ bH as the angle between the horizontal direction of the camera and the semi-minor axis of the Sun' image, it describes the direction of the length contraction or semi-minor axis of the Sun's image, and can be calculated by Upon neglecting atmospheric refraction effects or only considering the Einstein's special relativity theory, and regarding the Sun's visual shape as an ellipse, its shape parameter k SRT and the angle ψ bH can be evaluated by Eqs. (70) and (73), respectively. In Fig. 18, the variation of the shape parameter k SRT is given, where the horizontal axis is the date in the whole year of 2018, and the vertical axis is k SRT . During the simulation, the observation time is set to 12:00:00 Beijing time every day, and the observation position is Dalian. Figure 18 shows that the Sun's visual shape is the roundest at the winter solstice, and the flattest at the spring and autumnal equinoxes.
In Fig. 19, the variation of the angle ψ bH is given, where the horizontal axis is the Beijing time of a day in mid-December 2018, the vertical axis is the angle between the semi-minor axis of the Sun's ellipse image and the horizontal direction of the camera, which describes the direction of the length contraction in the Sun's image. Figure 19 shows that the angle ranges from 47 o to −47 o , and the extremes occur at sunrise and sunset. The two directions of the length contraction occurring at sunrise and sunset basically satisfy the mirror-image relation.

Numerical examples of the atmospheric refraction effects
Eq. (57) shows that the light in the atmosphere propagates in a plane, and the plane is determined by the Earth's center, the incidence point and incidence direction of light, as shown in Fig. 10 or Fig. 14. In the following numerical examples, the propagation path is simulated and analyzed by applying Eqs. (47) and (60).
When an observer on the Earth sees the Sun, the rays from the Sun pass through the thick atmosphere, and then focus on the observer's eyes. The Sun's visual shape is determined by the rays from the upper and lower edges of the Sun; therefore, during the analysis, only the rays from the edges of the Sun are simulated and analyzed. Meanwhile, due to the reversibility of optical path, it is assumed that the rays travel from Dalian to the Sun. In addition, when the altitude in the atmosphere is above 50 km, its refractive index approaches 1; therefore, the atmosphere below the altitude of 50 km is only considered. If the Sun's rays propagate in straight lines, the observer on the Earth has an angle of view of 0.5334 o . As a result, when simulating the ray of the upper edge of the Sun, the angle constraint τ Sun is set to +0.2667 o , and when simulating the ray of the lower edge of the Sun, τ Sun is set to −0.2667 o . This means that when the rays reach the atmosphere at the altitude of 50 km, the angles between the rays and the x axis are restricted to be ±0.2667 o .
Dalian is taken as the observation position. The coordinates o f D a l i a n a r e P ! 0 OAcos φ DL ð Þ; OAsin φ DL ð Þ ½ , OA = 6371 km, and the observation times are in the morning, at noon and in the afternoon on one day in mid-December 2018. In order to determine the initial angleθ 0 between the ray from Dalian and the radial direction of the Earth, the method of bisection is performed by setting the initial range  Fig. 20, Dalian is set as the origin, the horizontal and vertical directions are the same as the X and Y axes in Fig. 10. From the figure, it can be found that the horizon is basically tangent to the ray from the lower edge of the Sun in the morning. Figure 21 illustrates bending deformations of the rays in the atmosphere. In Fig. 21, the horizontal axis is the horizontal distance between the observer and the Sun, and the vertical axis is the angle between the propagation direction of the ray and the line connecting the Sun's and the Earth's centers. Compared with the morning, the propagation distance of the rays from the edge of the Sun becomes short in the 50-kmthick atmosphere. Meanwhile, the angles between the rays from the Sun and the Earth's radial direction passing through Dalian (point A) also become small. As a result, bending deformations of the rays in the atmosphere are weakened, which can easily be observed in Fig. 23. When the rays reach Dalian at noon, an observer's angle of view becomes 0.5327 o , bigger than 0.4457 o (in the morning) and a little less than 0.5334 o (the angular diameter of the Sun). According to Eq. (74), the shape parameter k Ref = 0.9986, which means that the Sun's image is rounder than that in the morning. In order to further investigate the effects of atmospheric refraction on the Sun's visual shape in the whole year, a variation of the shape parameter is given in Fig. 26, where the horizontal axis is the date in the whole year of 2018, and the vertical axis is the shape parameter k Ref , which is calculated by Eq. (74). During the simulation, the observation time is set as 12:00:00 Beijing time every day, the observation position is Dalian. Fig. 26 shows that the Sun's visual shape is the roundest at the summer solstice, and the flattest at the winter solstice.

Comparison and verification
When researching the effect of the Einstein's theory of special relativity on the Sun's visual shape, it is proved that the contraction ratio occurring in the Sun's image varies slowly and little with time. The shape parameter k SRT ranged from 0.99934 to 0.99948 in the whole year of 2018, as shown in Fig. 18. However, the contraction direction or the semi-minor axis direction in the Sun's ellipse image varies rapidly and greatly with time. The angle ψ bH between the horizontal direction of the image and the semi-minor axis of the Sun's ellipse image ranged from 47 o to −47 o on one day of mid-December 2018, as shown in Fig. 19.
When researching the effect of atmospheric refraction on the Sun's visual shape, it was shown that the shape parameter k Ref of the Sun's ellipse image changes greatly. k Ref ranged from 0.8357 to 0.9986 on one day of mid-December 2018, as shown in Figs. 21, 23 and 25. However, the contraction direction of the Sun's ellipse image caused by the atmospheric refraction remains unchanged. Because Eq. (57) proves that the propagation path of light in the atmosphere is a planar curve, the semi-minor axis of the Sun's ellipse image is always perpendicular to the horizontal direction of the image.
Assuming that the Sun is a circle, by applying the elementary transformation, it can be shown that after the Einstein's theory of special relativity and atmospheric refraction effects are applied to the circle, the circle will become an ellipse, and the shape parameter k of the ellipse can be calculated by Based on Eq. (70), the shape parameter k SRT can be figured out, which describes the effect of the Einstein's theory of special relativity on the Sun's visual shape; based on Eq. (74), the shape parameter k Ref can be evaluated, which describes the effects of atmospheric refraction on the Sun's visual shape; combining Eqs. (70), (73), (74) and (75), the shape parameter k can be calculated, which describes the effects of atmospheric refraction and the Einstein's theory of special relativity on the Sun's visual shape.
In order to show the reasons for the variation of the Sun's visual shape, the theoretical simulation results and the experimental results were compared, as shown in Figs. 27 and 28. The horizontal axis is Beijing time from sunrise to sunset, the vertical axis is the shape parameter, and the pink circles are the experimental results in Table 3. In Fig. 27, the blue line stands for the simulation results by only considering atmospheric refraction effects. In Fig. 28, the blue line stands for the simulation results by considering the combination of atmospheric refraction effects and the Einstein's theory of special relativity.
The relative error is defined as the difference between simulations and experimental results over the experimental results. Fig. 29 shows the relative errors of the simulated results, where the blue line is the errors calculated when only considering atmospheric refraction; the green dotted line is the errors calculated when considering the combination of atmospheric refraction and the Einstein's theory of special relativity. Fig. 29 shows that the Einstein's special relativity theory affects the Sun's visual shape much more at noon than at sunrise and sunset. For example, at 11:38 Beijing time, the shape parameter extracted from the Sun's photo is k Meas = 0.9995; the shape parameter caused by atmospheric refraction is k Ref = 0.9987 with a relative error of 0.8 × 10 −3 ; the shape 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:  parameter calculated by considering the combination of atmospheric refraction and the Einstein's theory of special relativity is k = 0.9993 with a relative error of 0.2 × 10 −3 . Although the relative error of k Ref is about four times than that of k, the error is still very small and can be neglected. Figs. 27, 28 and 29 illustrate that the simulated results are coincident with the experimental ones in Table 3.
The comparison study showed that the theoretical analyses in the paper are corrected; the reason for the variation of the Sun's visual shape is mainly due to atmospheric refraction effects, while the length contraction effect of the Einstein's special relativity theory also contributes a little except at sunrise and sunset. Therefore, the Einstein's special relativity theory can explain why the Sun's visual shape always appears to be an ellipse instead of a perfect circle.
In order to further analyze the contributions of atmospheric refraction and the Einstein's special relativity to the variation of the Sun's visual shape, three variations of the shape parameter were simulated in Fig. 30, where the horizontal axis is the date in the whole year of 2018, and the vertical axis is the shape parameter. The blue line stands for k Ref caused by atmospheric refraction effects, similar to the curve in Fig. 26; the pink dotted line stands for k SRT caused by the length contraction effect in the Einstein's special relativity theory, which is the same with the curve in Fig. 18; the black line is the shape parameter calculated when considering the combination of atmospheric refraction and the Einstein's special relativity theory. During the simulation, the observation time was set at 12:00:00 Beijing time every day, and the observation position was Dalian. Fig. 30 shows that in the whole year, the Sun's visual shape at Dalian is the roundest at the end of February, and is the flattest at the winter solstice. The conclusion differs from that drawn from only considering atmospheric refraction or the Einstein's special relativity, because both of them contribute to the variation of the Sun's visual shape.

Conclusions
From sunrise to sunset, the Sun's visual shape changes continuously from a flatter ellipse to almost a circle and then back to a flatter ellipse. In the paper, the experimental measurements and theoretical analyses were performed to investigate the reasons for the variation of the Sun's visual shape. Some meaningful conclusions were drawn: The method of image processing, the method of moments and the least-square method were combined to perform experimental measurements and calculations. The statistical error analyses showed that the relative measurements accuracy was about 0.023%, and the standard deviation of the fitting curve of ellipse was only 4 pixels. The experimental results showed that the Sun's visual shape can be approximated by an ellipse accurately. 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 -3 -2 -1 0 1 2 3 x 10 -3

Beijing Time errors of k between simulation results and measurement data
Atmospheric Refraction Effects Atmospheric Refraction Effects +Special Relativity Theory  Atmospheric refraction effects make the Sun's visual shape become an ellipse, and have a great influence on its shape during a day. Because the refraction index of the atmosphere was expressed as a function of altitude and wavelength of light, the trajectory of sunlight is a planar curve, when considering the atmospheric refraction effects only, the Sun's visual shape contracts in the zenith direction, resulting in an elliptic Sun. In the Sun's photo, the contraction direction remains unchanged, but the contraction ratio varies rapidly and greatly with time, especially at sunrise and sunset.
The length contraction effects in the Einstein's special relativity theory also transform the Sun's visual shape to an ellipse, and contribute a little to its shape variation on a day, especially at noon. Due to the relative movement of the Sun to the observer standing on the Earth, the Sun's visual shape contracts in a direction, resulting in an elliptic Sun. However, unlike atmospheric refraction effects, the contraction ratio varies slowly and little with time, and it is about 0.06%, while its contraction direction varies rapidly and greatly with time.
Comprehensively considering effects of atmospheric refraction and the Einstein's special relativity theory on the Sun's visual shape, although their contraction directions are different, it is shown that the Sun's visual shape is still an ellipse.
On a day, the Sun's visual shape is the roundest at noon and the flattest at sunrise and sunset; in the whole year, the Sun's visual shape for an observer at Dalian is the roundest at the end of February, and the flattest at the winter solstice. For observers at different geographic positions, if they see the Sun at the same moment, the Sun looks different. In theory, an observer at the equator could see the roundest Sun.
Comparing the theoretical simulations and experimental measurements, the relative errors were less than 0.3%. The results, thus, verified the theoretical analyses in the paper, including a set of formulas describing the relationship between the zenith angle of the Sun, the observation times and positions, the model of altitude and refraction index, the iterative algorithm for tracing rays of light in the atmosphere, and the calculation of feature parameters of the Sun's visual shape.