diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5fcfb04 --- /dev/null +++ b/.gitignore @@ -0,0 +1,29 @@ +# Files generated by invoking Julia with --code-coverage +*.jl.cov +*.jl.*.cov + +# Files generated by invoking Julia with --track-allocation +*.jl.mem + +# System-specific files and directories generated by the BinaryProvider and BinDeps packages +# They contain absolute paths specific to the host computer, and so should not be committed +deps/deps.jl +deps/build.log +deps/downloads/ +deps/usr/ +deps/src/ + +# Build artifacts for creating documentation generated by the Documenter package +docs/build/ +docs/site/ + +# File generated by Pkg, the package manager, based on a corresponding Project.toml +# It records a fixed state of all packages used by the project. As such, it should not be +# committed for packages, but should be committed for applications that require a static +# environment. +Manifest.toml + +*.swp + +*.vscode + diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..69552e5 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 roflmaostc and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..a348589 --- /dev/null +++ b/Project.toml @@ -0,0 +1,37 @@ +name = "WaveOpticsPropagation" +uuid = "c4c7a1f9-3adc-4a73-843a-d378b6c86436" +authors = ["Felix Wechsler (roflmaostc) "] +version = "0.1.0" + +[deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e" +IndexFunArrays = "613c443e-d742-454e-bfc6-1d7f8dd76566" +NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[compat] +CUDA = "5" +ChainRulesCore = "1" +EllipsisNotation = "1" +FFTW = "1" +FourierTools = "0.4" +IndexFunArrays = "0.2" +NDTools = "0.5" +Zygote = "0.6.60" +julia = "1.9" + +[extras] +ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" +FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +IndexFunArrays = "613c443e-d742-454e-bfc6-1d7f8dd76566" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[targets] +test = ["Test", "ChainRulesTestUtils", "Zygote", "IndexFunArrays", "Markdown", "InteractiveUtils", "FiniteDifferences"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..cceebf6 --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ + + +# WaveOpticsPropagation.jl + +| **Build Status** | **Code Coverage** | **Docs stable** | **Docs main** | +|:-----------------------------------------:|:-------------------------------:|:-----------------------:|:--------------------:| +| [![][CI-img]][CI-url] | [![][codecov-img]][codecov-url] | [![Documentation for stable version](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaPhysics.github.io/WaveOpticsPropagation.jl/stable) | [![Documentation for development version](https://img.shields.io/badge/docs-main-blue.svg)](https://JuliaPhysics.github.io/WaveOpticsPropagation.jl/dev) + + + +Propagate waves efficiently, optically, physically, differentiably with [Julia Lang](https://julialang.org/). +Those functions are fast and memory efficient implemented and hence are suited to be used in inverse problems. + +⚠️ Under development. Expect things to break. But feel free to try the examples, they should always work! + +## Installation +Not registered yet, hence install with: +```julia +julia> ]add https://github.com/JuliaPhysics/WaveOpticsPropagation.jl +``` + +## Examples +See [the rendered notebooks here](https://juliaphysics.github.io/WaveOpticsPropagation.jl/). +Otherwise, just look into the [examples folder](examples/). + +## Features +### Implemented +* Propagate (electrical) fields based on wave propagation +* Propagations + * [x] Angular Spectrum Method of Plane Waves (AS) + * [x] Fraunhofer Diffraction + * [x] [Scalable Angular Spectrum propagation](https://opg.optica.org/optica/viewmedia.cfm?uri=optica-10-11-1407&html=true) + * [x] Shifted Angular Spectrum propagation + * [ ] Fresnel Propagation with Scaling Behaviour (no priority yet, PR are welcome for that. In principle very similar to the other methods.) +* [x] CUDA support +* [x] Differentiable (mainly based on Zygote.jl and ChainRulesCore.jl) + +### Planned +Vectorial propagation in free space is just a propagation of each of the components. Right now, this is not a priority and is not implemented yet. +But of course, each vectorial component can be propagated separately. + +## Development +Contributions are very welcome! +File an [issue](https://github.com/JuliaPhysics/WaveOpticsPropagation.jl/issues) on [GitHub](https://github.com/JuliaPhysics/WaveOpticsPropagation.jl) if you encounter any problems. +Also file an issue if you want to discuss or propose features. + +## Related packages +There is the outdated [PhysicalOptics.jl](https://github.com/JuliaPhysics/PhysicalOptics.jl) which provided similar methods. +For geometrical ray tracing use [OpticSim.jl](https://github.com/brianguenter/OpticSim.jl). + +[CI-img]: https://github.com/JuliaPhysics/WaveOpticsPropagation.jl/actions/workflows/CI.yml/badge.svg +[CI-url]: https://github.com/JuliaPhysics/WaveOpticsPropagation.jl/actions/workflows/CI.yml + +[codecov-img]: https://codecov.io/gh/JuliaPhysics/WaveOpticsPropagation.jl/branch/main/graph/badge.svg?token=6XWI1M1MPB +[codecov-url]: https://codecov.io/gh/JuliaPhysics/WaveOpticsPropagation.jl diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..0fca31b --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +WaveOpticsPropagation = "c4c7a1f9-3adc-4a73-843a-d378b6c86436" diff --git a/docs/doc_server.sh b/docs/doc_server.sh new file mode 100755 index 0000000..7620444 --- /dev/null +++ b/docs/doc_server.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +python3 -m http.server --bind localhost diff --git a/docs/logo/logo-pluto.html b/docs/logo/logo-pluto.html new file mode 100644 index 0000000..f976d82 --- /dev/null +++ b/docs/logo/logo-pluto.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/docs/logo/logo-pluto.jl b/docs/logo/logo-pluto.jl new file mode 100644 index 0000000..717dedd --- /dev/null +++ b/docs/logo/logo-pluto.jl @@ -0,0 +1,665 @@ +### A Pluto.jl notebook ### +# v0.19.36 + +using Markdown +using InteractiveUtils + +# ╔═╡ 9d8bbf18-87b3-4ab3-8179-5019afe4fb53 +using Luxor + +# ╔═╡ e30e05ae-4ffa-4345-af5f-d835e4debdbf +cyl(r, ϕ) = Point(r * cos(deg2rad(ϕ)), r * sin(deg2rad(ϕ))) + +# ╔═╡ cc91bb32-9995-48ad-82dd-5fa6ec232cb3 +function ball(pos, α1, α2, radii=30:10:80) + circle(pos, 20, action=:fill) + for r in radii + arc(pos, r, deg2rad(α1), deg2rad(α2)) + strokepath() + end +end + +# ╔═╡ e0943f38-449d-4623-862a-f9c0e4f00f1d +@svg begin + #background("white") + sethue("brown3") + a = 200 + radii = 28:8:120 + Δt = 30 + ball(Point(0,0), -60 - Δt, Δt, radii) + sethue("mediumorchid3") + ball(Point(a,0), 180 - Δt, 240 + Δt, radii) + sethue("forestgreen") + ball(Point(a / 2, -a * √3 / 2), 60 - Δt, 180 - 60 + Δt, radii) +end 500 500 "logo.svg" + +# ╔═╡ 00000000-0000-0000-0000-000000000001 +PLUTO_PROJECT_TOML_CONTENTS = """ +[deps] +Luxor = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc" + +[compat] +Luxor = "~3.8.0" +""" + +# ╔═╡ 00000000-0000-0000-0000-000000000002 +PLUTO_MANIFEST_TOML_CONTENTS = """ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.10.0" +manifest_format = "2.0" +project_hash = "fcac3d4ae530321b609832702d21ffd4863fff29" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[deps.Cairo]] +deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"] +git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b" +uuid = "159f3aea-2a34-519c-b102-8c37f9878175" +version = "1.0.5" + +[[deps.Cairo_jll]] +deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2" +uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" +version = "1.16.1+1" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.4" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.10" + +[[deps.Compat]] +deps = ["UUIDs"] +git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.10.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.5+1" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "ac67408d9ddf207de5cfa9a97e114352430f01ed" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.16" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.5.0+0" + +[[deps.FFMPEG]] +deps = ["FFMPEG_jll"] +git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8" +uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" +version = "0.4.1" + +[[deps.FFMPEG_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] +git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e" +uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" +version = "4.4.4+1" + +[[deps.FileIO]] +deps = ["Pkg", "Requires", "UUIDs"] +git-tree-sha1 = "c5c28c245101bd59154f649e19b038d15901b5dc" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.16.2" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + +[[deps.Fontconfig_jll]] +deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03" +uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" +version = "2.13.93+0" + +[[deps.FreeType2_jll]] +deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"] +git-tree-sha1 = "d8db6a5a2fe1381c1ea4ef2cab7c69c2de7f9ea0" +uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" +version = "2.13.1+0" + +[[deps.FriBidi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91" +uuid = "559328eb-81f9-559d-9380-de523a88c83c" +version = "1.0.10+0" + +[[deps.Gettext_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] +git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" +uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" +version = "0.21.0+0" + +[[deps.Glib_jll]] +deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"] +git-tree-sha1 = "e94c92c7bf4819685eb80186d51c43e71d4afa17" +uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" +version = "2.76.5+0" + +[[deps.Graphics]] +deps = ["Colors", "LinearAlgebra", "NaNMath"] +git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd" +uuid = "a2bd30eb-e257-5431-a919-1863eab51364" +version = "1.1.2" + +[[deps.Graphite2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" +uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" +version = "1.3.14+0" + +[[deps.HarfBuzz_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] +git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3" +uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" +version = "2.8.1+1" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.5.0" + +[[deps.JpegTurbo_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "60b1194df0a3298f460063de985eae7b01bc011a" +uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" +version = "3.0.1+0" + +[[deps.Juno]] +deps = ["Base64", "Logging", "Media", "Profile"] +git-tree-sha1 = "07cb43290a840908a771552911a6274bc6c072c7" +uuid = "e5e0dc1b-0480-54bc-9374-aad01c23163d" +version = "0.8.4" + +[[deps.LAME_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" +uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" +version = "3.100.1+0" + +[[deps.LERC_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434" +uuid = "88015f11-f218-50d7-93a8-a6af411a945d" +version = "3.0.0+1" + +[[deps.LLVMOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "d986ce2d884d49126836ea94ed5bfb0f12679713" +uuid = "1d63c593-3942-5779-bab2-d838dc0a180e" +version = "15.0.7+0" + +[[deps.LZO_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" +uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" +version = "2.10.1+0" + +[[deps.LaTeXStrings]] +git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.1" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.4" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "8.4.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.11.0+1" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.Libffi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290" +uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" +version = "3.2.2+1" + +[[deps.Libgcrypt_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"] +git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae" +uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" +version = "1.8.7+0" + +[[deps.Libgpg_error_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9" +uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" +version = "1.42.0+0" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.17.0+0" + +[[deps.Libmount_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73" +uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" +version = "2.35.0+0" + +[[deps.Librsvg_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pango_jll", "Pkg", "gdk_pixbuf_jll"] +git-tree-sha1 = "ae0923dab7324e6bc980834f709c4cd83dd797ed" +uuid = "925c91fb-5dd6-59dd-8e8c-345e74382d89" +version = "2.54.5+0" + +[[deps.Libtiff_jll]] +deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "3eb79b0ca5764d4799c06699573fd8f533259713" +uuid = "89763e89-9b03-5906-acba-b20f662cd828" +version = "4.4.0+0" + +[[deps.Libuuid_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066" +uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" +version = "2.36.0+0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.Luxor]] +deps = ["Base64", "Cairo", "Colors", "DataStructures", "Dates", "FFMPEG", "FileIO", "Juno", "LaTeXStrings", "PrecompileTools", "Random", "Requires", "Rsvg"] +git-tree-sha1 = "aa3eb624552373a6204c19b00e95ce62ea932d32" +uuid = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc" +version = "3.8.0" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "b211c553c199c111d998ecdaf7623d1b89b69f93" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.12" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.2+1" + +[[deps.Media]] +deps = ["MacroTools", "Test"] +git-tree-sha1 = "75a54abd10709c01f1b86b84ec225d26e840ed58" +uuid = "e89f7d12-3494-54d1-8411-f7d8b9ae1f27" +version = "0.5.0" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2023.1.10" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.0.2" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.Ogg_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" +uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" +version = "1.3.5+1" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.23+2" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+2" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "3.0.12+0" + +[[deps.Opus_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720" +uuid = "91d4177d-7536-5919-b921-800302f37372" +version = "1.3.2+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.6.3" + +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.42.0+1" + +[[deps.Pango_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "4745216e94f71cb768d58330b059c9b76f32cb66" +uuid = "36c8627f-9965-5494-a995-c6b170f724f3" +version = "1.50.14+0" + +[[deps.Pixman_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] +git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87" +uuid = "30392449-352a-5448-841d-b1acce4e97dc" +version = "0.42.2+0" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.10.0" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.0" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.4.1" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.Rsvg]] +deps = ["Cairo", "Glib_jll", "Librsvg_jll"] +git-tree-sha1 = "3d3dc66eb46568fb3a5259034bfc752a0eb0c686" +uuid = "c4c386cf-5103-5370-be45-f3a111cca3b8" +version = "1.0.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.10.0" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.2.1+1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] +git-tree-sha1 = "801cbe47eae69adc50f36c3caec4758d2650741b" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.12.2+0" + +[[deps.XSLT_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] +git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" +uuid = "aed1982a-8fda-507f-9586-7b0439959a61" +version = "1.1.34+0" + +[[deps.Xorg_libX11_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] +git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495" +uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc" +version = "1.8.6+0" + +[[deps.Xorg_libXau_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8" +uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" +version = "1.0.11+0" + +[[deps.Xorg_libXdmcp_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "34d526d318358a859d7de23da945578e8e8727b7" +uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05" +version = "1.1.4+0" + +[[deps.Xorg_libXext_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3" +uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3" +version = "1.3.4+4" + +[[deps.Xorg_libXrender_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96" +uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa" +version = "0.9.10+4" + +[[deps.Xorg_libpthread_stubs_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "8fdda4c692503d44d04a0603d9ac0982054635f9" +uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74" +version = "0.1.1+0" + +[[deps.Xorg_libxcb_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"] +git-tree-sha1 = "b4bfde5d5b652e22b9c790ad00af08b6d042b97d" +uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b" +version = "1.15.0+0" + +[[deps.Xorg_xtrans_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "e92a1a012a10506618f10b7047e478403a046c77" +uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" +version = "1.5.0+0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+1" + +[[deps.Zstd_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "49ce682769cd5de6c72dcf1b94ed7790cd08974c" +uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" +version = "1.5.5+0" + +[[deps.gdk_pixbuf_jll]] +deps = ["Artifacts", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Xorg_libX11_jll", "libpng_jll"] +git-tree-sha1 = "e9190f9fb03f9c3b15b9fb0c380b0d57a3c8ea39" +uuid = "da03df04-f53b-5353-a52f-6a8b0620ced0" +version = "2.42.8+0" + +[[deps.libaom_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4" +uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" +version = "3.4.0+0" + +[[deps.libass_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" +uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" +version = "0.15.1+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.8.0+1" + +[[deps.libfdk_aac_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" +uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" +version = "2.0.2+0" + +[[deps.libpng_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] +git-tree-sha1 = "93284c28274d9e75218a416c65ec49d0e0fcdf3d" +uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" +version = "1.6.40+0" + +[[deps.libvorbis_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] +git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c" +uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" +version = "1.3.7+1" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.52.0+1" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+2" + +[[deps.x264_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" +uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" +version = "2021.5.5+0" + +[[deps.x265_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" +uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" +version = "3.5.0+0" +""" + +# ╔═╡ Cell order: +# ╠═9d8bbf18-87b3-4ab3-8179-5019afe4fb53 +# ╠═e30e05ae-4ffa-4345-af5f-d835e4debdbf +# ╠═cc91bb32-9995-48ad-82dd-5fa6ec232cb3 +# ╠═e0943f38-449d-4623-862a-f9c0e4f00f1d +# ╟─00000000-0000-0000-0000-000000000001 +# ╟─00000000-0000-0000-0000-000000000002 diff --git a/docs/logo/logo-pluto.plutostate b/docs/logo/logo-pluto.plutostate new file mode 100644 index 0000000..93f1229 Binary files /dev/null and b/docs/logo/logo-pluto.plutostate differ diff --git a/docs/logo/logo.png b/docs/logo/logo.png new file mode 100644 index 0000000..ec152e0 Binary files /dev/null and b/docs/logo/logo.png differ diff --git a/docs/logo/logo.svg b/docs/logo/logo.svg new file mode 100644 index 0000000..4fc8a15 --- /dev/null +++ b/docs/logo/logo.svg @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..d36c92b --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,13 @@ +using WaveOpticsPropagation, Documenter + +DocMeta.setdocmeta!(WaveOpticsPropagation, :DocTestSetup, :(using WaveOpticsPropagation); recursive=true) +makedocs(modules = [WaveOpticsPropagation], + sitename = "WaveOpticsPropagation.jl", + pages = Any[ + "WaveOpticsPropagation.jl" => "index.md", + "Function Docstrings" => "functions.md" + ], + warnonly=true, + ) + +deploydocs(repo = "github.com/JuliaPhysics/WaveOpticsPropagation.jl.git", devbranch="main") diff --git a/docs/src/functions.md b/docs/src/functions.md new file mode 100644 index 0000000..7d6e53c --- /dev/null +++ b/docs/src/functions.md @@ -0,0 +1,15 @@ +# Efficient Propagations +They create efficient functions which avoid recalculating some parts. +```@docs +AngularSpectrum +Fraunhofer +ScalableAngularSpectrum +ShiftedAngularSpectrum +``` + +# Propagation +Those were merely for testing purposes. +```@docs +fraunhofer +angular_spectrum +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..347bad3 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,46 @@ +```@raw html + +``` ⠀ + +# WaveOpticsPropagation.jl + + +Propagate waves efficiently, optically, physically, differentiably with [Julia Lang](https://julialang.org/). +Those functions are fast and memory efficient implemented and hence are suited to be used in inverse problems. + +⚠️ Under development. Expect things to break. But feel free to try the examples, they should always work! + +## Installation +Not registered yet, hence install with: +```julia +julia> ]add https://github.com/JuliaPhysics/WaveOpticsPropagation.jl +``` + +## Examples +See [the rendered notebooks here](https://juliaphysics.github.io/WaveOpticsPropagation.jl/). +Otherwise, just look into the [examples folder](examples/). + +## Features +### Implemented +* Propagate (electrical) fields based on wave propagation +* Propagations + * [x] Angular Spectrum Method of Plane Waves (AS) + * [x] Fraunhofer Diffraction + * [x] [Scalable Angular Spectrum propagation](https://opg.optica.org/optica/viewmedia.cfm?uri=optica-10-11-1407&html=true) + * [ ] Fresnel Propagation with Scaling Behaviour (no priority yet, PR are welcome for that. In principle very similar to the other methods.) +* [x] CUDA support +* [x] Differentiable (mainly based on Zygote.jl and ChainRulesCore.jl) + +### Planned +Vectorial propagation in free space is just a propagation of each of the components. Right now, this is not a priority and is not implemented yet. +But of course, each vectorial component can be propagated separately. + +## Development +Contributions are very welcome! +File an [issue](https://github.com/JuliaPhysics/WaveOpticsPropagation.jl/issues) on [GitHub](https://github.com/JuliaPhysics/WaveOpticsPropagation.jl) if you encounter any problems. +Also file an issue if you want to discuss or propose features. + +## Related packages +There is the outdated [PhysicalOptics.jl](https://github.com/JuliaPhysics/PhysicalOptics.jl) which provided similar methods. +For geometrical ray tracing use [OpticSim.jl](https://github.com/brianguenter/OpticSim.jl). + diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 0000000..12b1cf4 --- /dev/null +++ b/examples/Project.toml @@ -0,0 +1,22 @@ +[deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e" +ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" +ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31" +IndexFunArrays = "613c443e-d742-454e-bfc6-1d7f8dd76566" +NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d" +Napari = "ed46c112-542a-4e7e-ab0f-e39321fab6f0" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" +PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" +WaveOpticsPropagation = "c4c7a1f9-3adc-4a73-843a-d378b6c86436" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/examples/Scalable_Angular_Spectrum.html b/examples/Scalable_Angular_Spectrum.html new file mode 100644 index 0000000..1e5a3a8 --- /dev/null +++ b/examples/Scalable_Angular_Spectrum.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/Scalable_Angular_Spectrum.jl b/examples/Scalable_Angular_Spectrum.jl new file mode 100644 index 0000000..5556d09 --- /dev/null +++ b/examples/Scalable_Angular_Spectrum.jl @@ -0,0 +1,176 @@ +### A Pluto.jl notebook ### +# v0.19.36 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 1ea0729c-ad94-11ee-3b52-bb21d59e9509 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ aafab549-c977-48f4-9dba-9d693e439b44 +using WaveOpticsPropagation, CUDA, IndexFunArrays, NDTools, ImageShow, FileIO, Zygote, Optim, Plots, PlutoUI, FourierTools, NDTools + +# ╔═╡ 4aabd650-792c-407c-9590-e4960ac776d6 +md"# Setup Environment and packages" + +# ╔═╡ 3ed16b3f-08a8-4fc3-97cf-a884c452dff0 +TableOfContents() + +# ╔═╡ 20734e1c-faf7-4117-8652-528eca92c6a2 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ 7028f15f-7ec7-48b3-9b4d-66bd0c34bae6 +λ = 500f-9 + +# ╔═╡ f9869ff6-5b21-40a3-98cf-33e7192dfa8e +L = 128f-6 / 2 + +# ╔═╡ 1666eaa6-7857-415e-a43e-00849571db2a +N = 512 + +# ╔═╡ d418e87f-675f-4bc2-8eaf-70cf669aeb8d +y = fftpos(L, N, NDTools.CenterFT) + +# ╔═╡ 945ac45e-2773-4d34-878e-4f66a253afa1 +D_circ = N / 8 + +# ╔═╡ 4f7aca22-00d9-406c-90e1-f73394a718ee +U_circ = togoc(ComplexF32.(rr((N, N)) .< D_circ / 2) .* exp.(1im .* 2f0 * π ./ λ .* y .* sind(45f0)) .+ ComplexF32.(rr((N, N)) .< D_circ / 2) .* exp.(1im .* 2f0 * π ./ λ .* y' .* sind(-45f0))); + +# ╔═╡ 34c326ee-6663-432e-b884-a5419ce64827 +@bind M Slider(1:0.1:20, show_value=true, default=4) + +# ╔═╡ 4202033d-62e6-451f-a064-e61d420da5ff +z_circ = Float32(M / N / λ * L^2 * 2) + +# ╔═╡ 4dc0d972-88b7-4cb4-9ef3-925b154e45c0 +md"# First Example: Circular + +In the first example, one straight beam and one oblique beam are passing through a round aperture. + +The Fresnel number is $(round((D_circ / 2 * L / N)^2 / z_circ / λ, digits=3)) +" + +# ╔═╡ e5c50250-e3c0-4561-842a-d6c671999741 +simshow(Array(U_circ)) + +# ╔═╡ 484704ce-4497-4567-b387-aad06dee6a62 +@mytime SAS, _ = ScalableAngularSpectrum(U_circ, z_circ, λ, L) + +# ╔═╡ 9e3f7e53-3fdb-4638-89bc-27bced937193 +@time SAS_cpu, _ = ScalableAngularSpectrum(Array(U_circ), z_circ, λ, L) + +# ╔═╡ 9ec9d456-a33b-4605-9025-c82676eca7e2 +U_circ_cpu = Array(U_circ); + +# ╔═╡ cf8d99c6-120d-45a6-9fcc-9d1cb6e40a9b +@time SAS_cpu(U_circ_cpu); + +# ╔═╡ 342b000d-8865-4cd2-bb6b-2de0f6376e9c +@mytime U_p, t = SAS(U_circ) + +# ╔═╡ fc7dd2c9-7c7b-45a7-9fed-a319a62375cc +simshow(Array(abs2.(U_p)), γ=0.13, cmap=:inferno) + +# ╔═╡ 505bd6fa-2717-4297-b8b5-bf9452300655 +sum(abs2, U_p) + +# ╔═╡ bb1c90d3-62c0-4941-b8ac-00761856a3cd +sum(abs2, U_circ) + +# ╔═╡ 70afe821-c31c-4626-8b80-20695a46544f +L_box = 128f-6; + +# ╔═╡ cd2bdc4a-7abe-4184-a9c3-83b9820cf46d +N_box = 512; + +# ╔═╡ e227181f-017a-4601-914a-f63583755105 +y_box = fftpos(L_box, N_box, NDTools.CenterFT); + +# ╔═╡ f7926732-1e17-428b-9660-2bf3653a3f4a +D_box = L_box / 16f0 + +# ╔═╡ 4a30f342-cad6-46d5-a787-2144909c8765 +x_box = y_box'; + +# ╔═╡ f0cf4b3d-13fd-48d0-9890-1a8e21686a7e +U_box = togoc((x_box.^2 .<= (D_box / 2).^2) .* (y_box.^2 .<= (D_box / 2).^2) .* (exp.(1im .* 2f0 * π ./ λ .* y_box' .* sind(20f0)))); + +# ╔═╡ 11947eea-b2eb-4168-bc96-8e1fc7b74599 +@bind M_box Slider(1:0.5:20, show_value=true, default=8) + +# ╔═╡ aacc5f18-05b5-4a5e-9217-fa62966830e5 +z_box = M_box / N_box / λ * L_box^2 * 2 + +# ╔═╡ d325fa62-7895-44a9-87de-718e9d61f9bc +md"# Second Example: Quadratic + + +The Fresnel number is $(round((D_box)^2 / z_box / λ, digits=3)) +" + +# ╔═╡ 38bc91ff-3188-4129-aa6b-2af458cf59b1 +@mytime SAS2, _ = ScalableAngularSpectrum(U_box, z_box, λ, L_box, skip_final_phase=true) + +# ╔═╡ 2b11c87d-4746-4be6-81f6-f004d30beae4 +@mytime U_box_p, t_box = SAS2(U_box) + +# ╔═╡ 9e762dfc-4873-4fe3-b785-f6c4808417a6 +simshow(Array(abs2.(U_box_p)), γ=0.13, cmap=:inferno) + +# ╔═╡ Cell order: +# ╟─4aabd650-792c-407c-9590-e4960ac776d6 +# ╠═1ea0729c-ad94-11ee-3b52-bb21d59e9509 +# ╠═3ed16b3f-08a8-4fc3-97cf-a884c452dff0 +# ╠═aafab549-c977-48f4-9dba-9d693e439b44 +# ╠═20734e1c-faf7-4117-8652-528eca92c6a2 +# ╠═4dc0d972-88b7-4cb4-9ef3-925b154e45c0 +# ╠═7028f15f-7ec7-48b3-9b4d-66bd0c34bae6 +# ╠═f9869ff6-5b21-40a3-98cf-33e7192dfa8e +# ╠═1666eaa6-7857-415e-a43e-00849571db2a +# ╠═d418e87f-675f-4bc2-8eaf-70cf669aeb8d +# ╠═945ac45e-2773-4d34-878e-4f66a253afa1 +# ╠═4f7aca22-00d9-406c-90e1-f73394a718ee +# ╠═34c326ee-6663-432e-b884-a5419ce64827 +# ╠═4202033d-62e6-451f-a064-e61d420da5ff +# ╠═e5c50250-e3c0-4561-842a-d6c671999741 +# ╠═484704ce-4497-4567-b387-aad06dee6a62 +# ╠═9e3f7e53-3fdb-4638-89bc-27bced937193 +# ╠═9ec9d456-a33b-4605-9025-c82676eca7e2 +# ╠═cf8d99c6-120d-45a6-9fcc-9d1cb6e40a9b +# ╠═342b000d-8865-4cd2-bb6b-2de0f6376e9c +# ╠═fc7dd2c9-7c7b-45a7-9fed-a319a62375cc +# ╠═505bd6fa-2717-4297-b8b5-bf9452300655 +# ╠═bb1c90d3-62c0-4941-b8ac-00761856a3cd +# ╟─d325fa62-7895-44a9-87de-718e9d61f9bc +# ╠═70afe821-c31c-4626-8b80-20695a46544f +# ╠═cd2bdc4a-7abe-4184-a9c3-83b9820cf46d +# ╠═e227181f-017a-4601-914a-f63583755105 +# ╠═f7926732-1e17-428b-9660-2bf3653a3f4a +# ╠═4a30f342-cad6-46d5-a787-2144909c8765 +# ╠═f0cf4b3d-13fd-48d0-9890-1a8e21686a7e +# ╠═11947eea-b2eb-4168-bc96-8e1fc7b74599 +# ╠═aacc5f18-05b5-4a5e-9217-fa62966830e5 +# ╠═38bc91ff-3188-4129-aa6b-2af458cf59b1 +# ╠═2b11c87d-4746-4be6-81f6-f004d30beae4 +# ╠═9e762dfc-4873-4fe3-b785-f6c4808417a6 diff --git a/examples/Scalable_Angular_Spectrum.plutostate b/examples/Scalable_Angular_Spectrum.plutostate new file mode 100644 index 0000000..906bc53 Binary files /dev/null and b/examples/Scalable_Angular_Spectrum.plutostate differ diff --git a/examples/benchmark_CPU_CUDA.html b/examples/benchmark_CPU_CUDA.html new file mode 100644 index 0000000..90c14be --- /dev/null +++ b/examples/benchmark_CPU_CUDA.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/benchmark_CPU_CUDA.jl b/examples/benchmark_CPU_CUDA.jl new file mode 100644 index 0000000..227e8ee --- /dev/null +++ b/examples/benchmark_CPU_CUDA.jl @@ -0,0 +1,116 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +using Markdown +using InteractiveUtils + +# ╔═╡ 2fd337ac-8ad0-11ee-3739-459b5825a8c5 +begin + using Pkg + Pkg.activate(".") + using Revise +end + +# ╔═╡ 958b4e13-bb6c-4d4d-83f9-a922bbbfb842 +using WaveOpticsPropagation, FFTW, CUDA, Zygote + +# ╔═╡ 62c0385d-e712-43eb-adae-6731159f9f92 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ 8ee34fde-40a7-4cc1-a8ea-056a457901b0 +md"# Compared with FFT" + +# ╔═╡ a43a5d7e-5052-4ee1-a112-c8881e90b6a6 +sz = (2048, 2048) + +# ╔═╡ c33ed1de-feb3-4bab-80c5-cd422be96bb6 +array = randn(ComplexF32, sz); + +# ╔═╡ 5d7cf31a-7e08-45ea-80e4-a0e27a022258 +array_c = togoc(array); + +# ╔═╡ 4e651fc1-b6db-4c44-b830-f65c6f68f4a9 +p_cpu = plan_fft!(array, flags=FFTW.ESTIMATE); + +# ╔═╡ 5f5113ed-d335-4038-a230-b739242aafbc +p_cuda = plan_fft!(array_c); + +# ╔═╡ f897cd86-ddc1-4b84-a11c-d99161c1df3d +@mytime p_cpu * array; + +# ╔═╡ fcf4d2ab-512b-4a2c-a131-9d2519da1747 +@mytime p_cpu * array; + +# ╔═╡ f7633864-230a-4a06-a48e-8acf4adb58f1 +@mytime CUDA.@sync p_cuda * array_c; + +# ╔═╡ 0a7b6bc6-9a23-484b-ac26-6effb901ea3b +@mytime CUDA.@sync p_cuda * array_c; + +# ╔═╡ 181eb167-f944-480d-907a-e54a2ba113df +md"# Propagation" + +# ╔═╡ 864c37c5-6907-4b31-8a4e-8c3c9a684606 +AS = AngularSpectrum(array, 1f0, 1f0, 1f0)[1]; + +# ╔═╡ d9339c77-35cf-4575-b935-79d580badba5 +AS_c = AngularSpectrum(array_c, 1f0, 1f0, 1f0)[1]; + +# ╔═╡ ad07d942-4985-40dc-9371-d6579e5b92a3 +@mytime angular_spectrum(array, 1f0, 1f0, 1f0); + +# ╔═╡ 127e50d1-4538-4010-b7c3-716bffd804ed +@mytime AS(array); + +# ╔═╡ b0d00f6f-2a6c-42c0-b878-1791b2073353 +@mytime CUDA.@sync angular_spectrum(array_c, 1f0, 1f0, 1f0); + +# ╔═╡ 20935643-b699-46c0-8a44-d948551d9c44 +@mytime CUDA.@sync AS_c(array_c); + +# ╔═╡ a0637882-140c-466f-ab7d-af9cdb979d1f +md"# Gradient" + +# ╔═╡ ca2386b8-0882-4b31-90fa-7ab1df839457 +f(x) = sum(abs2, AS_c(x)[1]) + +# ╔═╡ 17679b22-61ac-42c2-959c-5e28fe880d15 +@mytime CUDA.@sync f(array_c) + +# ╔═╡ 02c3e595-b978-4ad4-961d-5d7398d3d320 +@mytime CUDA.@sync Zygote.gradient(f, array_c) + +# ╔═╡ 819f3128-3a42-425b-871d-a4913b17e94c +@mytime CUDA.@sync Zygote.withgradient(f, array_c) + +# ╔═╡ Cell order: +# ╠═2fd337ac-8ad0-11ee-3739-459b5825a8c5 +# ╠═958b4e13-bb6c-4d4d-83f9-a922bbbfb842 +# ╠═62c0385d-e712-43eb-adae-6731159f9f92 +# ╟─8ee34fde-40a7-4cc1-a8ea-056a457901b0 +# ╠═a43a5d7e-5052-4ee1-a112-c8881e90b6a6 +# ╠═c33ed1de-feb3-4bab-80c5-cd422be96bb6 +# ╠═5d7cf31a-7e08-45ea-80e4-a0e27a022258 +# ╠═4e651fc1-b6db-4c44-b830-f65c6f68f4a9 +# ╠═5f5113ed-d335-4038-a230-b739242aafbc +# ╠═f897cd86-ddc1-4b84-a11c-d99161c1df3d +# ╠═fcf4d2ab-512b-4a2c-a131-9d2519da1747 +# ╠═f7633864-230a-4a06-a48e-8acf4adb58f1 +# ╠═0a7b6bc6-9a23-484b-ac26-6effb901ea3b +# ╟─181eb167-f944-480d-907a-e54a2ba113df +# ╠═864c37c5-6907-4b31-8a4e-8c3c9a684606 +# ╠═d9339c77-35cf-4575-b935-79d580badba5 +# ╠═ad07d942-4985-40dc-9371-d6579e5b92a3 +# ╠═127e50d1-4538-4010-b7c3-716bffd804ed +# ╠═b0d00f6f-2a6c-42c0-b878-1791b2073353 +# ╠═20935643-b699-46c0-8a44-d948551d9c44 +# ╠═a0637882-140c-466f-ab7d-af9cdb979d1f +# ╠═ca2386b8-0882-4b31-90fa-7ab1df839457 +# ╠═17679b22-61ac-42c2-959c-5e28fe880d15 +# ╠═02c3e595-b978-4ad4-961d-5d7398d3d320 +# ╠═819f3128-3a42-425b-871d-a4913b17e94c diff --git a/examples/benchmark_CPU_CUDA.plutostate b/examples/benchmark_CPU_CUDA.plutostate new file mode 100644 index 0000000..d8a6443 Binary files /dev/null and b/examples/benchmark_CPU_CUDA.plutostate differ diff --git a/examples/double_slit.html b/examples/double_slit.html new file mode 100644 index 0000000..791c769 --- /dev/null +++ b/examples/double_slit.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/double_slit.jl b/examples/double_slit.jl new file mode 100644 index 0000000..e077f9c --- /dev/null +++ b/examples/double_slit.jl @@ -0,0 +1,262 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 75796f02-ac0d-11ee-3913-11fa1f818dea +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 9896f57c-0540-4be7-91fb-43f08f9e084c +using WaveOpticsPropagation, CUDA, IndexFunArrays, NDTools, ImageShow, FileIO, Zygote, Optim, Plots, PlutoUI + +# ╔═╡ 5a88b44b-ecce-4002-a4ba-87e8b9d55c81 +md"# Load and initialize packages" + +# ╔═╡ 0c95865d-7998-4191-a637-457c8d04d6d4 +TableOfContents() + +# ╔═╡ 6dc592f5-7fd3-4782-a72d-23b5dbfa80e2 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ aa9e7ca4-945a-4b8c-bc8b-29ab192350e3 +md"# Initialize physical parameters" + +# ╔═╡ a777cc7a-c79c-494c-9a70-83a90b2ee08d +L = 100f-3 + +# ╔═╡ 42ae161e-26b7-4d34-9e92-4c4af27fc447 +λ = 633f-9 + +# ╔═╡ 9facb2bf-f279-41de-8858-81dac7a032b5 +z = 1 + +# ╔═╡ 3cde3fbe-9641-4ffd-9630-9c7bf167ec34 +N = 256 + +# ╔═╡ 771642e4-c03e-41a2-a4de-07fdfd4bd535 +x = WaveOpticsPropagation.fftpos(L, N) + +# ╔═╡ 7f79a274-8f9c-43b5-bb7e-96bb495fb044 +L / N + +# ╔═╡ ac7c756a-b182-4458-b5ea-017e23fc1009 +2 / N * L + +# ╔═╡ 6acb2991-cba5-4cf5-8d82-099f9ad2335a +@bind ΔS Slider(0:N÷4, show_value=true, default=12) + +# ╔═╡ e5fc1e39-9300-4c96-b747-d532245cfd97 +S = x[1 + ΔS * 2] - x[1] + +# ╔═╡ 4f9c514a-5f24-4d49-8a92-88f250d9ef92 +@bind ΔW Slider(0:N÷2, show_value=true, default=2) + +# ╔═╡ d1734ddf-ec04-4134-82a8-187e697c4810 +begin + slit = zeros(ComplexF32, (N, N)) + + mid = N ÷ 2+ 1 + slit[:, mid-ΔS-ΔW:mid-ΔS+ΔW] .= 1 + slit[:, mid+ΔS-ΔW:mid+ΔS+ΔW] .= 1 +end; + +# ╔═╡ 3ebcb7d8-0ed6-497a-95e3-6e40bfadd8d0 +simshow(slit) + +# ╔═╡ ac72f69a-935e-4ba4-b8bd-e16f27367acd +W = x[1 + ΔW * 2 + 1] - x[1] + +# ╔═╡ 3330fef8-e6ed-412a-9321-9dc705ce41fe +md"# Propagate with Fraunhofer Diffraction" + +# ╔═╡ 7e797b20-77ee-4888-9b23-74bb2713912f +@time output, t = fraunhofer(slit, z, λ, L) + +# ╔═╡ 0c74a49c-e4de-4172-ac75-4c2692c505fb +# creating this function is more efficient! +efficient_fraunhofer, _ = Fraunhofer(slit, z, λ, L); + +# ╔═╡ 0e6f74f1-8723-4596-b754-973b253443a4 +@time output2, t2 = efficient_fraunhofer(slit) + +# ╔═╡ 5c0486d2-5fbd-4cd7-9fc5-583a8e188b2e +begin + intensity = abs2.(output) + intensity ./= maximum(intensity) +end; + +# ╔═╡ cd54d590-ea6f-4d3b-84ec-6f983a97c81c +simshow(intensity, γ=1) + +# ╔═╡ 5f3c7de8-a64e-42aa-b8a1-91f5869c9610 +md"# Compare to analytical solution" + +# ╔═╡ 4d836453-7e64-4f3b-9405-d3cdafc264dc +I_analytical(x) = sinc(W * x / λ / z)^2 * cos(π * S * x / λ / z)^2 + +# ╔═╡ 42e1a449-0e5c-4d67-a90c-10a54f479c8a +xpos_out = WaveOpticsPropagation.fftpos(t.L, N, NDTools.CenterFT) + +# ╔═╡ 5d36ba4d-e9fe-4b1b-b45d-c82b598dfe7a +begin + plot(xpos_out, intensity[(begin+end)÷2+1, :]) + plot!(xpos_out, I_analytical.(xpos_out)) +end + +# ╔═╡ ba6965bf-c7f1-4dfd-8381-e330582b3462 +all(≈(I_analytical.(xpos_out)[110:150] .+ 1, 1 .+ intensity[129, 110:150, :], rtol=1f-2)) + +# ╔═╡ 8193fd81-1580-4047-ac04-7c3837975c2d +md"# Optimization +We try to find the initial field from one diffraction pattern. + +We start with a initial guess of a double slit where distance and width are wrong. By using Fraunhofer as forward model, we are able to optimize for correct distance and width! +" + +# ╔═╡ 2c28b256-f60c-48b6-a34a-d2908c979e30 +intensity_cu = togoc(intensity); + +# ╔═╡ 3a19cbfe-6efc-45f7-9aa3-f319f34a534e +begin + rec0 = togoc(zeros(Float32, (N, N))); + + rec0[:, mid-40:mid-10] .= 1 + rec0[:, mid+10:mid+60] .= 1 +end; + +# ╔═╡ b19225e0-5536-4af1-bbbc-c6124e91536b +simshow(abs2.(Array(rec0))) + +# ╔═╡ e0aa0714-08e0-46f4-b6d4-8c5df044ddd9 +md"## Define Forward model and gradient +We use both methods: `fraunhofer` and `Fraunhofer`. +In principle `Fraunhofer` generates some buffers and stores the FFT plan and should be slightly faster. +" + +# ╔═╡ c8a3d8b0-f68d-4c38-ae1a-8081846bda23 +fwd = let z=z, L=L, λ=λ + x -> fraunhofer(x .+ 0im, z, λ, L) +end + +# ╔═╡ c9651b37-d911-4230-a2f5-f94e8f726fbc +fwd2 = let z=z, L=L, λ=λ, fr=Fraunhofer(rec0 .+ 0im, z, λ, L)[1] + x -> fr(x .+ 0im) +end + +# ╔═╡ e5c3f3dc-5eb3-421e-90ab-a16ea73e5621 +f = let intensity_cu=intensity_cu, fwd=fwd + f(x) = sum(abs2, intensity_cu .- abs2.(fwd(abs2.(x) .+ 0im)[1])) +end + +# ╔═╡ cb0b276f-d479-4024-8ca4-116a210fe2ed +f2 = let intensity_cu=intensity_cu, fwd2=fwd2 + f(x) = sum(abs2, intensity_cu .- abs2.(fwd2(abs2.(x) .+ 0im)[1])) +end + +# ╔═╡ 30555d35-04a8-44b7-afd7-345730d97a61 +@mytime f(rec0); + +# ╔═╡ c6c30873-9419-4a54-a9f5-6bf883856102 +@mytime f2(rec0); + +# ╔═╡ 8ab5cfd9-15e0-4be9-9936-cd499555fbec +@mytime gradient(f, rec0) + +# ╔═╡ dd8a5b39-b28a-4b79-857f-8c0a8923f2b3 +@mytime gradient(f2, rec0); + +# ╔═╡ c9a37a37-eef0-498e-8105-dd29284bfa4e +g!(G, x) = G .= Zygote.gradient(f, x)[1] + +# ╔═╡ 6795095b-bbbe-4dec-a79b-12d604873474 +g2!(G, x) = G .= Zygote.gradient(f2, x)[1] + +# ╔═╡ b75deb02-31e2-4b3d-adad-2f380b107e1b +md"## Optimize" + +# ╔═╡ c539d936-acda-4970-9f5d-92538fda2bc5 +res = optimize(f, g!, rec0, LBFGS(), Optim.Options(iterations=400)) + +# ╔═╡ 9563f911-b861-4cb1-87bc-fbbb3a3b77d2 +md"## Results +Up to a shift we obtain the correct result! The shift is not possible to find out since we only measure the intensity +" + +# ╔═╡ 8c53353c-7339-4060-8d2e-caf8f3025b48 +simshow(abs2.(Array(res.minimizer))) + +# ╔═╡ ae232c00-b4b7-41ff-8c49-a55382743956 +simshow(abs2.(Array(slit))) + +# ╔═╡ Cell order: +# ╟─5a88b44b-ecce-4002-a4ba-87e8b9d55c81 +# ╠═75796f02-ac0d-11ee-3913-11fa1f818dea +# ╟─0c95865d-7998-4191-a637-457c8d04d6d4 +# ╠═9896f57c-0540-4be7-91fb-43f08f9e084c +# ╠═6dc592f5-7fd3-4782-a72d-23b5dbfa80e2 +# ╟─aa9e7ca4-945a-4b8c-bc8b-29ab192350e3 +# ╠═a777cc7a-c79c-494c-9a70-83a90b2ee08d +# ╠═42ae161e-26b7-4d34-9e92-4c4af27fc447 +# ╠═9facb2bf-f279-41de-8858-81dac7a032b5 +# ╠═3cde3fbe-9641-4ffd-9630-9c7bf167ec34 +# ╠═d1734ddf-ec04-4134-82a8-187e697c4810 +# ╠═3ebcb7d8-0ed6-497a-95e3-6e40bfadd8d0 +# ╠═771642e4-c03e-41a2-a4de-07fdfd4bd535 +# ╠═ac72f69a-935e-4ba4-b8bd-e16f27367acd +# ╠═e5fc1e39-9300-4c96-b747-d532245cfd97 +# ╠═7f79a274-8f9c-43b5-bb7e-96bb495fb044 +# ╠═ac7c756a-b182-4458-b5ea-017e23fc1009 +# ╠═6acb2991-cba5-4cf5-8d82-099f9ad2335a +# ╠═4f9c514a-5f24-4d49-8a92-88f250d9ef92 +# ╟─3330fef8-e6ed-412a-9321-9dc705ce41fe +# ╠═7e797b20-77ee-4888-9b23-74bb2713912f +# ╠═0c74a49c-e4de-4172-ac75-4c2692c505fb +# ╠═0e6f74f1-8723-4596-b754-973b253443a4 +# ╠═5c0486d2-5fbd-4cd7-9fc5-583a8e188b2e +# ╠═cd54d590-ea6f-4d3b-84ec-6f983a97c81c +# ╟─5f3c7de8-a64e-42aa-b8a1-91f5869c9610 +# ╠═4d836453-7e64-4f3b-9405-d3cdafc264dc +# ╠═42e1a449-0e5c-4d67-a90c-10a54f479c8a +# ╠═5d36ba4d-e9fe-4b1b-b45d-c82b598dfe7a +# ╠═ba6965bf-c7f1-4dfd-8381-e330582b3462 +# ╟─8193fd81-1580-4047-ac04-7c3837975c2d +# ╠═2c28b256-f60c-48b6-a34a-d2908c979e30 +# ╠═3a19cbfe-6efc-45f7-9aa3-f319f34a534e +# ╠═b19225e0-5536-4af1-bbbc-c6124e91536b +# ╠═e0aa0714-08e0-46f4-b6d4-8c5df044ddd9 +# ╠═c8a3d8b0-f68d-4c38-ae1a-8081846bda23 +# ╠═c9651b37-d911-4230-a2f5-f94e8f726fbc +# ╠═e5c3f3dc-5eb3-421e-90ab-a16ea73e5621 +# ╠═cb0b276f-d479-4024-8ca4-116a210fe2ed +# ╠═30555d35-04a8-44b7-afd7-345730d97a61 +# ╠═c6c30873-9419-4a54-a9f5-6bf883856102 +# ╠═8ab5cfd9-15e0-4be9-9936-cd499555fbec +# ╠═dd8a5b39-b28a-4b79-857f-8c0a8923f2b3 +# ╠═c9a37a37-eef0-498e-8105-dd29284bfa4e +# ╠═6795095b-bbbe-4dec-a79b-12d604873474 +# ╟─b75deb02-31e2-4b3d-adad-2f380b107e1b +# ╠═c539d936-acda-4970-9f5d-92538fda2bc5 +# ╟─9563f911-b861-4cb1-87bc-fbbb3a3b77d2 +# ╠═8c53353c-7339-4060-8d2e-caf8f3025b48 +# ╠═ae232c00-b4b7-41ff-8c49-a55382743956 diff --git a/examples/double_slit.plutostate b/examples/double_slit.plutostate new file mode 100644 index 0000000..633dd94 Binary files /dev/null and b/examples/double_slit.plutostate differ diff --git a/examples/multi_slice_lens.html b/examples/multi_slice_lens.html new file mode 100644 index 0000000..66271ff --- /dev/null +++ b/examples/multi_slice_lens.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/multi_slice_lens.jl b/examples/multi_slice_lens.jl new file mode 100644 index 0000000..e1ae3a2 --- /dev/null +++ b/examples/multi_slice_lens.jl @@ -0,0 +1,243 @@ +### A Pluto.jl notebook ### +# v0.19.37 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ ac905778-bd13-11ee-2efc-c156980d778e +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ db5d6db1-0e20-48cf-90f9-f79a2f123197 +using WaveOpticsPropagation, CUDA, IndexFunArrays, NDTools, ImageShow, FileIO, Zygote, Optim, Plots, PlutoUI, FourierTools, NDTools + +# ╔═╡ af68a467-a2cb-4562-80b1-6c7bd57fe8ed +TableOfContents() + +# ╔═╡ 27266d4e-151e-4af7-becf-aecaa067e634 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ 286ef64e-275e-4ed5-8346-682962fb8c40 +ImageShow.simshow(x::CuArray) = ImageShow.simshow(Array(x)) + +# ╔═╡ ea0356ad-a47e-4895-adc4-ee9ce695e295 +Plotsheatmap(x::CuArray) = Plots.heatmap(Array(x)) + +# ╔═╡ 58b1397c-56b6-4f67-911c-f8d68802c405 +md"# 1. Define a ball lens" + +# ╔═╡ d059c790-9a91-41cd-8b54-f3628335010f +n1 = 1.0f0 + +# ╔═╡ 0228a79c-1518-4362-8268-8a2c39d1fec7 +n2 = 1.5f0 + +# ╔═╡ c2b77a7a-c359-4cf9-9459-db4a80cf1058 +sz = (512, 512) + +# ╔═╡ c88ba6a4-4b67-4ea8-826a-17e158158320 +radius = 15f-6 + +# ╔═╡ 2122864f-f77b-4aac-b347-a8209e5f4cad +md"# 2. Define field" + +# ╔═╡ ecf588f7-3bc4-4ccd-9837-4fef3c5e596f +L = 60f-6 + +# ╔═╡ 298962a7-a528-4322-ab4c-7e16a2cf2293 +x = fftpos(L, sz[2], CenterFT)' + +# ╔═╡ 797339e6-adbe-46df-a270-2737e550e3da +y = fftpos(L, sz[1], CenterFT) + +# ╔═╡ db47abea-4e8f-4c4e-a1bd-9a9796848638 +waist = 15f-6 + +# ╔═╡ 0f78e939-3b16-44c3-af17-ec6287ed6b23 +field = togoc(0im .+ exp.(.-(x.^2 .+ y.^2) ./ waist.^2)); + +# ╔═╡ 3e12fb1c-1067-498f-899f-56fbea5fe2ab +simshow(field) + +# ╔═╡ fe337699-a8ed-4691-be14-3d0f3f6847b8 +ztotal = L + +# ╔═╡ d92ecfaf-8a91-4e11-af67-0b58e3ee7a56 +Nz = 512 + +# ╔═╡ 24961b1c-0324-448c-a4a8-addd9add9bcd +z = radius / 3 .+ reshape(fftpos(ztotal, Nz, CenterFT), (1,1, Nz)); + +# ╔═╡ 1b6da0f1-267c-4df6-80d0-8454fb86611c +n = togoc(Float32.(n1 .+ (n2 - n1) .* (0 .* x .+ y.^2 .+ z.^2 .< radius .^2))); + +# ╔═╡ d20c628d-90fd-456e-8580-927b2bf6004a +Plots.heatmap(z[:], x[:], Array(n)[:,100,:]) + +# ╔═╡ 5484ba8e-cf44-4a6d-9fb4-0c787309c313 +Δz = ztotal / Nz + +# ╔═╡ eda88933-e0ef-4a5e-9b65-a52fe9f191c0 +λ = 632f-9 + +# ╔═╡ f184ba3c-84b5-4ad8-b8fd-a5ae39a11a66 +# for free space air propagation +AS = AngularSpectrum(field, Δz, λ / n1, L)[1] + +# ╔═╡ 0f7f0d36-5c3a-452e-a788-8f8e4713bddc +# used for WPM in medium n2 +AS2 = AngularSpectrum(field, Δz, λ / n2, L)[1] + +# ╔═╡ 67de273f-47c7-4978-a852-7ae06e104ff8 +md"# 3. Define Functions" + +# ╔═╡ c406fe60-f977-4a84-861e-8e75b4cbefa5 +# propagate in air -> phase shift -> propagate in air -> phase shift -> ... +function MSAS(field, n, AS, Nz, Δz, λ) + c = typeof(Δz)(2π) / λ * Δz + out = similar(field, (size(field)..., Nz)) + + out[:, :, 1] .= exp.(1im .* c .* (view(n, :, :, 1) .- n1)) .* view(field, :, :, 1) + + for i in 2:Nz + out[:, :, i] .= exp.(1im .* c .* (view(n, :, :, i) .- n1)) .* AS(@view out[:, :, i-1])[1] + end + + return out +end + +# ╔═╡ b80739d0-6417-47f4-b872-538b8d2c52d8 +function find_borders(n) + first = 0 + second = 0 + + for i in 2:size(n, 1) + if n[i] > n[i-1] + first = i + end + if n[i] < n[i-1] + second = i + end + end + return first, second +end + +# ╔═╡ e36f7a30-3a99-4f77-9e2a-c96ef8c37f21 +# similar to +# S. Schmidt, T. Tiess, S. Schröter, R. Hambach, M. Jäger, H. Bartelt, A. Tünnermann, and H. Gross, "Wave-optical modeling beyond the thin-element-approximation," Opt. Express 24, 30188-30200 (2016) +function WPM(field, n, AS1, AS2, Nz, Δz, λ) + c = typeof(Δz)(2π) / λ * Δz + out = similar(field, (size(field)..., Nz)) + + out[:, :, 1] .= view(field, :, :, 1) + + for i in 2:Nz + first, second = find_borders(Array(view(n, :, 127, i))) + + if first == 0 && second == 0 + out[:, :, i] .= AS1(@view out[:, :, i-1])[1] + else + res_n1 = AS1(@view out[:, :, i-1])[1] + res_n2 = AS2(@view out[:, :, i-1])[1] + + out[begin:first-1, :, i] .= res_n1[begin:first-1, :] + out[first:second-1, :, i] .= res_n2[first:second-1, :] + out[second:end, :, i] .= res_n1[second:end, :] + end + end + + return out +end + +# ╔═╡ 1990a28d-f87e-4b90-8027-4b27fe2f939b +@mytime out_WPM = WPM(field, n, AS, AS2, Nz, Δz, λ); + +# ╔═╡ 40547857-1a70-4194-a7ec-e90829d81aed +@mytime out_MAS = MSAS(field, n, AS, Nz, Δz, λ); + +# ╔═╡ 655f7c43-fa15-4882-b52a-8d7e2fb0b44d +md"# 4. Analyze results" + +# ╔═╡ f8825b79-85eb-4991-bf47-4ecaedd661e5 +begin + Plots.heatmap(z[:], x[:], Array(abs.(out_MAS[:,127,:])),title="Multi Slice Angular Spectrum") +end + +# ╔═╡ fb0e8b4b-ef66-4d40-8b6f-2232f7878303 +begin + Plots.heatmap(z[:], x[:], Array(abs.(out_WPM[:,127,:])), title="WPM") + #Plots.heatmap!(Array(0.0 .* n[:, 127, :] )) +end + +# ╔═╡ c386d509-e7ef-49ee-9d6d-7eb18d93a1ef +[simshow(out_WPM[:,127,:]) simshow(out_MAS[:,127,:])] + +# ╔═╡ 9ef44988-0a8d-44c1-b94d-155011182ba8 +@bind l Slider(1:4) + +# ╔═╡ 82341b70-eb05-4bca-9043-f90002ab6563 +simshow(Array(cat(n[:, 127, :], out_WPM[:, 127, :], out_MAS[:, 127, :], n[:, 127, :], dims=3)[:,:, l]), γ=0.4) + +# ╔═╡ d1f42c4b-0ccc-4d11-9ca3-dce49c98f3ad +simshow(Array(1 .* cispi.(n[:, 100, :]) .+ out_MAS[:,100,:]), γ=0.1) + +# ╔═╡ Cell order: +# ╠═ac905778-bd13-11ee-2efc-c156980d778e +# ╠═db5d6db1-0e20-48cf-90f9-f79a2f123197 +# ╟─af68a467-a2cb-4562-80b1-6c7bd57fe8ed +# ╠═27266d4e-151e-4af7-becf-aecaa067e634 +# ╠═286ef64e-275e-4ed5-8346-682962fb8c40 +# ╠═ea0356ad-a47e-4895-adc4-ee9ce695e295 +# ╟─58b1397c-56b6-4f67-911c-f8d68802c405 +# ╠═d059c790-9a91-41cd-8b54-f3628335010f +# ╠═0228a79c-1518-4362-8268-8a2c39d1fec7 +# ╠═c2b77a7a-c359-4cf9-9459-db4a80cf1058 +# ╠═298962a7-a528-4322-ab4c-7e16a2cf2293 +# ╠═797339e6-adbe-46df-a270-2737e550e3da +# ╠═24961b1c-0324-448c-a4a8-addd9add9bcd +# ╠═c88ba6a4-4b67-4ea8-826a-17e158158320 +# ╠═1b6da0f1-267c-4df6-80d0-8454fb86611c +# ╠═d20c628d-90fd-456e-8580-927b2bf6004a +# ╟─2122864f-f77b-4aac-b347-a8209e5f4cad +# ╠═0f78e939-3b16-44c3-af17-ec6287ed6b23 +# ╠═3e12fb1c-1067-498f-899f-56fbea5fe2ab +# ╠═ecf588f7-3bc4-4ccd-9837-4fef3c5e596f +# ╠═db47abea-4e8f-4c4e-a1bd-9a9796848638 +# ╠═fe337699-a8ed-4691-be14-3d0f3f6847b8 +# ╠═d92ecfaf-8a91-4e11-af67-0b58e3ee7a56 +# ╠═5484ba8e-cf44-4a6d-9fb4-0c787309c313 +# ╠═eda88933-e0ef-4a5e-9b65-a52fe9f191c0 +# ╠═f184ba3c-84b5-4ad8-b8fd-a5ae39a11a66 +# ╠═0f7f0d36-5c3a-452e-a788-8f8e4713bddc +# ╟─67de273f-47c7-4978-a852-7ae06e104ff8 +# ╠═c406fe60-f977-4a84-861e-8e75b4cbefa5 +# ╠═e36f7a30-3a99-4f77-9e2a-c96ef8c37f21 +# ╠═b80739d0-6417-47f4-b872-538b8d2c52d8 +# ╠═1990a28d-f87e-4b90-8027-4b27fe2f939b +# ╠═40547857-1a70-4194-a7ec-e90829d81aed +# ╟─655f7c43-fa15-4882-b52a-8d7e2fb0b44d +# ╠═f8825b79-85eb-4991-bf47-4ecaedd661e5 +# ╠═fb0e8b4b-ef66-4d40-8b6f-2232f7878303 +# ╠═c386d509-e7ef-49ee-9d6d-7eb18d93a1ef +# ╠═9ef44988-0a8d-44c1-b94d-155011182ba8 +# ╠═82341b70-eb05-4bca-9043-f90002ab6563 +# ╠═d1f42c4b-0ccc-4d11-9ca3-dce49c98f3ad diff --git a/examples/multi_slice_lens.plutostate b/examples/multi_slice_lens.plutostate new file mode 100644 index 0000000..abe32fb Binary files /dev/null and b/examples/multi_slice_lens.plutostate differ diff --git a/examples/old/Project.toml b/examples/old/Project.toml new file mode 100644 index 0000000..81648c0 --- /dev/null +++ b/examples/old/Project.toml @@ -0,0 +1 @@ +[deps] diff --git a/examples/old/gaussian_beam_SLM_propagation.html b/examples/old/gaussian_beam_SLM_propagation.html new file mode 100644 index 0000000..b4793a8 --- /dev/null +++ b/examples/old/gaussian_beam_SLM_propagation.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/old/gaussian_beam_SLM_propagation.jl b/examples/old/gaussian_beam_SLM_propagation.jl new file mode 100644 index 0000000..724f6ca --- /dev/null +++ b/examples/old/gaussian_beam_SLM_propagation.jl @@ -0,0 +1,184 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 97891d2b-facd-4f51-9774-5afde6a874e9 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 034e96b2-c88a-11ed-35f3-9fbfcefa62df +using WaveOpticsPropagation, ImageShow, FFTW, CUDA, FourierTools, NDTools, Plots, Colors + +# ╔═╡ 8082f7f3-3f47-4479-b3fa-20b2eb43ee78 +using LinearAlgebra + +# ╔═╡ b53bc181-42e6-410f-99fb-6fa28b9c4f09 +using IndexFunArrays + +# ╔═╡ f6d1a0dd-431f-4cbd-ae3c-3f031dc2beb5 +using PlutoUI + +# ╔═╡ febd1f55-b327-4e31-975b-afb635a82376 +using TestImages + +# ╔═╡ dfc85871-6697-44e1-a28a-86e7314a7d29 +FFTW.set_num_threads(4) + +# ╔═╡ 6e4cb613-6cec-4eb0-93a9-bc93b83bf755 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + + togoc(x) = use_CUDA[] ? CuArray(x) : x + toc(x) = x + toc(x::CuArray) = Array(x) + toc(x::LinearAlgebra.Adjoint{T, <:CuArray}) where T = Array(x) + +end + +# ╔═╡ 3dd7717c-0cbb-45f9-b03b-9316d37fd635 +ImageShow.simshow(x::CuArray; kwargs...) = simshow(toc(x); kwargs...) + +# ╔═╡ 15f89bfc-34c3-408e-9577-b79e88aea179 + + +# ╔═╡ 7e73e1e3-e881-4025-b88a-bb7b12f85e02 +N = 2000 + +# ╔═╡ e362673e-dbe9-4ad6-ba7a-d1f22da3f191 +field = togoc(zeros(ComplexF32, (N, N))) + +# ╔═╡ d6e11fcc-18ef-4c64-82b8-20abbb26d339 +function dmd_pattern() + dmd = togoc(zeros(ComplexF32, 768, 1024)); + + dmd[300:400, 300:400] .= 1 + dmd[390:410, 400:500] .= 1 + + dmd[450:451, 450:451] .= 1 + dmd[7 .+ (450:451), 450:451] .= 1 + + return dmd +end + +# ╔═╡ f0cd85d6-b8a4-4be5-8e99-0eebadc0d27b +simshow(dmd_pattern()) + +# ╔═╡ 5b586b1e-e785-4980-815f-0534f997ea41 +λ = 405f-9 / 1.5 + +# ╔═╡ 3ee91f4c-83a5-4e33-948f-88b36cf1828c +L = 14.0f-3 * 2000 / 1024 * 2 + +# ╔═╡ 25916008-4a44-4140-9558-6f055d01b101 +L / 2000 + +# ╔═╡ 49f045dd-b764-4e37-abc4-271a87db4677 +y = togoc(fftpos(L, N, CenterFT)) + +# ╔═╡ 2c301c7f-bc79-4141-8ae1-3022bfc54c66 +x = y'; + +# ╔═╡ b083214a-645f-499f-8aaa-e5511cbe64d8 +L ./ size(dmd) + +# ╔═╡ a11d219e-f6d2-4137-b0b7-be2ffa589893 +field_d = set_center!(copy(field), dmd_pattern()) .* cispi.(togoc(0.0000f0 .* rr2(field))); + +# ╔═╡ 5fc8ecab-6a7c-467f-a67d-d4ba7cd0b2e2 +simshow(field_d) + +# ╔═╡ 641e4caa-3d4a-46b1-af8f-a6ef8d246723 +@mytime field_p = angular_spectrum(field_d, 16f-3, λ, L)[1] + +# ╔═╡ b6db3c81-e960-4a29-b955-43721b3af769 +@bind iz PlutoUI.Slider(1:2, show_value=true) + +# ╔═╡ 3df00e4a-89f1-4ef3-a101-e02113fa471d +simshow(abs2.([toc(field_d);;; toc(field_p)][:, :, iz])) + +# ╔═╡ 6f31315a-f545-4989-91c9-b58db0e4ba18 +@view_image abs2.([toc(field_d);;; toc(field_p)]) + +# ╔═╡ 8e4c91eb-ed58-4d1f-9a43-48d2c262e4ea +begin + plot(abs.(resample(toc(field_d)[1000:1100, 938], (3000,)))) + plot!(abs2.(resample(toc(field_p)[1000:1100, 938], (3000,)))) +end; + +# ╔═╡ 661541a4-84f4-4f74-bbe7-3aab62fb52ce +begin + plot(abs2.(toc(field_d)[1000:1100, 938])) + plot!(abs2.(toc(field_p)[1000:1100, 938])) +end + +# ╔═╡ 3b90e7f5-daef-4181-81b7-db3d0a9a922b +z = range(100f-6, 10f-3, 5) + +# ╔═╡ 09ca6d0b-fbe0-4ff4-a900-2dd2179ab6bd +L / 1024 + +# ╔═╡ 8397198e-3b49-49db-af87-a932d4967a63 +heatmap(toc(x), toc(y), Float32.(simshow(abs2.(toc(field_p)), γ=0.4))) + +# ╔═╡ 3ba2a0a2-6535-464e-ae83-cddd68bff413 +toc(x) + +# ╔═╡ 0f232dcc-faf8-46d6-87bb-4fe318964968 +toc(x) + +# ╔═╡ 57cfdcf2-99c0-4db2-a740-ed74994116c9 +typeof(x) + +# ╔═╡ Cell order: +# ╠═97891d2b-facd-4f51-9774-5afde6a874e9 +# ╠═034e96b2-c88a-11ed-35f3-9fbfcefa62df +# ╠═8082f7f3-3f47-4479-b3fa-20b2eb43ee78 +# ╠═b53bc181-42e6-410f-99fb-6fa28b9c4f09 +# ╠═f6d1a0dd-431f-4cbd-ae3c-3f031dc2beb5 +# ╠═febd1f55-b327-4e31-975b-afb635a82376 +# ╠═3dd7717c-0cbb-45f9-b03b-9316d37fd635 +# ╠═dfc85871-6697-44e1-a28a-86e7314a7d29 +# ╠═6e4cb613-6cec-4eb0-93a9-bc93b83bf755 +# ╠═15f89bfc-34c3-408e-9577-b79e88aea179 +# ╠═7e73e1e3-e881-4025-b88a-bb7b12f85e02 +# ╠═e362673e-dbe9-4ad6-ba7a-d1f22da3f191 +# ╠═d6e11fcc-18ef-4c64-82b8-20abbb26d339 +# ╠═f0cd85d6-b8a4-4be5-8e99-0eebadc0d27b +# ╠═5b586b1e-e785-4980-815f-0534f997ea41 +# ╠═25916008-4a44-4140-9558-6f055d01b101 +# ╠═3ee91f4c-83a5-4e33-948f-88b36cf1828c +# ╠═49f045dd-b764-4e37-abc4-271a87db4677 +# ╠═2c301c7f-bc79-4141-8ae1-3022bfc54c66 +# ╠═b083214a-645f-499f-8aaa-e5511cbe64d8 +# ╠═a11d219e-f6d2-4137-b0b7-be2ffa589893 +# ╠═5fc8ecab-6a7c-467f-a67d-d4ba7cd0b2e2 +# ╠═641e4caa-3d4a-46b1-af8f-a6ef8d246723 +# ╠═b6db3c81-e960-4a29-b955-43721b3af769 +# ╠═3df00e4a-89f1-4ef3-a101-e02113fa471d +# ╠═6f31315a-f545-4989-91c9-b58db0e4ba18 +# ╠═8e4c91eb-ed58-4d1f-9a43-48d2c262e4ea +# ╠═661541a4-84f4-4f74-bbe7-3aab62fb52ce +# ╠═3b90e7f5-daef-4181-81b7-db3d0a9a922b +# ╠═09ca6d0b-fbe0-4ff4-a900-2dd2179ab6bd +# ╠═8397198e-3b49-49db-af87-a932d4967a63 +# ╠═3ba2a0a2-6535-464e-ae83-cddd68bff413 +# ╠═0f232dcc-faf8-46d6-87bb-4fe318964968 +# ╠═57cfdcf2-99c0-4db2-a740-ed74994116c9 diff --git a/examples/old/gaussian_beam_SLM_propagation.plutostate b/examples/old/gaussian_beam_SLM_propagation.plutostate new file mode 100644 index 0000000..07116f3 Binary files /dev/null and b/examples/old/gaussian_beam_SLM_propagation.plutostate differ diff --git a/examples/old/propagation_over_distance.html b/examples/old/propagation_over_distance.html new file mode 100644 index 0000000..2091cd0 --- /dev/null +++ b/examples/old/propagation_over_distance.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/old/propagation_over_distance.jl b/examples/old/propagation_over_distance.jl new file mode 100644 index 0000000..5124bd8 --- /dev/null +++ b/examples/old/propagation_over_distance.jl @@ -0,0 +1,196 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 6a9f04fe-ca20-11ed-1ed6-9967577ec81c +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ cba50443-4d13-4c05-b413-6a6f3d3ff1c2 +using WaveOpticsPropagation, ImageShow, FFTW, CUDA, FourierTools, NDTools, Plots, Colors, PlutoUI + +# ╔═╡ f07896d0-808a-470d-8dd8-ed5770a333db +using LinearAlgebra, IndexFunArrays + +# ╔═╡ 5a6fa509-7535-4188-9cb8-2ab0097d5bb0 +using TestImages + +# ╔═╡ ec354723-61c0-4966-a4f7-60efbc53b917 +Plots.plotly() + +# ╔═╡ 4d1cb1e3-327c-4fbc-a501-18efc2c86096 +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + + togoc(x) = use_CUDA[] ? CuArray(x) : x + toc(x) = Array(x) + toc(x::Array) = x + + ImageShow.simshow(x::CuArray; kwargs...) = simshow(toc(x); kwargs...) +end + +# ╔═╡ c3ba5b6b-33d2-4276-ab27-5009167d4d61 +N = 256 + +# ╔═╡ 7ef68ddb-474d-49ea-9510-91c0ab4d073d +field_img = togoc(ComplexF32.((rr2(Float32, (N, N), offset=(120, 200)) .< 10^2) .+ (rr2(Float32, (N, N), offset=(100, 100)) .< 40^2) .+ (rr2(Float32, (N, N), offset=(120, 170)) .< 3^2) )); + +# ╔═╡ e5fcbd81-c4f5-4bbc-8f22-dd248be8ec01 +L = 2f-3 + +# ╔═╡ 92d96248-06e5-446b-b856-163cbe8b4f45 +L / N + +# ╔═╡ da4d0740-46c4-460a-83ef-d78db575a3a4 +λ = 405f-9 / 1.5 + +# ╔═╡ 182fd5d8-857c-44d2-9451-8b6e8a9114e9 +zs = togoc(fftpos(16e-3, 256, CenterFT)); + +# ╔═╡ 38041742-1d46-4fa3-97b4-e5614f39ac33 +CUDA.@allowscalar zs[129] + +# ╔═╡ 851accf9-ea26-41f7-a2a5-4f0b9e53ad9c +@mytime field_img_p, tt = angular_spectrum(togoc(field_img), zs, λ, L); + +# ╔═╡ 62acd577-3402-4aee-9c62-4cc424b61758 +simshow(field_img) + +# ╔═╡ ddbe7f52-57dc-42d5-9edc-6377eb0bf31d +md"z=$(@bind iz PlutoUI.Slider(1:size(field_img_p, 3), show_value=true))" + +# ╔═╡ f9837f4a-e134-4d92-8984-51836bb89f74 +simshow((toc(abs2.(selectdim(field_img_p, 3, iz)))), γ=0.5) + +# ╔═╡ 4218847d-9a0c-4e4e-a8b6-6c76e0c3564e +begin + plot(fftpos(L, N, CenterFT), abs2.(toc(selectdim(field_img_p, 3, iz)))[120, :], title="z=$(1000 * toc(zs)[iz])mm", ticks=:native, marker=(:circle,1)) + plot!(fftpos(L, N, CenterFT), abs2.(toc(field_img)[120, :]), marker=(:circle,1)) +end + +# ╔═╡ 5fb47911-28b4-456c-b647-eed4af09f5e3 +sum(abs2.(field_img_p[:, :, 2])) + +# ╔═╡ a7325c95-d653-4e7e-ad9c-736a50ba9ea8 +sum(abs2.(field_img)) + +# ╔═╡ 03002431-9a33-4a6b-bb1a-d30f47c3bb45 + + +# ╔═╡ 5bb71ea0-cf7d-4be9-8df3-4d8bdab830d3 +simshow(togoc(selectdim(rs_kernel, 3, iz)), γ=0.1) + +# ╔═╡ 9cd1b7fd-d262-4ccc-8a2c-b82f0368ce35 +extrema(abs2.(rs_kernel)) + +# ╔═╡ 1b1144fe-b732-40bc-8382-727381f69b85 +size(field_img) + +# ╔═╡ 32ab279d-8c11-4753-9028-cb1983f1ed99 +function real_space_kernel(N, z, λ, L) + L = 2 * L + y = togoc(fftpos(L, N, CenterFT)) + x = y' + + r = sqrt.(x.^2 .+ y.^2 .+ z .^2 .+ 1f-7) + k = π / λ * 2 + kernel = z ./ r ./ 2 ./ π .* (1 ./ r .- 1im * k) .* exp.(1im .* k .* r) ./ r + return kernel +end + +# ╔═╡ dac63a43-73c2-44ab-bcaf-e3802b7e9714 +typeof(field_img) + +# ╔═╡ 5c7993ab-53bd-4991-900d-c8c64ebb76cb +field = CUDA.zeros(ComplexF32, 128, 128, 128); nothing + +# ╔═╡ 6ea1e637-112b-49f6-8526-b2a2ae20ef0c +sizeof(field) / 2^30 + +# ╔═╡ 31e0393f-9b1f-4dab-9a30-55b97d81dbac + + +# ╔═╡ 751c276f-6916-4db1-9a2f-5ec7379e45e3 +p = plan_fft!(field); nothing + +# ╔═╡ e61443e0-e3c6-4eff-98f7-e94d55f9f8fa +begin + GC.gc() + @mytime p * field; ; nothing +end + +# ╔═╡ 5ec4dbe5-9d88-45d4-9a47-3f2e888bb4c4 +begin + y = togoc(range(-10f0, 10f0, 100)) + x = togoc(range(-10f0, 10f0, 100)') + + z = togoc(reshape((range(0f0, 10f0, 100)), 1, 1, :)) +end + +# ╔═╡ 663ceff3-15be-4a61-8495-263af975de3c +@mytime exp.(1im .* z .* (x.^2 .+ y.^2)); + +# ╔═╡ d04ce01c-557f-441e-9968-c47df096349d +@mytime exp.(1im .* (x.^2 .+ y.^2)) .^ z; + +# ╔═╡ 81295271-bac4-44d7-8990-0741f0b8bb6e +NDTools.expand_dims + +# ╔═╡ 7feaec85-9604-4d3b-8f1d-c1da3df20a75 +@code_warntype WaveOpticsPropagation._prepare_angular_spectrum(field[:, :, 1], 1f0, 1f0, 1f0) + +# ╔═╡ Cell order: +# ╠═6a9f04fe-ca20-11ed-1ed6-9967577ec81c +# ╠═cba50443-4d13-4c05-b413-6a6f3d3ff1c2 +# ╠═f07896d0-808a-470d-8dd8-ed5770a333db +# ╠═ec354723-61c0-4966-a4f7-60efbc53b917 +# ╠═5a6fa509-7535-4188-9cb8-2ab0097d5bb0 +# ╠═4d1cb1e3-327c-4fbc-a501-18efc2c86096 +# ╠═c3ba5b6b-33d2-4276-ab27-5009167d4d61 +# ╠═7ef68ddb-474d-49ea-9510-91c0ab4d073d +# ╠═e5fcbd81-c4f5-4bbc-8f22-dd248be8ec01 +# ╠═92d96248-06e5-446b-b856-163cbe8b4f45 +# ╠═da4d0740-46c4-460a-83ef-d78db575a3a4 +# ╠═38041742-1d46-4fa3-97b4-e5614f39ac33 +# ╠═182fd5d8-857c-44d2-9451-8b6e8a9114e9 +# ╠═851accf9-ea26-41f7-a2a5-4f0b9e53ad9c +# ╠═62acd577-3402-4aee-9c62-4cc424b61758 +# ╠═f9837f4a-e134-4d92-8984-51836bb89f74 +# ╠═ddbe7f52-57dc-42d5-9edc-6377eb0bf31d +# ╠═4218847d-9a0c-4e4e-a8b6-6c76e0c3564e +# ╠═5fb47911-28b4-456c-b647-eed4af09f5e3 +# ╠═a7325c95-d653-4e7e-ad9c-736a50ba9ea8 +# ╠═03002431-9a33-4a6b-bb1a-d30f47c3bb45 +# ╠═5bb71ea0-cf7d-4be9-8df3-4d8bdab830d3 +# ╠═9cd1b7fd-d262-4ccc-8a2c-b82f0368ce35 +# ╠═1b1144fe-b732-40bc-8382-727381f69b85 +# ╠═32ab279d-8c11-4753-9028-cb1983f1ed99 +# ╠═dac63a43-73c2-44ab-bcaf-e3802b7e9714 +# ╠═6ea1e637-112b-49f6-8526-b2a2ae20ef0c +# ╠═5c7993ab-53bd-4991-900d-c8c64ebb76cb +# ╠═31e0393f-9b1f-4dab-9a30-55b97d81dbac +# ╠═751c276f-6916-4db1-9a2f-5ec7379e45e3 +# ╠═e61443e0-e3c6-4eff-98f7-e94d55f9f8fa +# ╠═5ec4dbe5-9d88-45d4-9a47-3f2e888bb4c4 +# ╠═663ceff3-15be-4a61-8495-263af975de3c +# ╠═d04ce01c-557f-441e-9968-c47df096349d +# ╠═81295271-bac4-44d7-8990-0741f0b8bb6e +# ╠═7feaec85-9604-4d3b-8f1d-c1da3df20a75 diff --git a/examples/old/propagation_over_distance.plutostate b/examples/old/propagation_over_distance.plutostate new file mode 100644 index 0000000..8e6b6e9 Binary files /dev/null and b/examples/old/propagation_over_distance.plutostate differ diff --git a/examples/optimize_phase.html b/examples/optimize_phase.html new file mode 100644 index 0000000..e996dcf --- /dev/null +++ b/examples/optimize_phase.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/optimize_phase.jl b/examples/optimize_phase.jl new file mode 100644 index 0000000..92e4352 --- /dev/null +++ b/examples/optimize_phase.jl @@ -0,0 +1,317 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 41c6fe0a-854c-11ee-0c0d-49cafb7df0d1 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 22b6b431-8bc1-4630-9a5e-998874831b84 +using Zygote, WaveOpticsPropagation, FFTW, ImageShow, FiniteDifferences, CUDA, Optim, IndexFunArrays, FourierTools, TestImages + +# ╔═╡ a715d3ef-b3ef-496f-b50f-7ef688df5efc +using PlutoUI + +# ╔═╡ dfacac83-565f-431b-a42c-1072c419a4be +using Plots + +# ╔═╡ d18059d4-9ff2-4afa-b40f-5c37db59b57f +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ a1c76c16-1eff-4cde-adaa-b544f6be45a7 +sz = (512, 512) + +# ╔═╡ 74c17e61-15d7-48cd-9fce-689d1c5a3148 +#field = zeros(ComplexF32, sz); + +# ╔═╡ fdd08423-020d-4dbd-9198-b7b03264a169 +#gield .= cis.(rr2(sz, scale=0.05)) .* gaussian(sz, sigma=50); + +field = (rr2(sz) .< 150 .^2)# .* conv((normal(Float32, (sz), sigma=2)), rand(Float32, sz) .* cispi.(2 .* ( 2 .* rand(Float32, (sz...))))) + +# ╔═╡ b2f8e5bf-fc0e-4dbc-9816-3cb0072132c2 +begin + field2 = box(Float32, sz, (100, 100)) + 1im .* box(Float32, sz, (100, 100), offset=(200, 300)) + field2 .= (0.1f0 .+ Float32.(Gray.(testimage("cameraman")))) .* cispi.(Float32.(Gray.(testimage("barbara_gray_512")))) +end + +# ╔═╡ c34da860-e131-483f-a2d7-ce2dcb1ba1e6 +simshow(field2) + +# ╔═╡ e661ac7e-63f9-40d7-ba9f-1784bce61e16 +field_c = togoc(cat(field2, field, dims=3)); + +# ╔═╡ a1ba22f5-8015-4e19-ab9c-bc82a4d38e44 +sum(abs2, field_c) + +# ╔═╡ 680c0d66-97af-4d21-b90c-97cb6d8dfbd9 +simshow(field) + +# ╔═╡ e67607af-d100-42f9-a830-4c5889584b3d +#sum(abs2, AS_c2(field_c)[1][:, :, 1]) + +# ╔═╡ 82a210f8-910b-42b8-bdd2-34b11011d34a +field_c + +# ╔═╡ 1beb527b-19f0-4352-839a-9ac7e91a3f15 +λ = 633f-9 + +# ╔═╡ c11c6881-3d50-49f6-8799-8e16f172ac9a +L = 5f-3 + +# ╔═╡ 3a6eb1f4-0763-4aa6-92c1-0abe491dd558 +z = 30f-3 + +# ╔═╡ 7be3d272-cf76-4950-b80b-f743ea200d4a +AS_c, _ = AngularSpectrum(field_c[:, :, 2], z, λ, L) + +# ╔═╡ f796abe0-00cc-4178-a945-96ae5a245227 +simshow(Array(AS_c(field_c[:, :, 2])[1]), γ=0.3) + +# ╔═╡ 9b03154d-bf00-4e95-a5a1-5ae5e7cb48a4 +simshow(Array(abs2.(AS_c(field_c)[1]))) + +# ╔═╡ 0180ea0b-87a8-4233-86fe-09308a70738d +N = 7 + +# ╔═╡ 1418e664-6207-4df4-aa79-f18633b63b2e +begin + diffuser = conv((normal(Float32, (sz), sigma=1)), cispi.(2 .* rand(Float32, (sz..., N)))) + diffuser_c = togoc(diffuser)# .* box(sz, (120, 150)) + + diffuser_c[:, :, 1] .= cis.(rr2(sz, scale=0.05)) .* gaussian(sz, sigma=100); + diffuser_c[:, :, 2] .= cis.(.-1 .* rr2(sz, scale=0.05)) .* gaussian(sz, sigma=100); + diffuser_c[:, :, 3] .= 1 + + diffuser_c[:, :, 4] .*= cis.(0.3 .* xx((sz))) + diffuser_c[:, :, 5] .*= cis.(.- 0.3 .* xx((sz))) + diffuser_c[:, :, 6] .*= cis.(0.3 .* yy((sz))) + diffuser_c[:, :, 7] .*= cis.(.- 0.3 .* yy((sz))) + +end + +# ╔═╡ 52f7af37-1a44-4a7c-8a9e-31df9bbcda9a +simshow(Array(diffuser_c[:, :, 5])) + +# ╔═╡ b436fdc0-ae8d-4dbe-94e6-8b7c26339a70 +#AS_c, _ = AngularSpectrum(field_c[:, :, 1], z, λ, L) + +# ╔═╡ c97f6332-82d9-4441-80ad-8a49c7739119 +AS_c2, _ = AngularSpectrum(togoc(zeros(ComplexF32, (sz..., N))), togoc(repeat([z], N)), λ, L) + +# ╔═╡ 0bd359cc-0bd4-4748-8905-361969811248 +p = plan_fft(diffuser_c, (1,2)) + +# ╔═╡ 23e64d50-8a6c-418f-86eb-1b4f9917fb32 +fwd(x) = begin + #first = AS_c(x)[1] + first = x[:, :, 2] + second = first .* diffuser_c + third = x[:, :, 1] .* AS_c2(second)[1] + #fourth = abs2.(fftshift(p * ifftshift(third, (1,2)), (1,2))) ./ size(x, 1).^2 + fourth = abs2.(AS_c2(third)[1]) + return fourth#, second, third +end + +# ╔═╡ a5593e55-5b77-4c99-b9a9-5e060748d6af +simshow(Array(fwd(field_c)[2][:, :, 8])) + +# ╔═╡ d00e3c52-a25d-4f13-9f76-f9fb0e255737 +simshow(Array(AS_c(fwd(field_c)[1][:, :, 3])[1])) + +# ╔═╡ 86c47180-3afe-4545-912e-12911155a02a +simshow(Array(diffuser_c)[:, :, 3]) + +# ╔═╡ 2c6ec01e-df45-4482-9d52-ce04ed343c5d +simshow(Array(AS_c2(diffuser_c)[1][:, :, 7])) + +# ╔═╡ df218193-d0f0-4c85-9439-9b70e2e55cf8 +simshow(Array(diffuser_c)[:, :, 9]) + +# ╔═╡ 2531c673-ab41-4374-85fe-e4e256f05c98 +field_c |> size + +# ╔═╡ 8516d9be-440c-45c5-8f20-cce74a9783e9 +function make_fg!(fwd, measured, loss=:L2) + L2 = let measured=measured, fwd=fwd + function L2(x) + return sum(abs2, fwd(x) .- measured) + end + end + + + f = let + if loss == :L2 + L2 + end + end + + g! = let f=f + function g!(G, rec) + if !isnothing(G) + return G .= Zygote.gradient(f, rec)[1] + end + end + end + return f, g! +end + +# ╔═╡ bdeb034b-8c34-4134-943c-eed63314cf20 +measurement_c = fwd(field_c); + +# ╔═╡ 525a1be2-286b-4626-b122-377084c6d456 +simshow(Array(measurement_c)[:, :, 1], γ=1) + +# ╔═╡ 18eb702e-ecad-4648-9e3a-b85cef60deb4 +begin + rec0 = togoc(ones(ComplexF32, (sz..., 2))) .* (rr2(sz) .< 150 .^2) + #rec0[:, :, 1] .= field_c[:, :, 1] + #rec0[:, :, 2] .= conv(CuArray(normal(sz, sigma=20)), field_c[:, :, 2]) + #rec0 .= cis.(rr2(sz, scale=0.04, offset=(200, 200))) .* box(sz, (80, 120)); +end + +# ╔═╡ c9776bdb-04ab-44e7-9c8e-9a7b8eb080f0 +simshow(Array(rec0)) + +# ╔═╡ 504074c0-9ff8-4b5e-8c6a-d196c8111a81 +f, g! = make_fg!(fwd, measurement_c) + +# ╔═╡ ad585dda-69c5-4bb5-a1bf-1e6001447669 +sum(abs2, measurement_c .* diffuser_c) + +# ╔═╡ 2d6dadd6-57a8-4c98-bf12-56276cb4cffe +@mytime res = Optim.optimize(f, g!, rec0, ConjugateGradient(), + Optim.Options(iterations = 300, + store_trace=true, + f_abstol=1e-10, + g_abstol=1e-10, )) + +# ╔═╡ 870e717a-577f-4097-a650-422e99fde221 +@time AS_c2(diffuser_c) + +# ╔═╡ a1bb8bb6-7b00-4478-84f3-1e350d530dbe +plot([t.iteration for t in res.trace], [t.value for t in res.trace], label="LBFGS", yaxis=:log) + +# ╔═╡ f252bf81-98d7-4383-9c29-fe4c2a034dc6 +simshow([Array(field_c[:, :, 1]) cispi(0.39) .* Array(res.minimizer)[:, :, 1]], γ=0.4) + +# ╔═╡ df93cc2e-4031-4011-a1d6-d76c62fd3000 +simshow(Array(field_c[:, :, 2])) + +# ╔═╡ 96c3bbfd-1b94-465e-b852-fadf7649ba57 +simshow(abs.([Array(field2) Array(res.minimizer)[:, :, 1]]), γ=1) + +# ╔═╡ 633dafe5-0f74-49a5-93bc-aaeec49d0254 +simshow(angle.([Array(field2) Array(res.minimizer)[:, :, 1]])) + +# ╔═╡ 13d98c90-296f-4fd4-bfc4-93ee81ae3484 +simshow(([Array(field2) Array(res.minimizer)[:, :, 2]])) + +# ╔═╡ afc890a3-25ae-4c4d-b886-a5433f317e73 +simshow(Array(res.minimizer)[:, :, 2]) + +# ╔═╡ 956c58a2-b182-48a9-b694-b8f6f9ba8077 +simshow(Array(field)) + +# ╔═╡ a9860e96-7888-4e68-9b9b-9adc0e260eef + + +# ╔═╡ 40c7f76c-2f40-403e-b160-d98ee9c54071 +@bind z_i Slider(1:size(measurement_c, 3)) + +# ╔═╡ 4823c477-38eb-47f0-89fa-6f02a0940d35 +#simshow(Array(fwd(res.minimizer)[:, :, z_i]) .- Array(measurement_c[:, :, z_i])) + +# ╔═╡ 00aaa61e-48a1-45dd-bdfe-5bd4879cd2d3 +simshow(Array(fwd(res.minimizer)[:, :, z_i]), γ=0.8) + +# ╔═╡ 6f57680b-7f7a-4d69-a62c-ccfce645151c +simshow(Array(measurement_c[:, :, z_i]), γ=0.8) + +# ╔═╡ 88239ba1-5528-46c3-9914-3ea62dfc0cfb + + +# ╔═╡ a3987c52-33d8-48f4-b373-4b6df1d6c17e +fwd(res.minimizer) + +# ╔═╡ Cell order: +# ╠═41c6fe0a-854c-11ee-0c0d-49cafb7df0d1 +# ╠═d18059d4-9ff2-4afa-b40f-5c37db59b57f +# ╠═22b6b431-8bc1-4630-9a5e-998874831b84 +# ╠═a715d3ef-b3ef-496f-b50f-7ef688df5efc +# ╠═dfacac83-565f-431b-a42c-1072c419a4be +# ╠═a1c76c16-1eff-4cde-adaa-b544f6be45a7 +# ╠═74c17e61-15d7-48cd-9fce-689d1c5a3148 +# ╠═fdd08423-020d-4dbd-9198-b7b03264a169 +# ╠═b2f8e5bf-fc0e-4dbc-9816-3cb0072132c2 +# ╠═c34da860-e131-483f-a2d7-ce2dcb1ba1e6 +# ╠═e661ac7e-63f9-40d7-ba9f-1784bce61e16 +# ╠═a1ba22f5-8015-4e19-ab9c-bc82a4d38e44 +# ╠═680c0d66-97af-4d21-b90c-97cb6d8dfbd9 +# ╠═7be3d272-cf76-4950-b80b-f743ea200d4a +# ╠═f796abe0-00cc-4178-a945-96ae5a245227 +# ╠═e67607af-d100-42f9-a830-4c5889584b3d +# ╠═82a210f8-910b-42b8-bdd2-34b11011d34a +# ╟─9b03154d-bf00-4e95-a5a1-5ae5e7cb48a4 +# ╠═1beb527b-19f0-4352-839a-9ac7e91a3f15 +# ╠═c11c6881-3d50-49f6-8799-8e16f172ac9a +# ╠═3a6eb1f4-0763-4aa6-92c1-0abe491dd558 +# ╠═0180ea0b-87a8-4233-86fe-09308a70738d +# ╠═1418e664-6207-4df4-aa79-f18633b63b2e +# ╠═52f7af37-1a44-4a7c-8a9e-31df9bbcda9a +# ╠═b436fdc0-ae8d-4dbe-94e6-8b7c26339a70 +# ╠═c97f6332-82d9-4441-80ad-8a49c7739119 +# ╠═0bd359cc-0bd4-4748-8905-361969811248 +# ╠═23e64d50-8a6c-418f-86eb-1b4f9917fb32 +# ╠═a5593e55-5b77-4c99-b9a9-5e060748d6af +# ╠═d00e3c52-a25d-4f13-9f76-f9fb0e255737 +# ╠═86c47180-3afe-4545-912e-12911155a02a +# ╠═2c6ec01e-df45-4482-9d52-ce04ed343c5d +# ╠═df218193-d0f0-4c85-9439-9b70e2e55cf8 +# ╠═2531c673-ab41-4374-85fe-e4e256f05c98 +# ╠═8516d9be-440c-45c5-8f20-cce74a9783e9 +# ╠═bdeb034b-8c34-4134-943c-eed63314cf20 +# ╠═525a1be2-286b-4626-b122-377084c6d456 +# ╠═18eb702e-ecad-4648-9e3a-b85cef60deb4 +# ╠═c9776bdb-04ab-44e7-9c8e-9a7b8eb080f0 +# ╠═504074c0-9ff8-4b5e-8c6a-d196c8111a81 +# ╠═ad585dda-69c5-4bb5-a1bf-1e6001447669 +# ╠═2d6dadd6-57a8-4c98-bf12-56276cb4cffe +# ╠═870e717a-577f-4097-a650-422e99fde221 +# ╠═a1bb8bb6-7b00-4478-84f3-1e350d530dbe +# ╠═f252bf81-98d7-4383-9c29-fe4c2a034dc6 +# ╠═df93cc2e-4031-4011-a1d6-d76c62fd3000 +# ╠═96c3bbfd-1b94-465e-b852-fadf7649ba57 +# ╠═633dafe5-0f74-49a5-93bc-aaeec49d0254 +# ╠═13d98c90-296f-4fd4-bfc4-93ee81ae3484 +# ╠═afc890a3-25ae-4c4d-b886-a5433f317e73 +# ╠═956c58a2-b182-48a9-b694-b8f6f9ba8077 +# ╠═a9860e96-7888-4e68-9b9b-9adc0e260eef +# ╠═40c7f76c-2f40-403e-b160-d98ee9c54071 +# ╠═4823c477-38eb-47f0-89fa-6f02a0940d35 +# ╠═00aaa61e-48a1-45dd-bdfe-5bd4879cd2d3 +# ╠═6f57680b-7f7a-4d69-a62c-ccfce645151c +# ╠═88239ba1-5528-46c3-9914-3ea62dfc0cfb +# ╠═a3987c52-33d8-48f4-b373-4b6df1d6c17e diff --git a/examples/optimize_phase.plutostate b/examples/optimize_phase.plutostate new file mode 100644 index 0000000..2d81cf8 Binary files /dev/null and b/examples/optimize_phase.plutostate differ diff --git a/examples/shifted_angular_spectrum.html b/examples/shifted_angular_spectrum.html new file mode 100644 index 0000000..6e96271 --- /dev/null +++ b/examples/shifted_angular_spectrum.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/examples/shifted_angular_spectrum.jl b/examples/shifted_angular_spectrum.jl new file mode 100644 index 0000000..36d103d --- /dev/null +++ b/examples/shifted_angular_spectrum.jl @@ -0,0 +1,139 @@ +### A Pluto.jl notebook ### +# v0.19.36 + +using Markdown +using InteractiveUtils + +# ╔═╡ b11e7be2-b315-11ee-27e7-abecfdbe64b6 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 3a5d9f20-a01d-481b-9858-b8e523ba7a20 +using WaveOpticsPropagation, CUDA, IndexFunArrays, NDTools, ImageShow, FileIO, Zygote, Optim, Plots, PlutoUI, FourierTools, NDTools + +# ╔═╡ dfc515b5-cfb5-4004-981f-a2262da47bab +begin + # use CUDA if functional + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ 517b00de-e25a-4688-ac2a-5ca067d7cef7 +md"# Define field with a ramp" + +# ╔═╡ 2f6871e8-7c11-49c0-ba9a-dc498e8eb39d +N = 256 + +# ╔═╡ 64b448ee-5ccc-4f87-8ee0-20d2d6a41a3b +sz = (N, N) + +# ╔═╡ 90286b89-aedd-4ece-b9d6-e5c26c6ad635 +α = deg2rad(10f0) + +# ╔═╡ dc01bc87-ffd7-400f-bbf2-3b00a3a84b78 +L = 50f-6 + +# ╔═╡ fdb36c00-57e6-4e3a-a9af-ed1282cf774a +y = fftpos(L[1], N, CenterFT) + +# ╔═╡ 89ec7708-f439-4881-9349-f46d0e75ea93 +λ = 405f-9 + +# ╔═╡ cfeb277b-3bc0-4371-b0ab-587ed626ea6c +field = box(Float32, sz, (20,20)) .* exp.(1im .* 2f0 * π ./ λ .* y .* sin(α)); + +# ╔═╡ ea02bb1c-7098-4c44-bc13-f9f62fcdce48 +z = 100f-6 + +# ╔═╡ 391ca41e-731d-4799-b09d-553c12b949d7 +simshow(field) + +# ╔═╡ 83982659-dccb-4691-bac1-53abcfc9a88b +md"# Compare AS and shifted AS" + +# ╔═╡ 2a0ef89b-d9ae-4186-9ccf-15d7785ff407 +res_AS = AngularSpectrum(field, z, λ, L)[1](field)[1]; + +# ╔═╡ 9cbafe25-3af6-4bcf-833c-8d3d7ca428a2 +res = ShiftedAngularSpectrum(field, z, λ, L, (α , 0), bandlimit=true)[1](field) + +# ╔═╡ 4efdc02b-4f69-4893-a410-6c6bbb765bab +shift = (z .* tan.(α) ./ L .* N)[1] + +# ╔═╡ 24b045a3-7828-4e25-a5bc-656f29cb8166 +shift2 = z .* tan.(α) ./ L + +# ╔═╡ 8afd2051-66dc-4b46-b7d8-13dd752b98da +md" +Left is shifted AS, middle is shifted AS but shifted in real space such that it roughly fits to AS. Right is AS +" + +# ╔═╡ cd6f41c7-532b-4681-98e9-fba8a05fb86b +[simshow(res[1][round(Int, shift)+1:end, :], γ=1) simshow(FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :], γ=1) simshow(res_AS[round(Int, shift)+1:end, :], γ=1)] + +# ╔═╡ 0c96fc1e-57ee-4f7f-8c0c-2cc8e65cc5b1 +simshow(FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :], γ=1) + +# ╔═╡ 6acb594c-2f4d-44a6-b795-786b71c9657b +simshow(res_AS[round(Int, shift)+1:end, :], γ=1) + +# ╔═╡ 7c9fd69b-9741-40b8-ae41-f4ce19d34594 + simshow(FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :] .- res_AS[round(Int, shift)+1:end, :], γ=1) + +# ╔═╡ 1c2b9b4b-6190-476b-8e6b-fab20147f0e9 +sum(abs2, res_AS[round(Int, shift)+1:end, :] .- FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :]) + +# ╔═╡ 365d1605-7cb7-43d3-a428-85621e56dcd6 +sum(abs2, res_AS[round(Int, shift)+1:end, :]) + +# ╔═╡ 154b48f2-58b8-4e4f-8648-8ea771125be5 +sum(abs2, FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :]) + +# ╔═╡ f3f7bc9e-a16f-49d3-920f-031326f5f1af +all(.≈(1 .+ FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :], 1 .+ res_AS[round(Int, shift)+1:end, :], rtol=5f-2)) + +# ╔═╡ d7eac41e-c5c0-46f8-9b2f-d2f116b65d95 + + +# ╔═╡ 4c3d7581-df26-4bec-9e0e-44279158d8b9 + + +# ╔═╡ ed3ec3b6-4825-4cd8-9d95-13d978d64584 + + +# ╔═╡ Cell order: +# ╠═b11e7be2-b315-11ee-27e7-abecfdbe64b6 +# ╠═3a5d9f20-a01d-481b-9858-b8e523ba7a20 +# ╠═dfc515b5-cfb5-4004-981f-a2262da47bab +# ╟─517b00de-e25a-4688-ac2a-5ca067d7cef7 +# ╠═2f6871e8-7c11-49c0-ba9a-dc498e8eb39d +# ╠═64b448ee-5ccc-4f87-8ee0-20d2d6a41a3b +# ╠═fdb36c00-57e6-4e3a-a9af-ed1282cf774a +# ╠═cfeb277b-3bc0-4371-b0ab-587ed626ea6c +# ╠═90286b89-aedd-4ece-b9d6-e5c26c6ad635 +# ╠═dc01bc87-ffd7-400f-bbf2-3b00a3a84b78 +# ╠═89ec7708-f439-4881-9349-f46d0e75ea93 +# ╠═ea02bb1c-7098-4c44-bc13-f9f62fcdce48 +# ╠═391ca41e-731d-4799-b09d-553c12b949d7 +# ╟─83982659-dccb-4691-bac1-53abcfc9a88b +# ╠═2a0ef89b-d9ae-4186-9ccf-15d7785ff407 +# ╠═9cbafe25-3af6-4bcf-833c-8d3d7ca428a2 +# ╠═4efdc02b-4f69-4893-a410-6c6bbb765bab +# ╠═24b045a3-7828-4e25-a5bc-656f29cb8166 +# ╟─8afd2051-66dc-4b46-b7d8-13dd752b98da +# ╠═cd6f41c7-532b-4681-98e9-fba8a05fb86b +# ╠═0c96fc1e-57ee-4f7f-8c0c-2cc8e65cc5b1 +# ╠═6acb594c-2f4d-44a6-b795-786b71c9657b +# ╠═7c9fd69b-9741-40b8-ae41-f4ce19d34594 +# ╠═1c2b9b4b-6190-476b-8e6b-fab20147f0e9 +# ╠═365d1605-7cb7-43d3-a428-85621e56dcd6 +# ╠═154b48f2-58b8-4e4f-8648-8ea771125be5 +# ╠═f3f7bc9e-a16f-49d3-920f-031326f5f1af +# ╠═d7eac41e-c5c0-46f8-9b2f-d2f116b65d95 +# ╠═4c3d7581-df26-4bec-9e0e-44279158d8b9 +# ╟─ed3ec3b6-4825-4cd8-9d95-13d978d64584 diff --git a/examples/shifted_angular_spectrum.plutostate b/examples/shifted_angular_spectrum.plutostate new file mode 100644 index 0000000..b69324e Binary files /dev/null and b/examples/shifted_angular_spectrum.plutostate differ diff --git a/index.html b/index.html new file mode 100644 index 0000000..78ff4d8 --- /dev/null +++ b/index.html @@ -0,0 +1,23 @@ + + + + + + + + + + + + +

