Modelling full-colour images of Earth: simulation of radiation brightness field of Earth’s atmosphere and underlying surface

ABSTRACT This article is dedicated to the problem of realistic colour rendering of space object images using the tools of computer graphics. In the form of a short essay, the authors describe the essence, sources and functionality of modern graphics applications. Particular attention is paid to the application of modern graphics in space science. The specific purpose of the study is the use of computer graphics in the field of remote sensing of the Earth’s surface. This article describes a method for synthesizing images to develop realistic 3D models of colour Earth images in the visible spectral range, observed from geostationary orbits. The method is based on the improved model of atmospheric radiation for arbitrary sighting conditions in an inhomogeneous spherical atmosphere. Physical models of horizontally inhomogeneous distributions of atmospheric density, temperature and albedo of the Earth were improved. All calculations were performed in accordance with the model of molecular scattering of radiation in a spherical atmosphere, taking into account sunlight forward-scattering and reflection from the planet’s surface. This allows us to obtain images of the Earth in its various phases, observed from arbitrary heights. The obtained theoretical colour images of the Earth were compared with black and white images from modern geostationary satellites.


Introduction
Methods of remote sensing of the Earth's surface, the use of which began in the 1960s with the launch of artificial space satellites, led to a real technical revolution in many disciplines of Earth sciences (Rodionov et al. 2020;Koteleva and Frenkel 2021). Their use made it possible to assess the planetary distribution of different natural objects and processes, as well as anthropogenic impact on the environment, in order to solve various fundamental and applied problems of natural sciences. The first and most widely used method of remote sensing is taking photographs of the Earth's surface and the cloud cover in the visible spectrum. At the same time, the beginning of the space era was marked by rapid development of computer graphics, stimulated by successful transition from mechanical calculators to electronic computers. A brief historical overview of the development of computer graphics in terms of Earth imaging is presented in Appendix A.

Studies of natural phenomena and computer graphics
Image synthesis for realistic 3D model elaboration is one of the most popular areas of relevant research studies (Zotti 2007;Chepyzhova, Pravdina, and Lepikhina 2019). In this way, we obtained colourful images of natural objects, such as mountains, trees, sea, clouds and the Earth images ( Figure 2) as well (Preetham, Shirley, and Smits 1999). The images are widely used in movies or commercial TV shows like CG Earth Images library (Shutterstock 2020). However, here they mostly strove to get realistic images, not physically relevant ones.
Image synthesis is used to simulate the structure views of building and urban area (Kaneda et al. 1991;HosseiniHaghighi et al. 2020) scenery under various atmospheric conditions, fog views in various sunlight and colours of sunset and dawn (Klassen 1987). Synthesis is also applicable for the calculation models comparison intended to predict daily variations of sky lighting for the energy saving tasks (HosseiniHaghighi et al. 2020;Littlefair 1994). These papers show the image synthesis considering atmospheric influence, but singlescattering approximation only.
Meanwhile, in our opinion, for perfect Earth image modelling, resembling the observations from meteorological satellites, and for their association with weather models, in this case, they can be used to study climate changes, as well as for space flight simulation, it is necessary to have physically based images that require multiple scattering. The subject was studied and described in few papers (Nishita et al. 1996;Tomoyuki et al. 1993;Nishita and Dobashi 2001). The authors of these articles showed detailed images of the landscape and the sea, but the scope of data on the sunlight diffused in the atmosphere and re-reflected by the Earth's surface (diffuse component) was not considered. Moreover, the non-homogeneous nature of temperature distribution and the Earth's albedo, affecting the geographical distribution of high-altitudinal atmospheric density changes, were not considered as well.
It should be noted that the synthesis of various natural phenomena and scene images is carried out mainly by digital artists; therefore, the atmospheric models' reality is rather neglected, with the consequential impact on the research results.
At the same time, the experts on atmospheric optics provide their results in graphical or tabular form, using numerous atmospheric and meteorological parameters that are poorly controlled in experiments. This is why it is difficult and sometimes even impossible to compare these data with a large scope of data obtained from modern space experiments, presented as images.

