Например, Бобцов

АЛГОРИТМ РАСЧЕТА ЗЕЙДЕЛЕВЫХ АБЕРРАЦИЙ ДЛЯ ОПТИЧЕСКОЙ СРЕДЫ С РАСПРЕДЕЛЕННЫМ ПОКАЗАТЕЛЕМ ПРЕЛОМЛЕНИЯ

64
УДК 535.317
А. Л. СУШКОВ
АЛГОРИТМ РАСЧЕТА ЗЕЙДЕЛЕВЫХ АБЕРРАЦИЙ ДЛЯ ОПТИЧЕСКОЙ СРЕДЫ
С РАСПРЕДЕЛЕННЫМ ПОКАЗАТЕЛЕМ ПРЕЛОМЛЕНИЯ
Рассмотрен аппарат расчета хода параксиальных лучей в осесимметричных неоднородных оптических средах с изменением показателя преломления по трем координатам. На примере линзы со сфероконцентрической неоднородностью показателя преломления показана практическая возможность расчета параксиальных отрезков и коэффициентов аберраций третьего порядка оптической системы с неоднородными линзами. Ключевые слова: неоднородный показатель преломления, коэффициенты аберраций, линза со сфероконцентрическим градиентом показателя преломления.
Известные численные методы расчета хода реальных лучей через осесимметричные неоднородные среды связаны с решением дифференциальных уравнений второго порядка [1]. При этом анализ влияния коэффициентов полинома, описывающего распределение показателя преломления (РПП), и толщин градиентных линз на аберрации оптической системы затруднен из-за большого объема вычислений.
Цель настоящей статьи — продемонстрировать возможность введения в инженерную практику универсального алгоритма расчета хода параксиальных лучей и коэффициентов аберраций третьего порядка систем, включающих неоднородные оптические элементы.
Теория аберраций третьего порядка неоднородного оптического элемента не позволяет точно оценить величину аберрации, но обеспечивает получение намного большего объема информации относительно влияния поверхностей линз и неоднородных оптических сред на аберрации третьего порядка.
Как известно, теория аберраций третьего порядка позволяет вычислить величины углов наклона и высоты первого и второго вспомогательных лучей. Для однородных систем формулы расчета коэффициентов аберраций третьего порядка широко применяются в программах анализа оптических систем [2, 3].
В работах [4—6] представлен аппарат расчета аберраций третьего порядка, построенный на вычислении параметров первого и второго вспомогательных лучей, проходящих через однородные и неоднородные осесимметричные оптические среды с произвольным законом распределения показателя преломления.
Условием расчета коэффициентов аберраций в системе координат OXYZ, где ось Z — оптическая ось, является наличие известной функциональной зависимости углов наклона и высот первого и второго вспомогательных лучей в неоднородных средах.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

Алгоритм расчета зейделевых аберраций оптической среды

65

В работе [5] функции высот и углов представлены в виде числовых рядов по координате z, которые обозначим через h(z), H(z) и α(z), β(z). В неоднородной среде угол наклона луча

является производной от функции высоты луча: α ( z) = h ( z) и β( z) = H ( z) .

Согласно [4, 6] параксиальная высота y(z) произвольного луча может быть записана в виде линейной функции через его высоту на предмете (T) и входном зрачке (Q):

y(z) = h(z)Q + H (z)T ,

(1)

где Q и T — нормированные координаты луча в плоскости входного зрачка и плоскости предмета.
Ход первого и второго вспомогательных лучей в системе координат, отнесенной к первой поверхности линзы, показан на рис. 1.

n1=1 Т

h1

Q n2= n00 –α1 β1=1
nz

n3=1

–sп –s H1

Рис. 1

В общем случае функция РПП оптической среды представляется в виде степенного ряда [5]

n(z, y) = n0 (z) + n1(z) y2 + n2 (z) y4 + ... ,

(2)

где коэффициенты n0 ( z ) , n1(z) , n2 (z) , в свою очередь, также представлены степенными ря-

дами:

n0 (z) n1 ( z )

= =

