function block_input = define_regular_block()
%
%USEAGE: block_filename = define_regular_block();
%
%INPUT ARGUMENTS: none;
%
% This function prompts the user for arguments to define a block parameter
% file.  It packages the information into the structures necessary to call
% the MatLab function blk_gen.m.
%
% This function makes many assumptions (below) about the block structure:
% a similar function, define_block.m, makes no assumptions
% about the block, but bugs the user for many more inputs.
%
% ASSUMPTIONS:
% The block has at most 1 radial layer with scintillator;
% There may be up to two layers of non-active material (e.g., housing,
% crystal wrap) before the scintiallator layer;
% There may be up to two layers of non-active material after the
% scintiallator layer;
% The scintillators are regularly spaced in the axial and transaxial
% dimensions (though the spacing and size can be different axially and
% transaxially);
% There can be an outer housing layer (inactive) around the scintillator
% array;
% There can be wrapping (inactive) around the scintillators in the
% transaxial and axial directions. This wrapping is the same thickness
% around each crystal, and the same thickness axially and transaxially.
%
% Robert Harrison, 2010-11-24

screen = 1;

% Request information for the file header (text)
fprintf(screen, '\n\nCreating a detector block parameter file for SimSET.\n')
fprintf(screen, '\n')

% Describe this program
fprintf(screen, 'This program creates regular scintillator arrays with\n')
fprintf(screen, 'block housing and crystal wrapping as the only additional\n')
fprintf(screen, 'options.  Another program, define_block.m, allows for more\n')
fprintf(screen, 'general arrays.  See the website for more info.\n')
fprintf(screen, '\n')

block_data.fn = input('Please enter the file name for the block parameters: ', 's');

fprintf(screen, 'Please enter the following information for the header.\n')
header.author = input('Author name: ', 's');

comment_num = 1;
fprintf(screen, 'Enter comment -\n')
fprintf(screen, 'terminate with 2 consecutive carriage returns:\n')
while (comment_num > 0)
    comment{comment_num} = input('', 's');
    if ( strcmp( comment{comment_num},'' ) )
        comment_num = 0;
    else
        comment_num = comment_num + 1;
    end
end

header.comments = comment;

% Request the information defining the block.
fprintf(screen, '\n')
fprintf(screen, 'Now define the block geometry and materials.\n')
fprintf(screen, '\n')

% Request the scintillator array dimensions.
fprintf(screen, 'How many scintillator elements (e.g. crystals) are there\n')
num_trans_xtal = input('        in the transaxial (y) direction: ');
while (num_trans_xtal < 1)
    fprintf(screen, '        There must be at least one scintillator element\n')
    num_trans_xtal = input('            in the transaxial (y) direction! : ');
end
num_axial_xtal = input('        in the axial (z) direction: ');
while (num_axial_xtal < 1)
    fprintf(screen, '        There must be at least one scintillator element\n')
    num_axial_xtal = input('            in the axial (z) direction! : ');
end

fprintf(screen, '\nWhat are the dimensions of these elements in the\n')
x_xtal_size = input('        in the radial (x) direction (cm): ');
while (x_xtal_size <= 0)
    fprintf(screen, '        Scintillator element size must be > 0\n')
    x_xtal_size = input('            in the radial (x) direction! : ');
end
y_xtal_size = input('        in the transaxial (y) direction (cm): ');
while (y_xtal_size <= 0)
    fprintf(screen, '        Scintillator element size must be > 0\n')
    y_xtal_size = input('            in the transaxial (y) direction! : ');
end
z_xtal_size = input('        in the axial (z) direction (cm): ');
while (z_xtal_size <= 0)
    fprintf(screen, '        Scintillator element size must be > 0\n')
    z_xtal_size = input('            in the axial (z) direction! : ');
end

% Request the wrap dimensions.
fprintf(screen, '\nThe transaxial and axial sides of the scintillator elements\n')
fprintf(screen, 'can have a layer of wrapping.(e.g., paint or tape).\n')
resp = input('    Are these elements wrapped? (y/n) [y]: ', 's');
if ( isempty(resp) || strcmp( resp,'y' ) || strcmp( resp,'Y' ) )
    is_wrapped = 1;
else
    is_wrapped = 0;
