function blk_out = blk_gen(blk,hdr)
%USEAGE: BLK = blk_gen(BLK,HDR);
%
%INPUT ARGUMENTS:
%
%BLK
%  Structure of block parameters. Fields include:
%      FIELD       DESCRIPTION
%      fn    Name of block parameter file
%
%      ref         3D coordinate of block reference point used to 
%                  position block in ring
%
%      min         3-element vector of lower bound on block volume (X,Y,Z)
%                  given with respect to "ref" point
%
%      max         3-element vector of upper bound on block volume (X,Y,Z)
%                  given with respect to "ref" point
%
%      x           Array of X coordinates for layer divisions (w.r.t. ref). 
%                  For Nx layers, "x" should be a (Nx-1)-element vector.
%                  Values should be monotonically increasing.
%
%      y           List of 1D-arrays (one array for each layer) of the form:
%                     {[...],[...],...}, where [...] is a 1D-array.
%                  The array for a given layer gives the y coordinates for 
%                  voxel boundaries/material changes. For (Ny+1) segments,
%                  there are Ny divisions.  Note that the divisions for
%                  each layer can be distinct.  
%
%      z           List of 1D-arrays (one array for each layer) of the form:
%                     {[...],[...],...}, where [...] is a 1D-array.
%                  The array for a given layer gives the z coordinates for 
%                  voxel boundaries/material changes. For Nz segments,
%                  there are (Nz-1) divisions.  Note that the divisions for
%                  each layer can be distinct.  
%
%      idx         List of 2D-arrays (one array for each layer) of the form:
%                     {[...],[...],...}, where [...] is a ((Ny+1)*(Nz+1)) 
%                  2D-array. The array for a given layer gives the material
%                  index for each of the ((Ny+1)*(Nz+1)) elements of that 
%                  layer voxel boundaries/material changes.
%
%      act         List of 2D-arrays (one array for each layer) of the form:
%                     {[...],[...],...}, where [...] is a ((Ny+1)*(Nz+1)) 
%                  2D-array.  The array for a given layer flags if each of 
%                  the ((Ny+1)*(Nz+1)) elements of that layer is active 
%                  (0=inactive, 1=active).
%
%
%HDR
%  Optional structure giving information for file header. The date 
%  and fn_blk are automatically included (no redundant input needed).
%  User input fields include:
%      FIELD       DESCRIPTION
%      author      string giving name of author
%
%      comments    list of strings giving comments (1 item per line) 
%
%OUTPUTS
%
%BLK
%  Updated structures of blk including array of azimuths and tilts
%
%William Hunter, 2008-06-04

blk.Nlayer = numel(blk.x)+1;

if isfield(blk,'y')
if blk.Nlayer~=numel(blk.y)
    error('number of layers for x and y must agree');
end
end

if isfield(blk,'z')
if blk.Nlayer~=numel(blk.z)
    error('number of layers for x and z must agree');
end
end

fid = fopen(blk.fn,'wt');

%Write Header
txt{1} = ['#  Block parameter file: ',blk.fn];
if nargin>1
    if isfield(hdr,'author')
        txt{end+1} = ['#  Author/Owner:         ',hdr.author];
    end
    if isfield(hdr,'comments')
        if iscell(hdr.comments)
            txt{end+1} = ['#  Comments:             ',hdr.comments{1}];
            for i = 2:numel(hdr.comments)
                txt{end+1} = ['#                        ',hdr.comments{i}];
            end
        elseif ischar(hdr.comments)
            txt{end+1} = ['#  Comments:             ',hdr.comments];
        else
            error('hdr.comments is not a string or list of strings');
        end
    end
end

max_len = 2;
for i = 1:numel(txt)
    max_len = max(max_len,length(txt{i}));
end

fprintf(fid,[repmat('#',1,max_len),'\n#\n']);
for i = 1:numel(txt), fprintf(fid,[txt{i},'\n']); end
fprintf(fid,['#\n',repmat('#',1,max_len),'\n']);
fprintf(fid,'\n');


