From 880cc5c579033add0bc891bcf739b42ef20d06eb Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sun, 18 Feb 2024 13:50:02 +0100 Subject: [PATCH 01/13] Add shifted SAS --- examples/Manifest.toml | 147 +++++++++--------- examples/shifted_scalable_angular_spectrum.jl | 121 ++++++++++++++ src/WaveOpticsPropagation.jl | 1 + src/shifted_SAS.jl | 101 ++++++++++++ 4 files changed, 299 insertions(+), 71 deletions(-) create mode 100644 examples/shifted_scalable_angular_spectrum.jl create mode 100644 src/shifted_SAS.jl diff --git a/examples/Manifest.toml b/examples/Manifest.toml index b85b1d6..9bf6230 100644 --- a/examples/Manifest.toml +++ b/examples/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.0" +julia_version = "1.10.1" manifest_format = "2.0" project_hash = "5adba9be67adc706c964b884ad91ae500b32e8e9" @@ -29,9 +29,9 @@ version = "1.2.3" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608" +git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.7.2" +version = "4.0.1" weakdeps = ["StaticArrays"] [deps.Adapt.extensions] @@ -54,9 +54,9 @@ version = "0.2.0" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa" +git-tree-sha1 = "c5aeb516a84459e0318a02507d2261edad97eb75" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.7.0" +version = "7.7.1" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" @@ -158,10 +158,10 @@ uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" version = "0.5.0" [[deps.CUDA]] -deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "Statistics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "95ac638373ac40e29c1b6d086a3698f5026ff6a6" +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"] +git-tree-sha1 = "baa8ea7a1ea63316fa3feb454635215773c9c845" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "5.1.2" +version = "5.2.0" weakdeps = ["ChainRulesCore", "SpecialFunctions"] [deps.CUDA.extensions] @@ -182,9 +182,9 @@ version = "0.2.3" [[deps.CUDA_Runtime_jll]] deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "9704e50c9158cf8896c2776b8dbc5edd136caf80" +git-tree-sha1 = "8e25c009d2bf16c2c31a70a6e9e8939f7325cc84" uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" -version = "0.10.1+0" +version = "0.11.1+0" [[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"] @@ -200,15 +200,15 @@ version = "0.2.2" [[deps.ChainRules]] deps = ["Adapt", "ChainRulesCore", "Compat", "Distributed", "GPUArraysCore", "IrrationalConstants", "LinearAlgebra", "Random", "RealDot", "SparseArrays", "SparseInverseSubset", "Statistics", "StructArrays", "SuiteSparse"] -git-tree-sha1 = "4cfc4916725289132746f730755262e1919cff38" +git-tree-sha1 = "4e42872be98fa3343c4f8458cbda8c5c6a6fa97c" uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2" -version = "1.59.1" +version = "1.63.0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "0d12ee16b3f62e4e33c3277773730a5b21a74152" +git-tree-sha1 = "ad25e7d21ce10e01de973cdc68ad0f850a953c52" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.20.0" +version = "1.21.1" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] @@ -228,15 +228,15 @@ version = "1.3.5" [[deps.CodecBzip2]] deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] -git-tree-sha1 = "c0ae2a86b162fb5d7acc65269b469ff5b8a73594" +git-tree-sha1 = "9b1ca1aa6ce3f71b3d1840c538a8210a043625eb" uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" -version = "0.8.1" +version = "0.8.2" [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8" +git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.3" +version = "0.7.4" [[deps.ColorSchemes]] deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] @@ -281,13 +281,13 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+1" +version = "1.1.0+0" [[deps.ComponentArrays]] deps = ["ArrayInterface", "ChainRulesCore", "ForwardDiff", "Functors", "LinearAlgebra", "PackageExtensionCompat", "StaticArrayInterface", "StaticArraysCore"] -git-tree-sha1 = "0bcacba787140924e954304216dcc09f968176ca" +git-tree-sha1 = "6e56489ecac859931b991c82dbcfcf65fb07f626" uuid = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -version = "0.15.7" +version = "0.15.9" [deps.ComponentArrays.extensions] ComponentArraysAdaptExt = "Adapt" @@ -327,9 +327,9 @@ version = "0.3.2" [[deps.ConcurrentUtilities]] deps = ["Serialization", "Sockets"] -git-tree-sha1 = "8cfa272e8bdedfa88b6aefbbca7c19f1befac519" +git-tree-sha1 = "9c4708e3ed2b799e6124b5673a712dda0b596a9b" uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" -version = "2.3.0" +version = "2.3.1" [[deps.Conda]] deps = ["Downloads", "JSON", "VersionParsing"] @@ -638,9 +638,9 @@ version = "1.0.10+0" [[deps.Functors]] deps = ["LinearAlgebra"] -git-tree-sha1 = "9a68d75d466ccc1218d0552a8e1631151c569545" +git-tree-sha1 = "166c544477f97bbadc7179ede1c1868e0e9b426b" uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196" -version = "0.4.5" +version = "0.4.7" [[deps.Future]] deps = ["Random"] @@ -660,15 +660,15 @@ version = "3.3.9+0" [[deps.GPUArrays]] deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"] -git-tree-sha1 = "85d7fb51afb3def5dcb85ad31c3707795c8bccc1" +git-tree-sha1 = "47e4686ec18a9620850bad110b79966132f14283" uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" -version = "9.1.0" +version = "10.0.2" [[deps.GPUArraysCore]] deps = ["Adapt"] -git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0" +git-tree-sha1 = "ec632f177c0d990e64d955ccc1b8c04c485a0950" uuid = "46192b85-c4d5-4398-a991-12ede77f4527" -version = "0.1.5" +version = "0.1.6" [[deps.GPUCompiler]] deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] @@ -725,9 +725,9 @@ version = "1.0.2" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "abbbb9ec3afd783a7cbd82ef01dcd088ea051398" +git-tree-sha1 = "ac7b73d562b8f4287c3b67b4c66a5395a19c1ae8" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.10.1" +version = "1.10.2" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] @@ -910,19 +910,20 @@ deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[deps.Interpolations]] -deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] -git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15" +deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "00a19d6ab0cbdea2978fc23c5a6482e02c192501" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.14.7" +version = "0.14.0" [[deps.IntervalSets]] -deps = ["Dates", "Random"] -git-tree-sha1 = "3d8866c029dd6b16e69e0d4a939c4dfcb98fac47" +git-tree-sha1 = "dba9ddf07f77f60450fe5d2e2beb9854d9a49bd0" uuid = "8197267c-284f-5f27-9208-e0e47529a953" -version = "0.7.8" -weakdeps = ["Statistics"] +version = "0.7.10" +weakdeps = ["Random", "RecipesBase", "Statistics"] [deps.IntervalSets.extensions] + IntervalSetsRandomExt = "Random" + IntervalSetsRecipesBaseExt = "RecipesBase" IntervalSetsStatisticsExt = "Statistics" [[deps.InvertedIndices]] @@ -947,9 +948,9 @@ version = "1.0.0" [[deps.JLD2]] deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "7c0008f0b7622c6c0ee5c65cbc667b69f8a65672" +git-tree-sha1 = "5ea6acdd53a51d897672edb694e3cc2912f3f8a7" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.45" +version = "0.4.46" [[deps.JLFzf]] deps = ["Pipe", "REPL", "Random", "fzf_jll"] @@ -1025,9 +1026,9 @@ version = "3.0.0+1" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"] -git-tree-sha1 = "cb4619f7353fc62a1a22ffa3d7ed9791cfb47ad8" +git-tree-sha1 = "9e70165cca7459d25406367f0c55e517a9a7bfe7" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.4.2" +version = "6.5.0" weakdeps = ["BFloat16s"] [deps.LLVM.extensions] @@ -1035,9 +1036,9 @@ weakdeps = ["BFloat16s"] [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "98eaee04d96d973e79c25d49167668c5c8fb50e2" +git-tree-sha1 = "114e3a48f13d4c18ddd7fd6a00107b4b96f60f9c" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.27+1" +version = "0.0.28+0" [[deps.LLVMLoopInfo]] git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea" @@ -1176,9 +1177,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" +git-tree-sha1 = "18144f3e9cbe9b15b070288eef858f71b291ce37" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.26" +version = "0.3.27" [deps.LogExpFunctions.extensions] LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -1244,9 +1245,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[deps.MathOptInterface]] deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] -git-tree-sha1 = "e2ae8cf5ac6daf5a3959f7f6ded9c2028b61d09d" +git-tree-sha1 = "569a003f93d7c64068d3afaab908d21f67a22cd5" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "1.25.1" +version = "1.25.3" [[deps.MbedTLS]] deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] @@ -1309,9 +1310,9 @@ version = "1.2.1" [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "806eea990fb41f9b36f1253e5697aa645bf6a9f8" +git-tree-sha1 = "302fd161eb1c439e4115b51ae456da4e9984f130" uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -version = "1.4.0" +version = "1.4.1" [[deps.NDTools]] deps = ["LinearAlgebra", "OffsetArrays", "PaddedViews", "Random", "Statistics"] @@ -1407,7 +1408,7 @@ version = "1.3.5+1" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+2" +version = "0.3.23+4" [[deps.OpenEXR]] deps = ["Colors", "FileIO", "OpenEXR_jll"] @@ -1446,9 +1447,9 @@ version = "0.5.5+0" [[deps.Optim]] deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "MathOptInterface", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "f55af9918e2a67dcadf5ec758a5ff25746c3819f" +git-tree-sha1 = "d024bfb56144d947d4fafcd9cb5cafbe3410b133" uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.8.0" +version = "1.9.2" [[deps.Opus_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1532,9 +1533,9 @@ version = "1.4.0" [[deps.Plots]] deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] -git-tree-sha1 = "38a748946dca52a622e79eea6ed35c6737499109" +git-tree-sha1 = "c4fa93d7d66acad8f6f4ff439576da9d2e890ee0" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.40.0" +version = "1.40.1" [deps.Plots.extensions] FileIOExt = "FileIO" @@ -1552,15 +1553,15 @@ version = "1.40.0" [[deps.Pluto]] deps = ["Base64", "Configurations", "Dates", "Downloads", "ExpressionExplorer", "FileWatching", "FuzzyCompletions", "HTTP", "HypertextLiteral", "InteractiveUtils", "Logging", "LoggingExtras", "MIMEs", "Malt", "Markdown", "MsgPack", "Pkg", "PrecompileSignatures", "PrecompileTools", "REPL", "RegistryInstances", "RelocatableFolders", "Scratch", "Sockets", "TOML", "Tables", "URIs", "UUIDs"] -git-tree-sha1 = "e1f12b2d765ff3f50e3732017126bbb5335b7911" +git-tree-sha1 = "449f468cbb80c3eec6e6d8443a0913d8bbad4d0d" uuid = "c3e4b0f8-55cb-11ea-2926-15256bba5781" -version = "0.19.37" +version = "0.19.38" [[deps.PlutoUI]] deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "FixedPointNumbers", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "MIMEs", "Markdown", "Random", "Reexport", "URIs", "UUIDs"] -git-tree-sha1 = "68723afdb616445c6caaef6255067a8339f91325" +git-tree-sha1 = "211cdf570992b0d977fda3745f72772e0d5423f2" uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8" -version = "0.7.55" +version = "0.7.56" [[deps.PooledArrays]] deps = ["DataAPI", "Future"] @@ -1642,9 +1643,9 @@ version = "5.15.3+2" [[deps.Quaternions]] deps = ["LinearAlgebra", "Random", "RealDot"] -git-tree-sha1 = "9a46862d248ea548e340e30e2894118749dc7f51" +git-tree-sha1 = "994cc27cdacca10e68feb291673ec3a76aa2fae9" uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" -version = "0.7.5" +version = "0.7.6" [[deps.REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] @@ -1742,9 +1743,13 @@ version = "1.4.2" [[deps.Rotations]] deps = ["LinearAlgebra", "Quaternions", "Random", "StaticArrays"] -git-tree-sha1 = "792d8fd4ad770b6d517a13ebb8dadfcac79405b8" +git-tree-sha1 = "2a0a5d8569f481ff8840e3b7c84bbf188db6a3fe" uuid = "6038ab10-8711-5258-84ad-4b1120ba62dc" -version = "1.6.1" +version = "1.7.0" +weakdeps = ["RecipesBase"] + + [deps.Rotations.extensions] + RotationsRecipesBaseExt = "RecipesBase" [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -1859,9 +1864,9 @@ version = "0.1.1" [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" +git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.8" +version = "0.8.10" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] @@ -1876,9 +1881,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5" +git-tree-sha1 = "7b0e9c14c624e435076d19aea1e5cbdec2b9ca37" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.1" +version = "1.9.2" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1991,9 +1996,9 @@ uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" version = "0.5.23" [[deps.TranscodingStreams]] -git-tree-sha1 = "1fbeaaca45801b4ba17c251dd8603ef24801dd84" +git-tree-sha1 = "54194d92959d8ebaa8e26227dbe3cdefcdcd594f" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.10.2" +version = "0.10.3" weakdeps = ["Random", "Test"] [deps.TranscodingStreams.extensions] @@ -2314,9 +2319,9 @@ version = "2.0.2+0" [[deps.libpng_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "93284c28274d9e75218a416c65ec49d0e0fcdf3d" +git-tree-sha1 = "873b4f805771d3e4bafe63af759a26ea8ca84d14" uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" -version = "1.6.40+0" +version = "1.6.42+0" [[deps.libsixel_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"] diff --git a/examples/shifted_scalable_angular_spectrum.jl b/examples/shifted_scalable_angular_spectrum.jl new file mode 100644 index 0000000..4d1582c --- /dev/null +++ b/examples/shifted_scalable_angular_spectrum.jl @@ -0,0 +1,121 @@ +### A Pluto.jl notebook ### +# v0.19.37 + +using Markdown +using InteractiveUtils + +# ╔═╡ 59413ae0-c2ce-11ee-31f3-c52d9d263afb +begin + using Pkg + Pkg.activate("../WaveOpticsPropagation.jl/examples/") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 87a307ea-506d-49d1-af79-ae1f3c52912a +using WaveOpticsPropagation, CUDA, IndexFunArrays, NDTools, ImageShow, FileIO, Zygote, Optim, Plots, PlutoUI, FourierTools, NDTools + +# ╔═╡ ea14a26f-9273-4f13-8edf-fc0691a1930f +TableOfContents() + +# ╔═╡ 0ee0e517-5408-42d7-9bb3-d0adef6df1e5 +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 + +# ╔═╡ d297ace0-deec-417e-88af-6fb7bf5d3358 + + +# ╔═╡ 53832390-bd97-4394-8f57-a253da354030 +N = 256 + +# ╔═╡ f1389ab1-44b9-485c-a2c1-ec8a7945c1e2 +sz = (N, N) + +# ╔═╡ a9ab5744-6a83-4928-99cc-19c7058724c1 +α = deg2rad(10f0) + +# ╔═╡ c487f254-5f32-4e8b-94d7-f18c45a2a3ca +L = 50f-6 + +# ╔═╡ 652adc09-d4f5-47b9-b5ac-c16e47ee03cc +y = fftpos(L[1], N, CenterFT) + +# ╔═╡ 5a06d1f8-ce32-4510-adc2-862a6c48a479 +λ = 405f-9 + +# ╔═╡ 4825d17b-eaff-4ff1-b492-b695666014f1 +field = box(Float32, sz, (20,20)) .* exp.(1im .* 2f0 * π ./ λ .* y .* sin(α)); + +# ╔═╡ 9c9560ca-d6c4-4fe8-9ac9-4de73e6af038 +z = 200f-6 + +# ╔═╡ 3cd421df-5f7b-4e48-b13e-35480200d4a7 +simshow(field) + +# ╔═╡ 8dd05266-5383-408b-a747-b50a0f4fba66 +res_AS = AngularSpectrum(field, z, λ, L)[1](field)[1]; + +# ╔═╡ 50f75305-9107-4495-8ae9-68cd11498e16 +res = ShiftedAngularSpectrum(field, z, λ, L, (α , 0), bandlimit=true)[1](field) + +# ╔═╡ c95c1b1e-7645-4488-883e-e45be15717d8 +res2 = WaveOpticsPropagation.ShiftedScalableAngularSpectrum(field, z, λ, L, (α , 0))[1](field) + +# ╔═╡ 2e41d692-4705-47f1-81f9-8e388d265847 +res3 = WaveOpticsPropagation.ScalableAngularSpectrum(field, z, λ, L)[1](field) + +# ╔═╡ 1be83926-2d7a-49e9-8b2b-9055d67c87ea +z * λ / L^2 * size(field,1) + +# ╔═╡ 016cde31-2cf7-4f56-b7d1-e4f709b2183e +Revise.errors() + +# ╔═╡ 50055786-7728-4580-81e5-4ebb2e05626b +shift = (z .* tan.(α) ./ L .* N)[1] + +# ╔═╡ 4a49a6e5-b50d-4d65-ad1a-75eda9b1a280 +shift2 = z .* tan.(α) ./ L + +# ╔═╡ 888cff79-9df6-4bb6-bdef-bf42ad0385cb +simshow(res2[1]) + +# ╔═╡ 3e012a80-466f-4ffe-ba80-8bd99e5fe230 + + +# ╔═╡ ba3721fa-9eae-4dc4-ae36-765ce2b4bd5a +simshow(res3[1]) + +# ╔═╡ 0dbc516f-34f3-4047-b3db-d96147b22b13 +[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)] + +# ╔═╡ Cell order: +# ╠═59413ae0-c2ce-11ee-31f3-c52d9d263afb +# ╠═ea14a26f-9273-4f13-8edf-fc0691a1930f +# ╠═87a307ea-506d-49d1-af79-ae1f3c52912a +# ╠═0ee0e517-5408-42d7-9bb3-d0adef6df1e5 +# ╠═d297ace0-deec-417e-88af-6fb7bf5d3358 +# ╠═53832390-bd97-4394-8f57-a253da354030 +# ╠═f1389ab1-44b9-485c-a2c1-ec8a7945c1e2 +# ╠═652adc09-d4f5-47b9-b5ac-c16e47ee03cc +# ╠═4825d17b-eaff-4ff1-b492-b695666014f1 +# ╠═a9ab5744-6a83-4928-99cc-19c7058724c1 +# ╠═c487f254-5f32-4e8b-94d7-f18c45a2a3ca +# ╠═5a06d1f8-ce32-4510-adc2-862a6c48a479 +# ╠═9c9560ca-d6c4-4fe8-9ac9-4de73e6af038 +# ╠═3cd421df-5f7b-4e48-b13e-35480200d4a7 +# ╠═8dd05266-5383-408b-a747-b50a0f4fba66 +# ╠═50f75305-9107-4495-8ae9-68cd11498e16 +# ╠═c95c1b1e-7645-4488-883e-e45be15717d8 +# ╠═2e41d692-4705-47f1-81f9-8e388d265847 +# ╠═1be83926-2d7a-49e9-8b2b-9055d67c87ea +# ╠═016cde31-2cf7-4f56-b7d1-e4f709b2183e +# ╠═50055786-7728-4580-81e5-4ebb2e05626b +# ╠═4a49a6e5-b50d-4d65-ad1a-75eda9b1a280 +# ╠═888cff79-9df6-4bb6-bdef-bf42ad0385cb +# ╠═3e012a80-466f-4ffe-ba80-8bd99e5fe230 +# ╠═ba3721fa-9eae-4dc4-ae36-765ce2b4bd5a +# ╠═0dbc516f-34f3-4047-b3db-d96147b22b13 diff --git a/src/WaveOpticsPropagation.jl b/src/WaveOpticsPropagation.jl index 1959d53..bfdd417 100644 --- a/src/WaveOpticsPropagation.jl +++ b/src/WaveOpticsPropagation.jl @@ -17,5 +17,6 @@ include("scalable_angular_spectrum.jl") include("fraunhofer.jl") include("beams.jl") include("conv.jl") +include("shifted_SAS.jl") end diff --git a/src/shifted_SAS.jl b/src/shifted_SAS.jl new file mode 100644 index 0000000..69cc542 --- /dev/null +++ b/src/shifted_SAS.jl @@ -0,0 +1,101 @@ +export ShiftedScalableAngularSpectrum + + # highly optimized version with pre-planning +struct ShiftedScalableAngularSpectrum{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 + +function ShiftedScalableAngularSpectrum(ψ₀::AbstractArray{CT}, z, λ, L, α; + skip_final_phase=true) where {CT} + @assert size(ψ₀, 1) == size(ψ₀, 2) "Restricted to auadratic fields." + + pad_factor = 2 + set_pad_zero = false + + + N = size(ψ₀, 1) + 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 + + W = 1#.*(smooth_f.(cx.^2 .* (1 .+ tx.^2) ./ tx.^2 .+ cy.^2, limits...), + # smooth_f.(cy.^2 .* (1 .+ ty.^2) ./ ty.^2 .+ cx.^2, limits...)) + + sinxy = .+ sin.(α) + tanxy = .+ tan.(α) + + # ΔH is the core part of Fresnel and AS + H_AS = sqrt.(0im .+ 1 .- abs2.(f_x .* λ .+ sinxy[2]) .- abs2.(f_y .* λ .+ sinxy[1])) + H_AS = (sqrt.(CT(1) .- abs2.(f_x .* λ .+ sinxy[2]) .- abs2.(f_y .* λ .+ sinxy[1]) + .+ λ .* (tanxy[2] .* f_x .+ tanxy[1] .* 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(CT), typeof(FFTplan)}( + ΔH, H₁, H₂, buffer, buffer2, Q / 2, FFTplan, pad_factor), (; L=Q / 2) +end + + + + +function (sas::ShiftedScalableAngularSpectrum)(ψ::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 + From 5317ae37849c56f3c87aa76cd61ef0626154d036 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sun, 18 Feb 2024 14:59:53 +0100 Subject: [PATCH 02/13] add more code --- examples/shifted_scalable_angular_spectrum.jl | 4 ++-- src/shifted_SAS.jl | 23 +++++++++++-------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/examples/shifted_scalable_angular_spectrum.jl b/examples/shifted_scalable_angular_spectrum.jl index 4d1582c..edf23fd 100644 --- a/examples/shifted_scalable_angular_spectrum.jl +++ b/examples/shifted_scalable_angular_spectrum.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.37 +# v0.19.38 using Markdown using InteractiveUtils @@ -7,7 +7,7 @@ using InteractiveUtils # ╔═╡ 59413ae0-c2ce-11ee-31f3-c52d9d263afb begin using Pkg - Pkg.activate("../WaveOpticsPropagation.jl/examples/") + Pkg.activate(".") Pkg.instantiate() using Revise end diff --git a/src/shifted_SAS.jl b/src/shifted_SAS.jl index 69cc542..874e57a 100644 --- a/src/shifted_SAS.jl +++ b/src/shifted_SAS.jl @@ -1,13 +1,15 @@ export ShiftedScalableAngularSpectrum # highly optimized version with pre-planning -struct ShiftedScalableAngularSpectrum{AP, AP2, AB, T, P} +struct ShiftedScalableAngularSpectrum{AP, AP2, AB, T, P, AV} ΔH::AP H₁::AP H₂::AP2 # could be nothing! buffer::AB buffer2::AB L::T + xout::AV + yout::AV FFTplan::P pad_factor::Int end @@ -28,7 +30,10 @@ function ShiftedScalableAngularSpectrum(ψ₀::AbstractArray{CT}, z, λ, L, α; ψ_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 - + + + + W = 1#.*(smooth_f.(cx.^2 .* (1 .+ tx.^2) ./ tx.^2 .+ cy.^2, limits...), # smooth_f.(cy.^2 .* (1 .+ ty.^2) ./ ty.^2 .+ cx.^2, limits...)) @@ -48,12 +53,12 @@ function ShiftedScalableAngularSpectrum(ψ₀::AbstractArray{CT}, z, λ, L, α; # 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 = similar(ψ_p, real(CT), (pad_factor * N,1)) + q_x = similar(ψ_p, real(CT), (1,pad_factor * N)) + q_yx_shift = tanxy .* z q_y .= fftpos(dq * pad_factor * N, pad_factor * N, CenterFT) - q_y = ifftshift(q_y) - q_x = q_y' + q_y = q_yx_shift[1] .+ ifftshift(q_y) + q_x .= q_yx_shift[2] .+ q_y' # calculate phases of Fresnel H₁ = exp.(1im .* k ./ (2 .* z) .* (x .^ 2 .+ y .^ 2)) @@ -70,9 +75,7 @@ function ShiftedScalableAngularSpectrum(ψ₀::AbstractArray{CT}, z, λ, L, α; buffer = similar(ψ_p) buffer2 = similar(buffer) - return ScalableAngularSpectrum{typeof(ΔH), typeof(H₂), typeof(buffer), - real(CT), typeof(FFTplan)}( - ΔH, H₁, H₂, buffer, buffer2, Q / 2, FFTplan, pad_factor), (; L=Q / 2) + return ShiftedScalableAngularSpectrum(ΔH, H₁, H₂, buffer, buffer2, L, q_x, q_y, FFTplan, pad_factor) end From 52e2a7c32e8d5475a917d76b12e7c45835f21aef Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sun, 18 Feb 2024 16:11:47 +0100 Subject: [PATCH 03/13] Change API --- src/angular_spectrum.jl | 21 ++++++++++++------ src/shifted_angular_spectrum.jl | 7 +++--- test/angular_spectrum.jl | 22 +++++++++--------- test/double_slit.jl | 4 ++-- test/gaussian_beam.jl | 2 +- test/shifted_angular_spectrum.jl | 38 ++++++++++++++++---------------- 6 files changed, 50 insertions(+), 44 deletions(-) diff --git a/src/angular_spectrum.jl b/src/angular_spectrum.jl index 70b2f7a..cb2541c 100644 --- a/src/angular_spectrum.jl +++ b/src/angular_spectrum.jl @@ -1,4 +1,3 @@ -export angular_spectrum export AngularSpectrum @@ -115,7 +114,7 @@ function angular_spectrum(field::AbstractArray{CT, 2}, z, λ, L; 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) + return field_out_cropped end @@ -145,10 +144,18 @@ See [`angular_spectrum`](@ref) for the full documentation. julia> field = zeros(ComplexF32, (4,4)); field[3,3] = 1 1 -julia> as, t = AngularSpectrum(field, 100e-9, 632e-9, 10e-6); +julia> as = AngularSpectrum(field, 100e-9, 632e-9, 10e-6) +WaveOpticsPropagation.AngularSpectrum3{Matrix{ComplexF32}, Float64, FFTW.cFFTWPlan{ComplexF32, -1, true, 2, Tuple{Int64, Int64}}}(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], ComplexF32[3.0f-45 + 6.0f-45im 6.3f-44 + 7.0f-44im … 3.46f-43 + 2.42f-43im 7.1f-44 + 7.1f-44im; 8.0f-45 + 1.1f-44im 7.3f-44 + 7.6f-44im … 3.48f-43 + 2.4f-43im 3.69f-43 + 6.9f-44im; … ; 4.9f-44 + 5.9f-44im 1.11f-43 + 1.15f-43im … 7.7f-44 + 7.7f-44im -327424.25f0 + 4.579f-41im; 1.3f-44 + 6.2f-44im 1.16f-43 + 1.23f-43im … 3.66f-43 + 7.4f-44im -328013.5f0 + 4.579f-41im], ComplexF32[3.0f-45 + 6.0f-45im 6.3f-44 + 7.0f-44im … 3.46f-43 + 2.42f-43im 7.1f-44 + 7.1f-44im; 8.0f-45 + 1.1f-44im 7.3f-44 + 7.6f-44im … 3.48f-43 + 2.4f-43im 3.69f-43 + 6.9f-44im; … ; 4.9f-44 + 5.9f-44im 1.11f-43 + 1.15f-43im … 7.7f-44 + 7.7f-44im -327424.25f0 + 4.579f-41im; 1.3f-44 + 6.2f-44im 1.16f-43 + 1.23f-43im … 3.66f-43 + 7.4f-44im -328013.5f0 + 4.579f-41im], 1.0e-5, FFTW in-place forward plan for 8×8 array of ComplexF32 +(dft-rank>=2/1 + (dft-direct-8-x8 "n1fv_8_avx2_128") + (dft-direct-8-x8 "n1fv_8_avx2_128")), true, 2) 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,)) +4×4 Matrix{ComplexF32}: + 6.51926f-8-3.53903f-7im -2.50973f-7+1.19667f-6im 0.000417544-0.000277363im -2.50973f-7+1.19667f-6im + -2.5137f-7+1.19668f-6im 8.54086f-7-4.05125f-6im -0.00142378+0.000938406im 8.54086f-7-4.05125f-6im + 0.000417545-0.00027737im -0.00142377+0.000938403im 0.549778+0.835303im -0.00142377+0.000938403im + -2.5137f-7+1.19668f-6im 8.54086f-7-4.05125f-6im -0.00142378+0.000938406im 8.54086f-7-4.05125f-6im ``` """ function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; @@ -172,7 +179,7 @@ function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; H .= H .* W HW = H - return AngularSpectrum3{typeof(H), typeof(L), typeof(p)}(HW, buffer, buffer2, L, p, padding, pad_factor), L + return AngularSpectrum3{typeof(H), typeof(L), typeof(p)}(HW, buffer, buffer2, L, p, padding, pad_factor) end """ @@ -188,7 +195,7 @@ function (as::AngularSpectrum3)(field) 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) + return field_out_cropped end @@ -198,7 +205,7 @@ function ChainRulesCore.rrule(as::AngularSpectrum3, field) 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] + y2 = ȳ fill!(as.buffer2, 0) field_new = as.padding ? set_center!(as.buffer2, y2, broadcast=true) : y2 diff --git a/src/shifted_angular_spectrum.jl b/src/shifted_angular_spectrum.jl index f58f637..796106e 100644 --- a/src/shifted_angular_spectrum.jl +++ b/src/shifted_angular_spectrum.jl @@ -1,6 +1,5 @@ export ShiftedAngularSpectrum - function _prepare_shifted_angular_spectrum(field::AbstractArray{CT}, z, λ, L, α; padding=true, pad_factor=2, bandlimit=true, @@ -159,7 +158,7 @@ function ShiftedAngularSpectrum(field::AbstractArray{CT, N}, z::Number, λ, L, 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) + buffer, buffer2, L, shift, p, padding, pad_factor, ramp_before, ramp_after) end @@ -185,7 +184,7 @@ function (as::ShiftedAngularSpectrum{A, T, T2, P, R})(field) where {A,T,T2,P,R} 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) + return field_out_cropped end @@ -193,7 +192,7 @@ function ChainRulesCore.rrule(as::ShiftedAngularSpectrum{A, T, T2, P, R}, field) field_and_tuple = as(field) function as_pullback(ȳ) f̄ = NoTangent() - y2 = ȳ.backing[1] + y2 = ȳ fill!(as.buffer2, 0) field_new = as.padding ? set_center!(as.buffer2, y2, broadcast=true) : y2 diff --git a/test/angular_spectrum.jl b/test/angular_spectrum.jl index 5180cdb..92b12cf 100644 --- a/test/angular_spectrum.jl +++ b/test/angular_spectrum.jl @@ -1,18 +1,18 @@ @testset "Angular Spectrum" begin @testset "Test gradient with Finite Differences" begin - field = zeros(ComplexF64, (24, 24)) - field[14:16, 14:16] .= 1 + field = zeros(ComplexF64, (14, 14)) + field[8:12, 8:12] .= 1 - gg(x) = sum(abs2.(x .- angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6)[1])) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6))) 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) + AS = AngularSpectrum(field, 100e-6, 633e-9, 100e-6) - f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) out3 = gradient(f_AS, field)[1] @@ -21,12 +21,12 @@ 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])) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6))) 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])) + AS = AngularSpectrum(field, 100e-6, 633e-9, 100e-6) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) out3 = gradient(f_AS, field)[1] @test out3 ≈ out1 @test out2 ≈ out1 @@ -34,12 +34,12 @@ 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])) + gg(x) = sum(abs2.(x .- angular_spectrum(cis.(x), [100e-6, 200e-6], 633e-9, 100e-6))) 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])) + AS = AngularSpectrum(field, [100e-6, 200e-6], 633e-9, 100e-6) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) out3 = gradient(f_AS, field)[1] @test out3 ≈ out1 @test out2 ≈ out1 diff --git a/test/double_slit.jl b/test/double_slit.jl index 98e6a8c..5329c1d 100644 --- a/test/double_slit.jl +++ b/test/double_slit.jl @@ -40,8 +40,8 @@ function double_slit(N, d, b, offset, z, λ, L) 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] + slit_prop = (angular_spectrum(slit_init, z, λ, L)) + slit_prop2 = (AngularSpectrum(slit_init, z, λ, L)(slit_init)) @test slit_prop2 ≈ slit_prop slit_prop = abs2.(slit_prop) diff --git a/test/gaussian_beam.jl b/test/gaussian_beam.jl index 9b40e70..c838d9f 100644 --- a/test/gaussian_beam.jl +++ b/test/gaussian_beam.jl @@ -30,7 +30,7 @@ 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_as = angular_spectrum(field, z, λ, L); field_z = gauss_beam.(y, x, z_init + z, λ, w_0); diff --git a/test/shifted_angular_spectrum.jl b/test/shifted_angular_spectrum.jl index 12f579e..52fbcad 100644 --- a/test/shifted_angular_spectrum.jl +++ b/test/shifted_angular_spectrum.jl @@ -7,28 +7,28 @@ 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_AS = AngularSpectrum(field, z, λ, L)(field); 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)) + res2 = ShiftedAngularSpectrum(field, z, λ, L, (α, 0), bandlimit=true)(field) + @test all(.≈(1 .+ FourierTools.shift(res, (shift, 0))[round(Int, shift)+1:end, :], 1 .+ res_AS[round(Int, shift)+1:end, :], rtol=5f-2)) + @test all(.≈(1 .+ FourierTools.shift(res2, (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_AS = AngularSpectrum(field, z, λ, L)(field); 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)) + res2 = ShiftedAngularSpectrum(field, z, λ, L, (0, α), bandlimit=true)(field) + @test all(.≈(1 .+ FourierTools.shift(res, (0, shift))[:, round(Int, shift)+1:end], 1 .+ res_AS[:, round(Int, shift)+1:end], rtol=5f-2)) + @test all(.≈(1 .+ FourierTools.shift(res2, (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] + res_AS = AngularSpectrum(field, z, λ, L)(field); + res = WaveOpticsPropagation.shifted_angular_spectrum(field, z, λ, L, (0, 0), bandlimit=true) + res2 = ShiftedAngularSpectrum(field, z, λ, L, (0, 0), bandlimit=true)(field) @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)) @@ -42,15 +42,15 @@ 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])) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.shifted_angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6, (α, 0)))) 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)) + AS = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0)) - f_AS(x) = sum(abs2.(x .- AS(cis.(x))[1])) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) out3 = gradient(f_AS, field)[1] @@ -59,20 +59,20 @@ 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])) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.shifted_angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6, (α, 0)))) 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])) + AS = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0)) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) 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])) + AS = ShiftedAngularSpectrum(field, 100e-6, 633e-9, 100e-6, (α, 0); extract_ramp=false) + f_AS(x) = sum(abs2.(x .- AS(cis.(x)))) out2 = FiniteDifferences.grad(central_fdm(5, 1), f_AS, field)[1] out3 = gradient(f_AS, field)[1] @test out3 ≈ out2 From 4b90d43227be3f740de1cd37700f4e0cc2fdbcc0 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sun, 18 Feb 2024 16:13:03 +0100 Subject: [PATCH 04/13] Change return statetment shifted AS --- src/shifted_angular_spectrum.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shifted_angular_spectrum.jl b/src/shifted_angular_spectrum.jl index 796106e..514510b 100644 --- a/src/shifted_angular_spectrum.jl +++ b/src/shifted_angular_spectrum.jl @@ -78,7 +78,7 @@ function shifted_angular_spectrum(field::AbstractArray{CT, 2}, z, λ, L, α; 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) + return field_out_cropped end From 67a9f610279e384b825c273dc9253ce30432b4b2 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sun, 18 Feb 2024 22:51:31 +0100 Subject: [PATCH 05/13] Update API --- src/WaveOpticsPropagation.jl | 12 ++++++++++++ src/angular_spectrum.jl | 38 +++++++++++++++++++----------------- src/fraunhofer.jl | 36 +++++++++++++++++++--------------- test/angular_spectrum.jl | 2 +- test/double_slit.jl | 2 +- test/fraunhofer.jl | 21 +++++++++----------- test/gaussian_beam.jl | 2 +- 7 files changed, 64 insertions(+), 49 deletions(-) diff --git a/src/WaveOpticsPropagation.jl b/src/WaveOpticsPropagation.jl index bfdd417..6af50dd 100644 --- a/src/WaveOpticsPropagation.jl +++ b/src/WaveOpticsPropagation.jl @@ -9,6 +9,16 @@ using FourierTools using IndexFunArrays using CUDA + +struct Params{M, M2} + y::M + x::M + yp::M + xp::M + L::M2 + Lp::M2 +end + include("utils.jl") include("propagation.jl") include("angular_spectrum.jl") @@ -19,4 +29,6 @@ include("beams.jl") include("conv.jl") include("shifted_SAS.jl") + + end diff --git a/src/angular_spectrum.jl b/src/angular_spectrum.jl index cb2541c..6d4d79b 100644 --- a/src/angular_spectrum.jl +++ b/src/angular_spectrum.jl @@ -17,7 +17,7 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; L = T.(L isa Number ? (L, L) : L) bandlimit_border = real(CT).(bandlimit_border) - L_new = padding ? pad_factor .* L : L + Lp = padding ? pad_factor .* L : L # applies zero padding if ndims(field) == 2 @@ -25,10 +25,10 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; 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 + fieldp = padding ? pad(field, pad_factor2) : field # helpful propagation variables - (; k, f_x, f_y) = Zygote.@ignore _propagation_variables(field_new, λ, L_new) + (; k, f_x, f_y, x, y) = Zygote.@ignore _propagation_variables(fieldp, λ, Lp) # transfer function kernel of angular spectrum H = exp.(1im .* k .* abs.(z) .* sqrt.(CT(1) .- abs2.(f_x .* λ) .- abs2.(f_y .* λ))) @@ -39,9 +39,11 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; # bandlimit according to Matsushima # as addition we introduce a smooth bandlimit with a Hann window # and fuzzy logic - Δu = 1 ./ L_new + Δu = 1 ./ Lp u_limit = Zygote.@ignore 1 ./ (sqrt.((2 .* Δu .* z).^2 .+ 1) .* λ) - + + params = Params(y, x, y, x, L, Lp) + W = let if bandlimit # bandlimit filter @@ -59,7 +61,7 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; end end - return (;field_new, H, W, fftdims) + return (;fieldp, H, W, fftdims, params) end @@ -105,11 +107,11 @@ function angular_spectrum(field::AbstractArray{CT, 2}, z, λ, L; @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, + (; fieldp, 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 = fftshift(ifft(fft(ifftshift(fieldp, fftdims), fftdims) .* H .* W, fftdims), fftdims) field_out_cropped = padding ? crop_center(field_out, size(field)) : field_out @@ -122,11 +124,11 @@ angular_spectrum(field::AbstractArray, z, λ, L; kwargs...) = throw(ArgumentErro # highly optimized version with pre-planning -struct AngularSpectrum3{A, T, P} +struct AngularSpectrum3{A, P, M, M2} HW::A buffer::A buffer2::A - L::T + params::Params{M, M2} p::P padding::Bool pad_factor::Int @@ -163,15 +165,15 @@ function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; bandlimit=true, bandlimit_border=(0.8, 1)) where {CT, N} - (; field_new, H, W, fftdims) = _prepare_angular_spectrum(field, z, λ, L; padding, + (; fieldp, H, W, fftdims, params) = _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))) + buffer2 = similar(field, complex(eltype(fieldp)), (size(fieldp, 1), size(fieldp, 2), length(z))) buffer = copy(buffer2) else - buffer2 = similar(field, complex(eltype(field_new)), (size(field_new, 1), size(field_new, 2))) + buffer2 = similar(field, complex(eltype(fieldp)), (size(fieldp, 1), size(fieldp, 2))) buffer = copy(buffer2) end @@ -179,7 +181,7 @@ function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; H .= H .* W HW = H - return AngularSpectrum3{typeof(H), typeof(L), typeof(p)}(HW, buffer, buffer2, L, p, padding, pad_factor) + return AngularSpectrum3(HW, buffer, buffer2, params, p, padding, pad_factor) end """ @@ -190,8 +192,8 @@ 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)) + fieldp = set_center!(as.buffer2, field, broadcast=true) + field_imd = as.p * ifftshift!(as.buffer, fieldp, (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 @@ -208,8 +210,8 @@ function ChainRulesCore.rrule(as::AngularSpectrum3, field) y2 = ȳ 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)) + fieldp = as.padding ? set_center!(as.buffer2, y2, broadcast=true) : y2 + field_imd = as.p * ifftshift!(as.buffer, fieldp, (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 diff --git a/src/fraunhofer.jl b/src/fraunhofer.jl index d105ba5..3c9cb00 100644 --- a/src/fraunhofer.jl +++ b/src/fraunhofer.jl @@ -56,7 +56,7 @@ function fraunhofer(U, z, λ, L; skip_final_phase=true) out = phasefactor .* fftshift(p * ifftshift(U)) ./ √(size(U, 1) * size(U, 2)) end - return out, (; L=L_new) + return out end """ @@ -85,33 +85,37 @@ julia> 4f-3 * 632f-9 * 256 / (100f-6)^2 """ function Fraunhofer(U, z, λ, L; skip_final_phase=true) @assert size(U, 1) == size(U, 2) - L_new = λ * z / L * size(U, 1) + L_new = λ * z / L Ns = size(U)[1:2] - + + 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))' + yp = similar(U, real(eltype(U)), (Ns[1], 1)) + yp .= (fftpos(L_new, Ns[1], CenterFT)) + xp = similar(U, real(eltype(U)), (1, Ns[2])) + xp .= (fftpos(L_new, Ns[2], CenterFT))' + 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) + params = Params(y, x, yp, xp, L, L_new) + return FraunhoferOp(buffer, phasefactor, params, FFTplan) end -struct FraunhoferOp{T, B, PF, P} +struct FraunhoferOp{B, PF, P, M, M2} buffer::B phasefactor::PF - L::T + params::Params{M, M2} FFTplan::P end @@ -124,7 +128,7 @@ function (fraunhofer::FraunhoferOp)(field) if !isnothing(fraunhofer.phasefactor) out .*= fraunhofer.phasefactor end - return out, (; fraunhofer.L) + return out end @@ -133,7 +137,7 @@ function ChainRulesCore.rrule(f::FraunhoferOp, U) function f_pullback(ȳ) buffer = f.buffer - y2 = ȳ.backing[1] + y2 = ȳ if !isnothing(f.phasefactor) y2 = y2 .* conj.((f.phasefactor)) end diff --git a/test/angular_spectrum.jl b/test/angular_spectrum.jl index 92b12cf..9d291b0 100644 --- a/test/angular_spectrum.jl +++ b/test/angular_spectrum.jl @@ -34,7 +34,7 @@ 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))) + gg(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), [100e-6, 200e-6], 633e-9, 100e-6))) out1 = gradient(gg, field)[1] out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] @test out1 .+ cis(1) ≈ out2 .+ cis(1) diff --git a/test/double_slit.jl b/test/double_slit.jl index 5329c1d..51a2993 100644 --- a/test/double_slit.jl +++ b/test/double_slit.jl @@ -40,7 +40,7 @@ function double_slit(N, d, b, offset, z, λ, L) 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)) + slit_prop = (WaveOpticsPropagation.angular_spectrum(slit_init, z, λ, L)) slit_prop2 = (AngularSpectrum(slit_init, z, λ, L)(slit_init)) @test slit_prop2 ≈ slit_prop diff --git a/test/fraunhofer.jl b/test/fraunhofer.jl index 68e04ce..cfdc830 100644 --- a/test/fraunhofer.jl +++ b/test/fraunhofer.jl @@ -1,4 +1,4 @@ -@testset "Angular Spectrum" begin +@testset "Fraunhofer" begin L = 100f-3 λ = 633f-9 z = 1 @@ -15,33 +15,30 @@ 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 + output = fraunhofer(slit, z, λ, L) - efficient_fraunhofer, t = Fraunhofer(slit, z, λ, L); - output2, L_new2 = efficient_fraunhofer(slit) + efficient_fraunhofer = Fraunhofer(slit, z, λ, L); + output2 = 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) + xpos_out = WaveOpticsPropagation.fftpos(L_new_ref, 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]) + fr = Fraunhofer(arr, z, λ, L) + f(x) = sum(abs2, arr .- fr(x)) + f2(x) = sum(abs2, arr .- fraunhofer(x, z, λ, L)) @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) + fr = 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] diff --git a/test/gaussian_beam.jl b/test/gaussian_beam.jl index c838d9f..0eb3e06 100644 --- a/test/gaussian_beam.jl +++ b/test/gaussian_beam.jl @@ -30,7 +30,7 @@ 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); + field_as = WaveOpticsPropagation.angular_spectrum(field, z, λ, L); field_z = gauss_beam.(y, x, z_init + z, λ, w_0); From 26c0c24ac66c6946aedf4334f70b7faaf2a2551f Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Sat, 24 Feb 2024 14:56:55 +0100 Subject: [PATCH 06/13] Adopt new API --- src/WaveOpticsPropagation.jl | 15 ++++++++++ src/scalable_angular_spectrum.jl | 46 ++++++++++++++++--------------- test/scalable_angular_spectrum.jl | 10 +++---- 3 files changed, 44 insertions(+), 27 deletions(-) diff --git a/src/WaveOpticsPropagation.jl b/src/WaveOpticsPropagation.jl index 6af50dd..c467d3e 100644 --- a/src/WaveOpticsPropagation.jl +++ b/src/WaveOpticsPropagation.jl @@ -9,7 +9,22 @@ using FourierTools using IndexFunArrays using CUDA +""" + Params{M, M2} +Calls such as `AS = Angular_Spectrum(field, z, λ, L)` will store a `Params` object in `AS.params` that +contains the physical parameters of the field and the propagated field. +This is useful to keep track of the physical parameters of the field and the propagated field. + +Has fields to store physical parameters +- `y`: is the y-axis, so the first dimension of the field +- `x`: is the x-axis, so the second dimension of the field +- `yp`: is the y-axis of the propagated field. +- `xp`: is the x-axis of the propagated field. +- `L`: is the physical size of the initial field. +- `Lp`: is the physical size of the propagated field. + +""" struct Params{M, M2} y::M x::M diff --git a/src/scalable_angular_spectrum.jl b/src/scalable_angular_spectrum.jl index 14a6e5e..5c4c8d6 100644 --- a/src/scalable_angular_spectrum.jl +++ b/src/scalable_angular_spectrum.jl @@ -1,15 +1,14 @@ export ScalableAngularSpectrum # highly optimized version with pre-planning -struct ScalableAngularSpectrum{AP, AP2, AB, T, P} +struct ScalableAngularSpectrumOp{AP, AP2, AB, T, P} ΔH::AP H₁::AP H₂::AP2 # could be nothing! buffer::AB buffer2::AB - L::T + params::T FFTplan::P - pad_factor::Int end """ @@ -39,10 +38,11 @@ function _SAS_propagation_variables(field::AbstractArray{T, M}, z, λ, L) where # y and x positions in real space #y = ifftshift(range(-L/2, L/2, length=N)) - y = similar(field, real(eltype(field)), (N,)) + y = similar(field, real(eltype(field)), (N, 1)) + x = similar(field, real(eltype(field)), (1, N)) y .= fftpos(L, N, CenterFT) y = ifftshift(y) - x = y' + x .= y' return (; k, dx, df, f_x, f_y, x, y) end @@ -158,13 +158,13 @@ function ScalableAngularSpectrum(ψ₀::AbstractArray{T}, z, λ, L ; # 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 = similar(ψ_p, (pad_factor * N, 1)) + q_x = similar(ψ_p, (1, pad_factor * N)) q_y .= fftpos(dq * pad_factor * N, pad_factor * N, CenterFT) q_y = ifftshift(q_y) - q_x = q_y' + q_x .= q_y' + # calculate phases of Fresnel H₁ = exp.(1im .* k ./ (2 .* z) .* (x .^ 2 .+ y .^ 2)) @@ -177,16 +177,21 @@ function ScalableAngularSpectrum(ψ₀::AbstractArray{T}, z, λ, L ; end FFTplan = plan_fft!(ψ_p, (1,2)) - + # output field size + yp = similar(ψ_p, real(T), (N, 1)) + xp = similar(ψ_p, real(T), (1, N)) + yp .= fftpos(dq * N, N, CenterFT) + yp = ifftshift(yp) + xp .= yp' + + params = Params(y, x, yp, xp, L, Q/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) + return ScalableAngularSpectrumOp(ΔH, H₁, H₂, buffer, buffer2, params, FFTplan) end -function (sas::ScalableAngularSpectrum)(ψ::AbstractArray{T}; return_view=false) where T +function (sas::ScalableAngularSpectrumOp)(ψ::AbstractArray{T}; return_view=false) where T p = sas.FFTplan fill!(sas.buffer2, 0) ψ_p = set_center!(sas.buffer2, ψ) @@ -204,20 +209,18 @@ function (sas::ScalableAngularSpectrum)(ψ::AbstractArray{T}; return_view=false) fftshift!(sas.buffer2, ψ_p_final) ψ_final = crop_center(sas.buffer2, size(ψ); return_view) - return ψ_final, (; sas.L) + return ψ_final end -function ChainRulesCore.rrule(sas::ScalableAngularSpectrum, ψ::AbstractArray{T}) where T - field_and_tuple = sas(ψ) +function ChainRulesCore.rrule(sas::ScalableAngularSpectrumOp, ψ::AbstractArray{T}) where T + field = sas(ψ) function sas_pullback(ȳ) p = sas.FFTplan - y2 = ȳ.backing[1] - fill!(sas.buffer2, 0) - set_center!(sas.buffer2, y2) + set_center!(sas.buffer2, ȳ) ifftshift!(sas.buffer, sas.buffer2) if !isnothing(sas.H₂) @@ -234,8 +237,7 @@ function ChainRulesCore.rrule(sas::ScalableAngularSpectrum, ψ::AbstractArray{T} final = crop_center(sas.buffer2, size(ψ)) - return NoTangent(), final end - return field_and_tuple, sas_pullback + return field, sas_pullback end diff --git a/test/scalable_angular_spectrum.jl b/test/scalable_angular_spectrum.jl index 0ffd8a0..c07564e 100644 --- a/test/scalable_angular_spectrum.jl +++ b/test/scalable_angular_spectrum.jl @@ -3,7 +3,7 @@ @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) + 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] @@ -11,10 +11,10 @@ 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] + ss = ScalableAngularSpectrum(cis.(field), 100e-6, 633e-9, 100e-6, skip_final_phase=false) + gg2(x) = sum(abs2.(x .- ss(cis.(x))[1])) + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg2, field)[1] + out1 = gradient(gg2, field)[1] @test out1 .+ cis(1) ≈ out2 .+ cis(1) end From b47ce9d85a8db79a23a44cf292864c9d82f225b9 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 18:57:18 +0100 Subject: [PATCH 07/13] Fix angular spectrum --- src/angular_spectrum.jl | 51 ++++++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 23 deletions(-) diff --git a/src/angular_spectrum.jl b/src/angular_spectrum.jl index 6d4d79b..74ed3dd 100644 --- a/src/angular_spectrum.jl +++ b/src/angular_spectrum.jl @@ -5,7 +5,7 @@ _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; +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} @@ -14,7 +14,7 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; T = real(CT) λ = T(λ) z = _transform_z(T, z) - L = T.(L isa Number ? (L, L) : L) + L = T.(_L isa Number ? (_L, _L) : _L) bandlimit_border = real(CT).(bandlimit_border) Lp = padding ? pad_factor .* L : L @@ -42,7 +42,13 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, L; Δu = 1 ./ Lp u_limit = Zygote.@ignore 1 ./ (sqrt.((2 .* Δu .* z).^2 .+ 1) .* λ) - params = Params(y, x, y, x, L, Lp) + # y and x positions in real space, use correct spacing -> fftpos + y1 = similar(field, real(eltype(field)), (size(field, 1), 1)) + y1 .= (fftpos(L[1], size(field, 1), CenterFT)) + x1 = similar(field, real(eltype(field)), (1, size(field, 2))) + x1 .= (fftpos(L[2], size(field, 2), CenterFT))' + + params = Params(y1, x1, y1, x1, L, L) W = let if bandlimit @@ -68,25 +74,6 @@ 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 @@ -138,6 +125,25 @@ end AngularSpectrum(field, z, λ, L; kwargs...) Returns a function for efficient reuse of pre-calculated kernels. +This function 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 +* `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. + + See [`angular_spectrum`](@ref) for the full documentation. @@ -168,7 +174,6 @@ function AngularSpectrum(field::AbstractArray{CT, N}, z, λ, L; (; fieldp, H, W, fftdims, params) = _prepare_angular_spectrum(field, z, λ, L; padding, pad_factor, bandlimit, bandlimit_border) - if z isa AbstractVector buffer2 = similar(field, complex(eltype(fieldp)), (size(fieldp, 1), size(fieldp, 2), length(z))) buffer = copy(buffer2) From 8a78f59e04f65d6c9270d9f703e6f40b348161b5 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 18:57:25 +0100 Subject: [PATCH 08/13] Fix fraunhofer --- src/fraunhofer.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/fraunhofer.jl b/src/fraunhofer.jl index 3c9cb00..62dc154 100644 --- a/src/fraunhofer.jl +++ b/src/fraunhofer.jl @@ -83,21 +83,23 @@ 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 +function Fraunhofer(U::AbstractArray{CT}, _z::Number, _λ, _L; skip_final_phase=true) where CT + λ = real(CT)(_λ) + z = real(CT)(_z) + L = real(CT).(_L isa Number ? (_L, _L) : _L) + L_new = λ .* z ./ L Ns = size(U)[1:2] k = eltype(U)(2π) / λ # output coordinates y = similar(U, real(eltype(U)), (Ns[1], 1)) - y .= (fftpos(L, Ns[1], CenterFT)) + y .= (fftpos(L[1], Ns[1], CenterFT)) x = similar(U, real(eltype(U)), (1, Ns[2])) - x .= (fftpos(L, Ns[2], CenterFT))' + x .= (fftpos(L[2], Ns[2], CenterFT))' yp = similar(U, real(eltype(U)), (Ns[1], 1)) - yp .= (fftpos(L_new, Ns[1], CenterFT)) + yp .= (fftpos(L_new[1], Ns[1], CenterFT)) xp = similar(U, real(eltype(U)), (1, Ns[2])) - xp .= (fftpos(L_new, Ns[2], CenterFT))' + xp .= (fftpos(L_new[2], Ns[2], CenterFT))' if skip_final_phase phasefactor = nothing From c33fe19c6a0bf832ab743ef4562b99befc358b9f Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 18:57:33 +0100 Subject: [PATCH 09/13] Fix SAS --- src/scalable_angular_spectrum.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/scalable_angular_spectrum.jl b/src/scalable_angular_spectrum.jl index 5c4c8d6..86d34d9 100644 --- a/src/scalable_angular_spectrum.jl +++ b/src/scalable_angular_spectrum.jl @@ -181,7 +181,6 @@ function ScalableAngularSpectrum(ψ₀::AbstractArray{T}, z, λ, L ; yp = similar(ψ_p, real(T), (N, 1)) xp = similar(ψ_p, real(T), (1, N)) yp .= fftpos(dq * N, N, CenterFT) - yp = ifftshift(yp) xp .= yp' params = Params(y, x, yp, xp, L, Q/2) From 7eb9b087b09d8c9f2de8a8b71a4e7336d787b304 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 18:57:48 +0100 Subject: [PATCH 10/13] Fix examples --- examples/1_simple_aperture.jl | 217 ++++++++++++++++++++++++++++++++++ examples/double_slit.jl | 16 +-- 2 files changed, 225 insertions(+), 8 deletions(-) create mode 100644 examples/1_simple_aperture.jl diff --git a/examples/1_simple_aperture.jl b/examples/1_simple_aperture.jl new file mode 100644 index 0000000..62f0dd6 --- /dev/null +++ b/examples/1_simple_aperture.jl @@ -0,0 +1,217 @@ +### 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 + +# ╔═╡ 5ad8803f-4df4-4d8b-8ca5-a417f15cee58 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ 1e3670db-7a12-4187-a4b9-aff613748f84 +# this is our main package +using WaveOpticsPropagation + +# ╔═╡ 752bacd0-5915-4983-86c7-5ef945f96130 +using ImageShow, ImageIO, PlutoUI, IndexFunArrays, Plots, NDTools + +# ╔═╡ cd1b749d-8f50-428a-9d07-f4f0b797fb8d +using CUDA + +# ╔═╡ eb37d998-ef91-4f6e-8ccd-498693c6fd26 +md"# 0. Load packages +On the first run, Julia is going to install some packages automatically. So start this notebook and give it some minutes (5-10min) to install all packages. +No worries, any future runs will be much faster to start! +" + +# ╔═╡ fae7df2a-d9ab-4400-b1b5-da4ec9f42bd5 +TableOfContents() + +# ╔═╡ 997e13f0-5892-4dc5-9948-0bb8cecaa041 +use_CUDA = Ref(true && CUDA.functional()) + +# ╔═╡ 27603905-781f-47f9-92d4-ebe71cff55ee +md" ## CUDA +CUDA accelerates the pattern generation easily by 5-20 times! +Otherwise most of the code will be multithreaded on your CPU but we strongly recommended the usage of CUDA for large scale 3D pattern generation. + +Your CUDA is functional: **$(use_CUDA[])** +" + +# ╔═╡ dd520fa8-f149-478a-87a5-948d8e02da4d +var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + +# ╔═╡ 812f7288-3815-4e87-8a5d-57f20a41d3ba +togoc(x) = use_CUDA[] ? CuArray(x) : x + +# ╔═╡ 326d9c32-53db-4e61-84c9-1568beb37590 +md"# 1. Generature Aperture" + +# ╔═╡ c8dcc43d-9274-4619-a820-dee0a11552f6 +L = 1f-3 + +# ╔═╡ 7c3812e8-24d1-4408-86cc-2332c04dfc79 +λ = 633f-9 + +# ╔═╡ 7c95e8e2-0dc6-4d3e-b57a-3e8443847ba7 +N = 256 + +# ╔═╡ 1feb58bc-74cd-4dc6-93bd-caf0f318746c +x = WaveOpticsPropagation.fftpos(L, N, CenterFT) + +# ╔═╡ 85dbabf2-b9d5-4413-ac21-090625e63532 +@bind ΔS Slider(0:N÷4, show_value=true, default=12) + +# ╔═╡ cd0a4108-5d97-4baa-b3d7-38b6d7078240 +S = x[1 + ΔS * 2] - x[1] + +# ╔═╡ d9b13717-e79b-40ce-8fda-6d9d2cfdba8d +@bind ΔW Slider(0:N÷2, show_value=true, default=2) + +# ╔═╡ d6cbc0f5-8572-4654-a270-f4aa4c1e5866 +W = x[1 + ΔW * 2 + 1] - x[1] + +# ╔═╡ e624fa70-e85b-4107-841d-021e4a699235 +begin + aperture = zeros(ComplexF32, (N, N)) + + mid = N ÷ 2+ 1 + aperture[:, mid-ΔS-ΔW:mid-ΔS+ΔW] .= 1 + aperture[:, mid+ΔS-ΔW:mid+ΔS+ΔW] .= 1 +end; + +# ╔═╡ cad80e46-1699-4472-aae7-1a8f186f5d04 +simshow(aperture) + +# ╔═╡ 966b102f-5559-4914-a9ed-85515c249ca9 +md"# 2. Propagate with Angular Spectrum" + +# ╔═╡ 6d2941d4-810a-4edf-b909-159a8d7fb3c5 +z = 30f-3 + +# ╔═╡ 229bbcb4-b4e6-4dbb-a53d-d6c2b8b6b1e5 +AS = AngularSpectrum(aperture, z, λ, (L, L)) + +# ╔═╡ 0b2891c4-a0cd-41cc-ad1f-e650f5629c6d +@mytime I = abs2.(AS(aperture)); + +# ╔═╡ c7876f3a-12d5-4775-820e-61dc47adf189 +heatmap(AS.params.xp[:], AS.params.yp[:], I) + +# ╔═╡ 4db95ee0-eaf1-4e18-b47a-78e1309d0263 +simshow(I) + +# ╔═╡ 4514eec2-f2f9-4700-9ba5-7bccd36dd0ea +md"# 3. Propagate with Fraunhofer" + +# ╔═╡ d5ea29f8-d643-459b-85af-16a6155d4d46 +FR = Fraunhofer(aperture, z, λ, (L, L)) + +# ╔═╡ d02afd96-beb7-4693-97b6-1b08f610cbea +Revise.retry() + +# ╔═╡ 4cee7136-7f37-48a5-966c-aed9d674ae48 +xout = FR.params.xp + +# ╔═╡ 3fc36591-730e-43ad-8bad-f3c23e3c40b9 +yout = FR.params.yp; + +# ╔═╡ 28cf3957-8a7b-45e4-88f8-e5bcd5d56e3a +I_fr = abs2.(FR(aperture)); + +# ╔═╡ 4761a635-0fed-4407-be3f-cfdbf4fc827e +heatmap(xout[:], yout[:], I_fr) + +# ╔═╡ cabbdb5d-1f93-4eeb-8325-a0348a616024 +FR.params.Lp + +# ╔═╡ d0f4a46e-4245-4c59-9180-9571bf0cdebf +md"# 4. Propagate with Scalable Angular Spectrum" + +# ╔═╡ ca57e236-f4b9-401a-b996-d5aaa46d3d71 +SAS = ScalableAngularSpectrum(aperture, z, λ, L) + +# ╔═╡ a22b1ecb-894b-494c-8888-e1b3f1a48f97 +xsas = SAS.params.xp + +# ╔═╡ fedcd534-120c-41aa-bd60-19734d32597a +ysas = SAS.params.yp; + +# ╔═╡ b59ebf87-4b86-4678-baad-791eaa9c28d6 +I_sas = abs2.(SAS(aperture)); + +# ╔═╡ e4ef052b-8e17-49b1-913a-33e5c32e7436 +M = SAS.params.Lp / SAS.params.L + +# ╔═╡ fd4d69e4-828e-4c95-898b-c70fa003adf4 +xsas + +# ╔═╡ cecb84d0-fbcf-41ed-898e-c2b603ee0f04 +heatmap(xsas[:], ysas[:], I_sas) + +# ╔═╡ 88200bb4-06e4-49aa-8a31-390344f54da6 + + +# ╔═╡ a8ea0c36-9fa9-4933-be4d-551565758f59 + + +# ╔═╡ Cell order: +# ╠═5ad8803f-4df4-4d8b-8ca5-a417f15cee58 +# ╠═1e3670db-7a12-4187-a4b9-aff613748f84 +# ╠═752bacd0-5915-4983-86c7-5ef945f96130 +# ╠═cd1b749d-8f50-428a-9d07-f4f0b797fb8d +# ╠═eb37d998-ef91-4f6e-8ccd-498693c6fd26 +# ╠═fae7df2a-d9ab-4400-b1b5-da4ec9f42bd5 +# ╠═997e13f0-5892-4dc5-9948-0bb8cecaa041 +# ╠═27603905-781f-47f9-92d4-ebe71cff55ee +# ╠═dd520fa8-f149-478a-87a5-948d8e02da4d +# ╠═812f7288-3815-4e87-8a5d-57f20a41d3ba +# ╠═326d9c32-53db-4e61-84c9-1568beb37590 +# ╠═c8dcc43d-9274-4619-a820-dee0a11552f6 +# ╠═7c3812e8-24d1-4408-86cc-2332c04dfc79 +# ╠═7c95e8e2-0dc6-4d3e-b57a-3e8443847ba7 +# ╠═d6cbc0f5-8572-4654-a270-f4aa4c1e5866 +# ╠═cd0a4108-5d97-4baa-b3d7-38b6d7078240 +# ╠═1feb58bc-74cd-4dc6-93bd-caf0f318746c +# ╟─85dbabf2-b9d5-4413-ac21-090625e63532 +# ╟─d9b13717-e79b-40ce-8fda-6d9d2cfdba8d +# ╠═e624fa70-e85b-4107-841d-021e4a699235 +# ╠═cad80e46-1699-4472-aae7-1a8f186f5d04 +# ╟─966b102f-5559-4914-a9ed-85515c249ca9 +# ╠═6d2941d4-810a-4edf-b909-159a8d7fb3c5 +# ╠═229bbcb4-b4e6-4dbb-a53d-d6c2b8b6b1e5 +# ╠═0b2891c4-a0cd-41cc-ad1f-e650f5629c6d +# ╠═c7876f3a-12d5-4775-820e-61dc47adf189 +# ╠═4db95ee0-eaf1-4e18-b47a-78e1309d0263 +# ╟─4514eec2-f2f9-4700-9ba5-7bccd36dd0ea +# ╠═d5ea29f8-d643-459b-85af-16a6155d4d46 +# ╠═d02afd96-beb7-4693-97b6-1b08f610cbea +# ╠═4cee7136-7f37-48a5-966c-aed9d674ae48 +# ╠═3fc36591-730e-43ad-8bad-f3c23e3c40b9 +# ╠═28cf3957-8a7b-45e4-88f8-e5bcd5d56e3a +# ╠═4761a635-0fed-4407-be3f-cfdbf4fc827e +# ╠═cabbdb5d-1f93-4eeb-8325-a0348a616024 +# ╟─d0f4a46e-4245-4c59-9180-9571bf0cdebf +# ╠═ca57e236-f4b9-401a-b996-d5aaa46d3d71 +# ╠═a22b1ecb-894b-494c-8888-e1b3f1a48f97 +# ╠═fedcd534-120c-41aa-bd60-19734d32597a +# ╠═b59ebf87-4b86-4678-baad-791eaa9c28d6 +# ╠═e4ef052b-8e17-49b1-913a-33e5c32e7436 +# ╠═fd4d69e4-828e-4c95-898b-c70fa003adf4 +# ╠═cecb84d0-fbcf-41ed-898e-c2b603ee0f04 +# ╠═88200bb4-06e4-49aa-8a31-390344f54da6 +# ╠═a8ea0c36-9fa9-4933-be4d-551565758f59 diff --git a/examples/double_slit.jl b/examples/double_slit.jl index e077f9c..7c7f4e2 100644 --- a/examples/double_slit.jl +++ b/examples/double_slit.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.30 +# v0.19.37 using Markdown using InteractiveUtils @@ -91,14 +91,14 @@ W = x[1 + ΔW * 2 + 1] - x[1] md"# Propagate with Fraunhofer Diffraction" # ╔═╡ 7e797b20-77ee-4888-9b23-74bb2713912f -@time output, t = fraunhofer(slit, z, λ, L) +@time output = fraunhofer(slit, z, λ, L); # ╔═╡ 0c74a49c-e4de-4172-ac75-4c2692c505fb # creating this function is more efficient! -efficient_fraunhofer, _ = Fraunhofer(slit, z, λ, L); +efficient_fraunhofer = Fraunhofer(slit, z, λ, L); # ╔═╡ 0e6f74f1-8723-4596-b754-973b253443a4 -@time output2, t2 = efficient_fraunhofer(slit) +@time output2 = efficient_fraunhofer(slit); # ╔═╡ 5c0486d2-5fbd-4cd7-9fc5-583a8e188b2e begin @@ -160,18 +160,18 @@ fwd = let z=z, L=L, λ=λ end # ╔═╡ c9651b37-d911-4230-a2f5-f94e8f726fbc -fwd2 = let z=z, L=L, λ=λ, fr=Fraunhofer(rec0 .+ 0im, z, λ, L)[1] +fwd2 = let z=z, L=L, λ=λ, fr=Fraunhofer(rec0 .+ 0im, z, λ, L) 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])) + f(x) = sum(abs2, intensity_cu .- abs2.(fwd(abs2.(x) .+ 0im))) 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])) + f(x) = sum(abs2, intensity_cu .- abs2.(fwd2(abs2.(x) .+ 0im))) end # ╔═╡ 30555d35-04a8-44b7-afd7-345730d97a61 @@ -244,7 +244,7 @@ simshow(abs2.(Array(slit))) # ╠═2c28b256-f60c-48b6-a34a-d2908c979e30 # ╠═3a19cbfe-6efc-45f7-9aa3-f319f34a534e # ╠═b19225e0-5536-4af1-bbbc-c6124e91536b -# ╠═e0aa0714-08e0-46f4-b6d4-8c5df044ddd9 +# ╟─e0aa0714-08e0-46f4-b6d4-8c5df044ddd9 # ╠═c8a3d8b0-f68d-4c38-ae1a-8081846bda23 # ╠═c9651b37-d911-4230-a2f5-f94e8f726fbc # ╠═e5c3f3dc-5eb3-421e-90ab-a16ea73e5621 From 4404ce27517235bce4460fc5ca1a6133a81368f0 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 19:14:55 +0100 Subject: [PATCH 11/13] Update failing test --- src/angular_spectrum.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/angular_spectrum.jl b/src/angular_spectrum.jl index 74ed3dd..6862f50 100644 --- a/src/angular_spectrum.jl +++ b/src/angular_spectrum.jl @@ -44,9 +44,9 @@ function _prepare_angular_spectrum(field::AbstractArray{CT}, z, λ, _L; # y and x positions in real space, use correct spacing -> fftpos y1 = similar(field, real(eltype(field)), (size(field, 1), 1)) - y1 .= (fftpos(L[1], size(field, 1), CenterFT)) + Zygote.@ignore y1 .= (fftpos(L[1], size(field, 1), CenterFT)) x1 = similar(field, real(eltype(field)), (1, size(field, 2))) - x1 .= (fftpos(L[2], size(field, 2), CenterFT))' + Zygote.@ignore x1 .= (fftpos(L[2], size(field, 2), CenterFT))' params = Params(y1, x1, y1, x1, L, L) From 512ff1aca4a47d8d23c14652109d36516611cd51 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 20:02:27 +0100 Subject: [PATCH 12/13] Remove old code --- src/propagation.jl | 25 ------------------------- 1 file changed, 25 deletions(-) diff --git a/src/propagation.jl b/src/propagation.jl index 8ad5916..590d0fa 100644 --- a/src/propagation.jl +++ b/src/propagation.jl @@ -37,28 +37,3 @@ 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 From 3d9bac36bdfb351de894b651a960dbacd2789664 Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Fri, 1 Mar 2024 20:02:47 +0100 Subject: [PATCH 13/13] Add more tests --- test/angular_spectrum.jl | 31 +++++++++++++++++++++---------- test/fraunhofer.jl | 14 ++++++++++++++ 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/test/angular_spectrum.jl b/test/angular_spectrum.jl index 9d291b0..8eea059 100644 --- a/test/angular_spectrum.jl +++ b/test/angular_spectrum.jl @@ -21,31 +21,42 @@ field = zeros(ComplexF64, (15, 15)) field[5:6, 3:8] .= 1 - gg(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6))) - out1 = gradient(gg, field)[1] - out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + gg3(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), 100e-6, 633e-9, 100e-6))) + out1 = gradient(gg3, field)[1] + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg3, 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)))) - out3 = gradient(f_AS, field)[1] + f_AS3(x) = sum(abs2.(x .- AS(cis.(x)))) + out3 = gradient(f_AS3, field)[1] @test out3 ≈ out1 @test out2 ≈ out1 field = zeros(ComplexF64, (15, 15)) field[5:6, 3:8] .= 1 - gg(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), [100e-6, 200e-6], 633e-9, 100e-6))) - out1 = gradient(gg, field)[1] - out2 = FiniteDifferences.grad(central_fdm(5, 1), gg, field)[1] + gg2(x) = sum(abs2.(x .- WaveOpticsPropagation.angular_spectrum(cis.(x), [100e-6, 200e-6], 633e-9, 100e-6))) + out1 = gradient(gg2, field)[1] + out2 = FiniteDifferences.grad(central_fdm(5, 1), gg2, 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)))) - out3 = gradient(f_AS, field)[1] + f_AS2(x) = sum(abs2.(x .- AS(cis.(x)))) + out3 = gradient(f_AS2, field)[1] @test out3 ≈ out1 @test out2 ≈ out1 end + @testset "Test symmetry" begin + arr = randn(ComplexF32, (4,4)) + arr_ = permutedims(arr, (2,1)) + @test AngularSpectrum(arr, 100e-6, 633e-9, (100e-6, 10e-6))(arr)[:] ≈ permutedims(AngularSpectrum(arr_, 100e-6, 633e-9, (10e-6, 100e-6))(arr_), (2,1))[:] + + arr = randn(ComplexF32, (4,2)) + arr_ = permutedims(arr, (2,1)) + @test AngularSpectrum(arr, 100e-6, 633e-9, (100e-6, 10e-6))(arr)[:] ≈ permutedims(AngularSpectrum(arr_, 100e-6, 633e-9, (10e-6, 100e-6))(arr_), (2,1))[:] + end + + @testset "double slit" begin include("double_slit.jl") end diff --git a/test/fraunhofer.jl b/test/fraunhofer.jl index cfdc830..ab122ab 100644 --- a/test/fraunhofer.jl +++ b/test/fraunhofer.jl @@ -43,4 +43,18 @@ 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] + + @testset "Test symmetry" begin + arr = randn(ComplexF32, (4,4)) + arr_ = permutedims(arr, (2,1)) + @test Fraunhofer(arr, 100e-6, 633e-9, (100e-6, 10e-6))(arr)[:] ≈ permutedims(Fraunhofer(arr_, 100e-6, 633e-9, (10e-6, 100e-6))(arr_), (2,1))[:] + + arr = randn(ComplexF32, (4,2)) + arr_ = permutedims(arr, (2,1)) + @test Fraunhofer(arr, 100e-6, 633e-9, (100e-6, 10e-6))(arr)[:] ≈ permutedims(Fraunhofer(arr_, 100e-6, 633e-9, (10e-6, 100e-6))(arr_), (2,1))[:] + end end + + + +