end
if is_wrapped
    wrap_thickness = input('    How thick is the wrap (note this will be doubled between crystals) (cm)? ');
    while (wrap_thickness <= 0)
        wrap_thickness = input('            Wrap thickness must be > 0! : ');
    end
else
    wrap_thickness = 0;
end

% Request the transaxial/axial housing dimensions.
fprintf(screen, '\nThere can be housing material on the transaxial and axial\n')
fprintf(screen, 'sides of the scintillator block.\n')
resp = input('    Is there housing material on these sides? (y/n) [y]: ', 's');
if ( isempty(resp) || strcmp( resp,'y' ) || strcmp( resp,'Y' ) )
    is_housing = 1;
else
    is_housing = 0;
end
if is_housing
    housing_thickness = input('    How thick is the housing on the axial and transaxial sides (cm)? ');
    while (housing_thickness <= 0)
        housing_thickness = input('            housing thickness must be > 0! : ');
    end
else
    housing_thickness = 0;
end

% Request the radial layer dimensions.
fprintf(screen, '\nRadially there can be up to two layers of uniform material in front of\n')
fprintf(screen, 'and behind the scintillator array (e.g. module housing).\n')

num_layers_front = input('    How many layers of material are there before the scintillator (0, 1, or 2)? ');
while ( (num_layers_front ~= 0) && (num_layers_front ~= 1) && (num_layers_front ~= 2) )
    fprintf(screen, '\n    Number of layers must be 0, 1, or 2.\n')
    num_layers_front = input('        How many layers of material are there before the scintillator (0, 1, or 2)? ');
end
if (num_layers_front > 0)
    layer_thickness(1) = input('        How thick is radial layer 1 (cm)? ');
    while (layer_thickness(1) <= 0)
        layer_thickness(1) = input('            Layer thickness must be > 0! : ');
    end
else
    layer_thickness(1) = 0;
end
if (num_layers_front == 2)
    layer_thickness(2) = input('        How thick is radial layer 2 (cm)? ');
    while (layer_thickness(2) <= 0)
        layer_thickness(2) = input('            Layer thickness must be > 0! : ');
    end
else
    layer_thickness(2) = 0;
end

layer_thickness(3) = x_xtal_size;

num_layers_back = input('    How many layers of material are there after the scintillator (0, 1, or 2)? ');
while ( (num_layers_back ~= 0) && (num_layers_back ~= 1) && (num_layers_back ~= 2) )
    fprintf(screen, '\n    Number of layers must be 0, 1, or 2.\n')
    num_layers_back = input('        How many layers of material are there after the scintillator (0, 1, or 2)? ');
end
if (num_layers_back > 0)
    layer_thickness(4) = input('        How thick is radial layer 1 (cm)? ');
    while (layer_thickness(4) <= 0)
        layer_thickness(4) = input('            Layer thickness must be > 0! : ');
    end
else
    layer_thickness(4) = 0;
end
if (num_layers_back == 2)
    layer_thickness(5) = input('        How thick is radial layer 2 (cm)? ');
    while (layer_thickness(5) <= 0)
        layer_thickness(5) = input('            Layer thickness must be > 0! : ');
    end
else
    layer_thickness(5) = 0;
end

% Get materials
fprintf(screen, '\nNow define the materials used.\n')

% Display the current material indexes for SimSET.
fprintf(screen, '\nDefault material indices from phg.data/phg_att_table:\n')
fprintf(screen, 'air                    0      con_tissue               15\n')
fprintf(screen, 'water                  1      copper                   16\n')
fprintf(screen, 'blood                  2      perfect_absorber         17\n')
fprintf(screen, 'bone                   3      LSO                      18\n')
fprintf(screen, 'brain                  4      GSO                      19\n')
fprintf(screen, 'heart                  5      aluminum                 20\n')
fprintf(screen, 'lung                   6      tungsten                 21\n')
fprintf(screen, 'muscle                 7      liver                    22\n')
fprintf(screen, 'lead                   8      fat                      23\n')
fprintf(screen, 'NaI                    9      LaBr3                    24\n')
fprintf(screen, 'BGO                    10     polycarbonateLowVisc     25\n')
fprintf(screen, 'iron                   11     polyethyleneNEMA         26\n')
fprintf(screen, 'graphite               12     polymethylMethacryl      27\n')
fprintf(screen, 'tin                    13     polystyreneFibers        28\n')
fprintf(screen, 'GI_tract               14\n')

fprintf(screen, '\nWhen using the default phg_att_table, you can use these material\n')
fprintf(screen, 'indices to answer the following.  If, however, you are using a\n')
fprintf(screen, 'custom table you should use the indices appropriate to that table.\n')

xtal_material = input('What is the scintillator material: ');
if is_wrapped
    wrap_material = input('What is the wrap material: ');
end
if (is_housing)
    housing_material = input('What is the housing material on the axial and transaxial sides: ');
end

if ( layer_thickness(1) > 0 )
    layer_material(1) = input('What is the material for the first radial layer in front of the scintillator: ');
end
if ( layer_thickness(2) > 0 )
    layer_material(2) = input('What is the material for the second radial layer in front of the scintillator: ');
end
layer_material(3) = xtal_material;
if ( layer_thickness(4) > 0 )
    layer_material(4) = input('What is the material for the first radial layer behind the scintillator: ');
end
if ( layer_thickness(5) > 0 )
    layer_material(5) = input('What is the material for the second radial layer behind the scintillator: ');
end


% The block dimensions.
% Block dimensions are computed from the layer, crystal, wrap and housing
% information given above.  In the x (radial) direction, the block is
% defined to run from x = 0 to x = sum-of-the-layer thicknesses.  In the y
% and z directions, the min is set to -max, so the thicknesses of
% everything in the transaxial or axial direction is added up, then divided
% by 2.
block_data.min(1) = 0;
block_data.max(1) = sum(layer_thickness);

block_data.max(2) = y_xtal_size * num_trans_xtal;
if is_wrapped
    block_data.max(2) = block_data.max(2) + (2* wrap_thickness * num_trans_xtal);
end
if is_housing
    block_data.max(2) = block_data.max(2) + (2* housing_thickness);
end
block_data.max(2) = block_data.max(2) / 2;
block_data.min(2) = -block_data.max(2);

block_data.max(3) = z_xtal_size * num_axial_xtal;
if is_wrapped
    block_data.max(3) = block_data.max(3) + (2* wrap_thickness * num_axial_xtal);
end
if is_housing
    block_data.max(3) = block_data.max(3) + (2* housing_thickness);
end
block_data.max(3) = block_data.max(3) / 2;
block_data.min(3) = -block_data.max(3);

% The reference point.
% The reference point is set to the middle of the front face of the
% crystal, (0,0,0)
block_data.ref(1) = 0;
block_data.ref(2) = 0;
block_data.ref(3) = 0;

% Material changes
% These are computed from the same data as the min/max

% Material changes in the x (radial) direction (also referred to as
% layer changes).
numXChanges = num_layers_front + num_layers_back;

change = 0;

if (num_layers_front > 0)
    change = change+1;
    block_data.x(change) = layer_thickness(1);
    
    if (num_layers_front == 2)
        change = change+1;
        block_data.x(change) = block_data.x(change-1) + layer_thickness(2);
    end
    
end

change = change+1;

if (num_layers_back > 0)
    block_data.x(change) = block_data.x(change-1) + layer_thickness(3);
    change = change+1;
    
    if (num_layers_back == 2)
        block_data.x(change) = block_data.x(change-1) + layer_thickness(4);
    end
    
end

% Define the change boundaries and materials for each layer

layer = 0;

% First any layers before the crystal layer
if (num_layers_front > 0)
    layer = layer+1;
    block_data.y{layer} = []; % the outside layers before/after the crystal layer are uniform.
    block_data.z{layer} = [];
    
    % Assign the material and set is_active = false
    block_data.idx{layer}(1) = layer_material(1);
    block_data.act{layer}(1) = 0;
    
    if (num_layers_front == 2)
        layer = layer+1;
        if is_housing
            % the housing on the axial and transaxial faces also covers the
            % sides of the inner layers when the are 2 layers before/after
            % the crystal layer.
            block_data.y{layer}(1) = block_data.min(2) + housing_thickness;
            block_data.y{layer}(2) = block_data.max(2) - housing_thickness;
            block_data.z{layer}(1) = block_data.min(3) + housing_thickness;
            block_data.z{layer}(2) = block_data.max(3) - housing_thickness;

            % Assign the materials and set is_active = false
            for element = 1:4
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0;
            end
            block_data.idx{layer}(5) = layer_material(2);
            block_data.act{layer}(5) = 0;
            for element = 6:9
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0;
            end
            
        else
            
            block_data.y{layer} = [];
            block_data.z{layer} = [];

            % Assign the material and set is_active = false
            block_data.idx{layer}(1) = layer_material(2);
            block_data.act{layer}(1) = 0;
            
        end
    end
    
end

% Now the crystal layer - much more complicated!
    layer = layer+1;
    
    % compute the number of material changes in the transaxial and axial
    % directions.
    num_trans_changes = num_trans_xtal - 1;
    num_axial_changes = num_axial_xtal - 1;
    if is_wrapped
        num_trans_changes = num_trans_changes + num_trans_xtal + 1;
        num_axial_changes = num_axial_changes + num_axial_xtal + 1;
    end
    if is_housing
        num_trans_changes = num_trans_changes + 2;
        num_axial_changes = num_axial_changes + 2;
    end
    
    % Now compute the material boundaries in the transaxial and axial
    % directions.
    change_y = 0;
    change_z = 0;
    y_coord = block_data.min(2);
    z_coord = block_data.min(3);
    
    % a boundary after the housing, if any
    if is_housing
        change_y = change_y + 1;
        change_z = change_z + 1;
        y_coord = y_coord + housing_thickness;
        z_coord = z_coord + housing_thickness;
        block_data.y{layer}(change_y) = y_coord;
        block_data.z{layer}(change_z) = z_coord;
    end
    
    % a boundary for the first layer of wrapping, if any
    if is_wrapped
        change_y = change_y + 1;
        change_z = change_z + 1;
        y_coord = y_coord + wrap_thickness;
        z_coord = z_coord + wrap_thickness;
        block_data.y{layer}(change_y) = y_coord;
        block_data.z{layer}(change_z) = z_coord;
    end
    
    % now for the y boundaries for all crystals up to the penultimate
    % crystal, including a double layer of wrap, if any, after each
    crystal = 1;
    while (crystal < num_trans_xtal)
        change_y = change_y + 1;
        y_coord = y_coord + y_xtal_size;
        block_data.y{layer}(change_y) = y_coord;
        if is_wrapped
            change_y = change_y + 1;
            y_coord = y_coord + (2*wrap_thickness);
            block_data.y{layer}(change_y) = y_coord;
        end
        crystal = crystal + 1;
    end
    
    % if there is wrap or housing, then the boundary after the last crystal
    % gets added.
    if (is_wrapped || is_housing)
        change_y = change_y + 1;
        y_coord = y_coord + y_xtal_size;
        block_data.y{layer}(change_y) = y_coord;
        
        % and if both are there, lay in a boundary after the wrap as well
        if (is_wrapped && is_housing)
            change_y = change_y + 1;
            y_coord = y_coord + wrap_thickness;
            block_data.y{layer}(change_y) = y_coord;
        end
    end
        
    % now the same sequence for z
    % up to the penultimate crystal...
    crystal = 1;
    while (crystal < num_axial_xtal)
        change_z = change_z + 1;
        z_coord = z_coord + z_xtal_size;
        block_data.z{layer}(change_z) = z_coord;
        if is_wrapped
            change_z = change_z + 1;
            z_coord = z_coord + (2*wrap_thickness);
            block_data.z{layer}(change_z) = z_coord;
        end
        crystal = crystal + 1;
    end
    
    % if there is wrap or housing, then the boundary after the last crystal
    % gets added.
    if (is_wrapped || is_housing)
        change_z = change_z + 1;
        z_coord = z_coord + z_xtal_size;
        block_data.z{layer}(change_z) = z_coord;
        
        % and if both are there, lay in a boundary after the wrap as well
        if (is_wrapped && is_housing)
            change_z = change_z + 1;
            z_coord = z_coord + wrap_thickness;
            block_data.z{layer}(change_z) = z_coord;
        end
    end
        
    % now fill in all those elements with materials and active/inactive
    % flags.
    
    % if there is housing, lay in a row of housing as the first z-row
    element_z = 0;
    if is_housing
        element_z = element_z + 1;
        for element_y = 1:(num_trans_changes + 1)
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end
    end
    
    % if there is wrapping, lay in a row of wrap as the next z-row
    if is_wrapped
        element_z = element_z + 1;
        element_y = 0;
        
        % first a housing element, if any
        if is_housing
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end

        for wrap_element = 1:(2*num_trans_xtal + 1)
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = wrap_material;
            block_data.act{layer}(element) = 0; % not active
        end
        
         % end with a housing element, if any
        if is_housing
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end
       
    end
    
    % now add the crystal rows, with wrap layers if is_wrapped
    for crystal = 1:num_axial_xtal
        element_z = element_z + 1;
        element_y = 0;

        % begin row with housing, if any
        if is_housing
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end
        
        % a wrap element, if any, before the first crystal element
        if is_wrapped
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = wrap_material;
            block_data.act{layer}(element) = 0; % not active
        end
        
        % now we can alternate crystal and wrap (if any) elements
        for crystal = 1:num_trans_xtal
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = xtal_material;
            block_data.act{layer}(element) = 1; % active
            if is_wrapped
                element_y = element_y + 1;
                element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
                block_data.idx{layer}(element) = wrap_material;
                block_data.act{layer}(element) = 0; % not active
            end
        end
        
        % last element in row is housing, if any
        if is_housing
            element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end
           
        % if there is wrapping, lay in a row of wrap as the next z-row
        if is_wrapped
            element_z = element_z + 1;
            element_y = 0;

            % first a housing element, if any
            if is_housing
                element_y = element_y + 1;
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0; % not active
            end

            for wrap_element = 1:(2*num_trans_xtal + 1)
                element_y = element_y + 1;
                element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
                block_data.idx{layer}(element) = wrap_material;
                block_data.act{layer}(element) = 0; % not active
            end

             % end with a housing element, if any
            if is_housing
                element_y = element_y + 1;
                element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0; % not active
            end

        end
        
    end

    % if there is housing, the final row in the layer is housing
    if is_housing
        element_z = element_z + 1;
        for element_y = 1:(num_trans_changes + 1)
            element = ( (element_z - 1)*(num_trans_changes + 1) ) + element_y;
            block_data.idx{layer}(element) = housing_material;
            block_data.act{layer}(element) = 0; % not active
        end
    end
    
% Finally any layers after the crystal layer
if (num_layers_back == 1)
    layer = layer+1;
    block_data.y{layer} = []; % the outside layer before/after the crystal layer is uniform.
    block_data.z{layer} = [];
    
    % Assign the material and set is_active = false
    block_data.idx{layer}(1) = layer_material(4);
    block_data.act{layer}(1) = 0;
    
else
    if (num_layers_back == 2)
    
        layer = layer+1;
        if is_housing
            % the housing on the axial and transaxial faces also covers the
            % sides of the inner layers when the are 2 layers before/after
            % the crystal layer.
            block_data.y{layer}(1) = block_data.min(2) + housing_thickness;
            block_data.y{layer}(2) = block_data.max(2) - housing_thickness;
            block_data.z{layer}(1) = block_data.min(3) + housing_thickness;
            block_data.z{layer}(2) = block_data.max(3) - housing_thickness;

            % Assign the materials and set is_active = false
            for element = 1:4
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0;
            end
            block_data.idx{layer}(5) = layer_material(4);
            block_data.act{layer}(5) = 0;
            for element = 6:9
                block_data.idx{layer}(element) = housing_material;
                block_data.act{layer}(element) = 0;
            end

        else

            block_data.y{layer} = [];
            block_data.z{layer} = [];

            % Assign the material and set is_active = false
            block_data.idx{layer}(1) = layer_material(4);
            block_data.act{layer}(1) = 0;

        end

        layer = layer+1;
        block_data.y{layer} = []; % the outside layers before/after the crystal layer are uniform.
        block_data.z{layer} = [];

        % Assign the material and set is_active = false
        block_data.idx{layer}(1) = layer_material(5);
        block_data.act{layer}(1) = 0;
        
    end

end


block_input = { header , block_data };

blk_out = blk_gen(block_data,header);
