% function [fluence, nPhoton] = ReadHis( filenm, muas, freq )
%
% Read the History file generated by tMCimg and determine the
% number of photons received at each detector and calculate the
% fluence at each detector based on the given absorption
% coefficient of each tissue type and the modulation frequency
%
% filenm - The file name of the history file with out the
%          .his extention.  The .inp file of the same name must
%          also be present
%
% muas   - The mua of each tissue type.  Each tissue type is
%          a different column.  Rows can be added to perform
%          the fluence calculation simultaneously for different
%          sets of absorption coefficients. Units are per mm.
%        
% freq   - Modulation frequency of the source in MHz
%

function [fluence, nPhoton] = ReadHis( filenm, muas, freq )

  w = 2*pi*freq*1e6;
  clight = 2.99e11;
  
  file = sprintf( '%s.inp', filenm);
  fid = fopen(file, 'r');
  nPhotons = fscanf( fid, '%d', 1);
  seed = fscanf( fid, '%ld', 1);
  src = fscanf( fid, '%lf', 3);
  srcn = fscanf( fid, '%lf', 3);
  trange = fscanf( fid, '%lf', 3);
  imgfile = fscanf( fid, '%s', 1);
  xrange = fscanf( fid, '%lf', 4);
  yrange = fscanf( fid, '%lf', 4);
  zrange = fscanf( fid, '%lf', 4);
  nTissue = fscanf( fid, '%d', 1);
  for idx=1:nTissue
    layer(:,idx) = fscanf( fid, '%lf', 4);
    musps(idx) = layer(1,idx) * (1-layer(2,idx));
  end
  nDet = fscanf( fid, '%d', 2);
  fclose(fid);
  
  wv = w/clight/layer(1,4);
  
  file = sprintf( '%s.his', filenm);
  fid = fopen(file, 'rb');
  a = fread(fid,'float32');
  fclose(fid);
%  nPhotonsDet = size(a,1)/(nTissue+1);
  nPhotonsDet = floor(size(a,1)/(nTissue+1));

  a = reshape(a(1:nPhotonsDet*(nTissue+1)),[1+nTissue,nPhotonsDet])';

  for idxLoop = 1:size(muas,1)
    for idx=1:size(muas,2)
      b(:,idx+1) = a(:,idx+1) * (sqrt(-1)*wv  - muas(idxLoop,idx));
    end
  
    foo = b(:,2:size(a,2));
    foo = sum(foo,2);
    foo = exp(foo);

    for idx = 1:nDet(1)
      list = find((a(:,1)+1)==idx);
      nPhoton(idx) = length(list);

      if ~isempty(list)
	fluence(idxLoop, idx) = sum(foo(list));
      end
    end
  end

  fluence = fluence / (nPhotons * pi * nDet(2)^2);