Notebooks

+ + + + diff --git a/pluto_export.json b/pluto_export.json new file mode 100644 index 0000000..36d8109 --- /dev/null +++ b/pluto_export.json @@ -0,0 +1 @@ +{"notebooks":{"examples/shifted_angular_spectrum.jl":{"id":"examples/shifted_angular_spectrum.jl","hash":"A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo","html_path":"examples/shifted_angular_spectrum.html","statefile_path":"examples/shifted_angular_spectrum.plutostate","notebookfile_path":"examples/shifted_angular_spectrum.jl","current_hash":"A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo","desired_hash":"A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo","frontmatter":{"title":"shifted_angular_spectrum"}},"examples/old/propagation_over_distance.jl":{"id":"examples/old/propagation_over_distance.jl","hash":"-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ","html_path":"examples/old/propagation_over_distance.html","statefile_path":"examples/old/propagation_over_distance.plutostate","notebookfile_path":"examples/old/propagation_over_distance.jl","current_hash":"-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ","desired_hash":"-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ","frontmatter":{"title":"propagation_over_distance"}},"examples/multi_slice_lens.jl":{"id":"examples/multi_slice_lens.jl","hash":"9wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM","html_path":"examples/multi_slice_lens.html","statefile_path":"examples/multi_slice_lens.plutostate","notebookfile_path":"examples/multi_slice_lens.jl","current_hash":"9wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM","desired_hash":"9wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM","frontmatter":{"title":"multi_slice_lens"}},"examples/old/gaussian_beam_SLM_propagation.jl":{"id":"examples/old/gaussian_beam_SLM_propagation.jl","hash":"Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18","html_path":"examples/old/gaussian_beam_SLM_propagation.html","statefile_path":"examples/old/gaussian_beam_SLM_propagation.plutostate","notebookfile_path":"examples/old/gaussian_beam_SLM_propagation.jl","current_hash":"Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18","desired_hash":"Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18","frontmatter":{"title":"gaussian_beam_SLM_propagation"}},"examples/Scalable_Angular_Spectrum.jl":{"id":"examples/Scalable_Angular_Spectrum.jl","hash":"rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE","html_path":"examples/Scalable_Angular_Spectrum.html","statefile_path":"examples/Scalable_Angular_Spectrum.plutostate","notebookfile_path":"examples/Scalable_Angular_Spectrum.jl","current_hash":"rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE","desired_hash":"rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE","frontmatter":{"title":"Scalable_Angular_Spectrum"}},"test/double_slit.jl":{"id":"test/double_slit.jl","hash":"pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc","html_path":"test/double_slit.html","statefile_path":"test/double_slit.plutostate","notebookfile_path":"test/double_slit.jl","current_hash":"pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc","desired_hash":"pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc","frontmatter":{"title":"double_slit"}},"examples/double_slit.jl":{"id":"examples/double_slit.jl","hash":"Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0","html_path":"examples/double_slit.html","statefile_path":"examples/double_slit.plutostate","notebookfile_path":"examples/double_slit.jl","current_hash":"Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0","desired_hash":"Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0","frontmatter":{"title":"double_slit"}},"examples/optimize_phase.jl":{"id":"examples/optimize_phase.jl","hash":"nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA","html_path":"examples/optimize_phase.html","statefile_path":"examples/optimize_phase.plutostate","notebookfile_path":"examples/optimize_phase.jl","current_hash":"nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA","desired_hash":"nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA","frontmatter":{"title":"optimize_phase"}},"examples/benchmark_CPU_CUDA.jl":{"id":"examples/benchmark_CPU_CUDA.jl","hash":"GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA","html_path":"examples/benchmark_CPU_CUDA.html","statefile_path":"examples/benchmark_CPU_CUDA.plutostate","notebookfile_path":"examples/benchmark_CPU_CUDA.jl","current_hash":"GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA","desired_hash":"GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA","frontmatter":{"title":"benchmark_CPU_CUDA"}},"docs/logo/logo-pluto.jl":{"id":"docs/logo/logo-pluto.jl","hash":"ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U","html_path":"docs/logo/logo-pluto.html","statefile_path":"docs/logo/logo-pluto.plutostate","notebookfile_path":"docs/logo/logo-pluto.jl","current_hash":"ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U","desired_hash":"ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U","frontmatter":{"title":"logo-pluto"}},"test/gaussian_beam.jl":{"id":"test/gaussian_beam.jl","hash":"XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw","html_path":"test/gaussian_beam.html","statefile_path":"test/gaussian_beam.plutostate","notebookfile_path":"test/gaussian_beam.jl","current_hash":"XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw","desired_hash":"XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw","frontmatter":{"title":"gaussian_beam"}}},"pluto_version":"0.19.38","julia_version":"1.10.0","format_version":"1","title":null,"description":null,"collections":null,"binder_url":"https://mybinder.org/v2/gh/fonsp/pluto-on-binder/v0.19.38","slider_server_url":null} \ No newline at end of file diff --git a/pluto_state_cache/0_19_37-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate b/pluto_state_cache/0_19_37-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate new file mode 100644 index 0000000..68d3268 Binary files /dev/null and b/pluto_state_cache/0_19_37-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate differ diff --git a/pluto_state_cache/0_19_379wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate b/pluto_state_cache/0_19_379wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate new file mode 100644 index 0000000..21a65ce Binary files /dev/null and b/pluto_state_cache/0_19_379wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate differ diff --git a/pluto_state_cache/0_19_37A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate b/pluto_state_cache/0_19_37A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate new file mode 100644 index 0000000..0456170 Binary files /dev/null and b/pluto_state_cache/0_19_37A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate differ diff --git a/pluto_state_cache/0_19_37Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate b/pluto_state_cache/0_19_37Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate new file mode 100644 index 0000000..b6aa75f Binary files /dev/null and b/pluto_state_cache/0_19_37Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate differ diff --git a/pluto_state_cache/0_19_37Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate b/pluto_state_cache/0_19_37Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate new file mode 100644 index 0000000..ce5ca86 Binary files /dev/null and b/pluto_state_cache/0_19_37Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate differ diff --git a/pluto_state_cache/0_19_37GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate b/pluto_state_cache/0_19_37GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate new file mode 100644 index 0000000..af0ad86 Binary files /dev/null and b/pluto_state_cache/0_19_37GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate differ diff --git a/pluto_state_cache/0_19_37ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate b/pluto_state_cache/0_19_37ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate new file mode 100644 index 0000000..808f6be Binary files /dev/null and b/pluto_state_cache/0_19_37ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate differ diff --git a/pluto_state_cache/0_19_37XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate b/pluto_state_cache/0_19_37XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate new file mode 100644 index 0000000..4490a6a Binary files /dev/null and b/pluto_state_cache/0_19_37XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate differ diff --git a/pluto_state_cache/0_19_37nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate b/pluto_state_cache/0_19_37nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate new file mode 100644 index 0000000..a1bd651 Binary files /dev/null and b/pluto_state_cache/0_19_37nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate differ diff --git a/pluto_state_cache/0_19_37pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate b/pluto_state_cache/0_19_37pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate new file mode 100644 index 0000000..e6104b9 Binary files /dev/null and b/pluto_state_cache/0_19_37pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate differ diff --git a/pluto_state_cache/0_19_37rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate b/pluto_state_cache/0_19_37rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate new file mode 100644 index 0000000..f59048d Binary files /dev/null and b/pluto_state_cache/0_19_37rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate differ diff --git a/pluto_state_cache/0_19_38-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate b/pluto_state_cache/0_19_38-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate new file mode 100644 index 0000000..8e6b6e9 Binary files /dev/null and b/pluto_state_cache/0_19_38-ka3zcuzxh-YOaK8EUIr0SmXuIVxHctKYKenKPon4YQ.plutostate differ diff --git a/pluto_state_cache/0_19_389wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate b/pluto_state_cache/0_19_389wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate new file mode 100644 index 0000000..abe32fb Binary files /dev/null and b/pluto_state_cache/0_19_389wj2PZIHXBtrj3IZ8wMw4gDpQoaR-nCKujLD52Jo_mM.plutostate differ diff --git a/pluto_state_cache/0_19_38A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate b/pluto_state_cache/0_19_38A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate new file mode 100644 index 0000000..b69324e Binary files /dev/null and b/pluto_state_cache/0_19_38A6wzFHTBFXiWFyJAFyZri6foFsqWLJ8XPS-ZeC2Y0jo.plutostate differ diff --git a/pluto_state_cache/0_19_38Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate b/pluto_state_cache/0_19_38Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate new file mode 100644 index 0000000..07116f3 Binary files /dev/null and b/pluto_state_cache/0_19_38Fd2W88vhOghN87HsQH7OvTQn7KbdzTwkiyeB26S2R18.plutostate differ diff --git a/pluto_state_cache/0_19_38Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate b/pluto_state_cache/0_19_38Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate new file mode 100644 index 0000000..633dd94 Binary files /dev/null and b/pluto_state_cache/0_19_38Fp9Dq94fAoMu6CrQ6SovXXiSLbiC4FscPXG6qoQpRL0.plutostate differ diff --git a/pluto_state_cache/0_19_38GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate b/pluto_state_cache/0_19_38GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate new file mode 100644 index 0000000..d8a6443 Binary files /dev/null and b/pluto_state_cache/0_19_38GAhTvM70TPCYaUQci0yaHVWzmssu7ayYBoiqtuaJusA.plutostate differ diff --git a/pluto_state_cache/0_19_38ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate b/pluto_state_cache/0_19_38ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate new file mode 100644 index 0000000..93f1229 Binary files /dev/null and b/pluto_state_cache/0_19_38ME3lQ5efIT7pTqHjVPYKBS9LYdnYAUhbr8H5gAbhY6U.plutostate differ diff --git a/pluto_state_cache/0_19_38XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate b/pluto_state_cache/0_19_38XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate new file mode 100644 index 0000000..be95974 Binary files /dev/null and b/pluto_state_cache/0_19_38XxVlr_jDDoP_tENKAHP_AUBLZUnYnWGZH1FEzRhoMiw.plutostate differ diff --git a/pluto_state_cache/0_19_38nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate b/pluto_state_cache/0_19_38nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate new file mode 100644 index 0000000..2d81cf8 Binary files /dev/null and b/pluto_state_cache/0_19_38nEusZ7JE6zKHopY-eIcnMuO7TI2s0-vjx284Ia_eDlA.plutostate differ diff --git a/pluto_state_cache/0_19_38pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate b/pluto_state_cache/0_19_38pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate new file mode 100644 index 0000000..a142d2f Binary files /dev/null and b/pluto_state_cache/0_19_38pXsuKzsLET9WEGPWCRlcQP6JiJRkdndh8u925JYP4vc.plutostate differ diff --git a/pluto_state_cache/0_19_38rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate b/pluto_state_cache/0_19_38rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate new file mode 100644 index 0000000..906bc53 Binary files /dev/null and b/pluto_state_cache/0_19_38rwlZhE1IioD_S1IzI7WG_7waZ4LHQEjCv27bH8vQdzE.plutostate differ diff --git a/src/WaveOpticsPropagation.jl b/src/WaveOpticsPropagation.jl new file mode 100644 index 0000000..1959d53 --- /dev/null +++ b/src/WaveOpticsPropagation.jl @@ -0,0 +1,21 @@ +module WaveOpticsPropagation + +using EllipsisNotation +using FFTW +using ChainRulesCore +using Zygote +using NDTools +using FourierTools +using IndexFunArrays +using CUDA + +include("utils.jl") +include("propagation.jl") +include("angular_spectrum.jl") +include("shifted_angular_spectrum.jl") +include("scalable_angular_spectrum.jl") +include("fraunhofer.jl") +include("beams.jl") +include("conv.jl") + +end diff --git a/src/angular_spectrum.jl b/src/angular_spectrum.jl new file mode 100644 index 0000000..70b2f7a --- /dev/null +++ b/src/angular_spectrum.jl @@ -0,0 +1,218 @@ +export angular_spectrum +export AngularSpectrum + + +_transform_z(::Type{T}, z::Number) where {T<:Number} = T(z) +_transform_z(::Type{T}, z::AbstractArray{T}) where T = reshape(z, 1, 1, :) +_transform_z(::Type{T}, z::AbstractArray{T2}) where {T, T2} = reshape(T.(z), 1, 1, :) + +function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1.0)) where {CT<:Complex} + + fftdims = (1,2) + T = real(CT) + λ = T(λ) + z = _transform_z(T, z) + L = T.(L isa Number ? (L, L) : L) + + bandlimit_border = real(CT).(bandlimit_border) + L_new = padding ? pad_factor .* L : L + + # applies zero padding + if ndims(field) == 2 + pad_factor2 = pad_factor + elseif ndims(field) == 3 + pad_factor2 = (pad_factor * size(field, 1), pad_factor * size(field, 2), size(field, 3)) + end + field_new = padding ? pad(field, pad_factor2) : field + + # helpful propagation variables + (; k, f_x, f_y) = Zygote.@ignore _propagation_variables(field_new, λ, L_new) + + # transfer function kernel of angular spectrum + H = exp.(1im .* k .* abs.(z) .* sqrt.(CT(1) .- abs2.(f_x .* λ) .- abs2.(f_y .* λ))) + + # take complex conjugate, for negative zs + H = real.(H) .+ sign.(z) .* 1im .* imag(H) + + # bandlimit according to Matsushima + # as addition we introduce a smooth bandlimit with a Hann window + # and fuzzy logic + Δu = 1 ./ L_new + u_limit = Zygote.@ignore 1 ./ (sqrt.((2 .* Δu .* z).^2 .+ 1) .* λ) + + W = let + if bandlimit + # bandlimit filter + u_limit1 = u_limit isa AbstractArray ? u_limit[1:1, ..] : u_limit[1] + u_limit2 = u_limit isa AbstractArray ? u_limit[2:2, ..] : u_limit[2] + + W = .*(hann.(scale.(abs2.(f_y) ./ u_limit1 .^2 .+ abs2.(f_x) * λ^2, + bandlimit_border[1], bandlimit_border[2])), + hann.(scale.(abs2.(f_x) ./ u_limit2 .^2 .+ abs2.(f_y) * λ^2, + bandlimit_border[1], bandlimit_border[2]))) + else + # use an array here too, to avoid type instabilities + W = similar(field, real(eltype(field)), size(f_y, 1), size(f_x, 1)) + W .= 1 + end + end + + return (;field_new, H, W, fftdims) +end + + +""" + angular_spectrum(field, z, λ, L; kwargs...) + +Returns the electrical field with physical length `L` and wavelength `λ` propagated with the angular spectrum +method of plane waves (AS) by the propagation distance `z`. + +This method is efficient but to avoid recalculating some arrays (such as the phase kernel), see [`AngularSpectrum`](@ref). + +# Arguments +* `field`: Input field +* `z`: propagation distance. Can be a single number or a vector of `z`s (Or `CuVector`). In this case the returning array has one dimension more. +* `λ`: wavelength of field +* `L`: field size (can be a scalar or a tuple) indicating field size + + +# Keyword Arguments +* `padding=true`: applies padding to avoid convolution wraparound +* `pad_factor=2`: padding of 2. Larger numbers are not recommended since they don't provide better results. +* `bandlimit=true`: applies the bandlimit to avoid circular wraparound due to undersampling + of the complex propagation kernel [1] +* `bandlimit_border=(0.8, 1)`: applies a smooth bandlimit cut-off instead of hard-edge. + + +# Examples +```jldoctest +julia> field = zeros(ComplexF32, (4,4)); field[3,3] = 1 +1 +julia> as, t = angular_spectrum(field, 100e-9, 632e-9, 10e-6) +(ComplexF32[6.519258f-8 - 3.5390258f-7im -2.5097302f-7 + 1.1966712f-6im 0.00041754358f0 - 0.00027736276f0im -2.5097302f-7 + 1.1966712f-6im; -2.5136978f-7 + 1.1966767f-6im 8.540863f-7 - 4.0512546f-6im -0.001423778f0 + 0.00093840604f0im 8.540863f-7 - 4.0512546f-6im; 0.00041754544f0 - 0.0002773702f0im -0.0014237723f0 + 0.0009384034f0im 0.5497784f0 + 0.8353027f0im -0.0014237723f0 + 0.0009384034f0im; -2.5136978f-7 + 1.1966767f-6im 8.540863f-7 - 4.0512546f-6im -0.001423778f0 + 0.00093840604f0im 8.540863f-7 - 4.0512546f-6im], (H = ComplexF32[0.54519475f0 + 0.8383094f0im 0.5456109f0 + 0.8380386f0im … 0.54685986f0 + 0.8372242f0im 0.5456109f0 + 0.8380386f0im; 0.5456109f0 + 0.8380386f0im 0.5460271f0 + 0.83776754f0im … 0.54727626f0 + 0.83695203f0im 0.5460271f0 + 0.83776754f0im; … ; 0.54685986f0 + 0.8372242f0im 0.54727626f0 + 0.83695203f0im … 0.548526f0 + 0.8361335f0im 0.54727626f0 + 0.83695203f0im; 0.5456109f0 + 0.8380386f0im 0.5460271f0 + 0.83776754f0im … 0.54727626f0 + 0.83695203f0im 0.5460271f0 + 0.83776754f0im], L = 1.0e-5)) +``` + + +# References +* Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673. +""" +function angular_spectrum(field::AbstractArray{CT, 2}, z, λ, L; + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1.0)) where {CT<:Complex} + + @assert size(field, 1) == size(field, 2) "input field needs to be quadradically shaped and not $(size(field, 1)), $(size(field, 2))" + + (; field_new, H, W, fftdims) = _prepare_angular_spectrum(field, z, λ, L; padding, + pad_factor, bandlimit, bandlimit_border) + + # propagate field + field_out = fftshift(ifft(fft(ifftshift(field_new, fftdims), fftdims) .* H .* W, fftdims), fftdims) + + field_out_cropped = padding ? crop_center(field_out, size(field)) : field_out + + # return final field and some other variables + return field_out_cropped, (; H, L) +end + + +angular_spectrum(field::AbstractArray, z, λ, L; kwargs...) = throw(ArgumentError("Provided field needs to have a complex elementype")) + + + # highly optimized version with pre-planning +struct AngularSpectrum3{A, T, P} + HW::A + buffer::A + buffer2::A + L::T + p::P + padding::Bool + pad_factor::Int +end + +""" + AngularSpectrum(field, z, λ, L; kwargs...) + +Returns a function for efficient reuse of pre-calculated kernels. + +See [`angular_spectrum`](@ref) for the full documentation. + +# Example +```jldoctest +julia> field = zeros(ComplexF32, (4,4)); field[3,3] = 1 +1 + +julia> as, t = AngularSpectrum(field, 100e-9, 632e-9, 10e-6); + +julia> as(field) +(ComplexF32[6.519258f-8 - 3.5390258f-7im -2.5097302f-7 + 1.1966712f-6im 0.00041754358f0 - 0.00027736276f0im -2.5097302f-7 + 1.1966712f-6im; -2.5136978f-7 + 1.1966767f-6im 8.540863f-7 - 4.0512546f-6im -0.001423778f0 + 0.00093840604f0im 8.540863f-7 - 4.0512546f-6im; 0.00041754544f0 - 0.0002773702f0im -0.0014237723f0 + 0.0009384034f0im 0.5497784f0 + 0.8353027f0im -0.0014237723f0 + 0.0009384034f0im; -2.5136978f-7 + 1.1966767f-6im 8.540863f-7 - 4.0512546f-6im -0.001423778f0 + 0.00093840604f0im 8.540863f-7 - 4.0512546f-6im], (L = 1.0e-5,)) +``` +""" +function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1)) where {CT, N} + + (; field_new, H, W, fftdims) = _prepare_angular_spectrum(field, z, λ, L; padding, + pad_factor, bandlimit, bandlimit_border) + + + if z isa AbstractVector + buffer2 = similar(field, complex(eltype(field_new)), (size(field_new, 1), size(field_new, 2), length(z))) + buffer = copy(buffer2) + else + buffer2 = similar(field, complex(eltype(field_new)), (size(field_new, 1), size(field_new, 2))) + buffer = copy(buffer2) + end + + p = plan_fft!(buffer, (1, 2)) + H .= H .* W + HW = H + + return AngularSpectrum3{typeof(H), typeof(L), typeof(p)}(HW, buffer, buffer2, L, p, padding, pad_factor), L + end + +""" + (as:AngularSpectrum3)(field) + +Uses the struct to efficiently store some pre-calculated objects. +Propagate the field. +""" +function (as::AngularSpectrum3)(field) + fill!(as.buffer2, 0) + field_new = set_center!(as.buffer2, field, broadcast=true) + field_imd = as.p * ifftshift!(as.buffer, field_new, (1, 2)) + field_imd .*= as.HW + field_out = fftshift!(as.buffer2, inv(as.p) * field_imd, (1, 2)) + field_out_cropped = as.padding ? crop_center(field_out, size(field), return_view=true) : field_out + return field_out_cropped, (; as.L) +end + + +function ChainRulesCore.rrule(as::AngularSpectrum3, field) + field_and_tuple = as(field) + function as_pullback(ȳ) + f̄ = NoTangent() + # i tried to fix this once, but we somehow the Tangent type is missing the dimensionality + # which we need for set_center! and crop_center + y2 = ȳ.backing[1] + + fill!(as.buffer2, 0) + field_new = as.padding ? set_center!(as.buffer2, y2, broadcast=true) : y2 + field_imd = as.p * ifftshift!(as.buffer, field_new, (1, 2)) + field_imd .*= conj.(as.HW) + field_out = fftshift!(as.buffer2, inv(as.p) * field_imd, (1, 2)) + # that means z is a vector and we do plane to volume propagation + if size(as.buffer, 3) > 1 + sum!(view(as.buffer, :, :, 1), field_out) + field_out_cropped = as.padding ? crop_center(view(as.buffer, :, :, 1), size(field), return_view=true) : view(as.buffer, :, :, 1) + else + field_out_cropped = as.padding ? crop_center(field_out, size(field), return_view=true) : field_out + end + return f̄, field_out_cropped + end + return field_and_tuple, as_pullback +end diff --git a/src/beams.jl b/src/beams.jl new file mode 100644 index 0000000..21f2cb7 --- /dev/null +++ b/src/beams.jl @@ -0,0 +1,18 @@ +export gauss_beam + + +gauss_beam(y, x, z, λ, w_0) = throw(ArgumentError("All arguments need to have the same datatype (Float32, Float64, ...).")) + +function gauss_beam(y::T, x::T, z::T, λ::T, w_0::T) where T + k = π / λ * 2 + z_R = π * w_0^2 / λ + r² = x ^ 2 + y ^ 2 + # don't put exp(i * k * z) into the same exp, it causes some strange wraps + return w_0 / gauss_w(z, z_R, w_0) * exp(-r² / gauss_w(z, z_R, w_0)^2) * + exp.(1im * k * z) * + exp(1im * (k * r² / 2 / gauss_R(z, z_R) - gauss_ψ(z, z_R))) +end + +@inline gauss_R(z, z_R) = z * (1 + (z_R /z)^2) +@inline gauss_ψ(z, z_R) = atan(z, z_R) +@inline gauss_w(z, z_R, w_0) = w_0 * sqrt(1 + (z / z_R)^2) diff --git a/src/conv.jl b/src/conv.jl new file mode 100644 index 0000000..fdbf9b4 --- /dev/null +++ b/src/conv.jl @@ -0,0 +1,7 @@ + + +function conv(u_in, v, dims=(1,2)) + u = pad(u_in, size(v)[1:2]) + + return crop_center(ifft(fft(u, dims) .* fft(v, dims), dims), size(u_in)) +end diff --git a/src/fraunhofer.jl b/src/fraunhofer.jl new file mode 100644 index 0000000..d105ba5 --- /dev/null +++ b/src/fraunhofer.jl @@ -0,0 +1,146 @@ +export fraunhofer +export Fraunhofer + + +""" + fraunhofer(field, z, λ, L) + +Returns the electrical field with physical length `L` and wavelength `λ` +propagated with the Fraunhofer propagation (a single FFT) by the propagation distance `z`. +This is based on a far field approximation `z >> λ` + +This method is efficient but to save memory and avoiding recalculating some arrays (such as the phase kernel), see [`Fraunhofer`](@ref). + +# Arguments +* `field`: Input field +* `z`: propagation distance +* `λ`: wavelength of field +* `L`: field size indicating field size + +# Keyword Arguments +* `skip_final_phase=true` skip the final phase which is multiplied to the propagated field at the end + + +# Example +```jldoctest +julia> field = zeros(ComplexF32, (256,256)); field[130,130] = 1; + +julia> res, t = fraunhofer(field, 4f-3, 632f-9, 100f-6) +(ComplexF64[0.00390625 + 0.0im 0.003905073506757617 - 9.586417581886053e-5im … 0.003901544725522399 + 0.0001916706096380949im 0.003905073506757617 + 9.58640594035387e-5im; 0.003905073506757617 - 9.586417581886053e-5im 0.003901544725522399 - 0.0001916706096380949im … 0.003905073506757617 + 9.58640594035387e-5im 0.00390625 - 1.1641532182693481e-10im; … ; 0.003901544725522399 + 0.00019167049322277308im 0.003905073506757617 + 9.586406667949632e-5im … 0.003887440310791135 + 0.0003828795161098242im 0.0038956659846007824 + 0.00028736155945807695im; 0.003905073506757617 + 9.586417581886053e-5im 0.00390625 - 1.5902765215791703e-12im … 0.0038956659846007824 + 0.0002873614430427551im 0.003901544725522399 + 0.0001916706096380949im], (L = 0.006471681f0,)) + +julia> t.L / 100f-6 +64.71681f0 + +julia> 4f-3 * 632f-9 * 256 / (100f-6)^2 +64.71681f0 +``` + +""" +function fraunhofer(U, z, λ, L; skip_final_phase=true) + @assert size(U, 1) == size(U, 2) + L_new = λ * z / L * size(U, 1) + Ns = size(U)[1:2] + + p = Zygote.@ignore plan_fft(U, (1,2)) + + if skip_final_phase + out = fftshift(p * ifftshift(U)) ./ √(size(U, 1) * size(U, 2)) + else + k = eltype(U)(2π) / λ + # output coordinates + y = similar(U, real(eltype(U)), (Ns[1], 1)) + Zygote.@ignore y .= (fftpos(L, Ns[1], CenterFT)) + x = similar(U, real(eltype(U)), (1, Ns[2])) + Zygote.@ignore x .= (fftpos(L, Ns[2], CenterFT))' + phasefactor = (-1im) .* exp.(1im * k / (2 * z) .* (x.^2 .+ y.^2)) + out = phasefactor .* fftshift(p * ifftshift(U)) ./ √(size(U, 1) * size(U, 2)) + end + + return out, (; L=L_new) +end + +""" + Fraunhofer(U, z, λ, L; skip_final_phase=true) + + +This returns a function for efficient reuse of pre-calculated kernels. +See [`fraunhofer`](@ref) for the full documentation. + + +# Example +```jldoctest +julia> field = zeros(ComplexF32, (256,256)); field[130,130] = 1; + +julia> f, t = Fraunhofer(field, 4f-3, 632f-9, 100f-6); + +julia> f(field) +(ComplexF32[0.00390625f0 + 0.0f0im 0.0039050735f0 - 9.5864176f-5im … 0.0039015447f0 + 0.00019167061f0im 0.0039050735f0 + 9.586406f-5im; 0.0039050735f0 - 9.5864176f-5im 0.0039015447f0 - 0.00019167061f0im … 0.0039050735f0 + 9.586406f-5im 0.00390625f0 - 1.1641532f-10im; … ; 0.0039015447f0 + 0.0001916705f0im 0.0039050735f0 + 9.586407f-5im … 0.0038874403f0 + 0.00038287952f0im 0.003895666f0 + 0.00028736156f0im; 0.0039050735f0 + 9.5864176f-5im 0.00390625f0 - 1.5902765f-12im … 0.003895666f0 + 0.00028736144f0im 0.0039015447f0 + 0.00019167061f0im], (L = 0.006471681f0,)) + +julia> t.L / 100f-6 +64.71681f0 + +julia> 4f-3 * 632f-9 * 256 / (100f-6)^2 +64.71681f0 +``` +""" +function Fraunhofer(U, z, λ, L; skip_final_phase=true) + @assert size(U, 1) == size(U, 2) + L_new = λ * z / L * size(U, 1) + Ns = size(U)[1:2] + + if skip_final_phase + phasefactor = nothing + else + k = eltype(U)(2π) / λ + # output coordinates + y = similar(U, real(eltype(U)), (Ns[1], 1)) + y .= (fftpos(L, Ns[1], CenterFT)) + x = similar(U, real(eltype(U)), (1, Ns[2])) + x .= (fftpos(L, Ns[2], CenterFT))' + phasefactor = (-1im) .* exp.(1im * k / (2 * z) .* (x.^2 .+ y.^2)) + end + + buffer = zero.(U) + + FFTplan = plan_fft!(buffer, (1,2)) + + return FraunhoferOp{typeof(L), typeof(buffer), typeof(phasefactor), + typeof(FFTplan)}(buffer, phasefactor, L_new, FFTplan), (; L=L_new) +end + +struct FraunhoferOp{T, B, PF, P} + buffer::B + phasefactor::PF + L::T + FFTplan::P +end + +function (fraunhofer::FraunhoferOp)(field) + buffer = fraunhofer.buffer + ifftshift!(buffer, field) + fraunhofer.FFTplan * buffer + buffer ./= √(size(field, 1) * size(field, 2)) + out = fftshift(buffer) + if !isnothing(fraunhofer.phasefactor) + out .*= fraunhofer.phasefactor + end + return out, (; fraunhofer.L) +end + + +function ChainRulesCore.rrule(f::FraunhoferOp, U) + field_and_tuple = f(U) + + function f_pullback(ȳ) + buffer = f.buffer + y2 = ȳ.backing[1] + if !isnothing(f.phasefactor) + y2 = y2 .* conj.((f.phasefactor)) + end + ifftshift!(buffer, y2) + buffer .*= √(size(U, 1) * size(U, 2)) + res = fftshift(inv(f.FFTplan) * buffer) + return NoTangent(), res + end + return field_and_tuple, f_pullback +end diff --git a/src/propagation.jl b/src/propagation.jl new file mode 100644 index 0000000..8ad5916 --- /dev/null +++ b/src/propagation.jl @@ -0,0 +1,64 @@ +""" + _propagation_variables(field, λ, L) + +Internal method to create variables we need for propagation such as frequencies in Fourier space, etc.. +""" +function _propagation_variables(field::AbstractArray{T, M}, λ::Number, L::NTuple{2, T2}) where {T, T2, M} + # wave number + k = real(T)(2π) / λ + # number of samples + Ns = size(field)[1:2] + # sample spacing + dx = L ./ Ns + + # frequency spacing + df = 1 ./ L + + # total size in frequency space + Lf = Ns .* df + + # frequencies centered around first entry + # 1D vectors each + f_y = similar(field, real(eltype(field)), (Ns[1], 1)) + f_y .= fftfreq(Ns[1], Lf[1]) + f_x = similar(field, real(eltype(field)), (1, Ns[2])) + f_x .= fftfreq(Ns[2], Lf[2])' + + # y and x positions in real space, use correct spacing -> fftpos + y = similar(field, real(eltype(field)), (Ns[1], 1)) + y .= (fftpos(L[1], Ns[1], CenterFT)) + x = similar(field, real(eltype(field)), (1, Ns[2])) + x .= (fftpos(L[2], Ns[2], CenterFT))' + + return (; k, dx, df, f_x, f_y, x, y) +end + + +function _propagation_variables(field::AbstractArray{T, M}, λ, L::Number) where {T, M} + _propagation_variables(field, λ, (L, L)) +end + + + +function _real_space_kernel(field::AbstractArray{T, N}, z, λ, L) where {T, N} + L = L isa Number ? (L, L) : L + # wave number + k = T(2π) / λ + # number of samples + Ns = size(field)[1:2] + + # y and x positions in real space, use correct spacing -> fftpos + y = similar(field, real(eltype(field)), (Ns[1], 1)) + y .= (fftpos(L[1], Ns[1], CenterFT)) + x = similar(field, real(eltype(field)), (1, Ns[2])) + x .= (fftpos(L[2], Ns[2], CenterFT))' + + + #y = ifftshift(y, (1,)) + #x = ifftshift(x, (2,)) + + r = sqrt.(x.^2 .+ y.^2 .+ z .^2) + kernel = T(2 ./ π) .* z ./ r .* (1 ./ r .- 1im * k) .* exp.(1im .* k .* r) ./ r + kernel .*= CuArray(rr2(size(field)) .< 60 .^2) + return kernel +end diff --git a/src/scalable_angular_spectrum.jl b/src/scalable_angular_spectrum.jl new file mode 100644 index 0000000..14a6e5e --- /dev/null +++ b/src/scalable_angular_spectrum.jl @@ -0,0 +1,241 @@ +export ScalableAngularSpectrum + + # highly optimized version with pre-planning +struct ScalableAngularSpectrum{AP, AP2, AB, T, P} + ΔH::AP + H₁::AP + H₂::AP2 # could be nothing! + buffer::AB + buffer2::AB + L::T + FFTplan::P + pad_factor::Int +end + +""" + _SAS_propagation_variables(field, z, λ, L) + +Internal method to create variables we need for propagation such as frequencies in Fourier space, etc.. +""" +function _SAS_propagation_variables(field::AbstractArray{T, M}, z, λ, L) where {T, M} + @assert size(field, 1) == size(field, 2) "Quadratic fields only working currently" + + # wave number + k = T(2π) / λ + # number of samples + N = size(field, 1) + # sample spacing + dx = L / N + # frequency spacing + df = 1 / L + # total size in frequency space + Lf = N * df + + # frequencies centered around first entry + # 1D vectors each + f_y = similar(field, real(eltype(field)), (N,)) + f_y .= fftfreq(N, Lf) + f_x = f_y' + + # y and x positions in real space + #y = ifftshift(range(-L/2, L/2, length=N)) + y = similar(field, real(eltype(field)), (N,)) + y .= fftpos(L, N, CenterFT) + y = ifftshift(y) + x = y' + + return (; k, dx, df, f_x, f_y, x, y) +end + + +# helper function for smooth window +function find_width_window(ineq_1D::AbstractVector, bandlimit_border) + ineq_1D = Array(ineq_1D) + bs = ineq_1D .≤ 1 + ind_x_first = findfirst(bs) + ind_x_last = findlast(bs) + + if isnothing(ind_x_first) || isnothing(ind_x_last) + return (1,1) + end + + diff_b = round(Int, 0.5 * (1 - bandlimit_border[1]) * (Tuple(ind_x_last)[1] - Tuple(ind_x_first)[1])) + diff_e = round(Int, 0.5 * (1 - bandlimit_border[2]) * (Tuple(ind_x_last)[1] - + Tuple(ind_x_first)[1])) + + ineq_v_b = ineq_1D[ind_x_first + diff_b] + ineq_v_e = ineq_1D[ind_x_first + diff_e] + + return ineq_v_b, ineq_v_e +end + +""" + ScalableAngularSpectrum(field, z, λ, L; kwargs...) + +Returns the electrical field with physical length `L` and wavelength `λ` propagated with the scalable angular spectrum method of plane waves (SAS) by the propagation distance `z`. + + +# Arguments +* `field`: Input field +* `z`: propagation distance +* `λ`: wavelength of field +* `L`: field size (can be a scalar or a tuple) indicating field size + + +# Keyword Arguments +* `skip_final_phase=true`: avoid multiplying with final phase. This phase is also undersampled. +* `bandlimit_border=(0.8, 1)`: apply soft bandlimit instead of hard cutoff + + +# Example +```jldoctest +julia> field = zeros(ComplexF32, (256,256)); field[130,130] = 1; + +julia> sas, t = ScalableAngularSpectrum(field, 10f-3, 633f-9, 500f-6); + +julia> res, t2 = sas(field); + +julia> t2 +(L = 0.00162048f0,) + +julia> t +(L = 0.00162048f0,) + +julia> t2.L / 500f-6 # calculate magnification +3.24096f0 + +julia> 10f-3 * 633f-9 * 256 / 500f-6^2 / 2 # equation from paper +3.2409596f0 +``` + +# References +* [Rainer Heintzmann, Lars Loetgering, and Felix Wechsler, "Scalable angular spectrum propagation," Optica 10, 1407-1416 (2023)](https://opg.optica.org/optica/viewmedia.cfm?uri=optica-10-11-1407&html=true) +""" +function ScalableAngularSpectrum(ψ₀::AbstractArray{T}, z, λ, L ; + skip_final_phase=true, bandlimit_soft_px=20, + bandlimit_border=(0.8, 1)) where {T} + @assert bandlimit_soft_px ≥ 0 "bandlimit_soft_px must be ≥ 0" + @assert size(ψ₀, 1) == size(ψ₀, 2) "Restricted to auadratic fields." + + pad_factor = 2 + set_pad_zero = false + + N = size(ψ₀, 1) + z_limit = (- 4 * L * sqrt(8*L^2 / N^2 + λ^2) * sqrt(L^2 * inv(8 * L^2 + N^2 * λ^2)) / (λ * (-1+2 * sqrt(2) * sqrt(L^2 * inv(8 * L^2 + N^2 * λ^2))))) + + # vignetting limit + z > z_limit && @warn "Propagated field might be affected by vignetting" + L_new = pad_factor * L + + # applies zero padding + ψ_p = select_region(ψ₀, new_size=size(ψ₀) .* pad_factor) + k, dx, df, f_x, f_y, x, y = _SAS_propagation_variables(ψ_p, z, λ, L_new) + M = λ * z * N / L^2 / 2 + + # calculate anti_aliasing_filter for precompensation + cx = λ .* f_x + cy = λ .* f_y + tx = L_new / 2 / z .+ abs.(λ .* f_x) + ty = L_new / 2 / z .+ abs.(λ .* f_y) + + # smooth window function + smooth_f(x, α, β) = hann(scale(x, α, β)) + # find boundary for soft hann + ineq_x = fftshift(cx[1, :].^2 .* (1 .+ tx[1, :].^2) ./ tx[1, :].^2 .+ cy[1, :].^2) + limits = find_width_window(ineq_x, bandlimit_border) + + # bandlimit filter for precompensation + W = .*(smooth_f.(cx.^2 .* (1 .+ tx.^2) ./ tx.^2 .+ cy.^2, limits...), + smooth_f.(cy.^2 .* (1 .+ ty.^2) ./ ty.^2 .+ cx.^2, limits...)) + + # ΔH is the core part of Fresnel and AS + H_AS = sqrt.(0im .+ 1 .- abs2.(f_x .* λ) .- abs2.(f_y .* λ)) + H_Fr = 1 .- abs2.(f_x .* λ) / 2 .- abs2.(f_y .* λ) / 2 + # take the difference here, key part of the ScaledAS + ΔH = W .* exp.(1im .* k .* z .* (H_AS .- H_Fr)) + + + # new sample coordinates + dq = λ * z / L_new + Q = dq * N * pad_factor + q_y = similar(ψ_p, pad_factor * N) + # fftpos generates coordinates from -L/2 to L/2 but excluding the last + # final bit + q_y .= fftpos(dq * pad_factor * N, pad_factor * N, CenterFT) + q_y = ifftshift(q_y) + q_x = q_y' + + # calculate phases of Fresnel + H₁ = exp.(1im .* k ./ (2 .* z) .* (x .^ 2 .+ y .^ 2)) + + # skip final phase because often undersampled + if skip_final_phase + H₂ = nothing + else + H₂ = (exp.(1im .* k .* z) .* + exp.(1im .* k ./ (2 .* z) .* (q_x .^ 2 .+ q_y .^2))) + end + + FFTplan = plan_fft!(ψ_p, (1,2)) + + buffer = similar(ψ_p) + buffer2 = similar(buffer) + return ScalableAngularSpectrum{typeof(ΔH), typeof(H₂), typeof(buffer), + real(T), typeof(FFTplan)}( + ΔH, H₁, H₂, buffer, buffer2, Q / 2, FFTplan, pad_factor), (; L=Q / 2) +end + + +function (sas::ScalableAngularSpectrum)(ψ::AbstractArray{T}; return_view=false) where T + p = sas.FFTplan + fill!(sas.buffer2, 0) + ψ_p = set_center!(sas.buffer2, ψ) + ψ_p = ifftshift!(sas.buffer, ψ_p) + ψ_p_f = p * ψ_p + ψ_p_f .*= sas.ΔH + ψ_precomp = inv(p) * ψ_p_f + ψ_precomp .*= sas.H₁ + ψ_p_final = p * ψ_precomp + ψ_p_final .*= 1 / (1im * sqrt(T(size(ψ_precomp, 1) * size(ψ_precomp, 2)))) + + if !isnothing(sas.H₂) + ψ_p_final .*= sas.H₂ + end + fftshift!(sas.buffer2, ψ_p_final) + + ψ_final = crop_center(sas.buffer2, size(ψ); return_view) + return ψ_final, (; sas.L) +end + + + +function ChainRulesCore.rrule(sas::ScalableAngularSpectrum, ψ::AbstractArray{T}) where T + field_and_tuple = sas(ψ) + function sas_pullback(ȳ) + p = sas.FFTplan + y2 = ȳ.backing[1] + + fill!(sas.buffer2, 0) + + set_center!(sas.buffer2, y2) + ifftshift!(sas.buffer, sas.buffer2) + + if !isnothing(sas.H₂) + sas.buffer .*= conj.(sas.H₂) + end + sas.buffer .*= conj.(1 / (1im * sqrt(T(size(sas.buffer2, 1) * size(sas.buffer2, 2))))) + sas.buffer .*= (T(size(sas.buffer2, 1) * size(sas.buffer2, 2))) + inv(p) * sas.buffer + sas.buffer .*= conj.(sas.H₁) + (p) * sas.buffer + sas.buffer .*= conj.(sas.ΔH) + inv(p) * sas.buffer + fftshift!(sas.buffer2, sas.buffer) + + final = crop_center(sas.buffer2, size(ψ)) + + + return NoTangent(), final + end + return field_and_tuple, sas_pullback +end diff --git a/src/shifted_angular_spectrum.jl b/src/shifted_angular_spectrum.jl new file mode 100644 index 0000000..f58f637 --- /dev/null +++ b/src/shifted_angular_spectrum.jl @@ -0,0 +1,216 @@ +export ShiftedAngularSpectrum + + +function _prepare_shifted_angular_spectrum(field::AbstractArray{CT}, z, λ, L, α; + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1.0)) where {CT<:Complex} + + fftdims = (1,2) + T = real(CT) + λ = T(λ) + z = _transform_z(T, z) + L = T.(L isa Number ? (L, L) : L) + + sxy = .+ sin.(α) + txy = .+ tan.(α) + + bandlimit_border = real(CT).(bandlimit_border) + L_new = padding ? pad_factor .* L : L + + # applies zero padding + if ndims(field) == 2 + pad_factor2 = pad_factor + elseif ndims(field) == 3 + pad_factor2 = (pad_factor * size(field, 1), pad_factor * size(field, 2), size(field, 3)) + end + field_new = padding ? pad(field, pad_factor2) : field + + # helpful propagation variables + (; k, f_x, f_y, x, y) = Zygote.@ignore _propagation_variables(field_new, λ, L_new) + + + H = exp.(1im .* k .* z .* (sqrt.(CT(1) .- abs2.(f_x .* λ .+ sxy[2]) .- abs2.(f_y .* λ .+ sxy[1])) + .+ λ .* (txy[2] .* f_x .+ txy[1] .* f_y))) + + # bandlimit according to Matsushima + W = let + if bandlimit + # bandlimit filter + χ = max.(0, 1 / λ^2 .- abs2.(f_x .+ sxy[2] ./ λ) .- abs2.(f_y .+ sxy[1] ./ λ)) + + Ωx = z * (txy[2] .- (f_x .+ sxy[2] ./ λ) ./ (sqrt.(χ))) + Ωy = z * (txy[1] .- (f_y .+ sxy[1] ./ λ) ./ (sqrt.(χ))) + + W = ( (1 / L_new[2]) .<= abs.(1 ./ 2 ./ Ωx)) .* ( (1 / L_new[1]) .<= abs.(1 ./ 2 ./ Ωy)) + else + # use an array here too, to avoid type instabilities + W = similar(field, real(eltype(field)), size(f_y, 1), size(f_x, 2)) + W .= 1 + end + end + + shift = txy .* z + + ya = similar(field_new, real(eltype(field)), (size(field_new, 1), 1)) + Zygote.@ignore ya .= (fftpos(L_new[1], size(field_new, 1), CenterFT)) .+ shift[1] + xa = similar(field_new, real(eltype(field)), (1, size(field_new, 2))) + Zygote.@ignore xa .= (fftpos(L_new[2], size(field_new, 2), CenterFT))' .+ shift[2] + + ramp_before = ifftshift(exp.(1im .* 2 .* T(π) ./ λ .* (sxy[2] .* x .+ sxy[1] .* y)), (1,2)) + ramp_after = ifftshift(exp.(1im .* 2 .* T(π) ./ λ .* (sxy[2] .* xa .+ sxy[1] .* ya)), (1,2)) + return (;field_new, H, W, fftdims, ramp_before, ramp_after) +end + + +function shifted_angular_spectrum(field::AbstractArray{CT, 2}, z, λ, L, α; + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1.0)) where {CT<:Complex} + + @assert size(field, 1) == size(field, 2) "input field needs to be quadradically shaped and not $(size(field, 1)), $(size(field, 2))" + + (; field_new, H, W, fftdims, ramp_before, ramp_after) = _prepare_shifted_angular_spectrum(field, z, λ, L, real(CT).(α); padding, + pad_factor, bandlimit, bandlimit_border) + + # propagate field + field_new_is = ifftshift(field_new, fftdims) ./ (ramp_before) + field_out = fftshift(ramp_after .* ifft(fft(field_new_is, fftdims) .* H .* W, fftdims), fftdims) + field_out_cropped = padding ? crop_center(field_out, size(field)) : field_out + shift = z .* tan.(α) ./ L + # return final field and some other variables + return field_out_cropped, (; L, shift) +end + + + # highly optimized version with pre-planning +struct ShiftedAngularSpectrum{A, T, T2, P, R} + HW::A + buffer::A + buffer2::A + L::T + shift::T2 + p::P + padding::Bool + pad_factor::Int + ramp_before::R + ramp_after::R +end + + +""" + ShiftedAngularSpectrum(field, z, λ, L, α; kwargs...) + +Returns a method to propagate the electrical field with physical length `L` and wavelength `λ` with the shifted angular spectrum +method of plane waves (AS) by the propagation distance `z`. +`α` should be a tuple containing the offset angles with respect to the optical axis. + + +# Arguments +* `field`: Input field +* `z`: propagation distance. Can be a single number or a vector of `z`s (Or `CuVector`). In this case the returning array has one dimension more. +* `λ`: wavelength of field +* `L`: field size (can be a scalar or a tuple) indicating field size +* `α` is the tuple of shift angles with respect to the optical axis. + + +# Keyword Arguments +* `pad_factor=2`: padding of 2. Larger numbers are not recommended since they don't provide better results. +* `bandlimit=true`: applies the bandlimit to avoid circular wraparound due to undersampling + of the complex propagation kernel [1] +* `extract_ramp=true`: divides the field by phase ramp `exp.(1im * 2π / λ * (sin(α[2]) .* x .+ sin(α[1]) .* y))` and multiplies after + propagation the ramp (with new real space coordinates) back to the field + +# Examples +```jldoctest +julia> field = zeros(ComplexF32, (4,4)); field[3,3] = 1 + +julia> AS, _ = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (deg2rad(10), 0)); + +julia> AS(field) +(ComplexF32[1.5269792f-5 + 1.7594219f-5im -3.996831f-5 - 7.624799f-5im -0.0047351345f0 + 0.002100923f0im -3.996831f-5 - 7.624799f-5im; -8.294997f-5 - 1.8230454f-5im 0.00028230582f0 + 8.1745195f-5im 0.0051693693f0 - 0.016958509f0im 0.00028230582f0 + 8.1745195f-5im; 0.0029884572f0 + 0.0040671355f0im -0.009686601f0 - 0.014245203f0im -0.82990384f0 + 0.5566719f0im -0.009686601f0 - 0.014245203f0im; -1.4191573f-7 + 9.41665f-5im 2.111472f-5 - 0.00031620878f0im -0.017670793f0 - 0.0014212304f0im 2.111472f-5 - 0.00031620878f0im], (L = 0.0001, shift = (0.17632698070846498, 0.0))) +``` + +# References +* Matsushima, Kyoji. "Shifted angular spectrum method for off-axis numerical propagation." Optics Express 18.17 (2010): 18453-18463. +""" +function ShiftedAngularSpectrum(field::AbstractArray{CT, N}, z::Number, λ, L, α; + extract_ramp=true, + padding=true, pad_factor=2, + bandlimit=true, + bandlimit_border=(0.8, 1)) where {CT, N} + + (; field_new, H, W, fftdims, ramp_before, ramp_after) = + _prepare_shifted_angular_spectrum(field, z, λ, L, real(CT).(α); padding, + pad_factor, bandlimit, bandlimit_border) + + + buffer2 = similar(field, complex(eltype(field_new)), (size(field_new, 1), size(field_new, 2))) + buffer = copy(buffer2) + + p = plan_fft!(buffer, (1, 2)) + H .= H .* W + HW = H + shift = z .* tan.(α) ./ L + if !extract_ramp + ramp_before = nothing + ramp_after = nothing + end + + + return ShiftedAngularSpectrum{typeof(H), typeof(L), typeof(shift), typeof(p), typeof(ramp_before)}(HW, + buffer, buffer2, L, shift, p, padding, pad_factor, ramp_before, ramp_after), (;L, shift) + end + + + +""" + (shifted_as::ShiftedAngularSpectrum)(field) + +Uses the struct to efficiently store some pre-calculated objects. +Propagate the field. +""" +function (as::ShiftedAngularSpectrum{A, T, T2, P, R})(field) where {A,T,T2,P,R} + fill!(as.buffer2, 0) + field_new = set_center!(as.buffer2, field, broadcast=true) + ifftshift!(as.buffer, field_new, (1, 2)) + if !(R === Nothing) + as.buffer ./= as.ramp_before + end + field_imd = as.p * as.buffer + field_imd .*= as.HW + field_imd = inv(as.p) * field_imd + if !(R === Nothing) + field_imd .*= as.ramp_after + end + field_out = fftshift!(as.buffer2, field_imd, (1, 2)) + field_out_cropped = as.padding ? crop_center(field_out, size(field), return_view=true) : field_out + return field_out_cropped, (; as.L, as.shift) +end + + +function ChainRulesCore.rrule(as::ShiftedAngularSpectrum{A, T, T2, P, R}, field) where {A,T,T2,P,R} + field_and_tuple = as(field) + function as_pullback(ȳ) + f̄ = NoTangent() + y2 = ȳ.backing[1] + + fill!(as.buffer2, 0) + field_new = as.padding ? set_center!(as.buffer2, y2, broadcast=true) : y2 + ifftshift!(as.buffer, field_new, (1, 2)) + if !(R === Nothing) + as.buffer .*= conj.(as.ramp_after) + end + field_imd = as.p * as.buffer + field_imd .*= conj.(as.HW) + + field_imd = inv(as.p) * field_imd + if !(R === Nothing) + field_imd ./= conj.(as.ramp_before) + end + field_out = fftshift!(as.buffer2, field_imd, (1, 2)) + field_out_cropped = as.padding ? crop_center(field_out, size(field), return_view=true) : field_out + return f̄, field_out_cropped + end + return field_and_tuple, as_pullback +end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..48a92b2 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,237 @@ +export crop_center +export set_center! +export pad + + +hann(x) = sin(π * x / 2)^2 +scale(x, α, β) = clamp(1 - (x - α) / (β - α), 0, 1) +smooth_f(x, α, β) = hann(scale(x, α, β)) + +""" + get_indices_around_center(i_in, i_out) + +A function which provides two output indices `i1` and `i2` +where `i2 - i1 = i_out` + +The indices are chosen in a way that the range `i1:i2` +cuts the interval `1:i_in` in a way that the center frequency +stays at the center position. + +Works for both odd and even indices. +""" +function get_indices_around_center(i_in, i_out) + if (mod(i_in, 2) == 0 && mod(i_out, 2) == 0 + || mod(i_in, 2) == 1 && mod(i_out, 2) == 1) + x = (i_in - i_out) ÷ 2 + return 1 + x, i_in - x + elseif mod(i_in, 2) == 1 && mod(i_out, 2) == 0 + x = (i_in - 1 - i_out) ÷ 2 + return 1 + x, i_in - x - 1 + elseif mod(i_in, 2) == 0 && mod(i_out, 2) == 1 + x = (i_in - (i_out - 1)) ÷ 2 + return 1 + x, i_in - (x - 1) + end +end + +""" + pad(arr::AbstractArray{T, N}, new_size; value=zero(T)) + +Pads `arr` with `values` around such that the resulting array size is `new_size`. + +See also [`crop_center`](@ref), [`set_center!`](@ref). + +```jldoctest +julia> pad(ones((2,2)), 2) +4×4 Matrix{Float64}: + 0.0 0.0 0.0 0.0 + 0.0 1.0 1.0 0.0 + 0.0 1.0 1.0 0.0 + 0.0 0.0 0.0 0.0 +``` +""" +function pad(arr::AbstractArray{T, N}, new_size::NTuple; value=zero(T)) where {T, N} + n_arr = similar(arr, new_size) + fill!(n_arr, value) + # don't do broadcast. Does not fit with the meaning of pad + set_center!(n_arr, arr, broadcast=false) +end + + +""" + pad(arr::AbstractArray{T, N}, M; value=zero(T)) + +Pads `arr` with `values` around such that the resulting array size is `round(Int, M .* size(arr))`. +If `M isa Integer`, no rounding is needed. + + +!!! warn "Automatic Differentation" + If `M isa AbstractFloat`, automatic differentation might fail because of size failures + + +See also [`crop_center`](@ref), [`set_center!`](@ref). + +```jldoctest +julia> pad(ones((2,2)), 2) +4×4 Matrix{Float64}: + 0.0 0.0 0.0 0.0 + 0.0 1.0 1.0 0.0 + 0.0 1.0 1.0 0.0 + 0.0 0.0 0.0 0.0 +``` +""" +function pad(arr::AbstractArray{T, N}, M::Number; value=zero(T)) where {T, N} + pad(arr, round.(Int, size(arr) .* M), value=value) +end + + +""" + crop_center(arr, new_size) + +Extracts a center of an array. +`new_size_array` must be list of sizes indicating the output +size of each dimension. Centered means that a center frequency +stays at the center position. Works for even and uneven. +If `length(new_size_array) < length(ndims(arr))` the remaining dimensions +are untouched and copied. + + +See also [`pad`](@ref), [`set_center!`](@ref). + +# Examples +```jldoctest +julia> crop_center([1 2; 3 4], (1,)) +1×2 Matrix{Int64}: + 3 4 + +julia> crop_center([1 2; 3 4], (1,1)) +1×1 Matrix{Int64}: + 4 + +julia> crop_center([1 2; 3 4], (2,2)) +2×2 Matrix{Int64}: + 1 2 + 3 4 +``` +""" +function crop_center(arr, new_size::NTuple{N}; return_view=true) where {N} + M = ndims(arr) + @assert N ≤ M "Can't specify more dimensions than the array has." + @assert all(new_size .≤ size(arr)[1:N]) "You can't extract a larger array than the input array." + @assert Base.require_one_based_indexing(arr) "Require one based indexing arrays" + + out_indices = ntuple(i -> let + if i ≤ N + inds = get_indices_around_center(size(arr, i), + new_size[i]) + inds[1]:inds[2] + else + 1:size(arr, i) + end + end, + Val(M)) + + + # if return_view + # return @inbounds view(arr, out_indices...) + # else + return @inbounds arr[out_indices...] + # end +end + + +""" + set_center!(arr_large, arr_small; broadcast=false) + +Puts the `arr_small` central into `arr_large`. + +The convention, where the center is, is the same as the definition +as for FFT based centered. + +Function works both for even and uneven arrays. +See also [`crop_center`](@ref), [`pad`](@ref), [`set_center!`](@ref). + +# Keyword +* If `broadcast==false` then a lower dimensional `arr_small` will not be broadcasted +along the higher dimensions. +* If `broadcast==true` it will be broadcasted along higher dims. + + +See also [`crop_center`](@ref), [`pad`](@ref). + + +# Examples +```jldoctest +julia> set_center!([1, 1, 1, 1, 1, 1], [5, 5, 5]) +6-element Vector{Int64}: + 1 + 1 + 5 + 5 + 5 + 1 + +julia> set_center!([1, 1, 1, 1, 1, 1], [5, 5, 5, 5]) +6-element Vector{Int64}: + 1 + 5 + 5 + 5 + 5 + 1 + +julia> set_center!(ones((3,3)), [5]) +3×3 Matrix{Float64}: + 1.0 1.0 1.0 + 1.0 5.0 1.0 + 1.0 1.0 1.0 + +julia> set_center!(ones((3,3)), [5], broadcast=true) +3×3 Matrix{Float64}: + 1.0 1.0 1.0 + 5.0 5.0 5.0 + 1.0 1.0 1.0 +``` +""" +function set_center!(arr_large::AbstractArray{T, N}, arr_small::AbstractArray{T1, M}; + broadcast=false) where {T, T1, M, N} + @assert N ≥ M "Can't put a higher dimensional array in a lower dimensional one." + + if broadcast == false + inds = ntuple(i -> begin + a, b = get_indices_around_center(size(arr_large, i), size(arr_small, i)) + a:b + end, + Val(N)) + arr_large[inds..., ..] .= arr_small + else + inds = ntuple(i -> begin + a, b = get_indices_around_center(size(arr_large, i), size(arr_small, i)) + a:b + end, + Val(M)) + arr_large[inds..., ..] .= arr_small + end + + + return arr_large +end + +function ChainRulesCore.rrule(::typeof(crop_center), arr, new_size::NTuple{N, <:Int}) where N + y = crop_center(arr, new_size) + function crop_center_pullback(ȳ) + c̄ = pad(ȳ, size(arr), value=zero(eltype(arr))) + return NoTangent(), c̄, NoTangent() + end + return y, crop_center_pullback +end + + +function ChainRulesCore.rrule(::typeof(pad), arr, M::Union{NTuple, Number}; value=zero(eltype(arr))) + y = pad(arr, M; value=value) + function pad_pullback(ȳ) + @assert size(y) == size(ȳ) + p̄ = crop_center(ȳ, size(arr)) + return NoTangent(), p̄, NoTangent() + end + return y, pad_pullback +end diff --git a/test/angular_spectrum.jl b/test/angular_spectrum.jl new file mode 100644 index 0000000..5180cdb --- /dev/null +++ b/test/angular_spectrum.jl @@ -0,0 +1,56 @@ +@testset "Angular Spectrum" begin + + @testset "Test gradient with Finite Differences" begin + field = zeros(ComplexF64, (24, 24)) + field[14:16, 14:16] .= 1 + + gg(x) = sum(abs2.(x .- angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6)[1])) + + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + + out1 = gradient(gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + AS, _ = AngularSpectrum(field, 100e-6, 633e-9, 100e-6) + + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + + out3 = gradient(f_AS, field)[1] + + @test out3 ≈ out1 + + + field = zeros(ComplexF64, (15, 15)) + field[5:6, 3:8] .= 1 + gg(x) = sum(abs2.(x .- angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6)[1])) + out1 = gradient(gg, field)[1] + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + AS, _ = AngularSpectrum(field, 100e-6, 633e-9, 100e-6) + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + out3 = gradient(f_AS, field)[1] + @test out3 ≈ out1 + @test out2 ≈ out1 + + + field = zeros(ComplexF64, (15, 15)) + field[5:6, 3:8] .= 1 + gg(x) = sum(abs2.(x .- angular_spectrum(cis.(x), [100e-6, 200e-6], 633e-9, 100e-6)[1])) + out1 = gradient(gg, field)[1] + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + AS, _ = AngularSpectrum(field, [100e-6, 200e-6], 633e-9, 100e-6) + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + out3 = gradient(f_AS, field)[1] + @test out3 ≈ out1 + @test out2 ≈ out1 + + end + + @testset "double slit" begin + include("double_slit.jl") + end + + @testset "Gaussian beam" begin + include("gaussian_beam.jl") + end +end diff --git a/test/double_slit.html b/test/double_slit.html new file mode 100644 index 0000000..791c769 --- /dev/null +++ b/test/double_slit.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/test/double_slit.jl b/test/double_slit.jl new file mode 100644 index 0000000..98e6a8c --- /dev/null +++ b/test/double_slit.jl @@ -0,0 +1,89 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +using Markdown +using InteractiveUtils + +# ╔═╡ 9745c33c-c8d4-11ed-216d-cb75e0aaec28 +# ╠═╡ skip_as_script = true +#=╠═╡ +begin + using Pkg, Revise + Pkg.activate("../examples/.") +end + ╠═╡ =# + +# ╔═╡ 4acbbe98-f033-42c5-b6d0-60a5ddf2570e +using Test, WaveOpticsPropagation, IndexFunArrays, NDTools, FourierTools + +# ╔═╡ 2b5650e9-db71-459d-ae13-7b52ca8d2a85 +# ╠═╡ skip_as_script = true +#=╠═╡ +using Plots, ImageShow + ╠═╡ =# + +# ╔═╡ 83144a55-ac66-4b92-b9c7-36daacbab791 +""" + double_slit + +`d` is distance between slits in pixel +`b` is the width of the slits in pixel +""" +function double_slit(N, d, b, offset, z, λ, L) + + slit_init = (box(ComplexF32, (N, N), (N, b), offset=(11, offset + N÷2 + 1 - d ÷2 )) .+ + box(ComplexF32, (N, N), (N, b), offset=(11, offset + N÷2 + 1 + d ÷2 ))) + + xpos = fftpos(L, N, NDTools.CenterFT) .- offset .* L ./ N + sinθ = sin.(atan.(xpos, z)) + d_m = d * L / N + b_m = b * L / N + slit_ana = cos.(π .* d_m .* sinθ ./ λ).^2 .* sinc.(b_m .* sinθ ./ λ).^2 + + slit_prop = (angular_spectrum(slit_init, z, λ, L)[1]) + slit_prop2 = (AngularSpectrum(slit_init, z, λ, L)[1](slit_init))[1] + @test slit_prop2 ≈ slit_prop + + slit_prop = abs2.(slit_prop) + slit_prop2 = abs2.(slit_prop2) + + slit_prop ./= maximum(slit_prop[11, :]) + slit_prop2 ./= maximum(slit_prop2[11, :]) + + @test all(.≈(0.01f0 .+ slit_ana, 0.01f0 .+ slit_prop[11, :], rtol=0.3)) + @test findmax(slit_ana) == findmax(slit_prop[11, :]) + + return slit_init, slit_prop, slit_prop2, slit_ana +end + +# ╔═╡ c12b2aac-5e48-47b2-a0f8-e626a0020b20 +slit_init, slit_prop, slit_prop2, slit_ana = double_slit(100, 10, 3, 0, 25e-3, 1000e-9, 1e-3) + +# ╔═╡ 33a6b4e4-ee04-45fe-96b4-6fa8614b5c6c +#=╠═╡ +begin + #plot(abs2.(slit_init[11, :])) + plot(slit_prop[11, :]) + plot!(slit_ana) +end + ╠═╡ =# + +# ╔═╡ 6b2554ac-613a-4b06-b9a5-fb96a81816ee +double_slit(512, 10, 3, -30, 5e-3, 632e-9, 1e-3) + +# ╔═╡ 05587986-fe48-4ac8-91eb-5fedfbf6634e +double_slit(512, 10, 3, 30, 5e-3, 1000e-9, 1e-3) + +# ╔═╡ 8c9133cd-ae62-4a3d-81b8-4e30c795a04f + double_slit(100, 10, 3, 0, 25e-3, 1000e-9, 1e-3) + +# ╔═╡ Cell order: +# ╠═9745c33c-c8d4-11ed-216d-cb75e0aaec28 +# ╠═4acbbe98-f033-42c5-b6d0-60a5ddf2570e +# ╠═2b5650e9-db71-459d-ae13-7b52ca8d2a85 +# ╠═83144a55-ac66-4b92-b9c7-36daacbab791 +# ╠═c12b2aac-5e48-47b2-a0f8-e626a0020b20 +# ╠═33a6b4e4-ee04-45fe-96b4-6fa8614b5c6c +# ╠═6b2554ac-613a-4b06-b9a5-fb96a81816ee +# ╠═05587986-fe48-4ac8-91eb-5fedfbf6634e +# ╠═8c9133cd-ae62-4a3d-81b8-4e30c795a04f diff --git a/test/double_slit.plutostate b/test/double_slit.plutostate new file mode 100644 index 0000000..a142d2f Binary files /dev/null and b/test/double_slit.plutostate differ diff --git a/test/fraunhofer.jl b/test/fraunhofer.jl new file mode 100644 index 0000000..68e04ce --- /dev/null +++ b/test/fraunhofer.jl @@ -0,0 +1,49 @@ +@testset "Angular Spectrum" begin + L = 100f-3 + λ = 633f-9 + z = 1 + N = 256 + ΔS = 18 + ΔW = 2 + slit = zeros(ComplexF32, (N, N)) + + x = WaveOpticsPropagation.fftpos(L, N) + W = x[1 + ΔW * 2 + 1] - x[1] + S = x[1 + ΔS * 2] - x[1] + + L_new_ref = λ * z / L * N + mid = N ÷ 2+ 1 + slit[:, mid-ΔS-ΔW:mid-ΔS+ΔW] .= 1 + slit[:, mid+ΔS-ΔW:mid+ΔS+ΔW] .= 1 + output, t = fraunhofer(slit, z, λ, L) + L_new = t.L + + efficient_fraunhofer, t = Fraunhofer(slit, z, λ, L); + output2, L_new2 = efficient_fraunhofer(slit) + + + @test L_new ≈ L_new_ref + @test L_new2.L ≈ L_new_ref + + @test output2 ≈ output + + I_analytical(x) = sinc(W * x / λ / z)^2 * cos(π * S * x / λ / z)^2 + intensity = abs2.(output) + intensity ./= maximum(intensity) + xpos_out = WaveOpticsPropagation.fftpos(L_new, N, NDTools.CenterFT) + + @test all(≈(I_analytical.(xpos_out)[110:150] .+ 1, 1 .+ intensity[129, 110:150, :], rtol=1f-2)) + + arr = randn(ComplexF32, (N, N)) + fr, t = Fraunhofer(arr, z, λ, L) + f(x) = sum(abs2, arr .- fr(x)[1]) + f2(x) = sum(abs2, arr .- fraunhofer(x, z, λ, L)[1]) + @test Zygote.gradient(f, arr)[1] ≈ Zygote.gradient(f2, arr)[1] + + arr = randn(ComplexF32, (15, 15)) + fr, t = Fraunhofer(arr, z, λ, L, skip_final_phase=false) + f(x) = sum(abs2, arr .- fr(x)[1]) + f2(x) = sum(abs2, arr .- fraunhofer(x, z, λ, L, skip_final_phase=false)[1]) + @test Zygote.gradient(f, arr)[1] ≈ Zygote.gradient(f2, arr)[1] + +end diff --git a/test/gaussian_beam.html b/test/gaussian_beam.html new file mode 100644 index 0000000..35d2951 --- /dev/null +++ b/test/gaussian_beam.html @@ -0,0 +1,16 @@ + + + + + + +
\ No newline at end of file diff --git a/test/gaussian_beam.jl b/test/gaussian_beam.jl new file mode 100644 index 0000000..9b40e70 --- /dev/null +++ b/test/gaussian_beam.jl @@ -0,0 +1,75 @@ +### A Pluto.jl notebook ### +# v0.19.30 + +using Markdown +using InteractiveUtils + +# ╔═╡ 9902bc0e-c8dc-11ed-05fe-a9f16c85ba7b +# ╠═╡ skip_as_script = true +#=╠═╡ +begin + using Pkg, Revise + Pkg.activate("../examples/.") +end + ╠═╡ =# + +# ╔═╡ 748c98db-f175-4b8d-87d4-eae2ec6f26ab +# ╠═╡ skip_as_script = true +#=╠═╡ +using Plots, ImageShow, PlutoUI + ╠═╡ =# + +# ╔═╡ 38f7c686-48c2-4051-9e2e-8f2aad6e95ee +using Test, WaveOpticsPropagation, IndexFunArrays, NDTools, FourierTools + +# ╔═╡ 0299632f-e805-4e27-82e8-d67a756409f2 + + +# ╔═╡ 96986249-3edb-4fae-b4a5-ec7ae18fdb59 +function test_gauss_consistency(λ, L, N, z, z_init, w_0; do_test=true) + y = reshape(Float32.(fftpos(L[1], N[1], CenterFT)) ,:, 1) + x = reshape(Float32.(fftpos(L[2], N[2], CenterFT)), 1, :) + field = gauss_beam.(y, x, z_init, λ, w_0); + field_as = angular_spectrum(field, z, λ, L)[1]; + field_z = gauss_beam.(y, x, z_init + z, λ, w_0); + + + if do_test + @test ≈(field_as .+ maximum(abs.(field_as)) , maximum(abs.(field_as)) .+ field_z, + rtol=0.004) + @test all(≈(field_as .+ maximum(abs.(field_as)) , maximum(abs.(field_as)) .+ field_z, + rtol=0.1)) + @test sum(abs2.(field_z)) ≈ sum(abs2.(field_as)) ≈ sum(abs2.(field)) + end + + return field, field_as, field_z +end + +# ╔═╡ c732c1d7-6ae8-41ab-ae2e-d635fa9d66a3 +res1 = test_gauss_consistency(405f-9, (10f-4, 20.0f-4), (512, 512), 1f-2, 1f-10, 0.01f-3) + +# ╔═╡ a5e7878f-4ea1-4f75-87f8-fafd09f5ad82 +res2 = test_gauss_consistency(10f-9, (1f-4, 0.5f-4), (400, 400), 0.02f-2, 0.05f-2, 0.003f-3, do_test = true) + +# ╔═╡ 22e070fd-b27f-4cb2-8ad7-71f409e881bf +# ╠═╡ skip_as_script = true +#=╠═╡ +simshow([res1[2] res1[3]], γ=0.1) + ╠═╡ =# + +# ╔═╡ 1f026add-c513-42fa-a26c-5da40badf8a8 +# ╠═╡ skip_as_script = true +#=╠═╡ +simshow([res2[1] res2[2] res2[3]], γ=0.2) + ╠═╡ =# + +# ╔═╡ Cell order: +# ╠═9902bc0e-c8dc-11ed-05fe-a9f16c85ba7b +# ╠═748c98db-f175-4b8d-87d4-eae2ec6f26ab +# ╠═38f7c686-48c2-4051-9e2e-8f2aad6e95ee +# ╠═0299632f-e805-4e27-82e8-d67a756409f2 +# ╠═96986249-3edb-4fae-b4a5-ec7ae18fdb59 +# ╠═c732c1d7-6ae8-41ab-ae2e-d635fa9d66a3 +# ╠═a5e7878f-4ea1-4f75-87f8-fafd09f5ad82 +# ╠═22e070fd-b27f-4cb2-8ad7-71f409e881bf +# ╠═1f026add-c513-42fa-a26c-5da40badf8a8 diff --git a/test/gaussian_beam.plutostate b/test/gaussian_beam.plutostate new file mode 100644 index 0000000..be95974 Binary files /dev/null and b/test/gaussian_beam.plutostate differ diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..cb985f6 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,17 @@ +using WaveOpticsPropagation +using Test +using Zygote +using ChainRulesTestUtils +using IndexFunArrays, NDTools, FourierTools +using FiniteDifferences + +@testset "WaveOpticsPropagation.jl" begin + include("utils.jl") + include("angular_spectrum.jl") + include("shifted_angular_spectrum.jl") + include("scalable_angular_spectrum.jl") + include("fraunhofer.jl") +end + +return nothing + diff --git a/test/scalable_angular_spectrum.jl b/test/scalable_angular_spectrum.jl new file mode 100644 index 0000000..0ffd8a0 --- /dev/null +++ b/test/scalable_angular_spectrum.jl @@ -0,0 +1,23 @@ +@testset "Scalable Angular Spectrum" begin + + @testset "Test gradient with Finite Differences" begin + field = zeros(ComplexF64, (24, 24)) + field[14:16, 14:16] .= 1 + ss, _ = ScalableAngularSpectrum(cis.(field), 100e-6, 633e-9, 100e-6) + gg(x) = sum(abs2.(x .- ss(cis.(x))[1])) + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + out1 = gradient(gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + + field = zeros(ComplexF64, (19, 19)) + field[14:16, 14:16] .= 1 + ss, _ = ScalableAngularSpectrum(cis.(field), 100e-6, 633e-9, 100e-6, skip_final_phase=false) + gg(x) = sum(abs2.(x .- ss(cis.(x))[1])) + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + out1 = gradient(gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + + end + + +end diff --git a/test/shifted_angular_spectrum.jl b/test/shifted_angular_spectrum.jl new file mode 100644 index 0000000..12f579e --- /dev/null +++ b/test/shifted_angular_spectrum.jl @@ -0,0 +1,83 @@ +@testset "Shifted Angular Spectrum" begin + L = 50f-6 + N = 128 + α = deg2rad(10f0) + λ = 405f-9 + z = 100f-6 + sz = (N, N) + y = fftpos(L, N, CenterFT) + field = box(Float32, sz, (20,20)) .* exp.(1im .* 2f0 * π ./ λ .* y .* sin(α)); + res_AS = AngularSpectrum(field, z, λ, L)[1](field)[1]; + res = WaveOpticsPropagation.shifted_angular_spectrum(field, z, λ, L, (α , 0), bandlimit=true) + shift = z * tan(α) / L * N + shift2 = z * tan(α) / L + res2 = ShiftedAngularSpectrum(field, z, λ, L, (α, 0), bandlimit=true)[1](field) + @test all(.≈(1 .+ FourierTools.shift(res[1], (shift, 0))[round(Int, shift)+1:end, :], 1 .+ res_AS[round(Int, shift)+1:end, :], rtol=5f-2)) + @test all(.≈(1 .+ FourierTools.shift(res2[1], (shift, 0))[round(Int, shift)+1:end, :], 1 .+ res_AS[round(Int, shift)+1:end, :], rtol=5f-2)) + + + field = box(Float32, sz, (20,20)) .* exp.(1im .* 2f0 * π ./ λ .* y' .* sin(α)); + res_AS = AngularSpectrum(field, z, λ, L)[1](field)[1]; + res = WaveOpticsPropagation.shifted_angular_spectrum(field, z, λ, L, (0, α), bandlimit=true) + shift = z * tan(α) / L * N + shift2 = z * tan(α) / L + res2 = ShiftedAngularSpectrum(field, z, λ, L, (0, α), bandlimit=true)[1](field) + @test all(.≈(1 .+ FourierTools.shift(res[1], (0, shift))[:, round(Int, shift)+1:end], 1 .+ res_AS[:, round(Int, shift)+1:end], rtol=5f-2)) + @test all(.≈(1 .+ FourierTools.shift(res2[1], (0, shift))[:, round(Int, shift)+1:end], 1 .+ res_AS[:, round(Int, shift)+1:end], rtol=5f-2)) + + + res_AS = AngularSpectrum(field, z, λ, L)[1](field)[1]; + res = WaveOpticsPropagation.shifted_angular_spectrum(field, z, λ, L, (0, 0), bandlimit=true)[1] + res2 = ShiftedAngularSpectrum(field, z, λ, L, (0, 0), bandlimit=true)[1](field)[1] + @test all(.≈(1 .+ res, 1 .+ res_AS, rtol=1f-1)) + @test ≈(1 .+ res, 1 .+ res_AS, rtol=1f-2) + @test all(.≈(1 .+ res2, 1 .+ res_AS, rtol=1f-1)) + @test ≈(1 .+ res2, 1 .+ res_AS, rtol=1f-2) + + + + + @testset "Test gradient with Finite Differences" begin + field = zeros(ComplexF64, (24, 24)) + field[14:16, 14:16] .= 1 + + α = deg2rad(10) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.shifted_angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6, (α, 0))[1])) + + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + + out1 = gradient(gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + AS, _ = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0)) + + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + + out3 = gradient(f_AS, field)[1] + + @test out3 ≈ out1 + + + field = zeros(ComplexF64, (15, 15)) + field[5:6, 3:8] .= 1 + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.shifted_angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6, (α, 0))[1])) + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + out1 = gradient(gg, field)[1] + @test out1 .+ cis(1) ≈ out2 .+ cis(1) + AS, _ = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0)) + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + out3 = gradient(f_AS, field)[1] + @test out3 ≈ out1 + + + field = zeros(ComplexF64, (15, 15)) + field[5:6, 3:8] .= 1 + AS, _ = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0); extract_ramp=false) + f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + out2 = FiniteDifferences.grad(central_fdm(5, 1), f_AS, field)[1] + out3 = gradient(f_AS, field)[1] + @test out3 ≈ out2 + + + end + +end diff --git a/test/utils.jl b/test/utils.jl new file mode 100644 index 0000000..836cece --- /dev/null +++ b/test/utils.jl @@ -0,0 +1,43 @@ + +@testset "Pad" begin + @test pad(ones((2, 2)), 2) == [0.0 0.0 0.0 0.0; 0.0 1.0 1.0 0.0; 0.0 1.0 1.0 0.0; 0.0 0.0 0.0 0.0] + @test pad(ones((1,)), 2) == [0.0, 1.0] + @test pad(ones((1,)), 3) == [0.0, 1.0, 0.0] + @test pad(ones((1,)), 4) == [0.0, 0.0, 1.0, 0.0] + @test pad(ones((1,)), 5) == [0.0, 0.0, 1.0, 0.0, 0.0] + @test pad(ones((2, 2)), 3) == [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 1.0 0.0 0.0; 0.0 0.0 1.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0] + @test pad(ones((3, 3)), 3) == [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] + + @test pad(ones((2, 2)), (2, 2, 2)) == [0.0 0.0; 0.0 0.0;;; 1.0 1.0; 1.0 1.0] + @test pad(ones((1, 1)), (1, 3, 4)) == [0.0 0.0 0.0;;; 0.0 0.0 0.0;;; 0.0 1.0 0.0;;; 0.0 0.0 0.0] + @test pad(ones((1, 1)), (1, 3, 4), value = 10) == [10.0 10.0 10.0;;; 10.0 10.0 10.0;;; 10.0 1.0 10.0;;; 10.0 10.0 10.0] + @test pad(ones((2, 2)), 2.5) == [0.0 0.0 0.0 0.0 0.0; 0.0 1.0 1.0 0.0 0.0; 0.0 1.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0] +end + +@testset "set_center!" begin + @test set_center!([-1, 1, 1, 13, 1, 4], [5, 5, 5]) == [-1, 1, 5, 5, 5, 4] + @test set_center!(ones((3, 3)), [5]) == [1.0 1.0 1.0; 1.0 5.0 1.0; 1.0 1.0 1.0] + @test set_center!(ones((3, 3)), [5], broadcast = true) == [1.0 1.0 1.0; 5.0 5.0 5.0; 1.0 1.0 1.0] +end + +@testset "crop center" begin + @test crop_center([1 2; 3 4], (1,)) == [3 4] + @test crop_center([1 2 3; 3 4 5], (1,)) == [3 4 5] + @test crop_center([1 2 3; 3 4 5], (1, 2)) == [3 4] + @test crop_center([1 2 3; 3.3 4 5], (1, 2)) == [3.3 4.0] + @test crop_center([1 2 3; 3.3 4 5], (2, 2)) == [1.0 2.0; 3.3 4.0] +end + + +@testset "test rrule for pad and crop_center" begin + test_rrule(pad, randn((4,2)), 2, check_thunked_output_tangent=false) + test_rrule(pad, randn((2,8)), 2.5, check_thunked_output_tangent=false) + test_rrule(pad, randn((1,9)), 1, check_thunked_output_tangent=false) + test_rrule(pad, randn((1,9)), (2,20), check_thunked_output_tangent=false) + test_rrule(pad, randn((1,9)), (1,9), check_thunked_output_tangent=false) + + test_rrule(crop_center, randn((6,6)), (1,1), check_thunked_output_tangent=false) + test_rrule(crop_center, randn((6,4)), (2,3), check_thunked_output_tangent=false) + test_rrule(crop_center, randn((5,5)), (3,3), check_thunked_output_tangent=false) + test_rrule(crop_center, randn((5,5)), (5,5), check_thunked_output_tangent=false) +end