function status = DanMC(varargin)
% function status = DanMC(progressDMCfile [machine name])
% Daniel's Monte Carlo modelling function
% can call with a dmc file to resume a simulation (must end in .dmc
% or can call with a name of a text file containing the model description to begin a new analysis
% optionally can specify the machine name with the second parameter
%
% Daniel Fraser Septermber-Octomber 2000
% This is a wrapper function for the DanMCcore function
% This wrapper defines the model parameters, definitions
% saves the progress from the DanMCcore function
% And when complete saves the results file and plots graphs is requested.
%
% Has been modified to accept a progress file as input to enable to continuation
% of an interupted model

%---------------------------------------------------------------------------------
%Local definitions
PI              = pi;
D2P             = (PI/180);       %convert degrees to pi
P2D             = 1/D2P;

MODEL_COMPLETE           = 0;
MODEL_COMPLETE_NO_ACTION = 1;
MODEL_INCOMPLETE_ERROR   = 98;
FILE_DOES_NOT_EXIST      = 99;

%determine machine name
computername = computer;
if size(computername,2)>5
   computername = computername((size(computername,2)-4):size(computername,2));
end

if computername  == 'PCWIN'
   if length(varargin) == 2 %second parameter is the machine being launched on
      MACHINENAME = varargin{2}
   else
      MACHINENAME = machine;
      MACHINENAME = MACHINENAME(3:length(MACHINENAME));  %strips off preceeding //
      if isempty(MACHINENAME)
         MACHINENAME        = getenv('MACHINENAME');
      end
   end
      
elseif computername  == 'ALPHA'
   MACHINENAME = getenv('HOSTNAME');
   MACHINENAME = MACHINENAME(1:min(find(MACHINENAME=='.')-1));
elseif computername  == 'LNX86'
   MACHINENAME = getenv('HOSTNAME');
else
   error(['Don''t recongise your OS (',computername ,'), didn''t save progress']); 
end
   
if isempty(MACHINENAME)
   warning('MACHINENAME enviroment variable hasn''t been defined');
   MACHINENAME='';
end

%---------------------------------------------------------------------------------
%Definitions to share
def.I_WANT_TRAJECTORIES   = 1;   %1=true, 0=false
def.I_WANT_DATA_SHARING   = 1;   %regular dumping of data for pooling use
def.I_WANT_OUTPUT         = 1;
def.I_WANT_DUMP           = 1;   % To file!! :)
def.I_WANT_PLOTS          = 1;
def.I_WANT_SOUND          = 0;
def.I_WANT_DIAGNOSTICS      = 1;
def.I_WANT_TO_LOAD_PARAMS = 1;   %if want to check for the parameters file to load them first
def.SHARING_PER_MIN       = 5;   %how often the progress file is saved

def.DEAD       = 0;    			   % if photon is to be terminated 
def.ALIVE      = 1;   			   % if photon not yet terminated 
def.OUT        = 2;              % photon has left the sample :)
def.REFLECTED  = 3;              % photon has been reflected off wall

def.THRESHOLD  = 0.01;			   % used in roulette 
def.CHANCE     = 0.1;  			   % used in roulette 

def.ONE_MINUS_COSZERO = 1.0E-12;
def.JUST_INSIDE= 0.000001;       %depth inside boundary where reflection occurs

def.MAX_TRAJECTORYS    = 100;   %no more than this many trajectories can be stored.
def.MAX_INTERCEPTS     = 1000;  %no more than this many boundary reflections can be stored.
%---------------------------------------------------------------------------------

if length(varargin)==1
   progressDMCfile = varargin{1};
   if strcmp([progressDMCfile(size(progressDMCfile,2)-3:size(progressDMCfile,2))],'.dmc')
      if fileexists(progressDMCfile)
         [m,data] = examineDMC(progressDMCfile);
         if m.NumPhotons == data.NumPhotonsSoFar
            disp('Nothing to do, this model has finished!')
            status = MODEL_COMPLETE_NO_ACTION;
            return;
         end
         if isempty(findstr(progressDMCfile,MACHINENAME))
            warning('Resuming of a different machine...')
         end
         disp(['resuming ' progressDMCfile ' ...']);
         disp([num2str(data.NumPhotonsSoFar), ' Photons sent of ',num2str(m.NumPhotons)]);
      else
         disp('Error: File does not exist?!');
         status = FILE_DOES_NOT_EXIST;
         return
      end
   else %not a DMC file might be a model file
      m = ReadDMCModelTextFile(progressDMCfile);
   end
   
elseif length(varargin) >2
   disp('Error: Can''t resume more than one file')
   status = MODEL_INCOMPLETE_ERROR;
   return
   
else %must have called with no arguments
   %---------------------------------------------------------------------------------
   %Model Parameters
   if fileexists('DefaultDMCModel.txt') & def.I_WANT_TO_LOAD_PARAMS
      %  danmcdefaultparams'
      %handle = fopen('.danmcdefaultparams','rt');
      %m = readmodelparams(handle);
      m = ReadDMCModelTextFile('DefaultDMCModel.txt');
   else
      % sources include: 
      %1 = isotropic point, starts at origion and radiates evenly in all directions 
      %2 = laser z direction, starts at origion and radated in +ve z direction
      %3 = like 2 but displaced by radial_size in -ve z direction
      %4 = like 1 but limited to radiate only in NA 0.2
      %5 = like 1 but limited to radiate only in NA 0.8
      %6 = .5cm wide laser source
      %7 = .25cm wide laser source (Mand experiment)
      
      m.NumPhotons    = 200;      % set number of photons in simulation
      m.source        = 3;  

      %materal
      m.mua           = .01;        % cm^-1 (0.4)
      m.mus           = 0;      % cm^-1 (10)
      m.g             = 0;      % scattering anisotropy factor 0=isotropic 1=forward scattering
       

      %boundary
      m.BoundaryShape   = 1       % 1=sphere 2=box 3=cylinder
      m.radial_size   = 2.5;      % cm, total range over which bins extend
      m.BoundaryReflect = .9;     % Probability of internal reflection from boundary
      m.BoundaryRefType = 4;      % 1=reflects in direction of normal 2=cos spread, 3=isotropic scatter 4=spectral
      m.Iwantcutboundary= 0;      % 0=if intact boundary 1=if want cut else 2=one way door 3=partial boundary 4=no boundary
      m.CutSizeAngle    = 70*D2P;     % angle defining region of removed boundary 
      
      %datacollection
      m.NumBins         = 101;	            % set number of bins.
      m.dr              = m.radial_size/(m.NumBins-1);  % cm 
      
      m.GridBinsz        = 201;               %number of bins in grid
      m.GridRange1       = -m.radial_size;    %z postion extremity of grid
      m.GridRange2       = m.radial_size;                 %z postion extremity of grid
      m.dgz              = (m.GridRange2-m.GridRange1)/(m.GridBinsz-1);
      m.GridBinsx        = 201;               %number of bins in grid
      m.GridWidth        = 2*m.radial_size;     %width of the grid
      m.dgx              = m.GridWidth/(m.GridBinsx-1);
   
      %detector
      m.Diameter        = .2;        %Collecting Fibre 400micron = 0.4mm = 0.04cm
      m.AcceptanceAngle = 16.8*D2P;      %Collecting Fibre
   end %Model Parameters definitions
end %if resuming from a progress file or not

if ~exist('data','var')
   %--------------------------------------------------------------------------------
   %Data Initalisations
   data.NumPhotonsSoFar = 0;            %number of the photon currently tracked
   data.flops =0;
   data.CSph = zeros(m.NumBins,1);      %storage space for spherical shells from 0,0,0 or zaxis cylindrical shells
   data.CCylGrid = zeros(m.NumBins,m.GridBinsz); %storage for cylindical shells grid
   data.CFibre = zeros(m.NumBins,1);    %storage space for fibre sampling
   data.CFibreiso = zeros(m.NumBins,1); %storage space for isotropic colection fibre sampling
   data.CFibreiso2 = zeros(m.NumBins,1); %storage space for isotropic colection fibre sampling
   data.CGrid = zeros(m.GridBinsx,m.GridBinsz);
   data.SurfaceO=[];             %storage of the boundary intercepts when OUT
   data.SurfaceR=[];             %storage of the boundary intercepts when REFLECTED
   data.Field={};                %storage of all the trajectories
end

%local model based definiations
if def.I_WANT_DATA_SHARING
   %unique identifer for this combination of parameters
   ModelCode             = m.g *100;     
   ModelCode = ModelCode + m.CutSizeAngle*P2D *100;
   ModelCode = ModelCode + m.Iwantcutboundary *10000;
   ModelCode = ModelCode + m.BoundaryRefType *100000;
   ModelCode = ModelCode + m.source * 1000000;
else
   ModelCode     = 999999;   
end
   
%Initalisations for both with and without a file to resume
%--------------------------------------------------------------------------------
tic;
flops(data.flops);
rand('state',sum(100*clock));        %set the intall random number position to diff value
OldOutputFileName = [];  %The previously saved progress file name
OutputFileName = ['MC ' MACHINENAME ' ' num2str(ModelCode)];
%---------------------------------------------------------------------------------


%**************************************************************************************
while data.NumPhotonsSoFar < m.NumPhotons
   %**************************************************************************************
   %**************************************************************************************
   data = DanMCcore(m,def,data);   %Monte Carlo function.
   %**************************************************************************************
   %**************************************************************************************
   if data.NumPhotonsSoFar < m.NumPhotons
      %save progress
      ProgressP= round((data.NumPhotonsSoFar/m.NumPhotons)*100);
      NewOutputFileName = [OutputFileName ' new.dmc'];
      fwriteid = fopen(NewOutputFileName,'w');
      if (fwriteid == -1) 
         error('Error trying to open file for binary writing');
      else
         %*****************************************************************************
         [count,expectedcount] = WriteDMCFile(fwriteid,m,data);
         %*****************************************************************************
         status = fclose(fwriteid);
         if count < expectedcount
            error(['Error writing file ',num2str(count),' writen of ',expectedcount,' expected']);
         end
      end
      
      %check if old file exists if does remove it and rename new data to replace
      if ~isempty(OldOutputFileName) & fileexists(OldOutputFileName)
         if computername =='PCWIN'
            eval(['!del "' OldOutputFileName '"']);
            OldOutputFileName = [OutputFileName ' progress ' num2str(ProgressP) '%.dmc'];
            eval(['!rename "' NewOutputFileName '" "' OldOutputFileName '"']);
         elseif or(computername =='ALPHA',computername =='LNX86')
            eval(['!rm "' OldOutputFileName '"']);
            OldOutputFileName = [OutputFileName ' progress ' num2str(ProgressP) '%.dmc'];
            eval(['!cp "' NewOutputFileName '" "' OldOutputFileName '"']);
            eval(['!rm "' NewOutputFileName '"']);
         else
            error(['Don''t recongise your OS (',computername ,'), didn''t save progress']);   
         end
      else
         %must be the first progress save rename the new file to the old
         OldOutputFileName = [OutputFileName ' progress ' num2str(ProgressP) '%.dmc'];
         if computername =='PCWIN'
            if fileexists(OldOutputFileName) %an old file here with the same name!
               eval(['!del "' OldOutputFileName '"']);
            end
            eval(['!rename "' NewOutputFileName '" "' OldOutputFileName '"']);
         elseif or(computername =='ALPHA',computername =='LNX86')
            eval(['!cp "' NewOutputFileName '" "' OldOutputFileName '"']);
            eval(['!rm "' NewOutputFileName '"']);
         else
            error(['Don''t recongise your OS (',computername ,'), didn''t save progress']);   
         end
      end
   end
end   
%**************************************************************************************

%--------------------------------------------------------------------------------------
% output data into a solution file .dmc

str = datestr(now);
str = str([1:size(str,2)-6,size(str,2)-4:size(str,2)-3]);
OutputFileNameStr = [OutputFileName ' ' str '.dmc'];

if def.I_WANT_OUTPUT
  
   %file io
   if def.I_WANT_DUMP
      disp(['writing to file: ',OutputFileNameStr]);

      fwriteid = fopen(OutputFileNameStr,'w');
      if (fwriteid == -1) 
         error('Error trying to open file for binary writing');
      else
         %*****************************************************************************
         [count,expectedcount] = WriteDMCFile(fwriteid,m,data);
         %*****************************************************************************
         status = fclose(fwriteid);

         if count < expectedcount
            error(['Error writing file ',num2str(count),' writen of ',num2str(expectedcount)]);
         else
            %it wrote OK go ahead and delete the progress file
            if def.I_WANT_DATA_SHARING & ~isempty(OldOutputFileName)
               if computername =='PCWIN'
                  eval(['!del "' OldOutputFileName '"']);
               elseif or(computername =='ALPHA',computername =='LNX86')
                  eval(['!rm "' OldOutputFileName '"']);
               else
                  error(['Don''t recongise your OS (',computername ,'), didn''t save progress']);
               end
            end
            %graphical output
            if def.I_WANT_PLOTS
               [m,d]=examineDMC(OutputFileNameStr,1);
            end
         end
      end
   end
   
   %text output
   disp([num2str(toc),' seconds of calculations.']);
   disp([num2str(flops),' floating point operations']);

   %sound output
   if def.I_WANT_SOUND
      playjobsdone;
   end
end


      