The role of optical models in the visualization of natural phenomena
Further development of technologies for photo and video recording made it possible to replace black-and-white representations of captured objects with colour ones, photo and video films -with carriers of digital information, as well as significantly increase the resolving power of camera lenses. Colour images of the underlying surface acquired special importance in this process. They made it possible to assess distribution areas of surface water masses and zones of biological productivity in the World Ocean and to determine various details of terrestrial landscapes in a very wide colour spectrum. Now the largest image area of the underlying Earth's surface can be obtained from satellites in geostationary orbits.
Various computer programs help to increase the processing speed of colour digital information, but their algorithms are mostly based on principles (parallel computing and delays in individual tasks) that do not provide in the visible spectrum an objective picture of realistic colour rendering for natural objects и явлений. This is largely due to the shift in microprocessors from single-core to multi-core chips. This trend, most evident in CPU families, aims to run parallel workloads. This maximizes overall throughput although it may be necessary to sacrifice the sequential performance of a single task and the impact of the light transmission functions of Earth's atmosphere by photo and video cameras. This is largely due to the influence of light transmission of photo and video cameras through the Earth's atmosphere. These functions are characterized by strong fluctuations in the three-dimensional plane of the Earth's surface, as well as in the time range of their seasonal and synoptic-scale variations.
The MODTRAN (atmospheric transmittance/radiation with moderate resolution) software in different versions is currently especially popular. The first version of this program was published in 1987 (Berk, Bernstein, and Robertson MODTRAN 1987). Related subroutines were integrated into LOWTRAN 6. The spectral resolution of this new option exceeds 5 cm −1 (FWHM). The focus is on FASCOD (Fast Atmospheric Signature Code) calculations with the HITRAN (High-Resolution Transmission) molecular spectroscopy spectral database (Rothman et al. 2013). It presents a collection of spectroscopic parameters used to model and analyse the transmission and emission of light in the gas atmospheres of the planet.
The first LOWTRAN programs, which preceded the MODTRAN series and appeared simultaneously with the personal computers in the early seventies, were semi-empirical. A LOWTRAN 3 (see reference 3 in Selby and A 1975) computer program is described for calculating the transmittance of the atmosphere in the spectral region from 0.25 to 28.5 micrometres at a spectral resolution of 20 cm −1 . The program provides a choice of six atmospheric models covering seasonal and latitudinal variations from the sea level to 100 km, two haze models, and accounts for molecular absorption, molecular scattering and aerosol extinction. Refraction and earth curvature effects are also included. This program provides some modifications to the molecular absorption and aerosol extinction data provided in an earlier LOWTRAN 2 report. In addition, input modifications have been made, making the LOWTRAN 3 program considerably more flexible in terms of the input of meteorological data.
LOWTRAN 7 is a low-resolution spreading model (Kneizys, Eric Shettle, and Chetwynd 1988) and simultaneously a computer code for predicting atmospheric transmittance and background radiation from 0 to 50,000 cm −1 at 20 cm −1 resolution. The code is based on the LOWTRAN 6 model (Kneizys et al. 1983). Multiple scattered radiation has been added to the model, as well as new molecular band model parameters and new ozone and molecular oxygen absorption parameters for the UV radiation. Other modifications include a winddependent desert model and new models of cirrus clouds, clouds and rain.
The SENTRAN computer code, originally developed by researchers at the University of Pennsylvania, is worth noting. It allows users to quickly estimate transmittance and emission from LOWTRAN in response to perturbations of input atmospheric conditions. The code provides a convenient and efficient interface for creating multiple LOWTRAN input decks as well as a graphical representation of the output data in 2D and 3D formats. The main purpose of this paper is to update the SENTRAN source code for use with LOWTRAN7 (Longtin et al. 1991) as well as MODTRAN (Longtin and Hummel 1993).
The MODTRAN computer code series is now used worldwide by researchers and scientists in government agencies, commercial organizations and educational institutions to predict and analyse optical measurements through the atmosphere. The code is integrated into many operational and research sensor and data processing systems. They are commonly referred to as atmospheric correction in remotely sensed multispectral and hyperspectral imagery (MSI and HSI). The most advanced version is the MODTRAN6 code (Berk et al. 2014). This code provides a modern modular object-oriented software architecture, including an application programming interface, extended physics functions, a line-by-line algorithm, an additional set of physics tools and new documentation.
Unlike the other codes, this code contains a large number of input parameters. Therefore, taking into account interrelated physical processes, complex approximating operations are applied, which makes the use complicated and inconsistent in terms of selfconsistency of the results obtained. Moreover, involving a large number of arguments to a function describing any process can be a fitting.
Furthermore, since MODTRAN is the successor of LOWTRAN, it has also inherited a number of 'rudiments' and 'atavisms'. For example, it has six layer models of the Earth's atmosphere, the program code in FORTRAN, the problem of the reality of sphericity calculations and other less significant.
Additionally, despite the extensive apparatus of available local atmospheric correction capabilities, MODTRAN is unable to solve global problems.
This article examines the influence of the Earth's atmosphere on the possibility of determining the colour rendering of images of its outgoing radiation received from artificial Earth satellites (Khokhlov et al. 2005) and groundbased astronomical observatories (Khokhlov 2001a), in the process of observing space objects of the Solar system.
At present, experimental and theoretical material on optical and other parameters of planetary atmospheres allows us to set the task of creating a self-consistent model that will reflect the global distribution of atmospheric radiation. In this model, the atmosphere is regarded as a single entity, all the components of which are interconnected. The model should have a small number of adjustable input data that adequately describe the observation conditions, geophysical situation and the state of the atmosphere. Therefore, we have created a model of Rayleigh scattering in a locally isothermal spherical atmosphere with an inhomogeneous latitudinal distribution of temperature and albedo of the planet's surface.
Using a previously created physical model of Rayleigh scattering in the atmospheres of planets with inhomogeneous distribution of temperature and surface albedo (Khokhlov 2001a), we developed an algorithm for realistic calculation of images of the Earth observed from space, from an arbitrary point of observation. This model has a tightly limited set of input data, namely, spectral range, Julian date, geographic coordinates and altitude, which determine the position of the observer, as well as angular range covered by the observer of an object (or a part of it).
The following sections describe the algorithm of the proposed model for constructing realistic images of the Earth in detail.

