function [m,data] = ReadDMCFile(freadid);
%function [outcount,expectedcount = ReadDMCFile(freadid,m,data);

%definisions note must be same as in WriteDMCFile
BLOCKTYPEPARAM = 1001;
BLOCKTYPEFIBRE = 1002;
BLOCKTYPEFIBREISO = 1003;
BLOCKTYPEFIBREISO2 = 1004;
BLOCKTYPESPHSHELL = 1005;
BLOCKTYPEREFBOUNDARY = 1006;
BLOCKTYPEOUTBOUNDARY = 1007;
BLOCKTYPETRAJ = 1008;
BLOCKTYPEGRID = 1009;
BLOCKTYPECYLSHELLGRID = 1010;

TrajNum = 1;
fseek(freadid,0,'eof');
eof = ftell(freadid);
frewind(freadid);

data.NumPhotonsSoFar = fread(freadid,1,'double');
data.flops = fread(freadid,1,'double');

while ftell(freadid)<eof
   %disp(num2str(ftell(freadid)));
   blocktype =fread(freadid,1,'double');
   size =fread(freadid,1,'double');
   switch blocktype
   case BLOCKTYPEPARAM,
      m = readmodelparams(freadid,size);
   case BLOCKTYPEFIBRE,
      data.CFibre = fread(freadid,size,'double');
   case BLOCKTYPEFIBREISO,
      data.CFibreiso = fread(freadid,size,'double');
   case BLOCKTYPEFIBREISO2,
      data.CFibreiso2 = fread(freadid,size,'double');
   case BLOCKTYPESPHSHELL,
      data.CSph = fread(freadid,size,'double');
   case BLOCKTYPEREFBOUNDARY,
      data.SurfaceR = fread(freadid,[3,size/3],'double');
   case BLOCKTYPEOUTBOUNDARY,
      data.SurfaceO = fread(freadid,[6,size/6],'double');
   case BLOCKTYPETRAJ,
      %disp([num2str(size),' ',num2str(3),' ',num2str(size/3),' ']);
      Traj = fread(freadid,[3,size/3],'double');
      data.Field{TrajNum} =Traj;
      TrajNum=TrajNum+1;
   case BLOCKTYPEGRID,
      data.CGrid = fread(freadid,size,'double');
      data.CGrid=reshape(data.CGrid,m.GridBinsx,m.GridBinsz);
   case BLOCKTYPECYLSHELLGRID
      data.CCylGrid = fread(freadid,size,'double');
      data.CCylGrid=reshape(data.CCylGrid,m.NumBins,m.GridBinsz);
   otherwise
      disp('don''t recongnise this blocktype!');
   end
end


