Structure variable changes to non-structure within function

1 view (last 30 days)
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;
I would appreciate any help you can offer. Please let me know if you would like any more information.
Roman

Accepted Answer

Walter Roberson
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.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!