How do I get power from gaussian beam numerically?

Question

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
        r = radius_sunflower(k,nPoints,b)
        θ = 2*pi*k/ϕ^2
        screen[k,1] = w₀*r*cos(θ)
        screen[k,2] = w₀*r*sin(θ)
        screen[k,3] = 0
    end
    return screen
end

function radius_sunflower(k,n,b)
    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")



Radius = 50/k₀
θ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 ???

Answer 1

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))


Radius = 1000/k₀
θ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
    ΔA = (Radius^2)*sin(   (θ_range[θ]+θ_range[θ+1])/2  )*(θ_range[θ+1]-θ_range[θ])*(φ_range[φ+1]-φ_range[φ])
    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