Materials, methods, data structures and schema algorithms
Currently, a great amount of satellite data on the landscape, vegetation and clouds is provided as colour or black and white images (Strizhenok, Ivanov, and Suprun 2019;Kovyazin et al. 2020a;Kazantsev, Boikov, and Valkov 2020;Musta, Zhurov, and Belyaev 2019;Ignatiev and Kessel 2016;Kopylova and Starikov 2018;Kovyazin et al. 2020b). For its analysis, we need to know about the Earth's atmosphere impact and lighting used by satellite equipment to get the image of various natural objects. These data can be obtained if we compare the experimental data in the form of colour images with the results of model calculations for similar images of various landscapes, vegetation, clouds, snow and ice cover (Kozhaev et al. 2017). This paragraph describes the way of the Earth outgoing radiation modelling results transformation into a colour image. The colours show the outgoing radiation under different illuminations of the inhomogeneous Earth surface by direct and reflected solar radiation.

Realistic model of Rayleigh scattering in the Earth's atmosphere
As mentioned above, the purpose of this article is to visualize the calculation results of Earth's radiation brightness using a realistic model of a molecular, spherical, 3D inhomogeneous, terrestrial atmosphere with an inhomogeneous (integral) distribution of albedo and surface temperature. Such a realistic model was published in 2001 (Khokhlov 2001a) by one of the authors of this article. However, having considered that a mere reference to that study would not suffice, we included an overall description of its basic concepts and key solutions in the form of a block structure, presented in Figure 1.
Since the models of albedo and temperature distributions in (Khokhlov 2001a) were simplified, we found it useful to present them in an updated form (see Sections 2.2 and 2.4.2).
As for the description of the block structure diagram, it should be noted that all the parameters indicated in the upper block are mandatory and the absence of any  parameter deprives us of the possibility to compare visualization results for theoretical and experimental images.
The block 'Input parameters and data arrays' contains a set of 7 parameters and 2 arrays. The parameters include Julian date J d , geographic latitude φ, longitude λ and altitude h of the observer, zenith angle z and azimuth A of the observation direction, as well as spectral range ΔΛ. The full Julian date consists of the day d, month m, year y and time UT.
The set of parameters φ, λ, h, z and A characterize the position of the observation point hodograph relative to the Earth's surface. The arrays of surface temperature T (φ,λ) and albedo distribution Al(φ,λ) are calculated separately using dependent parameters and stored in binary files. The temperature and albedo simulation results are the input parameters we use to calculate the Earth's outgoing radiation as a planet.
In the block 'Astro', the Julian date J d is used to determine the position of the Earth in the ecliptic orbit (r e , ν e , φ e ); angular geocentric coordinates θ ax and ψ ax are determined using ecliptic pairs ν e , φ e and ν ax , φ ax for the Earth's rotation axis. All coordinates are in the standard spherical representation. For example, r e is the distance to the Sun, the first angular coordinate ν e is the colatitude and the second one φ e is the longitude. The subscript refers to the object, e.g. 'e' stands for Earth and 'ax' -for Earth's rotation axis. The ampersand (&) at the beginning indicates resulting parameters. Parameters without an ampersand are the initial known values that shape the result.
The block 'Transformation and modification of the coordinate systems' is designed to determine relationships between several systems, used for calculating components of scattered radiation in the Earth's atmosphere. The coordinate systems include the ecliptic one with the Sun at its centre, geocentric ones (geographic, horizontal and the resulting main system), as well as two others related to the modelling of optical processes in the Earth's atmosphere. Ecliptic, geographic and horizontal systems, due to the differences in their angular characteristics, are transformed into a single main coordinate system. In its own turn, the latter is modified into two convenient systems for calculating optical phenomena, including illuminated regions separated by a terminator.
The block 'Rayleigh scattering model' contains calculations of spectral brightness for the outgoing radiation of the atmosphere and the underlying Earth's surface (surface and cloud regions).
The tasks solved in this block include • Estimation of single-scattering brightness of direct solar radiation (direct component B 1 R ); • Estimation of scattering brightness of solar radiation, reflected from the underlying surface and attenuated by passing through the atmosphere (diffuse component B 2 Rs ); • Allowance for the contribution of attenuated radiation of the underlying surface, illuminated by the Sun (B E ).
While solving these problems, great attention is paid to the complex issues of atmospheric transmittance and phase separation of illuminated and shaded regions of the Earth.

