Skip to content

Commit

Permalink
fixed bugs with luminaire rotation
Browse files Browse the repository at this point in the history
  • Loading branch information
@Frudawski
Frudawski committed Dec 7, 2023
1 parent d12bf80 commit 2e4fe2f
Show file tree
Hide file tree
Showing 9 changed files with 304 additions and 130 deletions.
179 changes: 109 additions & 70 deletions LUMOS functions/Meshing.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,11 @@
% Date: 08.11.2017 - last updated: 05.10.2021


function mesh = Meshing(input_surface,density,plot_mode,luminaires,measurements)
function mesh = Meshing(input_surface,density,plot_mode,luminaires,measurements,MeshOptimization)

if ~exist('MeshOptimization','var')
MeshOptimization = 1;
end

surface = input_surface;
% initialize
Expand Down Expand Up @@ -54,7 +57,7 @@
else
rotax = cross(normal,[0 0 1]);
R = rotMatrixD(rotax,rad2deg(elevation));

%{
figure(2)
hold off
Expand All @@ -74,13 +77,13 @@
dummy = 1;
%}
end

% rotate
a = R*[surfx surfy surfz]';
px = a(1,:);
py = a(2,:);
pz = a(3,:);

else
R = eye(3);
px = surfx;
Expand Down Expand Up @@ -125,7 +128,7 @@
px = a(1,:);
py = a(2,:);
pz = a(3,:);
else
else
P = eye(3);
end

Expand Down Expand Up @@ -200,7 +203,7 @@
nborderz = [];
nbordery = [];
for i = 1:numel(pz)-1
m = round(3*density*norm([r(i+1)-r(i);pz(i+1)-pz(i)]))-2;
m = round(3*density*norm([r(i+1)-r(i);pz(i+1)-pz(i)]))-2;
if m<2
m = 2;
end
Expand Down Expand Up @@ -246,70 +249,106 @@
[y2,z2] = meshgrid(xgrid2,zgrid2);


% luminaire subgrid
lumrayx = [];
lumrayz = [];
if MeshOptimization

lumrep = {};
lim = 1/3;

% luminare equidist ray2plane grid points

% check surface type
if ~strcmp(surface.type,'luminaire')

% delta angle stepwidth in degree
%da = 30;
% point on surface plane
%Q = [px(1) py(1) pz(1)];
% surface normal (y-z plane)
n = [1 0 0];
for L = 1:size(luminaires,2)
% luminare emittance vectors
[alpha,theta] = meshgrid(15:27.5:345,15:27.5:345);
theta(:,2:4:end) = theta(:,2:4:end)+11.25;
theta(:,4:4:end) = theta(:,4:4:end)+11.25;

