Amplitude method for diagnosing the frequency of collisions in the Earth’s ionosphere using the bicharacteristics method

Толық мәтін

Введение

Вопросы дистанционного зондирования атмосферы Земли — это важные и актуальные направления развития современных научных исследований [1, 2]. Одной из главных задач является контроль и диагностика структуры ионосферной плазмы, поскольку ионосфера и магнитосфера Земли оказывают существенное влияние на работу радиосистем различного назначения [3–5]. Это системы радиосвязи, радиолокации, позиционирования и радионавигации. Особо актуальной является задача осуществления мониторинга параметров атмосферы для прогнозирования и предотвращения экстремальных природных явлений.

Настоящая работа посвящена определению эффективной частоты соударений электронов с нейтральными молекулами. Этот параметр определяет поглощение электромагнитных волн в ионосферных слоях и поэтому играет, совместно с распределением электронной концентрации, а также величиной и направлением магнитного поля Земли, важную роль при моделировании характеристик распространения радиоволн на ионосферных трассах [6].

Для определения эффективной частоты соударений электронов в настоящее время разработаны различные методы. Среди них следует отметить ракетные измерения плотности и температуры в нижних слоях ионосферы, по которым можно рассчитать эффективную частоту соударений, а также метод определения эффективной частоты столкновения в области нижней ионосферы по данным радиозатменных измерений [7, 8].

Классические и наиболее логичные подходы к определению эффективной частоты соударений электронов основаны на измерении поглощения радиоволн декаметрового диапазона при прохождении ионосферной плазмы. В монографии К. Дэвиса [9] указаны пять способов измерения поглощения: риометрический метод, импульсный метод, метод непрерывного излучения, минимальной частоты и метод развертывания по частоте [10, 11]. Развивая амплитудный подход к определению эффективной частоты соударений, в работе использована возможность прямого измерения излучения антенной системы и исследования её диаграммы направленности с помощью беспилотных летательных аппаратов [12]. Применение методов математического моделирования распространения декаметровых радиоволн в возмущенной ионосферной плазме [13], основанное на оперативных данных наклонного и вертикального зондирования [14–16], является перспективным методом решения данной задачи.

В работах [17–19] амплитудный подход рассмотрен в случае, когда методом получения данных об ионосфере Земли являлось вертикальное зондирование. В этой статье результаты, полученные [17–19], обобщены на случай слабонаклонного зондирования ионосферы. Актуальность работы связана с активно проводящимися в настоящее время экспериментальными исследованиями по мониторингу ионосферы Земли на слабонаклонных трассах [14].

Работа доложена на Всероссийской открытой научной конференции «Современные проблемы дистанционного зондирования, радиолокации, распространения и дифракции волн», проводимой в рамках Армандовских чтений, в г. Муроме в июне 2024 года [20].

1. Модель электронной концентрации ионосферы. Лучевые траектории

В работе рассмотрена двуслоевая модель электронной концентрации ионосферы N, определяемая формулой:

$\begin{array}{c}N\left(z\right)\approx 7.025426\cdot {10}^4\exp \left[\frac{1}{2}\left(1-\exp {\theta }_1+{\theta }_1\right)\right]+5.6182804\cdot {10}^5\times \\ \times \exp \left[\frac{1}{2}\left(1-\exp {\theta }_2+{\theta }_2\right)\right],\end{array}$ (1)

в которой z — высота,

${\theta }_1=\frac{1}{40}\left(196-z\right)$, ${\theta }_2=\frac{1}{55}\left(263-z\right)$. (2)

Как следует из (2), максимум слоя F1 расположен на высоте 196 км, а слоя F2 — на высоте 263 км.

Для определения лучевых траекторий в такой изотропной среде нами использована бихарактеристическая система Гамильтона-Лукина [21–24]:

$\frac{d\overrightarrow{r}}{dt}=2\overrightarrow{k}{c}^2{\left(\frac{\partial \left(\epsilon {\omega }^2\right)}{\partial \omega }\right)}^{-1}$, $\frac{d\overrightarrow{k}}{dt}=-{\omega }^2\frac{\partial \epsilon }{\partial \overrightarrow{r}}{\left(\frac{\partial \left(\epsilon {\omega }^2\right)}{\partial \omega }\right)}^{-1}$, (3)