%Write block reference point
fprintf(fid,'# Reference point on block used for ring position & orientation.\n');
fprintf(fid,'REAL	block_reference_x       =   %0.14f\n',blk.ref(1));
fprintf(fid,'REAL	block_reference_y       =   %0.14f\n',blk.ref(2));
fprintf(fid,'REAL	block_reference_z       =   %0.14f\n',blk.ref(3));
fprintf(fid,'\n');


%Write block bounding volume
fprintf(fid,'# BLOCK OUTER DIMENSIONS\n');
fprintf(fid,'# Provides the outer bounds of the block.  The maximum and minimum extent of the block\n');
fprintf(fid,'# are given with respect to the reference point.  The coordinates (X,Y,Z) represent the\n');
fprintf(fid,'# radial, transaxial, and axial directions, respectively.\n');
fprintf(fid,'REAL	block_x_minimum         =   %0.14f\n',blk.min(1));
fprintf(fid,'REAL	block_x_maximum         =   %0.14f\n',blk.max(1));
fprintf(fid,'REAL	block_y_minimum         =   %0.14f\n',blk.min(2));
fprintf(fid,'REAL	block_y_maximum         =   %0.14f\n',blk.max(2));
fprintf(fid,'REAL	block_z_minimum         =   %0.14f\n',blk.min(3));
fprintf(fid,'REAL	block_z_maximum         =   %0.14f\n',blk.max(3));
fprintf(fid,'\n');


%Write number of layers in block
fprintf(fid,'INT    block_num_layers        =   %d\n',blk.Nlayer);
fprintf(fid,'\n');


%Write blocks
for n = 1:blk.Nlayer
    fprintf(fid,'NUM_ELEMENTS_IN_LIST	block_layer_info_list = 5\n');
    fprintf(fid,'\n');
    
    if n==1
        fprintf(fid,'        REAL           block_layer_inner_x   = %f\n',blk.min(1));
    else
        fprintf(fid,'        REAL           block_layer_inner_x   = %f\n',blk.x(n-1));
    end
    fprintf(fid,'                       # (must be == to outer x of previous layer)\n');
    fprintf(fid,'\n');
    
    if n==blk.Nlayer
        fprintf(fid,'        REAL           block_layer_outer_x   = %f\n',blk.max(1));
    else
        fprintf(fid,'        REAL           block_layer_outer_x   = %f\n',blk.x(n));
    end
    fprintf(fid,'\n');
    
    
    if isfield(blk,'y')
        Ny = numel(blk.y{n});
    else
        Ny = 0;
    end
    fprintf(fid,'        NUM_ELEMENTS_IN_LIST	block_layer_num_y_changes = %d\n',Ny);
    for i=1:Ny
        fprintf(fid,'        REAL           block_layer_y_change  = %f\n',blk.y{n}(i));
    end
    fprintf(fid,'\n');
    
    
    if isfield(blk,'z')
        Nz = numel(blk.z{n});
    else
        Nz = 0;
    end
    fprintf(fid,'        NUM_ELEMENTS_IN_LIST	block_layer_num_z_changes = %d\n',Nz);
    for i=1:numel(blk.z{n})
        fprintf(fid,'        REAL           block_layer_z_change  = %f\n',blk.z{n}(i));
    end
    fprintf(fid,'\n');
    
    
    J=(Ny+1)*(Nz+1);
    fprintf(fid,'        NUM_ELEMENTS_IN_LIST	block_layer_material_elements  = %d\n',J);
    
    
    for j=1:J
        fprintf(fid,'            # Material element  %d\n',j-1);
        fprintf(fid,'            NUM_ELEMENTS_IN_LIST	block_material_info_list = 2\n');
        fprintf(fid,'                INT	block_material_index = %d\n',blk.idx{n}(j));
        if blk.act{n}(j)
            fprintf(fid,'                BOOL	block_material_is_active = TRUE\n');
        else
            fprintf(fid,'                BOOL	block_material_is_active = FALSE\n');
        end
    end
    fprintf(fid,'\n');

end

            
fclose(fid);

if nargout>0; blk_out = blk; end