function block_input = define_block()

%USEAGE: block_filename = define_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 no assumptions about the block structure:  a similar
% function, define_regular_block.m, bugs the user for less inputs by making
% some assumptions about the block, e.g., regularly spaced crystals of
% the same material, the same material in all gaps between the crystals.
%
% Robert Harrison, 2010-11-24

screen = 1;

% Request information for the file header (text)
fprintf(screen, 'Creating a detector block parameter file for SimSET.\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')

% The block dimensions.
fprintf(screen, 'Enter the minimum and maximum x, y, and z coordinates\n')
fprintf(screen, 'of the block in cm.  These coordinates define the\n')
fprintf(screen, 'interior of the block.  The coordinates are relative\n')
fprintf(screen, 'coordinates for the block, they do not give a position\n')
fprintf(screen, 'in the tomograph.  They are related to tomograph\n')
fprintf(screen, 'coordinates, however.\n')
fprintf(screen, 'In a the default block orientation (blocks can be\n')
fprintf(screen, 'rotated transaxially when positioned in the ring\n')
fprintf(screen, 'parameter file):\n')
fprintf(screen, '   x  increases in the radial direction;\n')
fprintf(screen, '	y increases in the clockwise direction;\n')
fprintf(screen, '	z is parallel to the tomograph z coordinate.\n')
fprintf(screen, '\n')
block_data.min(1) = input('Enter x coordinate of block front face, x minimum (cm): ');
block_data.max(1) = input('Enter x coordinate of block back face, x maximum (cm): ');
while ( block_data.min(1) >= block_data.max(1) )
    fprintf(screen, '\nInput error: x maximum must be greater than x minimum!\n')
    block_data.min(1) = input('Enter x minimum (cm): ');
    block_data.max(1) = input('Enter x maximum (cm): ');
end
fprintf(screen, '\n')
block_data.min(2) = input('Enter y coordinate of block face in counterclockwise direction, y minimum (cm): ');
block_data.max(2) = input('Enter y coordinate of block face in clockwise direction, y maximum (cm): ');
while ( block_data.min(2) >= block_data.max(2) )
    fprintf(screen, '\nInput error: y maximum must be greater than y minimum!\n')
    block_data.min(2) = input('Enter y minimum (cm): ');
    block_data.max(2) = input('Enter y maximum (cm): ');
end
fprintf(screen, '\n')
block_data.min(3) = input('Enter minimum axial coordinate for block, z minimum (cm): ');
block_data.max(3) = input('Enter maximum axial coordinate for block, z maximum (cm): ');
while ( block_data.min(3) >= block_data.max(3) )
    fprintf(screen, '\nInput error: z maximum must be greater than z minimum!\n')
    block_data.min(3) = input('Enter z minimum (cm): ');
    block_data.max(3) = input('Enter z maximum (cm): ');
end



% The reference point.
fprintf(screen, '\nA reference point within the block is used to position\n')
fprintf(screen, 'the block in the ring - in the ring parameter file,\n')
fprintf(screen, 'when a block position is specified, this is the point \n')
fprintf(screen, 'within the block placed at that position.\n\n')

block_data.ref(1) = input('Enter reference point x coordinate (cm): ');
while ( (block_data.min(1) > block_data.ref(1)) || (block_data.max(1) < block_data.ref(1)) )
    fprintf(screen, '\nInput error: x reference coordinate must be between\n')
    fprintf(screen, '        x minimum (%f) and x maximum (%f)!\n', block_data.min(1), block_data.max(1))
    block_data.ref(1) = input('Enter reference point x coordinate (cm): ');
end

block_data.ref(2) = input('Enter reference point y coordinate (cm): ');
while ( (block_data.min(2) > block_data.ref(2)) || (block_data.max(2) < block_data.ref(2)) )
    fprintf(screen, '\nInput error: y reference coordinate must be between\n')
    fprintf(screen, '        y minimum (%f) and y maximum (%f)!\n', block_data.min(2), block_data.max(2))
    block_data.ref(2) = input('Enter reference point y coordinate (cm): ');
end

block_data.ref(3) = input('Enter reference point z coordinate (cm): ');
while ( (block_data.min(3) > block_data.ref(3)) || (block_data.max(3) < block_data.ref(3)) )
    fprintf(screen, '\nInput error: z reference coordinate must be between\n')
    fprintf(screen, '        z minimum (%f) and z maximum (%f)!\n', block_data.min(3), block_data.max(3))
    block_data.ref(3) = input('Enter reference point z coordinate (cm): ');
end

% Material changes
fprintf(screen, '\nNow specify the x, y and z coordinates where there are\n')
fprintf(screen, 'material changes within the block.  We start with the x\n')
fprintf(screen, 'changes.  They must be input as a monotonically\n')
fprintf(screen, 'increasing sequence.\n\n')

% Material changes in the x (radial) direction (also referred to as
% layer changes).
numXChanges = input('Enter number of material changes in the x direction (i.e. number radial layers - 1): ');
while ( numXChanges < 0 )
    fprintf(screen, '\nInput error: the number of material changes cannot be less than 0!\n')
    numXChanges = input('Enter number of material changes in the x direction (number radial layers - 1): ');
end

if (numXChanges > 0)
    change = 1;
    block_data.x(change) = input('Enter first x material boundary/radial layer change (cm): ');
    
    while ( (block_data.min(1) >= block_data.x(change)) || (block_data.max(1) <= block_data.x(change)) )
        fprintf(screen, '\nInput error: first x material boundary must be between\n')
        fprintf(screen, '        block x minimum (%f) and x maximum (%f)!\n', block_data.min(1), block_data.max(1))
        block_data.x(change) = input('Enter first layer change (cm): ');
    end
    change = change + 1;
    
    while (change <= numXChanges)
        block_data.x(change) = input('Enter next x material boundary/radial layer change (cm): ');
    
        while ( (block_data.x(change-1) >= block_data.x(change)) || (block_data.max(1) <= block_data.x(change)) )
            fprintf(screen, '\nInput error: x material boundary must be between\n')
            fprintf(screen, '        previous x boundary (%f) and x maximum (%f)!\n', block_data.x(change-1), block_data.max(1))
            block_data.x(change) = input('Enter next layer change (cm): ');
        end
        change = change + 1;
    
    end
    
end

fprintf(screen, '\nNow each layer is defined:  first the material\n')
fprintf(screen, 'boundaries in the y znd y directions;  then the\n')
fprintf(screen, 'materials for each element these boundaries define.\n\n')

fprintf(screen, 'Within each layer, the y znd z material boundaries\n')
fprintf(screen, 'are assumed to extend across the entire layer.\n')
fprintf(screen, 'Boundaries for each dimension must be input\n')
fprintf(screen, 'as a monotonically increasing series.\n\n')