где $\epsilon$– эффективная диэлектрическая проницаемость ионосферной плазмы,  $\overrightarrow{k}$— волновой вектор,  $\omega$— круговая частота,  $\overrightarrow{r}$— координаты луча, t — время, c — скорость света.

Источник излучения находится в начале координат (t = 0, x = y = z = 0). Относительно источника излучения необходимо задать начальные условия, которые имеют вид:

$x\left(0\right)=y\left(0\right)=z\left(0\right)=0$, $k_x\left(0\right)=\frac{\omega }{c}\sqrt{\epsilon \left(0\right)}\cos \alpha \cos \beta$, $k_y\left(0\right)=\frac{\omega }{c}\sqrt{\epsilon \left(0\right)}\cos \alpha \sin \beta$, $k_z\left(0\right)=\frac{\omega }{c}\sqrt{\epsilon \left(0\right)}\sin \alpha$. (4)

В формулах (4) $\alpha$ и $\beta$ — это начальные углы выхода луча из источника.

Как известно, эффективная диэлектрическая проницаемость изотропной плазмы имеет вид [9]:

$\epsilon =1-X$, $X={\left(\frac{{\omega }_p}{\omega }\right)}^2$, ${\omega }_p^2=\frac{4\pi e^{2}N\left(z\right)}{m_e}$. (5)

В формулах (5) ${\omega }_p$ — это круговая плазменная частота, e — заряд электрона,  $m_e$— масса электрона.

Основным отличием амплитудного метода восстановления эффективной частоты соударений электронов на слабонаклонных трассах от такого же подхода в случае вертикального зондирования [17–19] является необходимость определения всех лучей, пришедших из передатчика в приёмник, то есть решения задачи «пристрелки». Результат решения этой задачи показан на рис. 1. Предполагалось, что угол $\beta =0$.

На рисунке приведены только те лучи, которые попали из излучателя в приёмник, расположенный справа на расстоянии 100 км. В построении рисунка участвовало 45 лучей. Нижний фиолетовый луч соответствует частоте 1 МГц; верхний красный луч — частоте 6.993 МГц.

Тот же набор лучей приведен на рис. 2 в координатах время–высота, то есть показана зависимость высоты луча от группового времени. Видно, что все лучи, показанные на рис. 1 и приходящие в один и тот же приёмник, попадают туда в разное время.

На рис. 3 показаны зависимости волнового вектора в абсолютных (рис. 3(а)) и относительных (рис. 3(б)) единицах от высоты. На рис. 3(а) кривые образуют аркаду, так как начальное значение компоненты волнового вектора уменьшается с уменьшением частоты.

 

Рис. 1. Траектории лучей, попадающих из источника в приёмник на разных частотах.

 

Рис. 2. Зависимости высоты лучей от группового времени.

 

На рис. 3(б) кривые образуют купол. В отличие от рис. 1, на рис. 3 лучи движутся слева на право. Частота красного, самого верхнего луча близка к частоте просачивания волны сквозь ионосферный слой. Её небольшое увеличение приведет к тому, что луч не придет в приёмник, а пройдёт сквозь ионосферу.

Зависимости вертикальной компоненты $k_z$ волнового вектора от группового времени приведены на рис. 4. Вертикальная компонента до входа в ионосферу остаётся постоянной, затем она начинает убывать и меняет знак. На рис. 3 и 4 введено обозначение: $k_0=\left|\overrightarrow{k}\left(0\right)\right|$.

 

Рис. 3. Зависимость вертикальной компоненты ${\mathbf{k}}_{\mathbf{z}}$ волнового вектора от высоты; а) абсолютные значения; б) относительные значения.

 

Рис. 4. График зависимости вертикальной компоненты волнового вектора от группового времени; а) абсолютные значения, б) относительные значения.

 

На рис. 5 показана фаза радиосигнала, вычисленная в точке приёма по формуле (6)

 $\Phi =\underset{0}{\overset{t}{\int }}{k}_{x}dx+k_{y}dy+k_{z}dz$(6)

в зависимости от группового времени прихода сигнала, а поскольку время прихода растет с ростом частоты, то рис. 5 отражает и эту зависимость. Видно, что с ростом t фаза сначала ускоренно растёт, потом возникает точка перегиба, и уже при приближении к критической частоте (частоте просачивания) рост замедляется.

