Structure variable changes to non-structure within function

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

Answer 1

[Walter Roberson]

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.