Visualizing Photometry with Octave and Radiance

Summary

This GNU Octave code creates grid of points on a large sphere centered on the origin and runs an illuminance calculation on the surface of the sphere, with the virtual sensors pointed at the origin. If you place a simulated Radiance luminaire at the origin, this will calculate luminaire photometry.

For example, to view the photometry of a luminaire stored in light.rad, load the photometer.m script into Octave, then use the following command:

photometer('light.rad')

Results will be placed in the DX folder and may be viewed with OpenDX.

Code

photometer.m:

function flux_total = flux_calc(name)
%version 1.00 - this includes max intensity direction, upward and downward component

unix('mkdir tmp');

unix('mkdir dx');

%this is the radius of the sphere, beware of size of fixture
%radius should be significantly larger!
r = 100000;

%angle
x = 360;
y = 180;

%steps in degrees, higher accuracy for 1 degree
%lower accuracy for higher values (i.e. 5 - 10 degrees)
step = 1;

i = 0:step:x;
j = 0:step:y;

%how many points
nx = size(i,2);
ny = size(j,2);

i = repmat(i',ny,1);
j = reshape(repmat(j,nx,1),(nx)*(ny),1);
R = repmat(r,(nx)*(ny),1);

vx = R.*(cos(i.*pi/180)).*cos((j-90).*pi/180);
vy = R.*sin(i.*pi/180).*cos((j-90).*pi/180);
vz = R.*sin((j-90).*pi/180);
vvx = -1*cos(i.*pi/180).*cos((j-90).*pi/180);
vvy = -1*sin(i.*pi/180).*cos((j-90).*pi/180);
vvz = -1*sin((j-90).*pi/180);
grid = [vx vy vz vvx vvy vvz];
angles = grid(:,1:3);
angles = [angles(:,1)./R angles(:,2)./R angles(:,3)./R];

save -ascii tmp/grid grid;
save -ascii tmp/angles  angles;

unix(sprintf('oconv %s > tmp/light.oct',name));
unix(sprintf('echo x = %d > dx/light.data',ny));
unix(sprintf('echo y = %d >> dx/light.data',nx));
unix(sprintf('echo z = 1 >> dx/light.data'));

unix(sprintf('rtrace -w- -h- -ov -ab 0 -aa 0 -I tmp/light.oct < tmp/grid | rcalc -e ''$1=$1*179*%f*%f'' > tmp/ill',r,r));
unix('rlam  tmp/angles  tmp/ill  >> dx/light.data');

i = 0:step:x;
j = 0:step:y;

low_end = [j-step/2]';
high_end = [j+step/2]';
low_end(1) = 0;
high_end(ny) = y;

solid_angle = (step)*(pi/y)*(cos(low_end*pi/180)-cos(high_end*pi/180));
solid_angle = repmat(solid_angle,1,nx);
solid_angle = reshape(solid_angle',nx*ny,1);

ill = load('tmp/ill');
flux = solid_angle.*ill;

flux_up = sum(flux(angles(:,3)>0,:));
flux_down = sum(flux(angles(:,3)<0,:));
flux_total = sum(flux);


[high_spot direction_index ]= max(ill);
direction = angles(direction_index,:);

unix('rm -r tmp');

%%%%%%%%%%%%%%%%%%
disp(sprintf('Completed. The total luminaire flux is: %.0f lm',flux_total))
disp(sprintf(' '))
disp(sprintf('The upward component is: %.0f lm',flux_up))
disp(sprintf('The downward component is: %.0f lm',flux_down))
disp(sprintf(' '))
disp(sprintf('The highest intensity is : %.0f cd',high_spot))
disp(sprintf('and the direction is: %.3f %.3f %.3f ',direction(1),direction(2),direction(3)))
disp(sprintf(' '))
disp('Now go to DX and open the result file')