alpha = deg2rad(alpha);
theta = deg2rad(theta);
[Lx,Ly,Lz] = sph2cart(alpha(:),theta(:),ones(size(alpha(:))));
Lxyz = (R*(P*[Lx Ly Lz]'))';
Lx = Lxyz(:,1);
Ly = Lxyz(:,2);
Lz = Lxyz(:,3);
% intersection calculation parameters
p = dot(luminaires{L}.coordinates,n);
vecr = dot([Lx Ly Lz],repmat(n,size(Lx,1),1),2);
para = -p./vecr;
S = (para.*[Lx,Ly,Lz]);
% intersection points
lumrayx = [lumrayx;S(:,2)+luminaires{L}.coordinates(2)];
lumrayz = [lumrayz;S(:,3)+luminaires{L}.coordinates(3)];

%{
for lum = 1:length(luminaires)
lumgrid = [luminaires{lum}.ldt.header.luminous_width luminaires{lum}.ldt.header.luminous_length 0]./1000;
N = round(lumgrid(1)./lim);
M = round(lumgrid(2)./lim);

if isequal(mod(N,2),0)
N = N+1;
end
if isequal(mod(M,2),0)
M = M+1;
end

[xgrid,ygrid] = DINgrid(lumgrid(1),lumgrid(2),0,'12464',[N M]);
lumrepcoord = [xgrid(:) ygrid(:) zeros(size(xgrid(:)))];
lumrepcoord(:,1) = lumrepcoord(:,1)-lumgrid(1)/2;
lumrepcoord(:,2) = lumrepcoord(:,2)-lumgrid(2)/2;
%lumreprotM = rotMatrix(luminaires{lum}.rotation);
a = luminaires{lum}.rotation;
M1 = rotMatrixD([1 0 0],a(1));
M2 = rotMatrixD([0 1 0],a(2));
M3 = rotMatrixD([0 0 1],a(3));
lumreprotM = M3*M2*M1;
lumrepcoord = (lumreprotM*lumrepcoord')';
for numb = 1:N*M
lumrep{end+1} = luminaires{lum};
%lumrep{end+1}.geometry{1} = [0 0 0 0 0 0;0 0 0 0 0 0];
lumrep{end}.coordinates = lumrep{end}.coordinates + lumrepcoord(numb,:);
end
end
luminaires = lumrep;

% luminare equidist ray2plane grid points

% check surface type
if ~strcmp(surface.type,'luminaire')

% delta angle stepwidth in degree
%da = 30;
% point on surface plane
%Q = [px(1) py(1) pz(1)];
% surface normal (y-z plane)
n = [1 0 0];
for L = 1:size(luminaires,2)
% luminare emittance vectors
[alpha,theta] = meshgrid(15:27.5:345,15:27.5:345);
theta(:,2:4:end) = theta(:,2:4:end)+11.25;
theta(:,4:4:end) = theta(:,4:4:end)+11.25;

alpha = deg2rad(alpha);
theta = deg2rad(theta);
[Lx,Ly,Lz] = sph2cart(alpha(:),theta(:),ones(size(alpha(:))));
Lxyz = (R*(P*[Lx Ly Lz]'))';
Lx = Lxyz(:,1);
Ly = Lxyz(:,2);
Lz = Lxyz(:,3);
% intersection calculation parameters
p = dot(luminaires{L}.coordinates,n);
vecr = dot([Lx Ly Lz],repmat(n,size(Lx,1),1),2);
para = -p./vecr;
S = (para.*[Lx,Ly,Lz]);
% intersection points
lumrayx = [lumrayx;S(:,2)+luminaires{L}.coordinates(2)];
lumrayz = [lumrayz;S(:,3)+luminaires{L}.coordinates(3)];

%{
% debuging
if dummy
hold on
plot3(luminaires{L}.coordinates(1),luminaires{L}.coordinates(2),luminaires{L}.coordinates(3),'r*')
plot3(Lx+luminaires{L}.coordinates(1),Ly+luminaires{L}.coordinates(2),Lz+luminaires{L}.coordinates(3),'r.')
end
%}
%}

end
end

ind = lumrayx>min(py) & lumrayx<max(py);
lumrayx = lumrayx(ind);
lumrayz = lumrayz(ind);
ind = lumrayz>min(pz) & lumrayz<max(pz);
lumrayx = lumrayx(ind);
lumrayz = lumrayz(ind);
ind = lumrayx>min(py) & lumrayx<max(py);
lumrayx = lumrayx(ind);
lumrayz = lumrayz(ind);
ind = lumrayz>min(pz) & lumrayz<max(pz);
lumrayx = lumrayx(ind);
lumrayz = lumrayz(ind);

end
end


% observer ray2plane grid points
orayx = [];
orayz = [];

% TODO: change to user defined resolution
% Tregenza table
% almucantar number, number of patches, almucantar center angle, azimuth increament, solid angle?
% almucantar number, number of patches, almucantar center angle, azimuth increament, solid angle?
TR = [1 30 6 12 12; 2 30 18 12 12; 3 24 30 15 12; 4 24 42 15 12; 5 18 54 20 12; 6 12 66 30 12; 7 6 78 60 12; 8 1 90 0 12];
% patch numbers and angles
% line 1: almucantars, line 2: azimuths, line 3: Patchnumber
Expand Down Expand Up @@ -402,8 +441,8 @@
x = x(~inblank);
z = z(~inblank);
% density correction factor
correction = 3;
correction = 3;

blankrgrid = [];
blankzgrid = [];
for v = 1:size(rotgap{b},1)-1
Expand Down Expand Up @@ -440,7 +479,7 @@
%surfedges = sort([chw [chw(2:end);chw(1)]],2);

for i = 1:numel(pz)-1
m = round(3*density*norm([r(i+1)-r(i);pz(i+1)-pz(i)]))-2;
m = round(3*density*norm([r(i+1)-r(i);pz(i+1)-pz(i)]))-2;
if m<2
m = 2;
end
Expand All @@ -449,7 +488,7 @@
chw = find(ismember(ltfround(DT.Points,12),ltfround([y' z'],12),'rows'));
surfedges = [chw(1:end-1) chw(2:end)];
if ~isempty(surfedges)
DT = force_edge(DT,surfedges);
DT = force_edge(DT,surfedges);
end
end

Expand Down Expand Up @@ -587,18 +626,18 @@ function plot_mesh_2D(dt,windows)

% window
for win = 1:size(windows,2)
%triplot(windows{win},'Color','b')
dt = [];
dt.Points = windows{win}.Points;
dt.ConnectivityList = windows{win}.ConnectivityList;
datax = [dt.Points(dt.ConnectivityList(:,1),1) dt.Points(dt.ConnectivityList(:,2),1) dt.Points(dt.ConnectivityList(:,3),1) dt.Points(dt.ConnectivityList(:,1),1)];
dataz = [dt.Points(dt.ConnectivityList(:,1),2) dt.Points(dt.ConnectivityList(:,2),2) dt.Points(dt.ConnectivityList(:,3),2) dt.Points(dt.ConnectivityList(:,1),2)];
%dataz = [dt.Points(dt.ConnectivityList(:,1),3) dt.Points(dt.ConnectivityList(:,2),3) dt.Points(dt.ConnectivityList(:,3),3) dt.Points(dt.ConnectivityList(:,1),3)];
for i = 1:size(datax,1)
fill(datax(i,:),dataz(i,:),[0 0.5 0.75],'Facecolor',[0 0.5267 0.6461])%,'Facealpha',0.5
end
%triplot(windows{win},'Color','b')

dt = [];
dt.Points = windows{win}.Points;
dt.ConnectivityList = windows{win}.ConnectivityList;

datax = [dt.Points(dt.ConnectivityList(:,1),1) dt.Points(dt.ConnectivityList(:,2),1) dt.Points(dt.ConnectivityList(:,3),1) dt.Points(dt.ConnectivityList(:,1),1)];
dataz = [dt.Points(dt.ConnectivityList(:,1),2) dt.Points(dt.ConnectivityList(:,2),2) dt.Points(dt.ConnectivityList(:,3),2) dt.Points(dt.ConnectivityList(:,1),2)];
%dataz = [dt.Points(dt.ConnectivityList(:,1),3) dt.Points(dt.ConnectivityList(:,2),3) dt.Points(dt.ConnectivityList(:,3),3) dt.Points(dt.ConnectivityList(:,1),3)];
for i = 1:size(datax,1)
fill(datax(i,:),dataz(i,:),[0 0.5 0.75],'Facecolor',[0 0.5267 0.6461])%,'Facealpha',0.5
end
end
axis equal

Expand All @@ -607,7 +646,7 @@ function plot_mesh_2D(dt,windows)


function plot_mesh_3D(wall)
% 3D coordinates
% 3D coordinates

dt.Points = wall.mesh.points;
dt.ConnectivityList = wall.mesh.list;
Expand All @@ -627,14 +666,14 @@ function plot_mesh_3D(wall)
% window(s)
try
for win = 1:size(wall.windows,2)

dt.Points = wall.windows{win}.mesh.points;
dt.ConnectivityList = wall.windows{win}.mesh.list;

wdatax = [dt.Points(dt.ConnectivityList(:,1),1) dt.Points(dt.ConnectivityList(:,2),1) dt.Points(dt.ConnectivityList(:,3),1) dt.Points(dt.ConnectivityList(:,1),1)];
wdatay = [dt.Points(dt.ConnectivityList(:,1),2) dt.Points(dt.ConnectivityList(:,2),2) dt.Points(dt.ConnectivityList(:,3),2) dt.Points(dt.ConnectivityList(:,1),2)];
wdataz = [dt.Points(dt.ConnectivityList(:,1),3) dt.Points(dt.ConnectivityList(:,2),3) dt.Points(dt.ConnectivityList(:,3),3) dt.Points(dt.ConnectivityList(:,1),3)];

for i = 1:size(wdatax,1)
fill3(wdatax(i,:),wdatay(i,:),wdataz(i,:),[0 0.5 0.75],'Facecolor',[0 0.5267 0.6461])
end
Expand Down
13 changes: 11 additions & 2 deletions LUMOS functions/copy_object.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -968,16 +968,25 @@ function update_Lumos_object(handles)
case 'luminaire'
luminaire_table(H,[],handles.Lumos)
refresh_2D(handles.Lumos,[],H)

try
refresh_2D_objects(handles.Lumos,[],H,H.data.object,'object')
catch
try
refresh_2D_objects(handles.Lumos,[],H,H.data.object,'luminaire')
catch
end
end
refresh_2D_objects(handles.Lumos,[],H,H.data.object,'luminaire')
try
refresh_3DObjects(handles.Lumos, [], H, H.data.object, 'object')
catch
try
refresh_3DObjects(handles.Lumos, [], H, H.data.object, 'luminaire')
catch
end
end
refresh_3DObjects(handles.Lumos, [], H, H.data.object, 'luminaire')

end


Expand Down
32 changes: 18 additions & 14 deletions LUMOS functions/hyperspecfisheye.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,6 @@

% loop over surfaces
for s = 1:numel(surfaces)

% other surfaces index
o = 1:numel(surfaces);
o = o(~(o==s));
Expand Down Expand Up @@ -130,21 +129,22 @@
end
end
% luminaire I and resulting E at observer
[Iv,Ev] = ldc2IE(luminaires{lumnum},observer);
%[Iv,Ev] = ldc2IE(luminaires{lumnum},observer);
% find all surfaces of current luminaire
lumsurf = [];
%lumsurf = [];
lumI = [];
lumA = [];
obsE = [];
ind1 = [];
ind2 = [];
for sn = 1:numel(surfaces)
if strcmp(luminaires{lum}.name,surfaces{sn}.name(1:end-2))
%[num2str(s),': ',surfaces{s}.name,' ', surfaces{sn}.name]
if strcmp(luminaires{lumnum}.name,surfaces{sn}.name(1:end-2))
% patch visibility matrix, incidence angles, emission angles, distance R
[vis,ina,ema,R,vec1] = pointVisibilityMatrix(P,surfaces{s},[]);
vis = pointVisibilityMatrix(P,surfaces{s},[]);
if sum(vis)>0
% if visible get patch areas and patch radiant
% intensities and at observer resulting E
% intensities and resulting E at observer
if sn == s
ind1 = numel(lumI)+1;
ind2 = numel(lumI)+numel(surfaces{s}.mesh.patch_area);
Expand All @@ -153,19 +153,18 @@
[I,E] = lum2IE(luminaires{lumnum},observer,surfaces{s});
lumI = [lumI; I];
obsE = [obsE; E];

end
end
end
if sum(vis)>0
% factors
f1 = Ev/sum(obsE);
f2 = f1.*lumI./lumA;
f2 = f2(ind1:ind2);
% radiance
lumL = lumI./lumA;
lumL = lumL(ind1:ind2)./size(surfaces{s}.mesh.list,1);
%lumL = repmat(mean(lumL(ind1:ind2))./numel(lumA(ind1:ind2)),numel(lumA(ind1:ind2)),1);%
lamidx = ismember(luminaires{lumnum}.spectrum.data(1,:),lambda);
L = ciespec2Y(lambda,luminaires{lumnum}.spectrum.data(2,lamidx));
f3 = f2./L;
L = f3.*repmat(luminaires{lumnum}.spectrum.data(2,lamidx),numel(f3),1);
specL = ciespec2Y(lambda,luminaires{lumnum}.spectrum.data(2,lamidx));
factor = lumL./specL;
L = factor.*repmat(luminaires{lumnum}.spectrum.data(2,lamidx),numel(factor),1);
else
continue
end
Expand All @@ -183,7 +182,12 @@
singularity = (-360:90:720)';
razd = [razd-360; razd; razd+360; singularity];
reld = [repmat(reld,3,1);repmat(90,numel(singularity),1)];
try
L = [repmat(L,3,1);repmat(L(idx,:),numel(singularity),1)];
catch
dummy = 1;
end


% surface polygon
N = 100;
Expand Down
Loading

0 comments on commit 2e4fe2f

Please sign in to comment.