レーザービーム品質(M²)の解説と高次ガウシアンモードのM²

using SpecialPolynomials
using Polynomials
using Latexify
using LaTeXStrings
using PyPlot
using CSV
using DataFrames
using QuadGK

# calculation beam size depemd on M2
function beam_waist_calcuation(M2, rin, z)
    # this calculation is considered as High order Gaussian beam
    # We can not use q-prameter because q-parameter is only for 1st order Gaussian beam

    # unit is um 
    λ = 1
    w = rin*sqrt(1+(M2*λ*z/(π*rin^2))^2) # beam size
    return w
end

z_list = range(0, 100, length=100)
rin = 2
w_1 = map(z -> beam_waist_calcuation(1, rin, z), z_list)
w_2 = map(z -> beam_waist_calcuation(1.5, rin, z) , z_list) 
figure()
plot(z_list, w_1, label="M2 = 1")
plot(z_list, w_2, label="M2 = 1.5")

text(100,w_1[end], L"\theta", fontsize=12)
text(100,w_2[end], L"M^2\theta", fontsize=12)
xlabel("z [um]")
ylabel("w [um]")
title("Beam size depend on M2 - when rin size is same")
ylim(0)
legend()
grid()
display(gcf())

z_list = range(0, 150, length=100)
w_1 = map(z -> beam_waist_calcuation(1, 2, z), z_list)
w_2 = map(z -> beam_waist_calcuation(1.5, 2*sqrt(1.5^2), z), z_list)
figure()
plot(z_list, w_1, label="M2 = 1")
plot(z_list, w_2, label="M2 = 1.5")
text(0,w_1[1], L"w", fontsize=12)
text(0,w_2[1], L"Mw", fontsize=12)
xlabel("z [um]")
ylabel("w [um]")
title("Beam size depend on M2 - when divergence angle is same")
ylim(0)
legend()
grid()
display(gcf())

function Hermite_gaussian(x , y, z, l, m)
    # parameters
    k = 1
    z₀ = 0
    ω₀ = 1
    η  = atan(-(z-z₀)/(k*ω₀^2))
    E₀ = 1
    ω_z = ω₀/cos(η) 

    Eₗₘ = ( 
            E₀ * (ω₀/ω_z) * basis(Hermite, l)(√2*x/ω_z) * basis(Hermite, m)(√2*y/ω_z) 
            * exp(complex( -(x^2 + y^2)/(ω_z^2), -k*z+(l+m+1)*η) ) 
        )
    return Eₗₘ
end

function Laguerre_gaussian(x,y,z,p,m)
    # cluculation cylinderical coordinates
    r = sqrt(x^2 + y^2)
    ϕ = atan(y, x)

    # parameters
    k = 1
    z₀ = 0
    ω₀ = 1
    η  = atan(-(z-z₀)/(k*ω₀^2))
    E₀ = 1
    ω_z = ω₀/cos(η)

    Eₚₘ = ( 
            E₀ * (ω₀/ω_z) * (√2r/ω_z)^abs(m) * basis(Laguerre{m}, p)(√2*r/ω_z)  
            * exp(complex( -r^2/(ω_z^2), -k*z + m*ϕ + (2p + abs(m) + 1)*η) ) 
        )
    return Eₚₘ
end