n00 n10

+ +

n01z n11z

+ +

n02 z2 n12 z2

+ ... ,⎫

+

...

,

⎪⎪ ⎬

n2 (z)

=

n20

+

n21z

+

n22 z2

+ ...

⎪ ⎭⎪

(3)

Контроль сходимости каждого из степенных рядов (3) производится путем сравнения

величин параксиального инварианта I для каждой поверхности оптической системы. Обоб-

щенный параксиальный инвариант I для k-й поверхности определяется как

I = nk ( Hk αk − hkβk ), k = 1, p .

(4)

Для первой и второй производных функции высоты луча y(z) (см. формулу (1)) справед-

ливы выражения

y(z) = h(z)Q + H (z)T , y(z) = h(z)Q + H (z)T .

(5)

Дифференциальное уравнение для хода луча [5] в меридиональной плоскости в системе

координат, отнесенной к предмету, имеет вид

∂n( z, ∂z

ξ)

(1+

y2

)

y

+

n(z,

ξ)

y



(1+

y2

)

∂n(z, ∂y

ξ)

=

0

,

(6)

где ξ = y2.

Подстановка уравнений (5) в выражение (6) с исключением величин высшего порядка

малости позволяет получить общую запись уравнения для обоих вспомогательных лучей.

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

66 А. Л. Сушков

Для первого вспомогательного луча при y(z)=h(z) выражение (6) преобразуется к виду

( ) ( )h(z) n00 + n01z + n02z2 + ... + h(z)(n01 + 2n01z + ...) − 2h(z) n10 + n11z + n12z2 + ... = 0 ⋅ (7)

С учетом того, что уравнение (7) должно удовлетворяться для всех Q и T, коэффициенты при них (в силу независимости величин Q и T) должны быть равны нулю [5].
Решение (7) получено при представлении функции высоты луча h(z) в виде ряда

h(z) =




Aν zν

=

A0 + A1z + A2 z2

+ ... ,

ν=0

(8)

где A0 и A1 определяются высотой луча и углом его наклона к оптической оси на входе в не-

однородную среду.

Первая и вторая производные функции высоты луча h(z) определяются как

h(z) h(z)

= =



ν=0


ν=0

(ν (ν

+1) Aν+1zν + 2)(ν +1)

= A1 + Aν+2 zν

2A2 z + 3A3z2 = 2A2 + 6A3z

+ ... , +12 A4z2

+

...

⎫ ⎪⎪ ⎬ ⎪ ⎭⎪

Подставив выражения (9) в уравнение (7), получим

(9)




(ν + 2)(ν +1) Aν+2 zν

⎡⎣n00

+ n01z + n02z2 ⎤⎦ +




(ν +1) Aν+1zν[n01 + 2n02z] −

ν=0 ν=0

−2




Aν zν ⎡⎣n10 + n11z + n12 z2 ⎤⎦ = 0 .

ν=0

(10)

Приравняв коэффициенты различных степеней в тождестве (10) к нулю, получим в об-

щем случае рекурсивное соотношение для Ат:

Am

=

⎧⎨m∑−1 ⎩ν=2

⎣⎡2

Aν−2

n1,m−ν



ν(ν

−1) Aνn0,m−ν







m

+ 1)(ν

−1) Aν−1n0,m−ν+1⎤⎦



−(m −1) Am−1n01 + Am−2n10} m(m −1)n00 .

(11)

После задания показателей преломления по формулам (3) вычисляются коэффициенты

Am и определяются степенные ряды для высоты любого параксиального луча. Угол преломления луча определяется по известным формулам параксиальной оптики

однородных сред с учетом кривизны поверхности и с использованием величин n00 и n0′0 в

качестве показателей преломления.

Для проверки сходимости рядов (8) и (9) для каждой поверхности по формуле (4) вы-

числяется параксиальный оптический инвариант. Если значения инвариантов идентичны, то

ряды сходятся и расчет продолжается. При несовпадении значений инвариантов необходимо

увеличить точность расчета путем введения большего числа членов в разложении h(z).

