Skip to content

Commit

Permalink
Update v0.4.0
Browse files Browse the repository at this point in the history 
Former-commit-id: 98cb3b6
  • Loading branch information
hochunlin committed Dec 17, 2023
1 parent a849307 commit d72169b
Show file tree
Hide file tree
Showing 27 changed files with 482 additions and 190 deletions.
68 changes: 40 additions & 28 deletions ...njugated_light_through_disorder/README.md → ...disorder_with_phase_conjugation/README.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Focusing Phase Conjugated Light Through Disorder
# Focusing Inside Disorder With Phase Conjugation

In this example, we show how to use mesti() to compute the field profile of a point source in a scattering disordered medium, do phase conjugation to determine an incident wavefront that can focus on the the disorder, and then use mesti2s() again to compute the field profile to show its focus.
In this example, we show how to use mesti() to project field generated from a point source inside disorder onto propagating channels through APF method, do phase conjugation to determine an incident wavefront that can focus inside the disorder, and then use mesti2s() again to compute the field profile to show its focus.

```julia
# Call necessary packages
Expand Down Expand Up @@ -43,7 +43,7 @@ build_epsilon_disorder(W, L, r_min, r_max, min_sep,
no_scatterer_center = true)
```

# Compute the field profile of a point source
# Projecting field generated from a point source onto propagating channels through APF

```julia
syst = Syst()
Expand All @@ -63,28 +63,48 @@ Bx.pos = [[m0_focus,l0_focus+pml_npixels+1,1,1]]
Bx.data = [ones(1,1)]

# put PML along z-direction
pml = get_optimal_PML(syst.wavelength/syst.dx)
pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)
pml.npixels = pml_npixels
pml.direction = "z"
syst.PML = [pml]

# field profile: input from a point source
Ex_field, _ = mesti(syst, [Bx])
# build channels for the low side (air)
channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)
N_prop_low = channels_low.N_prop # number of propagating channels on the low side

# build projection profiles C on the low side
Cx = Source_struct()
ny_Ex = Int(W/dx)
Cx.pos = [[1,pml_npixels+1,ny_Ex,1]]
# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx)
# here we project the result onto all propagating channels as a output basis
C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0
Cx.data = [C_low]

# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF
w, _ = mesti(syst, [Bx], [Cx])

# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles.
# It would generate the same output channel amplitude coefficients
# Ex_field, _ = mesti(syst, [Bx])
# C = transpose(C_low)
# w = C_low*Ex_field[:,pml_npixels+1]
```
```text:Output
===System size===
ny_Ex = 5400; nz_Ex = 1381 for Ex(y,z)
UPML on -z +z sides; ; yBC = periodic; zBC = PEC
Building B,C... elapsed time: 0.230 secs
Building A ... elapsed time: 9.228 secs
< Method: factorize_and_solve using MUMPS in single precision with AMD ordering >
Analyzing ... elapsed time: 3.977 secs
Factorizing ... elapsed time: 64.983 secs
Solving ... elapsed time: 4.065 secs
Total elapsed time: 76.477 secs
Building B,C... elapsed time: 1.664 secs
Building A ... elapsed time: 8.106 secs
< Method: APF using MUMPS in single precision with AMD ordering >
Building K ... elapsed time: 2.822 secs
false
Analyzing ... elapsed time: 5.076 secs
Factorizing ... elapsed time: 76.478 secs
Total elapsed time: 99.805 secs
```

# Compute the regular focusing and phase-conjugated focusing profile
# Compute the regular focusing and phase-conjugated focusing profiles

