# Structure variable changes to non-structure within function

1 view (last 30 days)
Roman Hochuli on 24 Nov 2011
I am writing a Monte Carlo model for photon transport. I attribute the photon's properties (position, direction, weight, etc.) using a structure called 'photon'. Using the debugging utility, I can see that 'photon' is a structure up to the point where I call a function that I have written. When I step into the function, 'photon' changes to a non-structure type. There's nothing in the code that suggests why this might be, and I call 'photon' in many other functions I have written without any trouble.
The script is shown below (note: the error appears on the first iteration of both the two outer for-loops and the inner while-loop):
% number of photons to be simulated per beam configuration
N = 2200;
% determine physical map size [mm]
physMapHeight = 50;
physMapWidth = 11;
% scaling factor for higher resolution
SF = 10;
% high resolution grid size
mapHeight = physMapHeight * SF;
mapWidth = physMapWidth * SF;
% minimum simulable photon weight
threshold = 0.001;
% generate simulation map at resolution factor SF
map = zeros(mapHeight,mapWidth,'single');
map(1:100,:) = 0;
mapSize = size(map);
% initialise Q map, direcMap & intMap
Q = zeros(mapSize);
direcMap = zeros(mapSize);
intMap = zeros(mapSize);
% generate mua map
muMap = map;
% initialise escaped photon counter and escaped weight
escape_top = 0;
escape_bottom = 0;
escape_corner = 0;
escape_sides = 0;
escWeight = 0;
%debug
% direcCheck = NaN(1000,2);
tic;
for inh = (1/(2*SF)):(1/(2*SF)):physMapWidth
% simn is ith photon being simulated
for simn = 1:1:N/numel((1/(2*SF)):(1/(2*SF)):physMapWidth)
% inject photon into map
ink = 0.1*SF;
% initialise photon in y-direction
mux = 0;
muy = 1;
% give photon attributes (initialise photon weight & position)
inWeight = 1;
photon.W = inWeight;
% show completion status
% status = ['photon number ',int2str(simn), ' out of ', int2str(N/physMapWidth), ' for line ' num2str(inh*SF)];
% disp(status);
% initialise interaction counter
count = 0;
while photon.W > 0
if count == 0
% uninitialise coordinates of first event & update coordinates
photon.pos = [ink, inh];
photon.direc = [muy, mux];
else
photon.pos = [newk, newh];
photon.direc = [newmuy, newmux];
end
% direcMap = direcWatch(photon, direcMap, SF);
% get medium and load values as properties
medData = getMedium(map(round(photon.pos(1)*SF),round(photon.pos(2)*SF)));
n = medData.n;
g = medData.g;
mus = medData.mus;
mua = medData.mua;
% interaction
[direcCosy, direcCosx, photon] = scatterPhoton(photon, map, SF, threshold);
% calculate step size
s = howFar(mua, mus);
% move photon to next event & update coordinates
[newk, newh, newmuy, newmux, W] = movePhoton(map, s, photon, direcCosy, direcCosx, SF, 'on', 'off', 'off', 'off');
The initial part of the function 'movePhoton' is shown here:
function [newk, newh, newmuy, newmux, W] = movePhoton(map, Q, s, photon, direcCosy, direcCosx, SF, QFang, stepAdjust, stepAdjustBC, BC)
% calculates new position of photon using direction cosines
% reflects/transmits photon if track intersects interface
% assumes interactions more frequent than border crossing >> see below
% DOES NOT CONSIDER DOUBLE BORDER CROSSING i.e. 0|1|0 looks like one
% interface if no interaction in region of 1
%%calculate photon position from direction cosines
% note: y in negative direction!!!!
% get map size
[mapHeight, mapWidth] = size(map);
physMapHeight = mapHeight ./ SF;
physMapWidth = mapWidth ./ SF;
% extract photon weight
W = photon.W;
% pass photon.pos as coordinates
k = photon.pos(1);
h = photon.pos(2);
% pass photon.direc as explicit direction cosines
muy = photon.direc(1);
mux = photon.direc(2);
% determine coordinates of scattered photon
newmuy = muy.*direcCosy - mux.*direcCosx;
RN = RNG(1);
if RN(1) > 0.5
newmux = sin(acos(newmuy));
else
newmux = -sin(acos(newmuy));
end
newk = k + newmuy.*s;
newh = h + newmux.*s;
Roman

Walter Roberson on 24 Nov 2011
The function declaration for movephoton has Q as the second argument between map and s, but your function call to movephoton does not have Q in its calling sequence.