# D4σの計算 at julia my cord
function D4σ_switch(z, l, m, x₀, y₀, mode)
    if mode == "Hermite"
        E = Hermite_gaussian
    elseif mode == "Laguerre"
        E = Laguerre_gaussian
    else
        error("mode must be 'Hermite' or 'Laguerre'")
    end
    I(x, y) = abs2(E(x, y, l, m, z))

    σ_x_numerator = quadgk(y -> 
        quadgk(x -> (x - x₀)^2 * I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_x_denominator = quadgk(y -> 
        quadgk(x -> I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_x = √( σ_x_numerator / σ_x_denominator )

    σ_y_numerator = quadgk(x -> 
        quadgk(y -> (y - y₀)^2 * I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_y_denominator = quadgk(x-> 
        quadgk(y -> I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_y = √( σ_y_numerator / σ_y_denominator )

    return 2 * σ_x, 2 * σ_y
end

x_max = y_max = 3
x_vals = range(-x_max, x_max, length = 100)
y_vals = range(-y_max, y_max, length = 100)

figure(figsize=(12, 12))

# Calculate Hermite_gaussian 0,0 to 3,3
for l in 0:3
    for m in 0:3
        Elm = zeros(Complex{Float64}, length(x_vals), length(y_vals))
        for i in 1:length(x_vals)
            for j in 1:length(y_vals)
                Elm[i,j] = Hermite_gaussian(x_vals[i], y_vals[j], 0, l, m)
            end
        end

        # Plotting
        subplot(4,4,l*4+m+1)
        title("l is $l, m is $m")
        imshow(abs2.(Elm), extent = (x_vals[1], x_vals[end], y_vals[1], y_vals[end]), cmap = "hot")
    end
end

display(gcf())

# Calculate Laguerre_gaussian 0,0 to 3,3
for l in 0:3
    for m in 0:3
        Elm = zeros(Complex{Float64}, length(x_vals), length(y_vals))
        for i in 1:length(x_vals)
            for j in 1:length(y_vals)
                Elm[i,j] = Laguerre_gaussian(x_vals[i], y_vals[j], 0, l, m)
            end
        end

        # Plotting
        subplot(4,4,l*4+m+1)
        title("p is $l, m is $m")
        imshow(abs2.(Elm), extent = (x_vals[1], x_vals[end], y_vals[1], y_vals[end]), cmap = "hot")
    end
end

display(gcf())

# 積分関連のjulia 参考文献：https://zenn.dev/ohno/articles/440234fbb2adec#%E3%83%91%E3%83%83%E3%82%B1%E3%83%BC%E3%82%B8
#using Pkg
#Pkg.add("QuadGK")

# calculation D4σ
using QuadGK
using SpecialPolynomials
using Polynomials

function Eₗₘ(x , y, z, l, m)
    k = 1
    z₀ = 0
    ω₀ = 1
    ϕ  = atan(-(z-z₀)/(k*ω₀^2))
    E₀ = 1
    ω_z = ω₀/cos(ϕ) 

    Eₗₘ = ( 
            E₀ * (ω₀/ω_z) * basis(Hermite, l)(√2*x/ω_z) * basis(Hermite, m)(√2*y/ω_z) 
            * exp(complex( -(x^2 + y^2)/(ω_z^2), -k*z+(l+m+1)*ϕ) ) 
        )
    return Eₗₘ
end

# D4σの計算 at julia my cord
function D4σ(z, l, m, x₀, y₀)
    I(x, y) = abs2(Eₗₘ(x, y, l, m, z))

    σ_x_numerator = quadgk(y -> 
        quadgk(x -> (x - x₀)^2 * I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_x_denominator = quadgk(y -> 
        quadgk(x -> I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_x = √( σ_x_numerator / σ_x_denominator )

    σ_y_numerator = quadgk(x -> 
        quadgk(y -> (y - y₀)^2 * I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_y_denominator = quadgk(x-> 
        quadgk(y -> I(x, y), -Inf, Inf; maxevals=10^3)[1],
        -Inf, Inf; maxevals=10^3)[1]

    σ_y = √( σ_y_numerator / σ_y_denominator )

    return 2 * σ_x, 2 * σ_y
end

function e²_width(x_array,array_1d)
    if length(x_array) != length(array_1d)
        error("x_array and array_1d must have the same length")
    end
    max_value = maximum(array_1d)
    threshold = max_value /exp(2)
    indices = findall(x -> x > threshold, array_1d)
    w = x_array[indices[end]] - x_array[indices[1]]
    return w
end

ω_D4σ = D4σ(0, 0, 0, 0, 0) # l=0, m=0
print("ω_D4σ = ", ω_D4σ[1], "\n")

x_array = range(-3, 3, length=100)
E_lm = map(x -> Eₗₘ(x, 0, 0, 0, 0), x_array)
I = abs2.(E_lm)
w = e²_width(x_array,I)
print("w = ", w/2, "\n")

using CSV
using DataFrames

# Base width
ωₓ_base, ωᵧ_base = D4σ_switch(0, 0, 0, 0, 0, "Hermite")

# 結果を格納するDataFrameを作成
results = DataFrame(Type = String[], l_or_p = Int[], m = Int[], M2_x = Float64[], M2_y = Float64[])

# Hermite-Gaussian
for l in 0:3
    for m in 0:3
        ωₓ, ωᵧ = D4σ_switch(0, l, m, 0, 0, "Hermite")
        M2_x = (ωₓ / ωₓ_base)
        M2_y = (ωᵧ / ωᵧ_base)
        push!(results, ("Hermite", l, m, M2_x, M2_y))
    end
end

# Laguerre-Gaussian
for p in 0:3
    for m in 0:3
        ωₓ, ωᵧ = D4σ_switch(0, p, m, 0, 0, "Laguerre")
        M2_x = (ωₓ / ωₓ_base)
        M2_y = (ωᵧ / ωᵧ_base)
        push!(results, ("Laguerre", p, m, M2_x, M2_y))
    end
end

# CSVファイルに保存
CSV.write("M2_parameters.csv", results)