% Define the materials for each layer
for layer = 1:(numXChanges + 1)
    
    % Material changes in the y (transaxial) direction.
    fprintf(screen, 'Define material boundaries in the y (transaxial) direction\n')
    fprintf(screen, 'for layer %d.\n\n', layer)
    
    numYChanges = input('\nEnter number of material changes in the y direction: ');
    while ( numYChanges < 0 )
        fprintf(screen, '\nInput error: the number of changes cannot be less than 0!\n')
        numYChanges = input('Enter number of material changes in the y direction: ');
    end

    if (numYChanges > 0)
        change = 1;
        block_data.y{layer}(change) = input('Enter first y/transaxial material boundary (cm): ');

        while ( (block_data.min(2) >= block_data.y{layer}(change)) || (block_data.max(2) <= block_data.y{layer}(change)) )
            fprintf(screen, '\nInput error: first y material boundary must be between\n')
            fprintf(screen, '        block y minimum (%f) and y maximum (%f)!\n', block_data.min(2), block_data.max(2))
            block_data.y{layer}(change) = input('Enter first y/transaxial change (cm): ');
        end
        change = change + 1;

        while (change <= numYChanges)
            block_data.y{layer}(change) = input('Enter next y material boundary (cm): ');

            while ( (block_data.y{layer}(change-1) >= block_data.y{layer}(change)) || (block_data.max(2) <= block_data.y{layer}(change)) )
                fprintf(screen, '\nInput error: y material boundary must be between\n')
                fprintf(screen, '        previous y boundary (%f) and y maximum (%f)!\n', block_data.y{layer}(change-1), block_data.max(2))
                block_data.y{layer}(change) = input('Enter next y/transaxial change (cm): ');
            end
            change = change + 1;

        end
    else
        block_data.y{layer} = [];
    end
    
    % Material changes in the z (transaxial) direction.
    fprintf(screen, 'Define material boundaries in the z (axial) direction\n')
    fprintf(screen, 'for layer %d.\n\n', layer)
    
    numZChanges = input('\nEnter number of material changes in the z direction: ');
    while ( numZChanges < 0 )
        fprintf(screen, '\nInput error: the number of changes cannot be less than 0!\n')
        numZChanges = input('Enter number of material changes in the z direction: ');
    end

    if (numZChanges > 0)
        change = 1;
        block_data.z{layer}(change) = input('Enter first z/axial material boundary (cm): ');

        while ( (block_data.min(3) >= block_data.z{layer}(change)) || (block_data.max(3) <= block_data.z{layer}(change)) )
            fprintf(screen, '\nInput error: first z material boundary must be between\n')
            fprintf(screen, '        block z minimum (%f) and z maximum (%f)!\n', block_data.min(3), block_data.max(3))
            block_data.z{layer}(change) = input('Enter first z/axial change (cm): ');
        end
        change = change + 1;

        while (change <= numZChanges)
            block_data.z{layer}(change) = input('Enter next z material boundary (cm): ');

            while ( (block_data.z{layer}(change-1) >= block_data.z{layer}(change)) || (block_data.max(3) <= block_data.z{layer}(change)) )
                fprintf(screen, '\nInput error: z material boundary must be between\n')
                fprintf(screen, '        previous z boundary (%f) and z maximum (%f)!\n', block_data.z{layer}(change-1), block_data.max(3))
                block_data.z{layer}(change) = input('Enter next z/axial change (cm): ');
            end
            change = change + 1;

        end
    else
        block_data.z{layer} = [];
    end
    
    % Define materials for each element in the layer.
    y_elements = size(block_data.y{layer});
    z_elements = size(block_data.z{layer});
    
    % 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')

    fprintf(screen, '\nEach rectangular element in the current layer\n')
    fprintf(screen, 'must now be assigned a material.  The material\n')
    fprintf(screen, 'must also be designated as active (e.g., scintillation\n')
    fprintf(screen, 'crystal) or not active (e.g., light shielding).\n')

    for z_element = 1:(z_elements(2) + 1)
        
        if z_element == 1
            fprintf(screen, '\n    Defining row %d from z = %f to\n', z_element, block_data.min(3))
        else
            fprintf(screen, '\n    Defining row %d from z = %f to\n', z_element, block_data.z{layer}(z_element - 1))
        end
        if z_element > z_elements(2)
            fprintf(screen, '                        z = %f\n', block_data.max(3))
        else
            fprintf(screen, '                        z = %f\n', block_data.z{layer}(z_element))
        end
        
        for y_element = 1:(y_elements(2) + 1)
            
            if y_element == 1
                fprintf(screen, '\n        Define element from y = %f to\n', block_data.min(2))
            else
                fprintf(screen, '\n        Define element from y = %f to\n', block_data.y{layer}(y_element - 1))
            end
            if y_element > y_elements(2)
                fprintf(screen, '                            y = %f\n', block_data.max(2))
            else
                fprintf(screen, '                            y = %f\n', block_data.y{layer}(y_element))
            end
            
            element = y_element + ((z_element - 1)*(y_elements(2) + 1));
            
            block_data.idx{layer}(element) = input('            Enter material index for this element: ');
            resp = input('            Is this material active (e.g. scintillating)? (y/n) [y]: ', 's');
            if ( isempty(resp) || strcmp( resp,'y' ) || strcmp( resp,'Y' ) )
                block_data.act{layer}(element) = 1;
            else
                block_data.act{layer}(element) = 0;
            end
            
        end
    end    
    
end

block_input = { header , block_data };

blk_out = blk_gen(block_data,header);