В системе координат, отнесенной к поверхности линзы, для коэффициентов аберраций

третьего порядка осесимметричной неоднородной оптической системы выражения известны [7].

Коэффициенты аберраций являются суммой двух составляющих, обусловленных преломле-

нием луча на поверхности ( Si,k ) и прохождением луча через неоднородную среду ( Si,k ):

Si = Si,k + Si,k ,

(12)

где Si , i = I…V, — коэффициент аберраций третьего порядка оптической системы: SI — сфе-
рическая аберрация, SII — кома, SIII — астигматизм, SIV — кривизна поля изображения, SV — дисторсия.

Формулы для вычисления составляющей Si,k имеют следующий вид:

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

Алгоритм расчета зейделевых аберраций оптической среды

67

( )p
SI = ∑

hk Pk + Kk hk4

;

k =1

⎫ ⎪ ⎪

SII

=

p


⎡ ⎢

hk

Pk

k =1 ⎣

⎛ ⎜ ⎝

δβk δαk

⎞ ⎟

+



Kk hk3Hk

⎤ ⎥

;



SIII

=

p⎡


k =1

⎢ ⎢⎣

hk

Pk

⎛ ⎜ ⎝

δβk δαk

⎞2 ⎟ ⎠

+



Kk

hk2

H

2 k

⎥ ⎦⎥

;

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪

SIV

=

p⎛ −∑ ⎜
k =1⎝

δµk rk

⎟⎞ ; ⎠

⎪⎪ ⎬ ⎪



SV

=

p

k =1

⎪⎨⎪⎩⎧⎣⎢⎢⎡hk

Pk

⎛ ⎜ ⎝

δαk δµk

⎞2 ⎟ ⎠

++

J2 hk

δ(αk

nk

)

⎤ ⎥

δβk

nk nk+1 ⎦⎥ δαk

+

Kk

hk

H

3 k

⎪⎫ ⎬

;

⎪⎭

⎪ ⎪ ⎪ ⎪

Kk

=

4δn1, k rk

+

δ n0,k rk2

δαk = α′k − αk , δβk

, Pk = β′k

=

⎛ ⎜ ⎝

δαk δµk

− βk ; J

⎞2 ⎟

δ(αk

µk

),



µk

1 = nk

;

= −n1(s1 − sп )α1β1, s1 ≠

∞;

J

=

− n1h1β1 ,

s1

=

⎪ ⎪ ⎪ ∞,⎭⎪⎪

(13)

где s1, sп — расстояние от первой поверхности до предмета и входного зрачка оптической

системы; αk , βk α′k , β′k — величины углов падения и преломления первого и второго вспо-

могательных лучей для k-й поверхности; hk, Hk — высоты первого и второго вспомогательных лучей на поверхности с радиусом кривизны rk; Pk Kk — коэффициенты, учитывающие

соответственно однородную и неоднородную природу оптической среды.

Составляющая Si,k определяется следующими интегральными выражениями, где интегрирование ведется в пределах от 0 до d (d — толщина линзы):

p⎡
( ) ∫ ( )SI = ∑ ⎢∆

n0,k hk α3k



k =1 ⎣⎢

8n2,k

+1hk4+1

+

4n1,k

+1hk2+1α

2 k +1



n0,k

+1α4k +1

⎤⎫

dz⎥ ;⎪

⎥⎦

⎪ ⎪

( ) ∫(SII

=

p

k =1

⎡⎣∆

n0,k hk α2kβk



8n2,k+1hk3+1Hk +1 +

⎪ ⎪ ⎪

( ) )+ 2n1,k+1hk+1αk+1

hk +1βk +1 + Hk +1αk +1

− n0,k +1α3k+1βk+1

dz

⎤ ⎦

;

⎪ ⎪ ⎪

( ) ∫ (SIII

=

p

k =1

⎡ ⎣



n0,k hk αkβ2k



8n2,k

+1hk2+1H

2 k +1

+

⎪ ⎪⎪ ⎬

)+ 4n1,k+1hk +1Hk+1αk+1βk+1 − n0,k+1α2k+1β2k+1 dz⎤⎦ ;

