function data = DanMCcore(m,d,olddata);
% function data = DanMC(m,d);
% Daniel's Monte Carlo modelling function
%
% Daniel Fraser Septermber-October 2000
% core code converted from code by Steven L. Jacques et al.
% Added features include:
% Boundary checking with reflection or trasmision
% different source types
% Different data storage regemes inc. collection by a fibre
% File semaphones to share the data problem (uncompleted)
%
% Basic function is to launch a number of photons into a sphere
% of homogeneous material with a defined boundary and measure the 
% regions where the photons are absorbed.
%
% This code has been writen to be instructive not for efficency.
% Optimation can be acheved by: 
% Writing in a compiled form eg using C
% Intoducing interum calculation steps that are reused
% Using an optimised random number generator instead of matlab's
%
% Most of the theory has been derived - have a look at my labbook...
%
%Parameters:
% data    Is the data structure that is returned contains, trajectories, 
%         boundary intercepts and photon distributions
% m       The model structure, contains all the prameters used in this model
% d       The definiations struction, holds the program option parameters
% olddata The data set to be added to for the output data set.
%
% Assumtions:
% timer must already be started
% 

%Initalise
%--------------------------------------------------------------------
data = olddata;
iPhoton = data.NumPhotonsSoFar;
ProgressUpdateFlag =0;
Trajectory=[];               %storage of trajectory of current Photon
tag = zeros(m.NumBins,1);    %variable taging a photon if fibre takes it.
tag2= zeros(m.NumBins,1);    %variable taging a photon if isofibre takes it.
                             %so it is not counted twice.
if isfield(m,'GridBinsx')                             
   tagGrid=zeros(m.GridBinsx,m.GridBinsz);
end
TempNow = toc;
ClockCounter = TempNow-mod(TempNow,60);           %the time when next to display progress
SharingCounter=ClockCounter+60*d.SHARING_PER_MIN; %the time when next to save progress
ProgressP =(iPhoton/m.NumPhotons)*100;   %percentage progress


%optimations
%--------------------------------------------------------------------
DiameterD2      = m.Diameter/2; %Collecting Fibre
albedo          = m.mus/(m.mus + m.mua);
mue             = m.mua+m.mus;
PI              = pi;
PIover2         = PI/2;
PI2             = 2*PI;


 
%Launch the Nphotons into the medium
%--------------------------------------------------------------------
while (iPhoton < m.NumPhotons) & (ProgressUpdateFlag ==0)
   
   %Display progress on screen
   if (toc>ClockCounter) & d.I_WANT_OUTPUT
      ProgressP = (iPhoton/m.NumPhotons)*100;
      disp([num2str(ClockCounter/60),' minutes elapsed, Photon ',num2str(iPhoton),...
            ' of ',num2str(m.NumPhotons),' = ',num2str(ProgressP),'%']);
      ClockCounter=ClockCounter+60;
   end
   
   
   iPhoton = iPhoton+1;
   Weight = 1;
   PhotonStatus = d.ALIVE;
   
   %source Definition:
   %--------------------------------------------------------------------
   if m.source==1
      Position = [0;0;0]; %set photon position to origin
      Direction = ones(3,1)-2*rand(3,1); %Randomly set photontragectory to yield isotropic
      Direction = Direction/norm(Direction); %normalise to unit length.
      theta  =  acos(Direction(3));
      phi = PI2*rand;
   elseif m.source==2
      Position = [0;0;0]; %set photon position to origin
      Direction = [0;0;1]; %point source "laser" in the z direction
      theta = 0;
      phi = PI2*rand;
   elseif m.source==3
      Position = [0;0;-m.radial_size]; %set photon position to edge of sphere
      Direction = [0;0;1]; %point source "laser" in the z direction
      Direction = Direction/norm(Direction);
      theta = 0;
      phi = PI2*rand;
   elseif m.source==4
      Position = [0;0;0]; %set photon position to origin
      Direction = [0.2;.2;1]-[0.4*rand;0.4*rand;0]; %Randomly set photontragectory to yield NA=0.2
      Direction = Direction/norm(Direction); %normalise to unit length.
      theta  =  acos(Direction(3));
      phi = PI2*rand;
   elseif m.source==5
      Position = [0;0;0]; %set photon position to origin
      Direction = [0.8;.8;1]-[1.6*rand;1.6*rand;0]; %Randomly set photontragectory to yield NA=0.8
      Direction = Direction/norm(Direction); %normalise to unit length.
      theta  =  acos(Direction(3));
      phi = PI2*rand;
   elseif m.source==6
	   width=0.5*rand;
	   angle=2*pi*rand;
   	x=width*(cos(angle));
      y=width*(sin(angle));
      Position = [x;y;-m.radial_size]; %set photon position to edge of sphere .5cm beam width      dz=m.radial_size-(sqrt((4*m.radial_size^2)-(4*(Position(1)^2+Position(2)^2)))/2)
      if m.BoundaryShape ==1 %a sphere model
       	dz=m.radial_size-(sqrt((4*m.radial_size^2)-(4*(Position(1)^2+Position(2)^2)))/2);
         Position(3)=Position(3)+dz;%start on edge of sphere
      end
      Direction = [0;0;1]; %point source "laser" in the z direction
      Direction = Direction/norm(Direction);
      theta = 0;
      phi = PI2*rand;
   elseif m.source==7
		width=0.25*rand;
	   angle=2*pi*rand;
   	x=width*(cos(angle));
	   y=width*(sin(angle));
      Position = [x;y;-m.radial_size]; %set photon position to edge of sphere .5cm beam width      dz=m.radial_size-(sqrt((4*m.radial_size^2)-(4*(Position(1)^2+Position(2)^2)))/2)
      if m.BoundaryShape ==1 %a sphere model
         dz=m.radial_size-(sqrt((4*m.radial_size^2)-(4*(Position(1)^2+Position(2)^2)))/2);
         Position(3)=Position(3)+dz;%start on edge of sphere
      end
      Direction = [0;0;1]; %point source "laser" in the z direction
      Direction = Direction/norm(Direction);
      theta = 0;
      phi = PI2*rand;
   end
   %store the start point
   if d.I_WANT_TRAJECTORIES
      Trajectory=[Trajectory;Position];
   end

   while (PhotonStatus == d.ALIVE)
      %propagate Absorbance and Scatter check boundary conditions
      %--------------------------------------------------------------------
      RandomNum = rand;
      while (RandomNum == 0)  %make sure random number is not 0
         RandomNum =rand;
      end
      StepSize = -log(RandomNum)/(mue);  %mus?
      %log is natural log (base e). Stepsize =distance till absorbed or scattered
      
      newP = Position+(StepSize*Direction); %predicted new Position
      if d.I_WANT_DIAGNOSTICS & InCylinder(newP,m.radial_size,m.GridRange1,m.GridRange2)==0
         %warning('warning stepped beyond cylinder!')
         if newP(1)^2+newP(2)^2>(m.radial_size)^2
            %warning('beyond cylinder edge');
         end
         
      end
           
      if (m.Iwantcutboundary<4)
         %Check if it hit the boundary
         if m.BoundaryShape==1
            Intercept = SphereIntercept(Position,newP,m.radial_size);
         elseif m.BoundaryShape==2
            Intercept = BoxIntercept(Position,newP,m.GridWidth,m.GridRange1,m.GridRange2);
         else %must be 3==cylinder
            Intercept = CylinderIntercept(Position,newP,m.radial_size,m.GridRange1,m.GridRange2);
            if d.I_WANT_DIAGNOSTICS & or((newP(3)<m.GridRange1),newP(1)^2+newP(2)^2>(m.radial_size)^2)
               warning('warning intercept computation wrong!')
            end

         end
         
         if Intercept==[0;0;0] %didn't hit boundary
            Position = newP;
         else %hit the boundary unless it is cut away
            %compute spherical coordinates of the intercept
            PTheta = acos(Intercept(3)/norm(Intercept)); %surface position Theta
            PPhi= atan2(Intercept(2),Intercept(1)); %surface position Phi
            if Intercept(2)<0
               PPhi= PI2+PPhi;
            end
            
            if or((m.Iwantcutboundary==1 & abs(PTheta)<(m.CutSizeAngle)),m.Iwantcutboundary==2)
               %then on the removed boundary so escapes!
               PhotonStatus = d.OUT;
               %Position = newP; %(newP will be outside sphere)
               Position = Intercept;
          
            elseif (m.Iwantcutboundary==3 & abs(PTheta)<(m.CutSizeAngle)) %have a partial ref. section and it was hit
               if (rand<(m.BoundaryReflect/3)) %internal reflect 33.3% of the given boundary internal reflection
                  PhotonStatus = d.REFLECTED;
                  if size(data.SurfaceR,2)<d.MAX_INTERCEPTS
                     data.SurfaceR=[data.SurfaceR,[Intercept]];%;PTheta;PPhi;Weight]];
                  end
                  %the generation of a new direction is done in the scatter section
                  %Position = Intercept-Normal*d.JUST_INSIDE; %the surface becomes the scatter site.
                  Position = Intercept;
               else %escapes boundary (transmitted)
                  PhotonStatus = d.OUT;
                  Position = Intercept;
               end
               
            else %photon didn't hit removed section 
               %Compute transmistion probability and new reflectance angle & position
               if (rand<m.BoundaryReflect) %internally reflected
                  PhotonStatus = d.REFLECTED;
                  if size(data.SurfaceR,2)<d.MAX_INTERCEPTS
                     data.SurfaceR=[data.SurfaceR,[Intercept]];%;PTheta;PPhi;Weight]];
                  end
                  %the generation of a new direction is done in the scatter section
                  %Position = Intercept-Normal*d.JUST_INSIDE; %the surface becomes the scatter site.
                  Position = Intercept;
               else %escapes boundary (transmitted)
                  PhotonStatus = d.OUT;
                  Position = Intercept;
               end
            end
         end
      else
         Position = newP;   
      end   
      %Drop photon weight into local bin.
      %--------------------------------------------------------------------
      if or(PhotonStatus == d.ALIVE,PhotonStatus == d.REFLECTED)
         WAbsorbed = Weight*(1-albedo); %albedo should be exp(mua*StepSize)?
         Weight = Weight - WAbsorbed;
         
         if m.BoundaryShape ~=3
            %Spherical
            r = norm(Position);
            ir = floor(r/m.dr)+1;
            if (ir >=m.NumBins); ir =m.NumBins;   end;
            data.CSph(ir) = data.CSph(ir) + WAbsorbed;
            
         else
            %cylindrical
            r = sqrt(Position(1)^2+Position(2)^2);
            ir = floor(r/m.dr)+1;

            if (ir >=m.NumBins); ir =m.NumBins;   end;
            data.CSph(ir) = data.CSph(ir) + WAbsorbed;
            
            %cylindrical grid
            %if ((Position(3) >= m.GridRange1) & (Position(3) <= m.GridRange2))
            iz = floor((Position(3)+-m.GridRange1)/(m.dgz))+1;
            if d.I_WANT_DIAGNOSTICS & or(iz<1, iz > m.GridBinsz)               
               warning('warning iz incoreect in reflection!')
            end
            data.CCylGrid(ir,iz) = data.CCylGrid(ir,iz)+WAbsorbed;
            %end
         end
                  
         %fibre
         if (Position(1) > -DiameterD2) & (Position(1) < DiameterD2)
            if (Position(2) > -DiameterD2) & (Position(2) < DiameterD2)
               %along fibre track
               if (Position(3) >= -m.radial_size) & (Position(3) <= m.radial_size)
                  if (theta >= 0) & (theta<(m.AcceptanceAngle))
                     %will be accepted by fibre
                     iz = floor((Position(3)+m.radial_size)/(2*m.dr))+1;
                     if iPhoton>tag(iz)
                        %haven't captureed this photon previously
                        data.CFibre(iz) = data.CFibre(iz)+WAbsorbed;
                        %tag(iz) = iPhoton;
                     end
                  end
               end
            end
         end
         
         %fibreiso
         if (Position(1) > -DiameterD2) & (Position(1) < DiameterD2)
            if (Position(2) > -DiameterD2) & (Position(2) < DiameterD2)
               %along fibre track
               if (Position(3) >= -m.radial_size) & (Position(3) <= m.radial_size)
                  %all angles are accepted by fibre
                  iz = floor((Position(3)+m.radial_size)/(2*m.dr))+1;
                  %if iPhoton>tag2(iz)
                     %haven't captureed this photon previously
                     data.CFibreiso(iz) = data.CFibreiso(iz)+WAbsorbed;
                     %tag2(iz) = iPhoton;
                  %end
               end
            end
         end
         
         %fibreiso2
         if (Position(1) > -DiameterD2) & (Position(1) < DiameterD2)
            if (Position(2) > -DiameterD2) & (Position(2) < DiameterD2)
               if (Position(3) >= -m.radial_size) & (Position(3) <= m.radial_size)%along fibre track
                  %all angles are accepted by fibre
                  iz = floor((Position(3)+m.radial_size)/(2*m.dr))+1;
                  if iPhoton>tag2(iz)
                     %haven't captureed this photon previously
                     data.CFibreiso2(iz) = data.CFibreiso2(iz)+WAbsorbed;
                     tag2(iz) = iPhoton;
                  end
               end
            end
         end

         %grid
         if isfield(data,'CGrid') %this model has a CGrid aray defined...
            if (Position(1) > -m.GridWidth/2) & (Position(1) < m.GridWidth/2)
               if (Position(2) > -m.dgx/2) & (Position(2) < m.dgx/2)
                  if (Position(3) >= m.GridRange1) & (Position(3) <= m.GridRange2)
                     iz = floor((Position(3)+-m.GridRange1)/(m.dgz))+1;
                     ix = round(((m.GridBinsx)/2)+(Position(1)/(m.dgx)));
                     if iPhoton>tagGrid(ix,iz)
                        %haven't captureed this photon previously
                        data.CGrid(ix,iz) = data.CGrid(ix,iz)+WAbsorbed;
                        %tagGrid(ix,iz) = iPhoton;
                     end
                  end
               end
            end
         end
                  
      else %Photon must have hit boundary or escaped (ie PhotonStatus=OUT)
         %dump the energy in the surface array
         if size(data.SurfaceO,2)<d.MAX_TRAJECTORYS
            data.SurfaceO=[data.SurfaceO,[Intercept;PTheta;PPhi;Weight]];
         end
      end
  
      %Scatter into new direction
      %--------------------------------------------------------------------
      if PhotonStatus == d.REFLECTED
         
         %find normal
         if m.BoundaryShape == 1 %sphere
            Normal = -Intercept;
         elseif m.BoundaryShape == 3 %cylinder
            if Intercept(3)==m.GridRange1
               Normal = [0;0;1];
            elseif Intercept(3)==m.GridRange2
                  Normal = [0;0;-1];
            else
                  Normal = [-Intercept(1);-Intercept(2);0];
            end
         else %2 = box
            Normal = [0;0;0];
            disp('warning normal not defined');
         end
         Position = Intercept+Normal*d.JUST_INSIDE;
         
         if m.BoundaryRefType~=4 %not spectral reflection
            %define the current direction as Normal.
            Direction = Normal/norm(Normal); %normal to surface.
            theta  =  PI-PTheta;
            if d.I_WANT_DIAGNOSTICS & imag(theta)~=0
               warning('warning theta imaginary in reflection!')
            end
            phi = PPhi+PI;
            if phi>=PI2
               phi = phi-PI2;
            end
         end
                           
         switch m.BoundaryRefType
         case 1, %reflected in direction of normal
            newCosTheta =1;
         case 2, %reflected with a cosine distribution in dir of normal
            newCosTheta = cos(acos((2*rand)-1)-(PIover2)); %deviation from normal 0:90 deg.
         case 3, %reflected back isotropicly 0 to 90 degrees
            newCosTheta = cos((PIover2)*rand); % 0:(pi/2)
         otherwise, %spectral reflection
            %k = magnitude of the -incommingdirection + outgoing direction vector
            k = sqrt(2-2*cos(PI-2*(acos(dot(-Direction,Normal/norm(Normal))))));
            if d.I_WANT_DIAGNOSTICS & imag(k)~=0
               warning('warning k imaginary in reflection!')
            end
            Direction = k*(Normal/norm(Normal))+Direction; %spectral reflection direction
            theta = acos(Direction(3)/norm(Direction)); 
            phi = atan2(Direction(2),Direction(1)); 
            if Direction(2)<0
               phi= PI2+phi;
            end
            newCosTheta =1;
         end
                    
         newSinTheta = sqrt(1-newCosTheta*newCosTheta); %sqrt faster than sin()
         PhotonStatus = d.ALIVE;
         
      else
         %sample for costheta using Henyey Greenstein
         if m.g==0
            newCosTheta=(2*rand)-1; % -1:1 so theta 0:pi
         else
            temp = (1-m.g*m.g)/(1-m.g+2*m.g*rand);
            newCosTheta = (1+m.g*m.g-temp*temp)/(2*m.g);
         end
         newSinTheta = sqrt(1-newCosTheta*newCosTheta); %sqrt faster than sin()
      end
      
      %sample for phi
      newPhi = PI2*rand;
      newCosPhi = cos(newPhi);
      if (newPhi < PI)
         newSinPhi = sqrt(1-newCosPhi*newCosPhi);
      else
         newSinPhi = -sqrt(1-newCosPhi*newCosPhi);
      end
      

      %new direction computation
      %--------------------------------------------------------------------
      if m.mus == 0  %no scatter occurs
         newDirection = Direction;
      elseif (1-Direction(3) <= d.ONE_MINUS_COSZERO)
         %effectly travinging parallel to the +ve z-axis
         newDirection = [newSinTheta*newCosPhi;newSinTheta*newSinPhi;newCosTheta*sign(Direction(3))];
      else
         %compute the new direction by transforming these new scatter angles with
         %respect to the current direction theta, phi.
         codeT = [ cos(theta)*cos(phi), -sin(phi), sin(theta)*cos(phi);
                   sin(phi)*cos(theta),  cos(phi) , sin(theta)*sin(phi);
                   -sin(theta)         ,  0         , cos(theta) ];
         newDirection = codeT*[newSinTheta*newCosPhi;newSinTheta*newSinPhi;newCosTheta];
      end;
      
      if d.I_WANT_DIAGNOSTICS & imag(newDirection)~=0
         warning('warning newDirection imaginary in scatter section!')
      end
      Direction = newDirection;
      %note doing an inefficient conversion to phi and theta?, can write witout doing this
      theta  =  acos(Direction(3));  %angle of vetor to z axis
      phi = atan2(Direction(2),Direction(1));
      if Direction(2)<0
         phi = PI2+phi;
      end

      %Check roulette
      %--------------------------------------------------------------------
      %If photon weight is below threshold, then terminated using roulette technique.
      if (Weight < d.THRESHOLD)
         if rand <= d.CHANCE
            Weight = Weight / d.CHANCE; %i.e. increase by 1/Chance
         else
            PhotonStatus = d.DEAD;
         end
      end
      
      if d.I_WANT_TRAJECTORIES
         Trajectory=[Trajectory,Position];
      end
      
   end %while alive
   
   if d.I_WANT_TRAJECTORIES
      if iPhoton == 1; 
         data.Field ={Trajectory}; 
      elseif size(data.Field,2) < d.MAX_TRAJECTORYS
         data.Field{iPhoton} =Trajectory; %stores the trajectory of this photon
      end
      Trajectory =[]; %clears the path for the next one.
   end
      
   %If multiple computers are solving problem, regularly output progress
   %--------------------------------------------------------------------
   if (toc>SharingCounter) & d.I_WANT_DATA_SHARING
      %write the new data buy terminating this routine (don't worry it'll get run again!)
      ProgressUpdateFlag = 1;
   end %Progress will be saved when routine terminates.

end %looping through photons

%need to update how many photons done so far
data.NumPhotonsSoFar = iPhoton;
data.flops = flops;