%modified ExamineDMC to correct the data for a rounding bug.
function [m,data] = examineDMC2(FileName,I_WANT_PLOTS);
% loads a DMC data file

freadid = fopen(FileName,'r');
if (freadid == -1) 
   error('Error trying to open file for binary reading');
else
   [m,data] = ReadDMCFile(freadid);
   status = fclose(freadid);
end

if nargin < 2
   I_WANT_PLOTS = 0;
end

%shift the data to approximate a correction for the rounding bug.
%f2;semilogy(data.CFibre,'-+');
data.CFibre(1)=data.CFibre(1)+data.CFibre(2)/2;
data.CFibreiso(1)=data.CFibreiso(1)+data.CFibreiso(2)/2;
data.CFibreiso2(1)=data.CFibreiso2(1)+data.CFibreiso2(2)/2;
for i = 2:(m.NumBins-2);
   data.CFibre(i)=(data.CFibre(i)/2)+data.CFibre(i+1)/2;
   data.CFibreiso(i)=(data.CFibreiso(i)/2)+data.CFibreiso(i+1)/2;
   data.CFibreiso2(i)=(data.CFibreiso2(i)/2)+data.CFibreiso2(i+1)/2;
end;
%can't use bin 101 because it may be abnormally large
slope =(data.CFibre(m.NumBins-2)-data.CFibre(m.NumBins-1))/2;
data.CFibre(m.NumBins-1)=(data.CFibre(m.NumBins-1))-(slope/2);

slope =(data.CFibreiso(m.NumBins-2)-data.CFibreiso(m.NumBins-1))/2;
data.CFibreiso(m.NumBins-1)=(data.CFibreiso(m.NumBins-1))-(slope/2);

slope =(data.CFibreiso2(m.NumBins-2)-data.CFibreiso2(m.NumBins-1))/2;
data.CFibreiso2(m.NumBins-1)=(data.CFibreiso2(m.NumBins-1))-(slope/2);
%corrections finish


if I_WANT_PLOTS
   FileName=FileName(4:(size(FileName,2)-4));
   figure
   if (size(data.Field,2)<1000)
      for i = 1:size(data.Field,2)
         hold on;
         plot3(data.Field{i}(1,:),data.Field{i}(2,:),data.Field{i}(3,:));
      end
   end
   if (size(data.SurfaceO,1)>0)
      plot3(data.SurfaceO(1,1:size(data.SurfaceO,2)),...
      data.SurfaceO(2,1:size(data.SurfaceO,2)),...
      data.SurfaceO(3,1:size(data.SurfaceO,2)),'ro');
   end
   if (size(data.SurfaceR,1)>0)
      hold on;
      plot3(data.SurfaceR(1,1:size(data.SurfaceR,2)),...
      data.SurfaceR(2,1:size(data.SurfaceR,2)),...
      data.SurfaceR(3,1:size(data.SurfaceR,2)),'go');
   end
   xlabel('X');
   ylabel('Y');
   zlabel('Z');
   title([FileName,' ',num2str(size(data.Field,2)),' Trajectories']);
   
   if ~isfield(m,'BoundaryShape')
      m.BoundaryShape=1;   %the default boundary for backward compatablity
   end
   
   if m.BoundaryShape==1 
      [X,Y,Z] = sphere(40);
      X = m.radial_size*X;
      Y = m.radial_size*Y;
      Z = m.radial_size*Z;
      plot3(X,Y,Z,'y.');
      
      if isfield(data,'CFibreiso')
         figure;
         semilogy((([1:m.NumBins]-1)*2*m.dr)-m.radial_size,data.CFibre);
         title([FileName,' ',num2str(m.NumPhotons),' Photons Fibre']);
      end
      if isfield(data,'CFibreiso')
         figure;
         semilogy((([1:m.NumBins]-1)*2*m.dr)-m.radial_size,data.CFibreiso);
         title([FileName,' ',num2str(m.NumPhotons),' Photons isoFibre']);
      end
      if isfield(data,'CFibreiso2')
         figure;
         semilogy((([1:m.NumBins]-1)*2*m.dr)-m.radial_size,data.CFibreiso2);
         title([FileName,' ',num2str(m.NumPhotons),' Photons isoFibre2']);
      end

   elseif m.BoundaryShape==2  %box
      Cube([-(m.GridWidth/2),-(m.GridWidth/2),GridRange1],[(m.GridWidth/2),(m.GridWidth/2),GridRange2]);
      
   else % m.BoundaryShape==3  cylinder
      [X,Y,Z] = cylinder([m.radial_size,m.radial_size],40);
      Z=(Z.*(m.GridRange2-m.GridRange1));
      Z = Z-m.GridRange1;
      plot3(X,Y,Z,'y.')
      
      figure;
      semilogy((([1:m.GridBinsz]-1)*m.dgz)-m.GridRange1,data.CGrid(round(m.GridBinsx/2),:));
      title([FileName,' ',num2str(m.NumPhotons),' Photons centre of grid']);
      
      figure;
      plot(([1:size(data.CSph,1)]-1)*m.dr,data.CSph);
      title([FileName,' ',num2str(m.NumPhotons),' Cylindrical shells about z axis']);
      
      Map=zeros(size(data.CCylGrid));
      for i=1:size(Map,1);
         for j=1:size(Map,2);
            r = (i*m.dr);
            prevr = ((i-1)*m.dr);
            A = pi*r^2;
            prevA = pi*prevr^2;
            Vol=abs((A-prevA)*m.dgz);
            Map(i,:)=data.CCylGrid(i,:)./Vol;
         end;
      end;
      Map=[flipud(Map);Map];
      figure;
      contour(log(Map));
      
   end
  
   if isfield(data,'CGrid')
      figure;
      map = reshape(log(data.CGrid),m.GridBinsx,m.GridBinsz);
      try
         contour(map);
      catch
         warning('Unable to produce contour map');
      end
      
      title([FileName,' ',num2str(m.NumPhotons),' Photons CGrid']);
   end
   

end