∑ ∫p
SIV = −2
k =1

n1,k +1 n02,k +1

dz;

⎪ ⎪ ⎪ ⎪ ⎪ ⎪

( ) ∫(SV

=

p

k =1

⎣⎡∆

n0,k hkβ3k



8n2,k

+1hk

+1H

3 k +1

+

⎪ ⎪ ⎪ ⎪

( ) )+ 2n1,k+1Hk+1βk+1 hk+1βk+1 + Hk+1αk+1 − n0,k+1αk+1β3k+1 dz ⎤⎦ ,

⎪ ⎪⎭

(14)

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

68 А. Л. Сушков
где ∆ — разность значений выражений в скобках, соответствующих параметру луча перед преломлением на второй поверхности и после преломления на первой поверхности.
При подстановке значений высот и углов наклона лучей (см. формулы (8), (9)) в систему уравнений (13) получаем значения составляющей Si,k , а при численном интегрировании
функций углов и высот — составляющей Si,k . В сферической системе координат линейное сфероконцентрическое распределение по-
казателя преломления согласно [8] имеет вид

n (r − ρ) = nρ0 + nρ1 (r − ρ) ,

(15)

где nρ0 — показатель преломления материала неоднородной среды в центре сферической системы координат, r — радиус сферического РПП, ρ — сферическая координата, nρ1— линейный коэффициент РПП.
После подстановки значений r в формулу (15) получаем

n(ρ) = nρ0,пр + nρ1,прρ .

В последующих формулах для упрощения записи под величинами nρ0 , nρ1 будем пони-

мать их приведенные значения.

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

вой поверхности линзы, согласно [8] имеем

n00 = nρ0 ,

n01 = nρ1 ,

n10

=

−nρ1

1 2r

,

n11

= nρ1

1 2r2

,

n20

=

nρ1

1 8r3

,

n21

=

nρ1

3 8r 4

.

Показатель преломления на оси z на выходе из неоднородной среды линзы определяется

по формуле

nz = n00 + n01d .

Неоднородный ПП задается полиномом (2), где функции n0( z), n1(z), n2(z) представляют собой линейные зависимости:

n0 ( z ) = n00 + n01z , n1 ( z) = n10 + n11z , n2 ( z ) = n20 + n21z .

Угол преломления луча на поверхности вычисляется по известной формуле для одно-

родных сред:

α′

=

h

1 R

(

n′

− n)
n′

+

αn

,

(16)

где R — радиус кривизны поверхности линзы.
Поскольку коэффициенты А0= h1, A1= –α2 (см. формулу (8)), для коэффициентов полиномов высоты луча и угла его наклона согласно (11) получим

A2

=

2 A0n10 − A1n01 2n00

,

A3

=

A1n10

+

A0 n11 3n00



A2 n01

,

A4

=

2 A2n10

− 2 A1n11 12n00

− 9 A3n01

,

A5 =

A3n10 − A2n11 − 8A4n01 10n00

и т. д.

Высота h2 луча на второй поверхности согласно формуле (8) и угол α2 на выходе из неоднородной среды согласно формуле (9) определяются полиномами

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

Алгоритм расчета зейделевых аберраций оптической среды

69

h2 = A0 + A1d + A2d 2 + A3d 3 + A4d 4 + A5d 5 ,

(17)

α2 = −( A1 + 2A2d + 3A3d 2 + 4A4d 3 + 5A5d 4 ) .

Угол осевого луча в пространстве изображений согласно выражению (16) вычисляем по

формуле

α3

=

h2

1 R2

(1 − nz )

+

α2nz

.

Заднее фокусное расстояние и задний фокальный отрезок линзы рассчитываются по из-

вестным зависимостям: f ′ = h1 / α3 и sF′ ′ = h2 / α3 . Расчет высоты второго вспомогательного луча осуществляется по формулам (17) с за-

меной коэффициентов Aν на Bν:

H(z) =




Bν zν