Calculation of time-space distribution of temperature over the ocean and land surface
A relevant experimental and theoretical stuff on the planetary atmospheres optical parameters enables us to set a challenge and create a self-consistent model of atmospheric radiation global distribution. The model should include a small number of easily controlled input data, describing the terms of observation, geophysical situation and atmosphere status in a proper way. In the course of this approach, the right choice of the atmosphere model structural parameters, which includes time and spatial variations, are very important (Makhov and Sytko 2018). The structural parameter model should also be selfconsistent and use the minimum number of inputs. Exactly this type of approach was chosen for the calculations of seasonal variation and latitudinal temperature distribution over the Earth surface.
The temperature was calculated according to (1) using the one-dimensional heat conductivity equation together with the surface boundary condition, the initial condition and function F(t,φ,λ) expression for the planet surface heating by solar radiation, The equations use the following symbols: t -time, z -height, T -temperature, , K and cdensity, conductivity and planet shell heat capacity.
The boundary condition at the interface with latitude φ and longitude λ is specified by the equation of heat balance between the heat inflow caused by solar radiation and its outflow due to the outgoing thermal radiation of the planet. In the expression for the boundary condition, we use the following parameters: σ and k s -Stefan-Boltzmann constant and planet greyness value. The exponential factor affecting the surface greyness value expðÀ τ T Þ was introduced to consider the attenuation of the thermal radiation emitted by the planet's surface. Exponential function, its form and parameters are determined empirically, proceeding from the condition of the best form of data delivery from (Tverskoy 1962) on the seasonal temperature variation at all latitudes. Herewith, the albedo A distribution was specified according to the map, Figure 2.
The solar zenith angle ζ, atmospheric transmission Tr and atmospheric mass m included in the expression for the heating function F were calculated by formula (2), The expression uses the following symbols: R -observation site radius (R 0 = 6371 км), δ and α -declination and right ascension of the Sun, Ω -the Earth angular velocity, S 0 and E s -sidereal time and the solar constant at a distance of 1 au, r e -the distance from the planet to the Sun and and kabsorption coefficient. Parameters α A = 2.5 and A e = 0.3, and μ equals to 1, if ζ < 90°, and 0 if ζ ≥ 90°. The atmosphere transmittance was calculated for the isothermal altitude distribution of atmospheric density with the height scale H.
The heat conductivity equation solution and daily averaged heating function F d were represented as Fourier series (3), and to represent the averaged heating function, n is the number of harmonics, in most cases, 1-2 harmonics are enough, but at short daylight hours (polar region) 4-5 and sometimes even 6 harmonics are required, Symbol ω is used for the mean velocity of the Earth's rotation around the Sun. The expansion coefficients of the heating function in Fourier series were determined from the expressions, When the temperature and heating function Fourier series are substituted into the heat conduction equation, the representation parameters are as follows: The calculated values of average January and July monthly temperatures latitudinal distribution were compared with the data from (Tverskoy 1962). The comparison showed that calculation provides reasonable results although it represents the surface temperature, not the atmosphere. The atmospheric temperature is specified by experimental data. The most significant discrepancies we observe are at polar latitudes, where atmospheric circulation strong impact was not considered in calculations.
Note that the proposed method of temperature calculation does not consider its daily distribution at all since at one of the stages, the heating function is averaged per day.
Global temperature distribution modelling, using the described way of integral albedo distribution, enables us to create a three-dimensionally inhomogeneous model of the atmospheric density distribution, based on the calculated temperature distribution without any other data involvement.

