# How do I get power from gaussian beam numerically?

I would like to get the power from a Gaussian beam given a set of points at which electric field is evaluated. Please follow my reasoning and tell me what assumption maybe are wrong

Power definition is:

$$P = \int |E|^2dA$$

We integrate over area, then Power is a quantity that does not depend on the geometry. If we compute the integral over a cross section of a Gaussian Beam or over the spherical surface, both results should be the same (I think ??)

To test so, I wrote a code to get the electric field on both regions (code below), and I'm here looking for help of how to integrate electric field on both scenarios.

## Gaussian Electric Field

I'm using expression (3.11) from this pdf and not changing the direction of propagation, it is a Gaussian beam of waist $$w_0$$ traveling at $$\hat{z}$$ direction

## Cross Section

Looking at Stackoverflow results I am using Sunflower Seed Arrangements with radius equal to the Gaussian waist

## Sphere

Nothing special, I have matrix of points uniformly distributed in $$\theta$$ and $$\phi$$ spherical coordinates.

Now, my source code wrote in Julia language :

using Plots
plotly()

function E_gaussian(x,y,z, k₀, w₀, E₀)
first_term = exp(im*k₀*z)/(1 + 2*im*z/(k₀*w₀^2))
second_term = exp(  - ((x^2 + y^2)/w₀^2)*(1/(1 + 2*im*z/(k₀*w₀^2)))  )
E = E₀*first_term*second_term
return E
end

function createGaussianSection(nPoints, α, w₀ )
b = round(α*sqrt(nPoints)) # number of boundary points
ϕ = (sqrt(5)+1)/2 # golden ratio
screen = zeros(nPoints,3)
for k=1:nPoints
θ = 2*pi*k/ϕ^2
screen[k,1] = w₀*r*cos(θ)
screen[k,2] = w₀*r*sin(θ)
screen[k,3] = 0
end
return screen
end

if k>n-b
r = 1 # put on the boundary
else
r = sqrt(k-1/2)/sqrt(n-(b+1)/2) # apply square root
end
return r
end

function getE_GassianSection(screen, k₀,  w₀, E₀)
nPoints = size(screen,1)
E_laser = zeros(Complex{Float64}, nPoints )

for k =1:nPoints
xₛ = screen[k,1]
yₛ = screen[k,2]
zₛ = screen[k,3]
E_laser[k] = E_gaussian(xₛ, yₛ, zₛ, k₀,  w₀, E₀)
end
return E_laser
end

function createSphere(Radius, θNintervals, φNintervals, w₀, k₀)
θ_range = range(0, stop=π, length=θNintervals)
φ_range = range(0, stop=2π, length=φNintervals)

xₛ = Radius.*[cos(φ)*sin(θ) for θ in θ_range, φ in φ_range]
yₛ = Radius.*[sin(φ)*sin(θ) for θ in θ_range, φ in φ_range]
zₛ = Radius.*[cos(θ) for θ in θ_range, φ in φ_range]

screen = zeros(size(xₛ,1),size(xₛ,2),3)
screen[:,:,1] = xₛ
screen[:,:,2] = yₛ
screen[:,:,3] = zₛ

return screen, θ_range, φ_range
end

function getE_Sphere(screen, k₀,  w₀, E₀)
nΘ = size(screen,1)
nΦ = size(screen,2)
E_sphere = zeros(Complex{Float64},nΘ, nΦ )

for θ in 1:nΘ, φ in 1:nΦ
xₛ = screen[θ,φ,1]
yₛ = screen[θ,φ,2]
zₛ = screen[θ,φ,3]
E_sphere[θ,φ] = E_gaussian(xₛ, yₛ, zₛ, k₀,  w₀, E₀)
end
return E_sphere
end

k₀ =  1
w₀ = 30/k₀
E₀ = 1

screenGaussin = createGaussianSection(5000, 0, w₀)
E_section = getE_GassianSection(screenGaussin, k₀,  w₀, E₀)
Intensity_section = real.( conj.(E_section).*E_section)
scatter(screenGaussin[:,1], screenGaussin[:,2], zcolor=Intensity_section, axis=:equal, label="", xlabel="x", ylabel="y")

θNintervals = φNintervals = 25
screen, θ_range, φ_range = createSphere(Radius, θNintervals, φNintervals, w₀, k₀)
E_sphere = getE_Sphere(screen, k₀,  w₀, E₀)
Intensity_sphere = real.( conj.(E_sphere).*E_sphere)
scatter3d(screen[:,:,1],screen[:,:,2],screen[:,:,3], zcolor =Intensity_sphere , label="")


The results are the following, the color is the local intensity:

From everything that I've shown, I only have Points. How do I find Areas and integrate them ???

Thank you for people on inbox to help me.

1. I don't need full sphere, only half of sphere contains the same energy of the cross section.

2. I can use my point grid and imagine a small area between them. I set a point and create a square based upon the nearest neighbors.

2.1 In spherical coordinates, an element of area is $$dA = r^2sin(\theta)d\theta d\phi$$. To covert this expression into discrete case, we use $$d\theta = \theta[i+1]-\theta[i]$$ and $$d\phi= \phi[i+1]-\phi[i]$$, and also $$\theta$$ is the average of the nearest neighbors $$\theta = (\theta[i+1]+\theta[i])/2$$

2.2 For simplicity, the cross section area can be just a square grid of points

3. This procedure works only for far field, that is, the radius of the sphere should be big if compared with the Gaussian waist ( I've used a wrong choice of values in my example)

For everybody that is interested, here the full working code:

function E_gaussian(x,y,z, k₀, w₀, E₀)
first_term = exp(im*k₀*z)/(1 + 2*im*z/(k₀*w₀^2))
second_term = exp(  - ((x^2 + y^2)/w₀^2)*(1/(1 + 2*im*z/(k₀*w₀^2)))  )
E = E₀*first_term*second_term
return E
end

function createGaussianSection(nPoints, w₀ )
x_range = range(-w₀, stop=w₀, length=nPoints)
y_range = range(-w₀, stop=w₀, length=nPoints)
screen = zeros(nPoints, nPoints, 3)
screen[:,:,1] = [x for x in x_range, y in y_range]
screen[:,:,2] = [y for x in x_range, y in y_range]
screen[:,:,3] .= 0
return screen, x_range, y_range
end

function createSphere(Radius, θNintervals, φNintervals, w₀, k₀)
θ_range = range(0, stop=π/2, length=θNintervals)
φ_range = range(0, stop=2π, length=φNintervals)

xₛ = Radius.*[cos(φ)*sin(θ) for θ in θ_range, φ in φ_range]
yₛ = Radius.*[sin(φ)*sin(θ) for θ in θ_range, φ in φ_range]
zₛ = Radius.*[cos(θ) for θ in θ_range, φ in φ_range]

screen = zeros(size(xₛ,1),size(xₛ,2),3)
screen[:,:,1] = xₛ
screen[:,:,2] = yₛ
screen[:,:,3] = zₛ

return screen, θ_range, φ_range
end

function getE(screen, k₀,  w₀, E₀)
nΘ = size(screen,1)
nΦ = size(screen,2)
E_sphere = zeros(Complex{Float64},nΘ, nΦ )

for θ in 1:nΘ, φ in 1:nΦ
xₛ = screen[θ,φ,1]
yₛ = screen[θ,φ,2]
zₛ = screen[θ,φ,3]
E_sphere[θ,φ] = E_gaussian(xₛ, yₛ, zₛ, k₀,  w₀, E₀)
end
return E_sphere
end

k₀ =  1
w₀ = 30/k₀
E₀ = 1

screenGaussin, x_range, y_range = createGaussianSection(150, 2*w₀)
E_section = getE(screenGaussin, k₀,  w₀, E₀)
Intensity_section = real.( conj.(E_section).*E_section)

I_section = 0
for x = 1:(length(x_range)-1), y = 1:(length(y_range)-1)
global I_section
E_mean = (E_section[x, y] + E_section[x+1, y] + E_section[x, y+1]+E_section[x+1, y+1])/4
ΔA = (x_range[x+1]-x_range[x])*(y_range[y+1]-y_range[y])
I_section += real(conj(E_mean)*E_mean)*ΔA
end
println("Total Power Section: ", round(I_section,digits=3))

θNintervals = φNintervals = 200
screen, θ_range, φ_range = createSphere(Radius, θNintervals, φNintervals, w₀, k₀)
E_sphere = getE(screen, k₀,  w₀, E₀)
Intensity_sphere = real.( conj.(E_sphere).*E_sphere)

I_sphere = 0
for θ = 1:(θNintervals-1), φ = 1:(φNintervals-1)
global I_sphere
E_medio = (E_sphere[θ, φ] + E_sphere[θ+1, φ] + E_sphere[θ, φ+1]+E_sphere[θ+1, φ+1])/4
I_sphere += real(conj(E_medio)*E_medio)*ΔA
end
println("Total Power Sphere: ", round(I_sphere,digits=3))

println("Relation:",round(I_section/I_sphere,digits=3))


Results:

Total Power Section: 1413.028
Total Power Sphere: 1413.567
Relation:1.0