olejorik/Zernike_single_tests.jl Secret

## Zernike_single_tests.jl
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 72c0ea10-1d40-11eb-00e5-df29414a7ef9
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Plots")
	Pkg.add("PlutoUI")
	using Plots
	using PlutoUI
end

# ╔═╡ 84d33a00-1d40-11eb-1056-9ba7e3cc85da
md"# Многочлены Цернике

Объяснение картинки"

# ╔═╡ 21a1ced2-1e88-11eb-3a1a-b9da1e88dfaf
md"Эта ячейка устанавливает пакет для рисования графиков и для создания интерактивных виджетов. Может занять некоторое время..."

# ╔═╡ c0cb5240-1e1c-11eb-0831-f74288ca016b
md"Графики в Julia рисовать просто"

# ╔═╡ 15b2c810-1d41-11eb-015a-456a66526fdf
plot(sin,-π,π)

# ╔═╡ 78a59792-1d41-11eb-037e-1324a4432289
plot(x->(sin(2x) +3cos(13x-.1)), -π,π, title="Strange function")


# ╔═╡ ad5732f0-1e1d-11eb-25a8-0d26ca7a4ca6
md"А контурные графики или плотностные можно с помощью соответсвующих команд"

# ╔═╡ c65cfdc0-1e1d-11eb-1a4c-4b1ed0858b3e
heatmap(rand(300,200))

# ╔═╡ de002b50-1e1d-11eb-3a63-5b0bd021b901
X = -1:0.01:1; Y= -1:0.01:1