Colour rendering of spectra in colour systems
There are colour models available for a long time and designed to reproduce realistic graphic computer images. In 1931 and 1964, the International Commission on Illumination (ICI) developed a standard observer with the corresponding colour matching functions. The CIE 1931 Standard Observer is also known as the CIE 1931 2° Standard Observer, and the chromaticity coordinates are based on tables (Wright 1929;Guild 1931). The alternative is CIE 1964 10° Standard Observer, which is derived from the work (Stiles and Burch 1959;Speranskaya 1959).
The CIE chromaticity diagram (Smith and Guild 1931-32;CIE Proceedings 1964) requires all colours visible to the average person to be included. Also, the CIE diagrams serve as a basis for chromaticity calculations in any representation system dealing with the realistic images. In practice, CIE XYZ chromaticity coordinates are the intermediaries between the computed colour spectrum and the RGB chroma render that drives the graphics monitor. There are well-defined routines available for such colour conversions.

Determination of trichromatic coordinates (XYZ)
The first of these procedures is the trichromatic coordinates (XYZ) definition: one should create a colour image of the target according to its radiation calculated spectrum and convert this spectrum to the three-colour XYZ parameters. To determine the XYZ values of the calculated spectral distribution E λ , the characteristic equations are used, where x, ỹ and z -colour matching functions, k -normalization factor (if k= 1, then X, Y and Z are W m −2 ·sr −1 ·nm −1 ) and integration is carried out on the entire visible range, typically 380 to 780 nm. The colour matching function ỹ is actually defined to be identical to the human eye visibility, that is max 683 lm W −1 at λ = 555 nm, i.e. it is possible to calculate the illumination L depending on the spectrum by substituting k= 683 lm W −1 into the equation for the Y-value (units E (λ) -W m −2 sr −1 nm −1 , therefore, the illumination units Llm m −2 sr −1 or cd m −2 .)

Transformation of XYZ colour coordinates to RGB colour space
The next stage is the transfer from the colour functions (XYZ) to the colour spaces (RGB). RGB and XYZ chromaticity diagrams are related by a linear transformation, which matrix form looks like where x r , y and z r are colour functions of red, green and blue phosphors at full output. The method of matrix component calculation for the transfer from colour functions (XYZ) to RGB functions involves the standard colour characteristic application. The most common are NTSC system ('NTSC'), EBU system ('EBU (PAL/SECAM)'), SMPTE system ('SMPTE'), HDTV system ('HDTV'), CIE system ('CIE') and Rec709 system ('CIE REC709').
The matrix element calculation is carried out by the sequential steps: • For the selected colour standard, we use the XYZ-to-RGB conversion matrix., • Then, we define the white spot vector for the emitter used by the standard.
• The preliminary inverse matrix is determined as follows:Ã • The white spot RGB functions are calculated as A rate setting for y w -white spot functions is carried out in accordance with the common convention.
• The final step is the preliminary inverse matrix elements rating in accordance with the expressions A RGB ¼ r x ¼~r x r w r y ¼~r y r w r z ¼~r z r w g x ¼g x g w g y ¼g y g w g z ¼g z The obtained results of conversion matrices from XYZ to RGB colour coordinates for different colour spaces are presented in