= B0 + B1z + B2 z2 + ...

ν=0

(18)

при следующих исходных данных: H1=sп, β1=1; коэффициенты B0, B1, B2 и т.д. вычисляются по формуле (18), при этом B0=H1, B1= –β2.
Значения параксиального инварианта для первой и второй поверхностей вычисляются

как

( )I1 = n1 ( H1α1 − h1β1 ) , I2 = nz H2α2 − h2β2 ,

Функции углов и высот представляются полиномами

α2 ( z ) = −( A1 + 2 A2 z + 3A3z2 + 4 A4 z3 + 5A5 z4 + ...), h2 ( z ) = A0 + A1z + A2 z2 + A3 z3 + A4 z4 + A5 z5 + ... ;

(19)

β( z ) = −(B1 + 2B2 z + 3B3z2 + 4B4 z3 + 5B5 z4 ), H ( z ) = B0 + B1z + B2 z2 + B3 z3 + B4 z4 + B5 z5 .

(20)

Подстановка формул (19) и (20) в интегральные выражения (14) и численное интегриро-

вание позволяют получить численные значения коэффициентов аберраций SI…SV.

Рассмотрим пример расчета коэффициентов аберраций SI и SII одиночной линзы с ли-

нейным сфероконцентрическим РПП. Расчет производился при следующих параметрах лин-

зы (рис. 2): r=12,792 мм, R1 =12,792 мм, R2= 197,706 мм, d=1 мм, nρ0=1,65, nρ1=0,031551 мм–1,

sп=0, f ′=20,00 мм, s′F′ =19,376 мм.

n2= n00

r ρ

R1 R2
Рис. 2
В результате расчета углов наклона и высот первого и второго вспомогательных лучей по формулам (19), (20) при исходных α1=0, β1=1, h1= f ′ получим:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

70 А. Л. Сушков

α2 =0,615964, β2 = 0,594224, h2=19,376600,

α3=1,000000, β3= 1,001287, H2=–0,600180. Значения параксиальных инвариантов для первой и второй поверхностей:
I1= –20,001824, I2 = –20,001825.
Расчет коэффициента SI. Для однородной и неоднородной сред коэффициент SI (см.
выражения (13)) рассчитывается соответственно как

SI′ = h1P1 + h2 P2 ; SI′′ = K1h14 + K2h24 , где

K1

=

4n10 R1

+

n01 R12

,

K2

=



4(n10 + n11z) R2



n01 R22

.

Согласно выражениям (14) для составляющей Si,k коэффициенты вычисляются сле-

дующим образом:

SI,1 = nz h2α32 − n2h1α32 ;

SI,2

=

d
∫ n0

( z ) α4

( z) dz ;

0

SI,3

=

d


−4n1

(z)

h2

( z ) α2

(

z)

dz

,

SI,4

=

d


−8n2

( z ) h4

( z) dz

.

00

По суммарному значению коэффициента SI, определяемому как

SI = SI′ + SI′′ + SI,1 + SI,2 + SI,3 + SI,4 ,

вычисляется величина продольной сферической аберрации для края зрачка диаметром 5 мм:

∆s′

=



1 2

tg2uSI

,

∆s′ =0,000 мм,

где u — апертурный угол линзы в пространстве изображений.

Расчет коэффициента SII. Для однородной среды коэффициент

для неоднородной —

SI′I

=

h1P1

∆β1 ∆α1

+

h2 P2

∆β2 ∆α2

,

SI′′I = K1h13H1 + K2h23H2 .

Для составляющей Si,k коэффициенты рассчитываются как

SII,1 = nz h2α22β2 − n2h1α22β2 ;

SII,2

=

d


n0

( z)α3

( z)β( z) dz

;

0

d
SII,3 = ∫ −2n1 ( z) h ( z) α ( z) ⎣⎡h ( z)β( z) + H ( z) α ( z)⎦⎤ dz , 0

SII,4

=

d


−8n2

( z) h3

(z)

H

( z) dz

.

0

Суммарная величина коэффициента SII:

SII = SI′I + SI′′I + SII,1 + SII,2 + SII,3 + SII,4 .

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

Алгоритм расчета зейделевых аберраций оптической среды

71

Коэффициент SI SI′ = 28,072 SI′′ = –27,182 SI,1 = 0,606
SI,2 = 0,254
SI,3 = 0,776
SI,4 = –2,527 SI = 0,000

Таблица 1 Коэффициент SII
SI′I = –0,744 SI′′I = –0,114 SII,1 =0,203
SII,2 =0,244
SII,3 = 0,366
SII,4 = 0,039 SII = –0,004

Анализ составляющих коэффициентов SI и SII, приведенных в табл. 1, показывает, что исправление сферической аберрации произошло за счет взаимной компенсации коэффициен-

тов SI′, SI′′I , имеющих разные знаки, а исправление комы — за счет компенсации коэффи-

циентов SI′I , SI′′I и SII,1, SII,2 , SII,3 , SII,4 . Суммарные величины коэффициентов SI и SII имеют

близкие к нулю значения, что свидетельствует об апланатической коррекции аберраций. Результаты расчета в программной среде OPAL продольной ( ∆s′ ) и поперечной ( ∆y′ ) сфериче-

ской аберраций, а также отклонения от условия изопланазии (η) неоднородной и однородной линз приведены в табл. 2 и 3.
Таблица 2

mзр, мм

∆s′,мм

∆y′, мм

η, %

2,500 2,165 1,767 1,250 0,000

–0,0014 –0,0007 –0,0003 –0,0000 0,0000

–0,0001 –0,0000 –0,0000 –0,0000 0,0000

–0,0072 –0,0041 –0,0018 –0,0005 0,0000

mзр, мм
2,500 2,165 1,767 1,250 0,000

∆s′,мм
–0,2306 –0,1724 –0,1145 –0,0571 0,0000

∆y′, мм
–0,0279 –0,0180 –0,0097 –0,0034 0,0000

Таблица 3
η, %
–0,0950 –0,0682 –0,0435 –0,0208 0,0000

Анализ табл. 2 и 3 показывает, что наличие линейной сфероконцентрической неоднородности показателя преломления линзы в области первой поверхности позволяет в 150 раз уменьшить сферическую аберрацию и более чем в 10 раз уменьшить отклонение от условия изопланазии.

СПИСОК ЛИТЕРАТУРЫ
1. Вычислительная оптика: Справочник / М. М. Русинов, А. П. Грамматин, П. Д. Иванов и др.; Под общ. ред. М. М. Русинова. Л.: Машиностроение, 1984.
2. Инструкция к пакету программ OPAL [Электронный ресурс]: .
3. Инструкция к программе Zemax [Электронный ресурс]: .
4. Buchdahl H. A. Optical Aberration Coefficients. N.Y., Douver, 1968. P. 1—18.
5. Moore D. T. Design of singlets with continuonsly varing indices of refraction // JOSA. 1971. Vol. 62, N 7. P. 886—894.
6. Sands P. J. Third-order aberrations of inhomogeneous lenses // JOSA. 1970. N 60. P. 1436.

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5

72 И. В. Гончар, А. С. Иванов, В. В. Манухов, А. Б. Федорцов

7. Сушков А. Л. Монохроматические аберрации граданов как базовых элементов жестких эндоскопов: Учеб. пособие. М.: Изд-во МГТУ им. Н. Э. Баумана, 2008.

8. Сушков А. Л. Параметры сфероконцентрического распределения показателя преломления в сферической и прямоугольной системах координат // Изв. вузов. Приборостроение. 2009. Т. 52, № 12. С. 54—60.

Александр Леонидович Сушков

Сведения об авторе — канд. техн. наук, доцент; Московский государственный технический
университет им. Н. Э. Баумана, кафедра оптико-электронных приборов научных исследований; E-mail: ale-sushkov@yandex.ru

Рекомендована кафедрой оптико-электронных приборов научных исследований

Поступила в редакцию 10.03.10 г.

ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2012. Т. 55, № 5