# ╔═╡ f3e55530-1e1d-11eb-0d95-6d1aa916be50
heatmap(X' .* Y)

# ╔═╡ 3a24e240-1e1e-11eb-3bb7-a135326f66d9
contour(X' .* Y, levels = -1:.25:1)

# ╔═╡ 2111ad20-1e22-11eb-3652-a98fa458015c
md"Или вот так задав функцию (выражение можно менять, график обновится)"

# ╔═╡ 40789de0-1e22-11eb-0fa2-43d0031d6a02
f(x,y) = cos(6.28x) * sin(5π*y)

# ╔═╡ 1ee4e20e-1e22-11eb-045d-37bc369eb991
heatmap(X, Y,f)

# ╔═╡ 5f094fb0-1e1e-11eb-165d-31c91a7f3b0a
md"# Теперь определим функцию для полиномов Цернике"

# ╔═╡ 89bd0350-1e1e-11eb-0cf4-f5e1f2660af8
md"Определим плотную сетку"

# ╔═╡ be2d3ce0-1e1e-11eb-2b4c-6d3adfab37c5
x = -1: 0.01 :1; y = -1: 0.01 :1

# ╔═╡ 94f7a850-1e1f-11eb-350d-29222d66ee20
md"Теперь можно построить много графиков разом или смотреть по одному (index можно менять, график обновится)"

# ╔═╡ e01e1000-1e22-11eb-35e0-d7f61d28e4b0
@bind index Slider(1:66,default=31, show_value=true)

# ╔═╡ 3fff2f80-1e1f-11eb-0234-c3024a523219
md"Можно менять разные параметры (например палитру) и добавлять контуры"

# ╔═╡ 26bae700-1e1d-11eb-1dca-19631509c896
md"### Definitions"

# ╔═╡ 3317e160-1e1d-11eb-054d-a7c53da8aac0
@doc raw"""
    zernike(x, y, maxord)

Calculate  unit-normalized Zernike polynomials and their x,y derivatives

Translation of the code from
[1] T. B. Andersen, “Efficient and robust recurrence relations for the Zernike circle polynomials
and their derivatives in Cartesian coordinates,” _Opt. Express_, vol. 26, no. 15, p. 18878, Jul. 2018.

     Numbering scheme:
     Within a radial order, sine terms come first
             ...
           sin((n-2m)*theta)   for m = 0,..., [(n+1)/2]-1
             ...
              1                for n even, m = n/2
             ...
           cos((n-2m)*theta)   for m = [n/2]+1,...,n
             ...

     INPUT:
     x, y normalized (x,y) coordinates in unit circle
     MaxOrd: Maximum Zernike radial order

     OUTPUT:
     dUdx[...]   array to receive each derivative dU/dx at (x,y)
     dUdy[...]   array to receive each derivative dU/dy at (x,y)

"""
function zernike(x, y, maxord::Int)

    ztotnumber = Int((maxord + 2) * (maxord + 1) / 2)
    # print(ztotnumber)
    zern = zeros(ztotnumber)
    dudx = zeros(ztotnumber)
    dudy = zeros(ztotnumber)

    zern[1] = 1
    dudx[1] = 0
    dudy[1] = 0

    zern[2] = y
    zern[3] = x
    dudx[2] = 0
    dudx[3] = 1
    dudy[2] = 1
    dudy[3] = 0

    kndx = 1                # index for term from 2 orders down
    jbeg = 2                # start index for current radial order
    jend = 3                # end index for current radial order
    jndx = 3                # running index for current Zern
    even = -1

    #  Outer loop in radial order index
    for nn = 2:maxord


        even   = -even          # parity of radial index
        jndx1  = jbeg           # index for 1st ascending series in x
        jndx2  = jend           # index for 1st descending series in y
        jndx11 = jndx1 - 1      # index for 2nd ascending series in x
        jndx21 = jndx2 + 1      # index for 2nd descending series in y
        jbeg   = jend + 1       # end of previous radial order +1
        nn2    = nn / 2
        nn1    = (nn - 1) / 2

      # Inner loop in azimuthal index
        for mm = 0:nn
            jndx += 1                  # increment running index for current Zern

            if mm == 0
                zern[jndx] = x * zern[jndx1] + y * zern[jndx2]
                dudx[jndx] = zern[jndx1] * nn
                dudy[jndx] = zern[jndx2] * nn

            elseif mm == nn
                zern[jndx] = x * zern[jndx11] - y * zern[jndx21]
                dudx[jndx] = zern[jndx11] * nn
                dudy[jndx] = -zern[jndx21] * nn

            elseif even > 0 && mm == nn2
                zern[jndx] = 2 * (x * zern[jndx1] + y * zern[jndx2]) - zern[kndx]
                dudx[jndx] = 2 * nn * zern[jndx1] + dudx[kndx]
                dudy[jndx] = 2 * nn * zern[jndx2] + dudy[kndx]
                kndx += 1

            elseif even < 0 && mm == nn1
                qval = zern[jndx2] - zern[jndx21]
                zern[jndx] = x * zern[jndx11] + y * qval - zern[kndx]
                dudx[jndx] = zern[jndx11] * nn + dudx[kndx]
                dudy[jndx] = qval * nn + dudy[kndx]
                kndx += 1

            elseif even < 0 && mm == nn1 + 1
                pval = zern[jndx1] + zern[jndx11]
                zern[jndx] = x * pval + y * zern[jndx2] - zern[kndx]
                dudx[jndx] = pval * nn + dudx[kndx]
                dudy[jndx] = zern[jndx2] * nn + dudy[kndx]
                kndx += 1

            else
                pval = zern[jndx1] + zern[jndx11]
                qval = zern[jndx2] - zern[jndx21]
                zern[jndx] = x * pval + y * qval - zern[kndx]
                dudx[jndx] = pval * nn + dudx[kndx]
                dudy[jndx] = qval * nn + dudy[kndx]
                kndx += 1

            end

            jndx11 = jndx1                   # update indices
            jndx1 += 1
            jndx21 = jndx2
            jndx2 += -1
        end

        jend = jndx
    end

    (z = zern, zx = dudx, zy = dudy)

end

# ╔═╡ 846ab550-1e1e-11eb-05dc-b73ecdcf87b4
 z(n)=((x,y) ->zernike(x,y,10)[:z][n] * (x^2 +y^2 > 1 ? NaN : 1))

# ╔═╡ a1f0bba0-1e1f-11eb-01b4-8142ce56f880
contour(x,y, z(index),
	c=:tab20, # цветовая палитра, см тут список названий http://docs.juliaplots.org/latest/generated/colorschemes/
	fill=true,
	levels=-1: 0.1:1, # линии уровня
	cb=false, # легенда
	aspect_ratio=:equal, axis=:none, grid=false, showaxis=false, ticks=false)

# ╔═╡ Cell order:
# ╟─84d33a00-1d40-11eb-1056-9ba7e3cc85da
# ╟─21a1ced2-1e88-11eb-3a1a-b9da1e88dfaf
# ╠═72c0ea10-1d40-11eb-00e5-df29414a7ef9
# ╟─c0cb5240-1e1c-11eb-0831-f74288ca016b
# ╠═15b2c810-1d41-11eb-015a-456a66526fdf
# ╠═78a59792-1d41-11eb-037e-1324a4432289
# ╟─ad5732f0-1e1d-11eb-25a8-0d26ca7a4ca6
# ╠═c65cfdc0-1e1d-11eb-1a4c-4b1ed0858b3e
# ╠═de002b50-1e1d-11eb-3a63-5b0bd021b901
# ╠═f3e55530-1e1d-11eb-0d95-6d1aa916be50
# ╠═3a24e240-1e1e-11eb-3bb7-a135326f66d9
# ╟─2111ad20-1e22-11eb-3652-a98fa458015c
# ╠═40789de0-1e22-11eb-0fa2-43d0031d6a02
# ╠═1ee4e20e-1e22-11eb-045d-37bc369eb991
# ╟─5f094fb0-1e1e-11eb-165d-31c91a7f3b0a
# ╠═846ab550-1e1e-11eb-05dc-b73ecdcf87b4
# ╟─89bd0350-1e1e-11eb-0cf4-f5e1f2660af8
# ╠═be2d3ce0-1e1e-11eb-2b4c-6d3adfab37c5
# ╟─94f7a850-1e1f-11eb-350d-29222d66ee20
# ╠═e01e1000-1e22-11eb-35e0-d7f61d28e4b0
# ╠═a1f0bba0-1e1f-11eb-01b4-8142ce56f880
# ╟─3fff2f80-1e1f-11eb-0234-c3024a523219
# ╟─26bae700-1e1d-11eb-1dca-19631509c896
# ╟─3317e160-1e1d-11eb-054d-a7c53da8aac0
	### A Pluto.jl notebook ###
	# v0.12.6

	using Markdown
	using InteractiveUtils

	# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
	macro bind(def, element)
	quote
	local el = $(esc(element))
	global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
	el
	end
	end

	# ╔═╡ 72c0ea10-1d40-11eb-00e5-df29414a7ef9
	begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Plots")
	Pkg.add("PlutoUI")
	using Plots
	using PlutoUI
	end

	# ╔═╡ 84d33a00-1d40-11eb-1056-9ba7e3cc85da
	md"# Многочлены Цернике

	Объяснение картинки"

	# ╔═╡ 21a1ced2-1e88-11eb-3a1a-b9da1e88dfaf
	md"Эта ячейка устанавливает пакет для рисования графиков и для создания интерактивных виджетов. Может занять некоторое время..."

	# ╔═╡ c0cb5240-1e1c-11eb-0831-f74288ca016b
	md"Графики в Julia рисовать просто"

	# ╔═╡ 15b2c810-1d41-11eb-015a-456a66526fdf
	plot(sin,-π,π)

	# ╔═╡ 78a59792-1d41-11eb-037e-1324a4432289
	plot(x->(sin(2x) +3cos(13x-.1)), -π,π, title="Strange function")


	# ╔═╡ ad5732f0-1e1d-11eb-25a8-0d26ca7a4ca6
	md"А контурные графики или плотностные можно с помощью соответсвующих команд"

	# ╔═╡ c65cfdc0-1e1d-11eb-1a4c-4b1ed0858b3e
	heatmap(rand(300,200))

	# ╔═╡ de002b50-1e1d-11eb-3a63-5b0bd021b901
	X = -1:0.01:1; Y= -1:0.01:1

	# ╔═╡ f3e55530-1e1d-11eb-0d95-6d1aa916be50
	heatmap(X' .* Y)

	# ╔═╡ 3a24e240-1e1e-11eb-3bb7-a135326f66d9
	contour(X' .* Y, levels = -1:.25:1)

	# ╔═╡ 2111ad20-1e22-11eb-3652-a98fa458015c
	md"Или вот так задав функцию (выражение можно менять, график обновится)"

	# ╔═╡ 40789de0-1e22-11eb-0fa2-43d0031d6a02
	f(x,y) = cos(6.28x) * sin(5π*y)

	# ╔═╡ 1ee4e20e-1e22-11eb-045d-37bc369eb991
	heatmap(X, Y,f)

	# ╔═╡ 5f094fb0-1e1e-11eb-165d-31c91a7f3b0a
	md"# Теперь определим функцию для полиномов Цернике"

	# ╔═╡ 89bd0350-1e1e-11eb-0cf4-f5e1f2660af8
	md"Определим плотную сетку"

	# ╔═╡ be2d3ce0-1e1e-11eb-2b4c-6d3adfab37c5
	x = -1: 0.01 :1; y = -1: 0.01 :1

	# ╔═╡ 94f7a850-1e1f-11eb-350d-29222d66ee20
	md"Теперь можно построить много графиков разом или смотреть по одному (index можно менять, график обновится)"

	# ╔═╡ e01e1000-1e22-11eb-35e0-d7f61d28e4b0
	@bind index Slider(1:66,default=31, show_value=true)

	# ╔═╡ 3fff2f80-1e1f-11eb-0234-c3024a523219
	md"Можно менять разные параметры (например палитру) и добавлять контуры"

	# ╔═╡ 26bae700-1e1d-11eb-1dca-19631509c896
	md"### Definitions"

	# ╔═╡ 3317e160-1e1d-11eb-054d-a7c53da8aac0
	@doc raw"""
	zernike(x, y, maxord)

	Calculate unit-normalized Zernike polynomials and their x,y derivatives

	Translation of the code from
	[1] T. B. Andersen, “Efficient and robust recurrence relations for the Zernike circle polynomials
	and their derivatives in Cartesian coordinates,” _Opt. Express_, vol. 26, no. 15, p. 18878, Jul. 2018.

	Numbering scheme:
	Within a radial order, sine terms come first
	...
	sin((n-2m)*theta) for m = 0,..., [(n+1)/2]-1
	...
	1 for n even, m = n/2
	...
	cos((n-2m)*theta) for m = [n/2]+1,...,n
	...

	INPUT:
	x, y normalized (x,y) coordinates in unit circle
	MaxOrd: Maximum Zernike radial order

	OUTPUT:
	dUdx[...] array to receive each derivative dU/dx at (x,y)
	dUdy[...] array to receive each derivative dU/dy at (x,y)

	"""
	function zernike(x, y, maxord::Int)

	ztotnumber = Int((maxord + 2) * (maxord + 1) / 2)
	# print(ztotnumber)
	zern = zeros(ztotnumber)
	dudx = zeros(ztotnumber)
	dudy = zeros(ztotnumber)

	zern[1] = 1
	dudx[1] = 0
	dudy[1] = 0

	zern[2] = y
	zern[3] = x
	dudx[2] = 0
	dudx[3] = 1
	dudy[2] = 1
	dudy[3] = 0

	kndx = 1 # index for term from 2 orders down
	jbeg = 2 # start index for current radial order
	jend = 3 # end index for current radial order
	jndx = 3 # running index for current Zern
	even = -1

	# Outer loop in radial order index
	for nn = 2:maxord


	even = -even # parity of radial index
	jndx1 = jbeg # index for 1st ascending series in x
	jndx2 = jend # index for 1st descending series in y
	jndx11 = jndx1 - 1 # index for 2nd ascending series in x
	jndx21 = jndx2 + 1 # index for 2nd descending series in y
	jbeg = jend + 1 # end of previous radial order +1
	nn2 = nn / 2
	nn1 = (nn - 1) / 2

	# Inner loop in azimuthal index
	for mm = 0:nn
	jndx += 1 # increment running index for current Zern

	if mm == 0
	zern[jndx] = x * zern[jndx1] + y * zern[jndx2]
	dudx[jndx] = zern[jndx1] * nn
	dudy[jndx] = zern[jndx2] * nn

	elseif mm == nn
	zern[jndx] = x * zern[jndx11] - y * zern[jndx21]
	dudx[jndx] = zern[jndx11] * nn
	dudy[jndx] = -zern[jndx21] * nn

	elseif even > 0 && mm == nn2
	zern[jndx] = 2 * (x * zern[jndx1] + y * zern[jndx2]) - zern[kndx]
	dudx[jndx] = 2 * nn * zern[jndx1] + dudx[kndx]
	dudy[jndx] = 2 * nn * zern[jndx2] + dudy[kndx]
	kndx += 1

	elseif even < 0 && mm == nn1
	qval = zern[jndx2] - zern[jndx21]
	zern[jndx] = x * zern[jndx11] + y * qval - zern[kndx]
	dudx[jndx] = zern[jndx11] * nn + dudx[kndx]
	dudy[jndx] = qval * nn + dudy[kndx]
	kndx += 1

	elseif even < 0 && mm == nn1 + 1
	pval = zern[jndx1] + zern[jndx11]
	zern[jndx] = x * pval + y * zern[jndx2] - zern[kndx]
	dudx[jndx] = pval * nn + dudx[kndx]
	dudy[jndx] = zern[jndx2] * nn + dudy[kndx]
	kndx += 1

	else
	pval = zern[jndx1] + zern[jndx11]
	qval = zern[jndx2] - zern[jndx21]
	zern[jndx] = x * pval + y * qval - zern[kndx]
	dudx[jndx] = pval * nn + dudx[kndx]
	dudy[jndx] = qval * nn + dudy[kndx]
	kndx += 1

	end

	jndx11 = jndx1 # update indices
	jndx1 += 1
	jndx21 = jndx2
	jndx2 += -1
	end

	jend = jndx
	end

	(z = zern, zx = dudx, zy = dudy)

	end

	# ╔═╡ 846ab550-1e1e-11eb-05dc-b73ecdcf87b4
	z(n)=((x,y) ->zernike(x,y,10)[:z][n] * (x^2 +y^2 > 1 ? NaN : 1))

	# ╔═╡ a1f0bba0-1e1f-11eb-01b4-8142ce56f880
	contour(x,y, z(index),
	c=:tab20, # цветовая палитра, см тут список названий http://docs.juliaplots.org/latest/generated/colorschemes/
	fill=true,
	levels=-1: 0.1:1, # линии уровня
	cb=false, # легенда
	aspect_ratio=:equal, axis=:none, grid=false, showaxis=false, ticks=false)

	# ╔═╡ Cell order:
	# ╟─84d33a00-1d40-11eb-1056-9ba7e3cc85da
	# ╟─21a1ced2-1e88-11eb-3a1a-b9da1e88dfaf
	# ╠═72c0ea10-1d40-11eb-00e5-df29414a7ef9
	# ╟─c0cb5240-1e1c-11eb-0831-f74288ca016b
	# ╠═15b2c810-1d41-11eb-015a-456a66526fdf
	# ╠═78a59792-1d41-11eb-037e-1324a4432289
	# ╟─ad5732f0-1e1d-11eb-25a8-0d26ca7a4ca6
	# ╠═c65cfdc0-1e1d-11eb-1a4c-4b1ed0858b3e
	# ╠═de002b50-1e1d-11eb-3a63-5b0bd021b901
	# ╠═f3e55530-1e1d-11eb-0d95-6d1aa916be50
	# ╠═3a24e240-1e1e-11eb-3bb7-a135326f66d9
	# ╟─2111ad20-1e22-11eb-3652-a98fa458015c
	# ╠═40789de0-1e22-11eb-0fa2-43d0031d6a02
	# ╠═1ee4e20e-1e22-11eb-045d-37bc369eb991
	# ╟─5f094fb0-1e1e-11eb-165d-31c91a7f3b0a
	# ╠═846ab550-1e1e-11eb-05dc-b73ecdcf87b4
	# ╟─89bd0350-1e1e-11eb-0cf4-f5e1f2660af8
	# ╠═be2d3ce0-1e1e-11eb-2b4c-6d3adfab37c5
	# ╟─94f7a850-1e1f-11eb-350d-29222d66ee20
	# ╠═e01e1000-1e22-11eb-35e0-d7f61d28e4b0
	# ╠═a1f0bba0-1e1f-11eb-01b4-8142ce56f880
	# ╟─3fff2f80-1e1f-11eb-0234-c3024a523219
	# ╟─26bae700-1e1d-11eb-1dca-19631509c896
	# ╟─3317e160-1e1d-11eb-054d-a7c53da8aac0