Точка отражения радиоволны от ионосферы (z_m) определяется из условия: $k_z=0$. Анализируя решения бихарактеристической системы, нетрудно найти зависимость высоты отражения луча z_m от ионосферного слоя как функцию времени прихода луча в точку отражения t_m: z_m = z(t_m). Эта зависимость показана на рис. 6.

Видно, что и высота отражения, и время прихода сигнала в точку отражения одновременно увеличиваются, но высота растет медленнее.

Для каждого луча, показанного на рис. 1, существует своя частота и своя точка отражения. Их можно связать функциональной зависимостью, представленной на рис. 7. Высота отражения плавно растет с частотой: сначала медленно, а после точки перегиба всё быстрее при приближении к максимуму ионосферного слоя.

 

Рис. 5. Зависимости фазы сигнала в точке приёма от группового времени.

 

Рис. 6. Высота отражения луча от ионосферы как функция от времени прихода луча в точку отражения.

 

Рис. 7. Высота отражения луча от ионосферы как функция частоты.

 

Комбинируя рис. 6 и 7, можно установить функциональную зависимость времени прихода луча t_m в точку отражения от частоты f, как это показано на рис. 8. По существу, это ионограмма наклонного зондирования. Видно, что сначала зависимость положительная и линейная, но она переходит в экспоненциальную при приближении к максимуму ионосферного слоя.

 

Рис. 8. Время прихода луча tm в точку отражения от ионосферы как функция частоты f.

 

2. Амплитуда, расходимость и поглощение радиосигнала

Рассмотрим амплитудные характеристики радиосигнала в точке приёма. Амплитуду A радиосигнала в точке приёма можно представить как

$A=E_{0}D{\mathrm{e}}^{-\psi }$. (7)

В выражении (7) E₀ — это амплитуда поля на некотором расстоянии r₀ от источника излучения желательно вне ионосферы, D — расходимость лучей, а  $\psi$— поглощение. Для того чтобы найти расходимость D, необходимо вычислить якобиан J:

$J=\left|\begin{array}{ccc}{x^{\text{'}}}_t& {y^{\text{'}}}_t& {z^{\text{'}}}_t\\ {x^{\text{'}}}_{\alpha }& {y^{\text{'}}}_{\alpha }& {z^{\text{'}}}_{\alpha }\\ {x^{\text{'}}}_{\beta }& {y^{\text{'}}}_{\beta }& {z^{\text{'}}}_{\beta }\end{array}\right|$ (8)

и найти отношение якобианов в начальной точке (r₀) и в конечной точке траектории:

$D=\sqrt{\left|\frac{J_0}{J}\right|}$, $J_0={\left.J\right|}_{r=r_0}$. (9)

Первую строчку в якобиане (8) можно найти из решения бихарактеристической системы Гамильтона–Лукина (3)–(4). Две нижние строчки можно определить из решения расширенной бихарактеристической системы Лукина [25, 26]

$\frac{d{\overrightarrow{k}}_{\alpha }}{dt}=\frac{\partial }{\partial \alpha }\left(\frac{\partial {\omega }^2\epsilon }{\partial \overrightarrow{r}}/\frac{\partial {\omega }^2\epsilon }{\partial \omega }\right)$, $\frac{d{\overrightarrow{r}}_{\alpha }}{dt}=\frac{\partial }{\partial \alpha }\left(\left(2c^2\overrightarrow{k}-\frac{\partial {\omega }^2\epsilon }{\partial \overrightarrow{k}}\right)/\frac{\partial {\omega }^2\epsilon }{\partial \omega }\right)$,

$\frac{d{\overrightarrow{k}}_{\beta }}{dt}=\frac{\partial }{\partial \beta }\left(\frac{\partial {\omega }^2\epsilon }{\partial \overrightarrow{r}}/\frac{\partial {\omega }^2\epsilon }{\partial \omega }\right)$, $\frac{d{\overrightarrow{r}}_{\beta }}{dt}=\frac{\partial }{\partial \beta }\left(\left(2c^2\overrightarrow{k}-\frac{\partial {\omega }^2\epsilon }{\partial \overrightarrow{k}}\right)/\frac{\partial {\omega }^2\epsilon }{\partial \omega }\right)$ (10)

с начальными условиями:

${k^{\text{'}}}_{x\alpha }\left(0\right)=-\frac{\omega }{c}\sqrt{{\epsilon }_0}\sin \alpha \cos \beta$, ${k^{\text{'}}}_{x\beta }\left(0\right)=-\frac{\omega }{c}\sqrt{{\epsilon }_0}\cos \alpha \sin \beta$,

${k^{\text{'}}}_{y\alpha }\left(0\right)=-\frac{\omega }{c}\sqrt{{\epsilon }_0}\sin \alpha \sin \beta$, ${k^{\text{'}}}_{y\beta }\left(0\right)=\frac{\omega }{c}\sqrt{{\epsilon }_0}\cos \alpha \cos \beta$,

${k^{\text{'}}}_{z\alpha }\left(0\right)=\frac{\omega }{c}\sqrt{{\epsilon }_0}\cos \alpha$, ${k^{\text{'}}}_{z\beta }\left(0\right)=0$, $\overrightarrow{r}=(x,y,z),{{\overrightarrow{r}}^{\text{'}}}_{\alpha }\left(0\right)=0$, ${{\overrightarrow{r}}^{\text{'}}}_{\beta }\left(0\right)=0$. (11)

Обычно ослабление сигнала за счет расходимости лучевого потока вычисляют в децибелах. В дальнейшем мы будем пользоваться формулой:

$A_m=-20\lg \sqrt{\left|\frac{J_0}{J}\right|}$. (12)

На рис. 9 показано вычисленное по этой формуле ослабление сигнала за счет расходимости. Рис. 9(а) — это зависимость ослабления от высоты, а рис. 9(б) — это зависимость ослабления от времени вдоль лучей.

При вычислении ослабления предполагалось, что r₀ = 1 км, в отличие от работы [20], где считалось, что r₀ = 1 м. Поэтому значения на графиках в этой работе и в [20] различаются на 60 дБ.

Сначала ослабление для всех частот одинаковое, и оно быстро возрастает с ростом расстояния или времени.

 

Рис. 9. Ослабление радиосигнала Am, обусловленное расходимостью: а) как функция высоты z, б) как функция группового времени t.

 

Потом, когда луч касается каустики, ослабление резко падает (для каждого луча на своей высоте, когда луч приближается к каустике), а далее ослабление опять увеличивается и на земле составляет от 48 до 58 дБ (рис. 9(а)). Из рис. 9(б) видно, что каустики возникают для каждой частоты при определенном значении группового времени. С ростом частоты это значение увеличивается.

На рис. 10 показана расходимость, то есть зависимость корня из отношения якобианов D (рис. 10(а))

$D=\sqrt{\left|\frac{J_0}{J}\right|}$, (13)

и ослабление A_m (рис. 10(б)) как функции частоты f, найденные в точке приёма сигнала.

 

Рис. 10. Зависимости от частоты f: а — расходимости D, б — ослабления радиосигнала A_m.

 

С ростом частоты ослабление сигнала растет, а расходимость D убывает.

Будем считать, что эффективная частота v_e соударений электронов известна и изначально определяется равенством [20]

$\begin{array}{l}{\nu }_e={10}^{b_e},b_e\approx -2.144+16.425\exp \left(-z/85\right)+0.024z-\\ -4.093\cdot {10}^{-5}{z}^2+2.053\cdot {10}^{-8}{z}^3\end{array}$, (14)

в котором высота z выражена в километрах. На рис. 11 показана функциональная зависимость десятичного логарифма эффективной частоты соударений (lgv_e) от высоты:

 

Рис. 11. Эффективная частота соударений электронов от высоты; а — полный профиль, б — фрагмент.

 

Так как плотность нейтральных молекул с высотой уменьшается, эффективная частота соударений электронов также быстро уменьшается с увеличением высоты z, несмотря на рост электронной концентрации (рис. 11(а)). В связи с этим для задач радиосвязи важными являются значения эффективной частоты соударений в нижней ионосфере (рис. 11(б)), которые и позволяет определить амплитудный метод.

Известно, что [17, 27]:

$\psi \left(t_k\right)=-\frac{1}{2}\omega \underset{0}{\overset{t_k}{\int }}{\epsilon }_m\left(z\right(t)dt$, (15)

где

${\epsilon }_m\cong -\frac{XZ}{1+Z^2}$, (16)

а

$Z=\frac{{\nu }_e\left(z\right(t\left)\right)}{\omega }$. (17)

Интеграл в формуле (15) вычисляется вдоль лучевой траектории.

На рис. 12 показана зависимость поглощения (в неперах) от частоты в точке приёма.

 

Рис. 12. Зависимость поглощения $\mathbf{\psi }$ от рабочей частоты f.

 

Видно, что поглощение с ростом частоты вплоть до точки просачивания уменьшается. Так, на частоте 1 МГц амплитуда за счет поглощения падает в 9 раз, а на частоте 6 МГц — только в 1.3 раза.

 

Рис. 13. Зависимость амплитуды A от рабочей частоты f.

 

На рис. 13 приведена зависимость амплитуды радиосигнала, рассчитанная по формуле (7). При этом предполагалось, что [28]

$E_0=\frac{\sqrt{30W}}{r_0}$, (18)

а мощность W = 1 кВт.

Амплитуда имеет максимум в окрестности 2.5 МГц, так как, с одной стороны, расходимость растет с частотой, но с другой стороны — поглощение уменьшается.

3. Определение эффективной частоты соударений по амплитудным данным

Выполним численный эксперимент. Предположим, что амплитуда сигнала, начальная амплитуда E₀ и расходимость известны. Найдем поглощение и восстановим эффективную частоту соударений электронов. Так как, с одной стороны, зная амплитуду, поглощение можно найти из формулы (7) (см. рис. 12)

$\psi =-\ln \frac{A}{D{E}_0}$, (19)

а с другой стороны, поглощение определяется по формуле (15) и в неё входит эффективная частота соударений через ${\epsilon }_m$, которая теперь считается неизвестной, можно получить интегральное уравнение:

$\frac{\omega }{2}\underset{0}{\overset{t_k}{\int }}{\epsilon }_{m}dt=\ln \frac{A}{D{E}_0}$. (20)

Предположим, что [17, 27]

${\epsilon }_m\approx -XZ$ (21)

и введём обозначения

$V_e=\frac{4\pi e^2}{m_e}$, $L=-\ln \frac{A}{D{E}_0}$. (22)

Тогда вместо уравнения (20) получаем:

$\frac{V_e}{{\omega }^2}\underset{0}{\overset{t_m}{\int }}{v}_{e}Ndt=L$. (23)

Это интегральное уравнение Вольтерра второго рода. Верхний предел интегрирования t_k = t_m — это время, за которое сигнал проходит расстояние от источника до точки отражения от ионосферы. Оно равно половине времени распространения сигнала из излучателя в приёмник и поэтому легко измеряется. Время t_m является функцией частоты. Величина L также является функцией частоты f и может быть найдена в точке приёма радиосигнала.

Эффективная частота соударений v_e и электронная концентрация N являются функциями высоты z. В свою очередь, высота z зависит от группового времени t вдоль лучевой траектории. Так как лучевые траектории для различных частот разные, то можно считать, что функция z_f(t) также зависит от частоты в (23) и вычисляется вдоль лучевой траектории.

Применим метод итераций [17, 18] для решения уравнения (23). Введем обозначение $\Theta ={\nu }_{e}N{V}_c$. Так как электронная концентрация равна нулю ниже определённой высоты, то можно предположить, что до некоторой частоты f₀, для которой время t_m равно t_m0, выполняется равенство:

$\underset{0}{\overset{t_{m0}}{\int }}\Theta dt=0$. (24)

Весь частотный интервал (f₀, f_n) разделим на n равных частей (f₀, f₁,…, f_j,…, f_n). Частоту f_n выберем чуть меньше частоты просачивания волны через ионосферный слой. Функцию $\Theta$ будем считать постоянной и равной ${\Theta }_j$ на каждом интервале. Для каждой частоты f_j найдём время прохождения волны от источника до точки отражения t_mj. В результате получим цепочку равенств:

${\Theta }_1=L_1\frac{{\omega }_1^2}{t_{m1}-t_{m0}}$, ${\Theta }_2=L_2\frac{{\omega }_2^2}{t_{m2}-t_{m1}}-{\Theta }_1\frac{t_{m1}-t_{m0}}{t_{m2}-t_{m1}}$, ${\Theta }_j=L_j\frac{{\omega }_j^2}{t_{m,j}-t_{m,j-1}}-\sum _{i=1}^{j-1}{\Theta }_i\frac{t_{m,i}-t_{m,i-1}}{t_{m,j}-t_{m,j-1}}$, $j\ge 2$. (25)

Решениями системы (25) являются величины ${\Theta }_j$. Будем считать, что функциональная зависимость электронной концентрации N от высоты z известна — например, восстановлена по зависимости задержки радиосигнала от частоты излучения. Тогда по формуле

${\nu }_{ej}=\frac{{\Theta }_j}{N\left(z_f\right(t_{mj}\left)\right)V_E}$ (26)

находим эффективную частоту соударений как функцию высоты.

Описанный нами алгоритм восстановления эффективной частоты соударений электронов, обобщает на случай наклонного зондирования подход, описанный в [17, 18] для случая вертикального зондирования.

4. Численные результаты восстановления высотного профиля эффективной частоты соударений

Результат восстановления высотного профиля эффективной частоты соударений электронов показан чёрной линией на рис. 14. Жёлтой широкой линией показан исходный профиль эффективной частоты соударений электронов. Видно, что на высотах более 100 км результат реконструкции показывает приемлемую точность.

 

Рис. 14. Сопоставление исходного и восстановленного высотных профилей эффективной частоты соударений.

 

Построим методом наименьших квадратов по восстановленным данным две аппроксимации: двухпараметрическую

$\lg {\nu }_e\approx 0.1731+\frac{487.069}{z}$, (27)

и четырёхпараметрическую

$\lg {\nu }_e\approx 2.64094-\frac{483.933}{z}+\frac{119511}{z^2}-\frac{4.60491\times {10}^6}{z^3}$. (28)

Аппроксимация, полученная по формуле (27), на рис. 15 показана зелёной линией, а аппроксимация, полученная по формуле (28), чёрной линией.

 

Рис. 15. Сопоставление высотных профилей эффективной частоты соударений, полученных по аппроксимационным формулам, с исходными данными.

 

Видно, что графики, построенные по аппроксимационным формулам, хорошо ложатся на жёлтую широкую кривую (рис. 15), причём аппроксимация (27) лучше работает в нижней части кривой, и аппроксимация (28) — в верхней.

Из анализа рис. 14 и 15 следует, что аппроксимационные формулы позволяют экстраполировать оценки эффективной частоты соударений на более низкие высоты, чем те, для которых получены оценки при непосредственном вычислении.

На рис. 16 выполнено сопоставление исходной эффективной частоты соударений (красная линия) и аппроксимации, полученной по формуле (28) (черная линия). Из рис. 16 следует, что аппроксимационная формула неплохо работает даже при больших высотах (до 450 км). Существенные отличия возникают на высотах за пределами главного максимума ионосферного слоя.

Заключение

Создан амплитудный метод восстановления эффективной частоты соударений электронов на основе данных об ослаблении радиосигналов при слабонаклонном зондировании ионосферы Земли. Координаты лучевых траекторий и волновые векторы рассчитаны методом бихарактеристик Гамильтона-Лукина. Для вычисления расходимости радиосигнала применена расширенная бихарактеристическая система Лукина. Для определения лучей, пришедших из источника в приёмник, решена задача «пристрелки». Построены графики лучевых траекторий вертикальной компоненты волновых векторов вдоль лучей как функции высоты. В точке приёма радиосигнала получены зависимости высот отражений лучей от времени прихода лучей в точки отражений, высот отражений лучей от частоты, а также зависимости времени прихода лучей в точки отражений от частоты. В точке приёма рассчитаны расходимость, поглощение и амплитуда радиосигнала как функции частоты. Установлено, что частотная зависимость амплитуды имеет максимум, так как расходимость с частотой растёт, а поглощение — падает.

 

Рис. 16. Сопоставление высотного профиля эффективной частоты соударений, полученного по четырёхпараметрической аппроксимационной формуле (28), с исходными данными.

 

В результате проведения численного эксперимента получено хорошее совпадение расчётных значений эффективной частоты соударений и исходных данных. Исследованы возможности экстраполяции результатов, полученных при восстановлении частотного профиля, для оценки эффективной частоты соударений электронов в ионосферной плазме.

Работа выполнена при поддержке гранта Российского научного фонда № 20-12-00299-П.

Авторлар туралы

A. Kryukovsky

Russian New University

Хат алмасуға жауапты Автор.
Email: kryukovsky56@yandex.ru

д. ф.- м. н., профессор

Ресей, Moscow

E. Mikhaleva

Russian New University

Email: kryukovsky56@yandex.ru
Ресей, Moscow

D. Rastyagaev

Russian New University

Email: kryukovsky56@yandex.ru

к. ф.- м. н., доцент

Ресей, Moscow

Әдебиет тізімі

Қосымша файлдар