A simplified method for colour rendering in different colour systems
In order to model colour images of the Earth observed from different altitudes, we used the calculation results of Earth's outgoing radiation according to the method described in Khokhlov (2001a). This method allows us to identify the specifications of radiation distribution over the disk of the planet, including the terminator region and the transition between the Earth's surface and the atmosphere. As the first step involves only the most general features of the global distribution, not a detailed description of its regional specifics, we did not take into account inhomogeneous nature of the albedo distribution in order to reduce the time for calculating colour images. The calculations were performed using a uniform distribution with albedo value of 0.3, which corresponds to the Earth's global albedo. However, even under these conditions, it took a PC with a clock frequency of 1 GHz about 12 hours to calculate a colour image of the Earth with a radius of 60 pixels (1 pixel corresponded to the Earth's surface area of ~106 × 106 km 2 ). The visualization results of the calculated data for eight colour spaces from Table 1 are presented in Figure 3. The images were calculated under the same conditions. The observation altitude was 376,000 km, and the nadir direction corresponded to a point on the Earth's surface with coordinates 45°N and 30°E. The observation time was 21:00 UT, 20 March 2000.
It can clearly be seen that allowance for radiation brightness adds less intensively illuminated regions to the image, as a result expanding it compared to the images that do not take brightness into account. It is also clearly visible how the use of different colour spaces affects the colour of the images. The most similarity is observed between the images in the CIE colour system.
Hence, the CIE colour system with a higher spatial resolution of 32 × 32 km and no allowance for image brightness (see Table 1) was used for calculating Earth images. The results are presented in Figure 4. The images represent the view of the Earth at 21:00 UT; the date and the observation conditions are the same as in Figure 3. The left image demonstrates a visualization for a wavelength of 550 nm according to the method from (Khokhlov 2001b), the middle one -a visualization of the Earth's radiation in the visible range (380-780 nm). The image on the left contains a calculated optical wedge, where the left side of the upper part corresponds to the minimum brightness, and the right side of the lower part -to the maximum brightness. For comparison, the same figure on the right shows an image of the Earth from (John 1996). The image is a 4° view of the Earth from the altitude of 240,000 km above the equator at 0° longitude, with the Sun overhead at 23°N and 110°E. This corresponds to a late evening (BST) around mid-summer in the Northern Hemisphere. The image only shows the light scattered by the atmosphere. The effect of the   Earth's shadow is clearly visible -the light reaches parts of the atmosphere above the surface points in the darkness. Unfortunately, the data set of the observation conditions is not complete, the image is grey and full comparison is not possible. The simplified method allows us to construct images of the Earth observed from different altitudes. Figure 5 demonstrates some of the calculated images. Each image has a size of 101 × 101 pixels, and it takes about 12-15 minutes to calculate one on a medium-powered PC. The first group (a) contains images of the Earth at 500 nm, calculated from different altitudes at a 94° zenith angle of the Sun. The second group (b) presents images of the Earth from an altitude of 36,000 km for wavelengths of 800 and 2200 nm. A thin streak of dark Earth is visible in some images. Figure 6 shows a 201 × 201 pixel image of the Earth from the near-lunar distance for the latitude of 40° and local midnight time. The atmosphere is clearly visible, especially its 'whiskers' over the dark Earth.
At this stage, we used the colour representation of brightness on a logarithmic scale to visualize the calculation results.

