9.5 Lampestyrke¶
In [2]:
using LinearAlgebra
using Plots
In [3]:
Lx = 20
Ly = 20
h = 4
nx = 200
ny = 200
# --- Grid over the floor ---
x_vals = range(0, Lx, length=nx)
y_vals = range(0, Ly, length=ny)
for i in y_vals[1:10]
println(i)
end
0.0 0.10050251256281408 0.20100502512562815 0.3015075376884422 0.4020100502512563 0.5025125628140703 0.6030150753768844 0.7035175879396985 0.8040201005025126 0.9045226130653267
In [28]:
lamp_xy = [
(0, 3.0),
(0, 16.0),
(9.0, 5.0),
(14, 10.0),
(14.0, 14),
(20, 3),
(20, 8),
(13, 20),
(20, 19),
(7, 12)
]
Out[28]:
10-element Vector{Tuple{Real, Real}}:
(0, 3.0)
(0, 16.0)
(9.0, 5.0)
(14, 10.0)
(14.0, 14)
(20, 3)
(20, 8)
(13, 20)
(20, 19)
(7, 12)
In [29]:
function lamp_basis(x, y, lamp; h)
lx, ly = lamp
r2 = (x - lx)^2 + (y - ly)^2 + h^2
return h / r2^(1.5)
end;
In [30]:
grid_points = vec([(x, y) for y in y_vals, x in x_vals]);
println(size(grid_points))
for i in grid_points[1:10]
println(i)
end
(40000,) (0.0, 0.0) (0.0, 0.10050251256281408) (0.0, 0.20100502512562815) (0.0, 0.3015075376884422) (0.0, 0.4020100502512563) (0.0, 0.5025125628140703) (0.0, 0.6030150753768844) (0.0, 0.7035175879396985) (0.0, 0.8040201005025126) (0.0, 0.9045226130653267)
In [31]:
# --- Column vector for one lamp ---
function lamp_column(lamp, grid_points; h=4.0)
[lamp_basis(x, y, lamp; h=h) for (x, y) in grid_points]
end;
In [32]:
# --- Build design matrix A by column stacking ---
columns = [lamp_column(lamp, grid_points; h=h) for lamp in lamp_xy]
println(size(columns))
A = hcat(columns...);
println(size(A))
A
(10,) (40000, 10)
Out[32]:
40000×10 Matrix{Float64}:
0.032 0.000891675 0.00296839 … 0.000184684 0.00132386
0.0331731 0.000907675 0.00300507 0.00018605 0.00134701
0.0343759 0.000924055 0.00304175 0.000187426 0.00137064
0.0356068 0.000940824 0.0030784 0.000188812 0.00139477
0.0368639 0.000957995 0.00311499 0.000190207 0.00141939
0.0381446 0.000975578 0.00315149 … 0.000191611 0.00144453
0.0394462 0.000993588 0.00318789 0.000193026 0.00147018
0.0407654 0.00101204 0.00322414 0.00019445 0.00149637
0.0420984 0.00103094 0.00326022 0.000195883 0.00152311
0.0434412 0.0010503 0.00329611 0.000197327 0.0015504
0.0447888 0.00107014 0.00333176 … 0.00019878 0.00157826
0.0461361 0.00109048 0.00336715 0.000200242 0.00160669
0.0474774 0.00111133 0.00340225 0.000201714 0.00163572
⋮ ⋱
0.000231354 0.000457542 0.00066708 0.0624348 0.00112806
0.000229701 0.000456586 0.000658672 0.0624999 0.00111797
0.000228056 0.000455601 0.00065038 … 0.0624466 0.0011079
0.000226422 0.000454588 0.000642203 0.0622756 0.00109783
0.000224796 0.000453546 0.00063414 0.0619885 0.00108778
0.00022318 0.000452476 0.000626189 0.0615878 0.00107775
0.000221574 0.000451378 0.000618348 0.0610774 0.00106773
0.000219977 0.000450253 0.000610618 … 0.0604618 0.00105773
0.00021839 0.000449101 0.000602995 0.0597465 0.00104776
0.000216813 0.000447922 0.000595479 0.0589377 0.00103782
0.000215245 0.000446717 0.000588068 0.0580421 0.00102791
0.000213686 0.000445486 0.000580761 0.0570672 0.00101803
In [9]:
strengths = ones(length(lamp_xy))
i_vec = A * strengths
Out[9]:
40000-element Vector{Float64}:
0.039699385135698206
0.040967572698400234
0.04226643087753457
0.043594329641325695
0.04494926844143329
0.04632884724096554
0.04773023934303014
0.04915016698251067
0.0505848807653037
0.05203014414677238
0.053481224225529096
0.054932890179514374
0.056379420677536464
⋮
0.08451574690607658
0.08434746887701766
0.08406278206201762
0.08366206605112116
0.08314681024853321
0.08251957923797287
0.08178395540813968
0.08094446106023306
0.08000646294101313
0.07897606269493243
0.07785997708331803
0.07666541196871449
In [10]:
I_2d = reshape(i_vec, ny, nx) # reshape back to 2D (x,y order
Out[10]:
200×200 Matrix{Float64}:
0.0396994 0.0397899 0.0398438 … 0.0450933 0.0450617 0.0449906
0.0409676 0.041059 0.0411116 0.0465391 0.0465093 0.0464376
0.0422664 0.0423588 0.0424098 0.0480228 0.0479949 0.0479225
0.0435943 0.0436876 0.0437369 0.0495429 0.0495171 0.049444
0.0449493 0.0450433 0.0450909 0.0510979 0.0510743 0.0510006
0.0463288 0.0464237 0.0464693 … 0.0526856 0.0526644 0.05259
0.0477302 0.0478259 0.0478694 0.0543037 0.054285 0.05421
0.0491502 0.0492466 0.0492878 0.0559493 0.0559332 0.0558576
0.0505849 0.0506821 0.0507209 0.0576191 0.0576056 0.0575295
0.0520301 0.052128 0.0521644 0.0593092 0.0592985 0.0592219
0.0534812 0.0535798 0.0536136 … 0.0610154 0.0610077 0.0609304
0.0549329 0.0550322 0.0550633 0.062733 0.0627282 0.0626503
0.0563794 0.0564794 0.0565078 0.0644568 0.064455 0.0643765
⋮ ⋱
0.0428107 0.042939 0.0430289 0.085209 0.0849146 0.0845157
0.0415318 0.0416588 0.0417498 0.0850314 0.0847418 0.0843475
0.0402851 0.0404108 0.0405027 … 0.0847381 0.0844528 0.0840628
0.0390716 0.0391961 0.0392888 0.0843296 0.0840479 0.0836621
0.037892 0.0380151 0.0381085 0.0838071 0.0835286 0.0831468
0.0367467 0.0368685 0.0369625 0.0831734 0.0828976 0.0825196
0.0356359 0.0357565 0.0358509 0.0824319 0.0821583 0.081784
0.0345597 0.0346789 0.0347737 … 0.081587 0.0813152 0.0809445
0.0335179 0.0336358 0.0337309 0.0806442 0.0803737 0.0800065
0.0325102 0.0326268 0.032722 0.0796094 0.0793399 0.0789761
0.0315362 0.0316515 0.0317467 0.0784893 0.0782206 0.07786
0.0305953 0.0307092 0.0308045 0.077291 0.0770228 0.0766654
In [12]:
# --- Normalize and plot ---
I_2d ./= maximum(I_2d)
heatmap(x_vals, y_vals, I_2d;
aspect_ratio = 1,
xlabel = "x [m]", ylabel = "y [m]",
xlims = (0, 20), ylims = (0, 20),
colorbar_title = "Relative illumination",
c = :viridis,
title = "Illumination with equal lamp strengths")
scatter!([p[1] for p in lamp_xy], [p[2] for p in lamp_xy];
m = :o, ms = 4, mc = :yellow, label = "lamps")
Out[12]:
Desired uniform illumination¶
In [13]:
# --- Desired uniform illumination ---
target = ones(nx*ny) # all ones → uniform brightness target
# --- Least squares solution for lamp strengths ---
strengths_opt = A \ target
# --- Compute optimized illumination map ---
i_vec_opt = A * strengths_opt
i_opt = reshape(i_vec_opt, nx, ny)
# --- Print lamp strengths ---
for (i, s) in enumerate(strengths_opt)
println("Lamp $i: strength = $(round(s, digits=2))")
end
Lamp 1: strength = 22.26 Lamp 2: strength = 22.62 Lamp 3: strength = 17.53 Lamp 4: strength = 6.61 Lamp 5: strength = 7.47 Lamp 6: strength = 18.37 Lamp 7: strength = 10.95 Lamp 8: strength = 15.85 Lamp 9: strength = 15.12 Lamp 10: strength = 12.89
In [14]:
# --- Normalize and plot ---
i_opt ./= maximum(i_opt)
heatmap(x_vals, y_vals, i_opt;
aspect_ratio = 1,
xlabel = "x [m]", ylabel = "y [m]",
xlims = (0, 20), ylims = (0, 20),
colorbar_title = "Relative illumination",
c = :viridis,
title = "Illumination with optimal lamp strengths")
scatter!([p[1] for p in lamp_xy], [p[2] for p in lamp_xy];
m = :o, ms = 4, mc = :yellow, label = "lamps")
Out[14]:
In [35]:
# Flatten the matrix with vec() and plot the histogram
histogram(vec(I_2d),
title = "Equal lamp strength",
xlabel = "Relative illumination",
ylabel = "Frequency",
legend = false,
bins = 5 # You can optionally control the number of bins
)
Out[35]:
In [36]:
# Flatten the matrix with vec() and plot the histogram
histogram(vec(i_opt),
title = "Optimized lamp strength",
xlabel = "Relative illumination",
ylabel = "Frequency",
legend = false,
bins = 5 # You can optionally control the number of bins
)
Out[36]:
In [ ]: