Improve the .bindos format and add documentation

Change the binary format to the one suggested by Randle Taylor, and add
information about the sparse matrix binary dose format in the DOSXYZnrc
documentation.

HEN_HOUSE/doc/src/pirs794-dosxyznrc/inputs/dosxyznrc.inp

@@ -842,7 +842,7 @@ Record SC1-20a and SC1-21a (required for sources 20 and 21)
  Record 13 (9 if NMED=0)
  NCASE, IWATCH, TIMMAX, INSEED1, INSEED2,
  BEAM_SIZE, ISMOOTH,IRESTART,IDAT,IREJECT,ESAVE_GLOBAL,NRCYCL,IPARALLEL,
- PARNUM,n_split,ihowfarless,i_phsp_out
+ PARNUM,n_split,ihowfarless,i_phsp_out,i_bindos
 
  NCASE Number of histories
  IWATCH 0=>no tracking output, 1=>list every interaction
@@ -1014,6 +1014,11 @@ Record SC1-20a and SC1-21a (required for sources 20 and 21)
  (synchronized sources) then you also have the option of
  storing frMU_indx in the file (see source inputs).
 
+ i_bindos For outputting sparse binary dose data instead of the
+ dense ASCII .3ddose format. If i_bindos=0, then the
+ dose will be output in .3ddose format. Any other value
+ will output a .bindos file instead.
+
 ----------------------------------------------------------------------------
 ----------------------------------------------------------------------------
  EGSnrc INPUTS

HEN_HOUSE/doc/src/pirs794-dosxyznrc/pirs794-dosxyznrc.tex

@@ -26,6 +26,7 @@
 % Dave Rogers, 2001
 %
 % Contributors: Frederic Tessier
+% Marc-Andre Renaud
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -3499,6 +3500,12 @@ \subsection{{\tt i\_phsp\_out}}
 
 For more information on the IAEA phase space file format, see the BEAMnrc Users Manual\cite{Ro09}.
 
+\subsection{{\tt i\_bindos}}
+\indexm{i\_bindos}
+If {\tt i\_bindos} is set to 1, DOSXYZnrc outputs a sparse binary dose file instead of a dense ASCII dose file.
+The resulting file typically has a much smaller memory footprint than the default {\tt .3ddose} format.
+See section~\ref{subsec:bindos_format} for a description of the format.
+
 \subsection{ECUTIN}
 \indexm{ECUTIN}
 \label{ECUTIN}
@@ -4316,11 +4323,13 @@ \section{Format of Dose Outputs}
 
 \indexm{.egslst}
 The dose distributions calculated by DOSXYZnrc can be found in the output
-files, ``{\tt .egslst}'', ``{\tt .3ddose}'' and ``{\tt .pardose}''.
+files, ``{\tt .egslst}'', ``{\tt .3ddose}'', ``{\tt .bindos}'' and ``{\tt .pardose}''.
 The file ``{\tt .egslst}''
 \indexm{.pardose}
 \indexm{.3ddose}
+\indexm{.bindos}
 \indexm{files!.3ddose}
+\indexm{files!.bindos}
 \indexm{files!.pardose}
 contains not only the dose (when asked for) and statistical data but also the information
 about simulation geometry, number of histories run, CPU time used, etc.
@@ -4423,6 +4432,35 @@ \subsection{A Sample {\tt .3ddose} File}
 
 \lfoot[]{}
 
+\subsection{{\tt .bindos} Files}
+\label{subsec:bindos_format}
+\indexm{.bindos}
+\indexm{files!.bindos}
+A ``{\tt .bindos}'' file is output instead of a ``{\tt .3ddose}'' file if the
+flag ``{\tt i\_bindos}'' is set to 1. The same dose content is contained in each file,
+but ``{\tt .bindos}'' is a sparse format, hence only information about non-zero dose
+voxels is included. The linearized start and end indices for blocks of contiguous non-zero
+doses are stored. Both integers and floats are stored as 4 bytes each.
+
+\begin{enumerate}
+\item Number of voxels in the x, y, z direction (e.g., $n_x$, $n_y$, $n_z$, 3 ints)
+\item Voxel boundaries (cm) in x direction($ n_x$ +1 floats)
+\item Voxel boundaries (cm) in y direction($ n_y$ +1 floats)
+\item Voxel boundaries (cm) in z direction($ n_z$ +1 floats)
+\item Number of non-zero dose voxels ($n_{nz}$, 1 int)
+\item Number of voxel index pairs ($n_{b}$, 1 int)
+\item Linearized voxel index pairs indicating the [start, end) of contiguous non-zero
+ voxel doses. ($2 \times n_{b}$ ints)
+\item Nonzero voxel doses ($n_{nz}$ floats)
+\item Nonzero voxel uncertainties ($n_{nz}$ floats)
+\end{enumerate}
+
+The linearization formula is as follows:
+\begin{equation}
+v_{lin} = (n_x \times n_y) \cdot v_z + (v_x) \cdot v_y + v_x
+\end{equation}
+where $v_x$, $v_y$, $v_z$ are the voxel indices in the x, y, z axes, respectively.
+
 \subsection{{\tt .pardose} Files}
 \indexm{files!.pardose}
 \indexm{.pardose}

HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran

@@ -290,7 +290,7 @@
 " Record 13 (9 if NMED=0)
 " NCASE, IWATCH, TIMMAX, INSEED1, INSEED2,
 " BEAM_SIZE, ISMOOTH,IRESTART,IDAT,IREJECT,ESAVE_GLOBAL,NRCYCL,IPARALLEL,
-" PARNUM,n_split,ihowfarless,i_phsp_out
+" PARNUM,n_split,ihowfarless,i_phsp_out,i_bindos
 "
 " NCASE Number of histories
 " IWATCH 0=>no tracking output, 1=>list every interaction
@@ -462,6 +462,11 @@
 " (synchronized sources) then you also have the option of
 " storing frMU_indx in the file (see source inputs).
 "
+" i_bindos For outputting sparse binary dose data instead of the
+" dense ASCII .3ddose format. If i_bindos=0, then the
+" dose will be output in .3ddose format. Any other value
+" will output a .bindos file instead.
+"
 "----------------------------------------------------------------------------
 "----------------------------------------------------------------------------
 " EGSnrc INPUTS
@@ -1173,7 +1178,7 @@ REPLACE {;COMIN/SCORE/;} WITH {
  "in phase space data (phsp or BEAM simulation source"
  "only)"
  i_unit_out, "unit no. for IAEA phsp output"
- i_bindos, "set to 1 to output sparse binary file"
+ i_bindos,  "set to 1 to output sparse binary file"
  IWATCH,mxnp;
  $SHORT_INT endep_last;
 }
@@ -4313,9 +4318,9 @@ subroutine write_dose(iimax,jjmax,kkmax,xbnd,ybnd,zbnd,dd,ddun,unitnum,
  writemax20,bindos);
 
 
-"Writes dose to file designated by unitnum. Currently the subroutine only
-"writes in the .3ddose format (unitnum=3) and replaces a portion of the
-"main code. However, this subroutine will eventually output the dose in
+"Writes dose to file designated by unitnum. Currently the subroutine writes
+"in the .3ddose format (unitnum=3) or the .bindos format and replaces a portion
+"of the main code. However, this subroutine will eventually output the dose in
 "other formats as well.
 
 "3 possible methods of specifying the array sizes for dd and ddun, the
@@ -4372,15 +4377,17 @@ INTEGER iimax, "max number of x cells"
  jjmax, "max number of y cells"
  kkmax, "max number of z cells"
  ijkmax, "iimax*jjmax*kkmax"
- irtmp, "tmp variable to hold ir"
+ irtmp,  "tmp variable to hold ir"
  writemax20, "1 for summary of max 20 doses"
- bindos, "1 for binary file"
+ bindos,  "1 for binary file"
  unitnum, "unit number in which to write data"
  ii,jj,kk,mm,nn, "indices"
  NUMFRAC, "number of voxels for av. rel. uncertainty"
  MAXI(20),MAXJ(20),MAXK(20),"indices for max 20 doses"
  num_nonzero, "number of nonzero voxels"
- voxel_indices($MAXDOSE), "indices of nonzero voxels"
+ num_blocks, "number of voxel blocks"
+ start_block, "keep track of contiguous nonzero voxels"
+ voxel_blocks($MAXDOSE), "indices of nonzero voxels"
  egs_open_file;
 REAL xbnd($IMAX+1), "voxel x boundaries"
  ybnd($JMAX+1), "voxel y boundaries"
@@ -4392,7 +4399,7 @@ REAL xbnd($IMAX+1), "voxel x boundaries"
  ERR50AVG, "average % error of doses > 50% of max dose"
  ERR50FRAC, "average error of doses > 50% of max"
  " dose as a fraction of max dose"
- DOSEFRAC, "fraction for calculating av. rel. uncertainty"
+ DOSEFRAC,  "fraction for calculating av. rel. uncertainty"
  DOSEmin, "DOSEFRAC*MAXDOS(1)"
  ARUFRAC; "average relative uncertainty on all voxels"
  "with dose > DOSEFRAC*Dmax"
@@ -4405,6 +4412,8 @@ REPLACE {$BINOUT#;} WITH {write(unitnum){P1};};
 REPLACE {$IRDWD(#,#,#)} WITH {({P1}+({P2}-1)*iimax+({P3}-1)*iimax*jjmax)};
 
 num_nonzero = 0;
+num_blocks = 0;
+start_block = -1;
 ijkmax=iimax*jjmax*kkmax;
 
 IF(bindos=0)[
@@ -4420,10 +4429,26 @@ IF(bindos=0)[
 ELSE[
  DO ii=1,ijkmax[
  IF(dd(ii)>0.0)[
+ IF(start_block = -1) [
+ start_block = ii;
+ ]
  num_nonzero = num_nonzero + 1;
- voxel_indices(num_nonzero) = ii;
+ ] ELSE [
+ IF(start_block ~= -1) [
+ "First 0 dose voxel after a block of nonzero voxels"
+ num_blocks = num_blocks + 1;
+ voxel_blocks(2*num_blocks) = start_block;
+ voxel_blocks(2*num_blocks+1) = ii;
+ start_block = -1;
+ ]
  ]
  ]
+ "Capture the final voxel block"
+ IF(start_block ~= -1) [
+ num_blocks = num_blocks + 1;
+ voxel_blocks(2*num_blocks) = start_block;
+ voxel_blocks(2*num_blocks+1) = ijkmax;
+ ]
 
  "open the .bindos file for output"
  unitnum=egs_open_file(unitnum,-1,1,'.bindos');
@@ -4433,21 +4458,27 @@ ELSE[
  $BINOUT (zbnd(kk),kk=1,kkmax+1);
 
  $BINOUT num_nonzero;
+ $BINOUT num_blocks;
 
- "Writing voxel indices as a zero-indexed array even though fortran arrays"
- "are one-indexed."
- DO ii=1,num_nonzero[
- $BINOUT voxel_indices(ii) - 1;
+ "Writing start/end indices of contiguous non-zero voxel doses, 0-indexed"
+ "To be interpreted as [start, end)"
+ DO ii=1,num_blocks [
+ $BINOUT voxel_blocks(2*ii) - 1;
+ $BINOUT voxel_blocks(2*ii+1) - 1;
  ]
 
- DO ii=1,num_nonzero[
- dtemp = dd(voxel_indices(ii));
- $BINOUT dtemp;
+ DO ii=1,num_blocks [
+ DO jj=voxel_blocks(2*ii),voxel_blocks(2*ii+1)-1 [
+ dtemp = dd(jj);
+ $BINOUT dtemp;
+ ]
  ]
 
- DO ii=1,num_nonzero[
- dtemp = ddun(voxel_indices(ii));
- $BINOUT dtemp;
+ DO ii=1,num_blocks [
+ DO jj=voxel_blocks(2*ii),voxel_blocks(2*ii+1)-1 [
+ dtemp = ddun(jj);
+ $BINOUT dtemp;
+ ]
  ]
 ]