The model of inhomogeneous distribution of the Earth's surface albedo
To design the colour images of the Earth, observed from geostationary orbits, we used the calculated values of the Earth outgoing radiation according to the model described in detail in (Khokhlov (2001b).
We consider a spectrally indifferent, diffusely reflecting, but spatially inhomogeneous Earth's surface, which we describe in terms of integral albedo global distribution, not a spectral one. We use a reference model of integral albedo global distribution, presented in Figure 2 with a resolution of 0.4° х 0.4°. This map is worked out based on physical maps and integral albedo data for various types of the Earth surface (Tverskoy 1962), presented in Table 2. To space the albedo values indicated in Table 1 over the surface we use the Earth physical map. Thus, for each point of the global surface, specified by geographic coordinates, we have the corresponding albedo value. For the result visual reference, we provide a map picture (Figure 2), not in black and white  (256 gradations), but in pseudo colour (768 gradations). This map does not contain any information on the Earth surface colours.
To convert albedo А values to RGB colours, they used parameter 0.3, which corresponds to the albedo value as Parameter ∑ value was converted into RGB colours according to the following rule: The map of the global albedo distribution shown in Figure 2 was converted into a binary file, which was further used for automatic albedo value calculation using geographic coordinates.

Results and discussion
(

1) Colour rendering of Earth images based on the model of radiation brightness of the spherical Earth's atmosphere, taking into account inhomogeneous temperature and surface characteristics
The model is literal and involves three-dimensionally inhomogeneous distribution of atmospheric density in a spherical atmosphere and not only direct solar radiation but also a surface reflected one. The remarkable thing is that the model enables us to calculate the atmosphere brightness without any restrictions of directions, heights and observation time. As the input parameters it uses only the spectral region data, place, direction and Julian date (combining calendar date and time) of observation. This molecular scattering model in association with the experimental data, the assumed substantiation, application areas and limitations are discussed in Khokhlov (2001aKhokhlov ( , 2001b). The described model enables us to reveal the features of the radiation distribution transformation over the planetary molecular atmosphere, including the terminator area and the 'Earthatmosphere' transfer. A typical example of the Earth image calculation with inhomogeneous surface temperature and albedo distribution is shown in Figure 7 in real colours, as it is perceived by human eyes.
In Figure 8, the experimental images are shown for comparison. The calculated images in real colours (top row) are shown together with black and white images from the geostationary satellites METEOSAT, GMS, GOES-W and GOES-E (The Satellite Data Services (SDS) group 2020) (bottom row).
The calculated images were elaborated for the same conditions as the experimental ones. With that, only the place and time of observation were considered. The observation altitude is 36,000 km, the nadir direction corresponds to the equatorial points of the Earth's surface with the longitudes indicated in Table 3 for each satellite. The date of satellite passage is 06.03.2006. The observation time is also indicated in Table 3.
The calculations of colour theoretical images were performed with an allowance for instrumental spectral characteristics of all satellites. Figure 9 shows spectral characteristics for geostationary satellites of the early 20th century generation. These characteristics, covering the visible range, allow us to convert black-and-white images into colour ones; however, they either do not provide a set of input data required for calculations (David Woods and O'Brien 1968), including geostationarity of the orbit, or are published in black-and-white, sacrificing chromaticity in order to increase image resolution in the visible range. In the studies of natural phenomena, pseudo-colour images are typically used in the ultraviolet, near-, mid-and farinfrared ranges, as well as in the thermal range. Numerous visual colour images from low-flying satellites or satellites in polar orbits, such as Landsat, NOAA, Spot, Terra and other regional-scale systems, are available but not considered in this article.  (Tverskoy 1962 It is possible to compare the images only with black and white ones since most of the available colour images are artistic and pseudo-colour or do not contain full information on the calendar date, time, place and orientation of the observation direction.
The calculated and experimental image comparison shows close matching in general Earth's brightness distribution, the continents outlines and random brightness distribution over the surface. Since cloudiness is not considered in the calculations, it is missed. The calculated images clearly reflect orangered colours in the terminator area, related to glow phenomena. The blue colour of the calculated images is associated with solar radiation scattered by molecular constituents of the atmosphere. The most intense colours are over the oceans, where the contribution of the Earth's surface radiation is less due to the small ocean waters' albedo.  Comparison of colour and black-and-white images can be carried out only by the location of the Earth's disk, phase and regional landscape of its surface. These characteristics are visually similar for all four satellites.

Comparison of theoretical and satellite colour images in the visible spectrum
To begin a discussion on the possibility of comparing the realistic colour rendering of theoretical models of the global state of the Earth's atmosphere with satellite images of the Earth view from geostationary altitudes, let us consider the images in Figure 10.
It is interesting to compare the theoretical image from the GMS satellite (Himawari-5) dated 2006/03/ 06 and the image from the modern Himawari 8 satellite (The Himawari-8 Real-time Web 2021) dated 2020/11/26. Figure 10 clearly shows the geometric and seasonal differences, explained by the interval of more than 14 years between the images.
At present, a number of satellites regularly take photographs of the full disk of the Earth. Three satellites, colour images of which are borrowed from West (2017) and shown in Figure 11, attract special attention. The first two satellites (Himawari-8 and GOES-16) are geostationary satellites that take images from an altitude of 35,785 km above the Earth's equator. The third one, DSCOVR (Deep Space Climate Observatory), is located in the orbit between the Earth and the Sun, at a distance of about 1,500,000 km from the Earth at the Lagrange point L 1 .   The geostationary satellite images in Figure 11 have high colour resolution and even greater resolution in black-and-white. The Earth Polychromatic Imaging Camera (EPIC) DSCOVR operates in low-resolution colour (EPIC is a 10-channel spectroradiometer (317-780 nm) onboard the NOAA's DSCOVR spacecraft).
Spectral characteristics or quantum efficiency of the matrices for all three satellites are presented in Figure 12.
Images from Natural Color Imagery (2016) were created using 10 narrow spectral bands from EPIC that are within the visual range. They have colour and brightness adjusted to represent what a conventional camera would produce.
Unfortunately, due to the remoteness of the Earth and incomplete information about the parameters of the mandatory set of input data, the calculation using the model of outgoing Earth radiation for the DSCOVR satellite is currently impossible. For the same reason, the set of input data for calculating images of geostationary satellites has not been formed either.

Conclusions
This article addresses the possibility of colour rendering for the images of luminous planetary objects with gas and dust envelopes. Presented results of performed calculations indicate that there is a significant dependence between the colour palette of the image and the colour system used to convert the calculated data into images. Nevertheless, almost all colour systems can realistically represent the distribution of colour depending on illumination of the Earth by the Sun.
A specific aspect of the proposed method for calculating colour images of the Earth is the use of a spherical model based on physical principles. The model is limited by the choice of Rayleigh scattering processes of solar radiation coming from the Earth's atmosphere. The model takes into account the following processes: • Single scattering of direct solar radiation in the atmosphere;  • Reflection of solar radiation from the inhomogeneous Earth's surface, transmitted and dispersed into space by the atmosphere.
A characteristic feature of the model is the representation of inhomogeneous distribution of the Earth's surface albedo, surface temperature and altitude profile of the density of the main gas components.
The calculated values, presented in this article with an allowance for instrumental spectral distributions, demonstrate a realistic picture of colour brightness distribution of the Earth image depending on its illumination by the Sun. The planet appears light blue in the brightly illuminated part of the Earth's disk and orange-red in the luminous segment of the terminator.
The authors indicate the difficulties associated with comparing theoretical and experimental images of the full disk of the Earth in the visible range in natural colours. The reasons behind those include the incompleteness of information about the observation conditions of satellite colour images and the preference for obtaining the results in black-and-white due to the fourfold increase in the resolution of Earth images. This is especially true for high-altitude satellites and spacecraft.
One application of the material in this article is the use of the Moon's ash light emission as a tool for studying the Earth's global albedo. Measurements of the Moon's ash light are a challenging task currently being undertaken by astronomers in the United States (Big Bear Lake Solar Observatory in California), Russia (Novolazarevskaya Station in Antarctica) and several other countries in the deployment phase.
To control the experimental data in our research on this subject, we use the results of theoretical calculations of images of the Earth as it is seen from the Moon and from a height of 36,000 km above the Novolazarevskaya station. For the same purpose, the experimental and theoretical images of the Moon at the time of the experiment are compared. CAD technologies), the Sketchpad complex, the Whirlwind computer, projects of art and digital installations, developed by long-term efforts of researchers, designers, engineers, as well as artists, inventors, software engineers and IT specialists, have become indispensable necessities in the life of the global community.
The very first CG applications took place in the fields of scientific research, computer video games and video films (William 1965). The first bulky prototypes of computers were used only for solving scientific and industrial problems. Modern diversified scientific CG makes it possible to perform computational experiments and to provide the visual representation of their results.
Nowadays, the most productive and well-known CG applications are associated with film and entertainment industries, where image synthesis is used to create visually appealing special effects for movies and video games. CG transforms the results of solved problems into visual representations of various information; it is widely used in all branches of science, technology, medicine, commercial and management activities. Being an indispensable element of CAD software, CG prevails in the work of design engineers, architects and inventors of new equipment. Design, construction and detailed drafting of workpieces are carried out using CG tools, which allow to obtain both flat and 3D images. Due to television and Internet, CG has also found its application in such spheres as art, illustration, animation and advertising.