```julia
# Specify the system for mesti2s() and mesti_build_channels()
Expand All @@ -101,19 +121,12 @@ syst.zPML = [pml]
# equivalent average epsilon for this disordered system
epsilon_ave = mean(epsilon)

# build channels for the equivalent average epsilon and low side (air)
# build channels for the equivalent average epsilon
channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave)
channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)
N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon
N_prop_low = channels_low.N_prop # number of propagating channels on the low side

dn = 0.5
# regular focus wavefront
wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn))

# build projection matrix C on the low side
C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop)))
proj_coefficient = C_low*Ex_field[:,pml_npixels+1]
wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0

# specify two input incident wavefronts:
# (1) regular focusing wavefront
Expand All @@ -123,12 +136,11 @@ input.v_low = zeros(ComplexF64, N_prop_low, 2)
input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)])

# for the phased conjugated input:
# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)
# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)
# perm() means permute a vector that switches one propagating channel with one
# having a complex-conjugated transverse profile
# for periodic boundary, this flips the sign of ky.
input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient)

# for the periodic boundary, this flips the sign of ky.
input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w)

# we will also get the field profile in the free spaces on the two sides, for
# plotting purpose.
Expand Down Expand Up @@ -252,4 +264,4 @@ println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(a
```
```text:Output
I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = 1800
```
```
0 ...hrough_disorder/build_epsilon_disorder.jl → ...ase_conjugation/build_epsilon_disorder.jl
View file
File renamed without changes.
48 changes: 30 additions & 18 deletions ...e_conjugated_light_through_disorder.ipynb → ...ide_disorder_with_phase_conjugation.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,15 @@
"id": "ae0e7aae-9838-44df-9330-8a1d8f9fc85d",
"metadata": {},
"source": [
"# Focusing Phase Conjugated Light Through Disorder"
"# Focusing Inside Disorder With Phase Conjugation"
]
},
{
"cell_type": "markdown",
"id": "96a5a5eb-94ba-49be-9f3f-3e6d71cea482",
"metadata": {},
"source": [
"In this example, we show how to use mesti() to compute the field profile of a point source in a scattering disordered medium, do phase conjugation to determine an incident wavefront that can focus on the the disorder, and then use mesti2s() again to compute the field profile to show its focus."
"In this example, we show how to use mesti() to project field generated from a point source inside disorder onto propagating channels through APF method, do phase conjugation to determine an incident wavefront that can focus inside the disorder, and then use mesti2s() again to compute the field profile to show its focus."
]
},
{
Expand Down Expand Up @@ -80,7 +80,7 @@
"id": "be175504-fd46-4922-9262-527fa6389963",
"metadata": {},
"source": [
"# Compute the field profile of a point source "
"# Projecting field generated from a point source onto propagating channels through APF"
]
},
{
Expand All @@ -107,21 +107,40 @@
"Bx.data = [ones(1,1)]\n",
"\n",
"# put PML along z-direction\n",
"pml = get_optimal_PML(syst.wavelength/syst.dx)\n",
"pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)\n",
"pml.npixels = pml_npixels\n",
"pml.direction = \"z\"\n",
"syst.PML = [pml]\n",
"\n",
"# field profile: input from a point source\n",
"Ex_field, _ = mesti(syst, [Bx]);"
"# build channels for the low side (air)\n",
"channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)\n",
"N_prop_low = channels_low.N_prop # number of propagating channels on the low side \n",
"\n",
"# build projection profiles C on the low side\n",
"Cx = Source_struct()\n",
"ny_Ex = Int(W/dx)\n",
"Cx.pos = [[1,pml_npixels+1,ny_Ex,1]]\n",
"# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx)\n",
"# here we project the result onto all propagating channels as a output basis \n",
"C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0\n",
"Cx.data = [C_low]\n",
"\n",
"# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF\n",
"w, _ = mesti(syst, [Bx], [Cx])\n",
"\n",
"# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles. \n",
"# It would generate the same output channel amplitude coefficients\n",
"# Ex_field, _ = mesti(syst, [Bx])\n",
"# C = transpose(C_low)\n",
"# w = C_low*Ex_field[:,pml_npixels+1]"
]
},
{
"cell_type": "markdown",
"id": "21c79756-b0be-4e62-b998-ce405190a866",
"metadata": {},
"source": [
"# Compute the field profile of a point source "
"# Compute the regular focusing and phase-conjugated focusing profiles"
]
},
{
Expand All @@ -145,19 +164,12 @@
"# equivalent average epsilon for this disordered system\n",
"epsilon_ave = mean(epsilon)\n",
"\n",
"# build channels for the equivalent average epsilon and low side (air)\n",
"# build channels for the equivalent average epsilon\n",
"channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave)\n",
"channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)\n",
"N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon \n",
"N_prop_low = channels_low.N_prop # number of propagating channels on the low side \n",
"\n",
"dn = 0.5\n",
"# regular focus wavefront\n",
"wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn))\n",
"\n",
"# build projection matrix C on the low side\n",
"C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop)))\n",
"proj_coefficient = C_low*Ex_field[:,pml_npixels+1]\n",
"wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0\n",
"\n",
"# specify two input incident wavefronts:\n",
"# (1) regular focusing wavefront\n",
Expand All @@ -167,11 +179,11 @@
"input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)])\n",
"\n",
"# for the phased conjugated input: \n",
"# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)\n",
"# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)\n",
"# perm() means permute a vector that switches one propagating channel with one\n",
"# having a complex-conjugated transverse profile\n",
"# for the periodic boundary, this flips the sign of ky. \n",
"input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient)\n",
"input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w)\n",
"\n",
"# we will also get the field profile in the free spaces on the two sides, for\n",
"# plotting purpose.\n",
Expand Down
60 changes: 36 additions & 24 deletions ...hase_conjugated_light_through_disorder.jl → ...inside_disorder_with_phase_conjugation.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
## Focusing Phase Conjugated Light Through Disorder
## Focusing Inside Disorder With Phase Conjugation

# In this example, we show how to use mesti() to compute the field profile of
# a point source in a scattering disordered medium, do phase conjugation to
# determine an incident wavefront that can focus on the the disorder, and then
# use mesti2s() again to compute the field profile to show its focus.
# In this example, we show how to use mesti() to project field generated from
# a point source inside disorder onto propagating channels through APF method,
# do phase conjugation to determine an incident wavefront that can focus inside
# the disorder, and then use mesti2s() again to compute the field profile to
# show its focus.

# Call necessary packages
using MESTI, GeometryPrimitives, LinearAlgebra, Statistics, Printf
Expand Down Expand Up @@ -42,7 +43,7 @@ build_epsilon_disorder(W, L, r_min, r_max, min_sep,
epsilon_scat, epsilon_bg, build_TM;
no_scatterer_center = true)

## Compute the field profile of a point source
## Projecting field generated from a point source onto propagating channels through APF
syst = Syst()
pml_npixels = 15
syst.length_unit = "lambda_0"
Expand All @@ -60,15 +61,34 @@ Bx.pos = [[m0_focus,l0_focus+pml_npixels+1,1,1]]
Bx.data = [ones(1,1)]

# put PML along z-direction
pml = get_optimal_PML(syst.wavelength/syst.dx)
pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)
pml.npixels = pml_npixels
pml.direction = "z"
syst.PML = [pml]

# field profile: input from a point source
Ex_field, _ = mesti(syst, [Bx])

## Compute the regular focusing and phase-conjugated focusing profile
# build channels for the low side (air)
channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)
N_prop_low = channels_low.N_prop # number of propagating channels on the low side

# build projection profiles C on the low side
Cx = Source_struct()
ny_Ex = Int(W/dx)
Cx.pos = [[1,pml_npixels+1,ny_Ex,1]]
# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx)
# here we project the result onto all propagating channels as a output basis
C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0
Cx.data = [C_low]

# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF
w, _ = mesti(syst, [Bx], [Cx])

# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles.
# It would generate the same output channel amplitude coefficients
# Ex_field, _ = mesti(syst, [Bx])
# C = transpose(C_low)
# w = C_low*Ex_field[:,pml_npixels+1]

## Compute the regular focusing and phase-conjugated focusing profiles
# Specify the system for mesti2s() and mesti_build_channels()
syst = Syst()
syst.epsilon_xx = epsilon
Expand All @@ -83,19 +103,12 @@ syst.zPML = [pml]
# equivalent average epsilon for this disordered system
epsilon_ave = mean(epsilon)

# build channels for the equivalent average epsilon and low side (air)
# build channels for the equivalent average epsilon
channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave)
channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)
N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon
N_prop_low = channels_low.N_prop # number of propagating channels on the low side

dn = 0.5
# regular focus wavefront
wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn))

# build projection matrix C on the low side
C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop)))
proj_coefficient = C_low*Ex_field[:,pml_npixels+1]
wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0

# specify two input incident wavefronts:
# (1) regular focusing wavefront
Expand All @@ -105,12 +118,11 @@ input.v_low = zeros(ComplexF64, N_prop_low, 2)
input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)])

# for the phased conjugated input:
# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)
# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)
# perm() means permute a vector that switches one propagating channel with one
# having a complex-conjugated transverse profile
# for the periodic boundary, this flips the sign of ky.
input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient)

input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w)

# we will also get the field profile in the free spaces on the two sides, for
# plotting purpose.
Expand Down Expand Up @@ -202,4 +214,4 @@ intensity_plot = plot(plt3, plt4, layout = @layout([a b]), size=(800, 400))
display(intensity_plot)

# compare the ratio of intensity on the focus point
println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(abs.(Ex[m0_focus,opts.nz_low+l0_focus,2]).^2/abs.(Ex[m0_focus,opts.nz_low+l0_focus,1]).^2, digits=-2)))
println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(abs.(Ex[m0_focus,opts.nz_low+l0_focus,2]).^2/abs.(Ex[m0_focus,opts.nz_low+l0_focus,1]).^2, digits=-2)))
0 ...through_disorder/intensity_comparison.png → ...hase_conjugation/intensity_comparison.png
View file
File renamed without changes
0 ...se_conjugated_focusing.gif.REMOVED.git-id → ...se_conjugated_focusing.gif.REMOVED.git-id
View file
File renamed without changes.
0 ...order/regular_focusing.gif.REMOVED.git-id → ...ation/regular_focusing.gif.REMOVED.git-id
View file
File renamed without changes.
Loading

0 comments on commit d72169b

Please sign in to comment.