Расчет интенсивности света, рассеянного системой экзосфера плюс поверхность,
основан на методе инвариантного погружения, применение которого подробно описано 
Mishchenko и Travis - Satellite retrieval of aerosol properties over ocean 
using polarization as well as intensity of reflected sunlight // 
J. Geophys. Res. 1997. V. 102 (D14). P. 16989–17013), 
а необходимые компьютерные коды можно найти на сайте 
https://www.giss.nasa.gov/staff/mmishchenko/brf/.

Программы, выложенные здесь - rt_refl_agr4.f и rt_refl_sph1.f - адаптированы для
вычисления спектров экзосферы заданной оптической толщины (на длине волны 0.55 мкм) 
над поверхностью с заданным спектральным ходом альбедо.
Они отличаются только "входным блоком", который читает файлы, содержащие необходимые 
для расчетов коэффициенты разложения характеристик однократного рассеяния 
агрегатов (4 варианта) или сферических частиц, соответственно (см. Single_Scat).

Для пересчета оптической толщины, заданной на одной длине волны, 
на каждую из длин волн, необходимо заранее подготовить файлы, содержащие 
значения эффективности рассеяния данных частиц в зависимости от длины волны. 
Файлы такого типа, "Qsca_F...dat", приведены в качестве примера. 
Они дополнительно содержат значения фазовой функции в направлении обратного рассеяния 
и альбедо однократного рассеяния на каждой из длин волн. 
 
Пример необходимой для расчетов спектральной зависимости альбедо поверхности
содержится в файле S_aver_int_91.dat. Принято, что поверхность рассеивает свет 
изотропно (закон Ламберта).

Результаты вычислений интенсивности света, рассеянного средой над поверхностью в 
каждой из длин волн, записываются в файл в виде зависимости от азимута и
косинусов зенитных углов падения и отражения света.
Примеры выходных файлов - Rho_ice-...dat и Rho_sil-...dat. 

Число значений азимутов в программе (NZ0) можно изменять произвольно (один раз в начале), 
а число зенитных углов (N0=NG) соответствует числу узлов в квадратурной формуле Гаусса.
Поэтому, если надо изменить число зенитных углов, необходимо внести изменения во все процедуры 
кода, где есть этот параметр, и скорректировать размерность соответствующих массивов.
========================================================================================
 
Поскольку используемый метод расчета переноса излучения работает для системы плоскопараллельных 
слоев, полученные значения интенсивности относятся к небольшим областям на поверхности объекта,
наблюдаемого интегрально. 
Если принять его форму шарообразной, то его яркость можно получить интегрированием рассчитанных
значений яркости по поверхности сферы с учетом заданного фазового угла наблюдений.
Эта процедура выполняется программой read_write_Rho1.for для каждой из длин волн, 
и получаемый спектр записывается в файл типа B0_sil-...dat, приведенный в качестве примера.
   
========================================================================================

Результаты вычислений интенсивности света, рассеянного экзосферой, состоящей из частиц определенного вида,
над поверхностью при разной геометрии наблюдений запакованы в соответсвующих папках "Rho..."

========================================================================================

Завершающий блок описанных выше кодов может быть изменен для вычисления 
сферического и геометрического альбедо астероида с экзосферой - программа rt_refl_alb_wl.f
Выходной файл этой программы содержит значения 
- сферического и геометрического альбедо астероида с экзосферой,
- геометрического альбедо астероида без экзосферы и
- альбедо поверхности в зависимости от длины волны 
Пример - файл Alb_ice-t05_S&72-300r12_C.dat

========================================================================================