function [dosecalcSetup patient_dir] = RDXTPS_dosecalcSetup(dosecalcSetup, patient_dir, Geometry)

if nargin < 1 
    dosecalcSetup = [];
elseif ~isempty(dosecalcSetup)
    inpAns = {...
        dosecalcSetup.treat_mode, ...
        num2str(dosecalcSetup.det_margin), ...
        num2str(dosecalcSetup.spread_parameters), ...
        num2str(dosecalcSetup.spacing_parameters(:)'), ...
        num2str(dosecalcSetup.aperture_size), ...
        num2str(dosecalcSetup.beam_angles)};
end

if nargin < 2
    patient_dir = [];
end

%%
prompt = {...
    sprintf('All unit in cm\n\nTreatment mode (SS or DET):'), ...
    sprintf('DET margin (default 0.1):'), ...
    sprintf('Beam width (1 Gaussian sigma), typically 0.3 to 1.0:'), ...
    sprintf('Spot spacing, typically <= 1 sigma:\n(3 numbers seperated by space, the last one along beam direction)'), ...
    sprintf('Cutoff distance (aperture size) of the dose calculation grid:'), ...
    sprintf('Gantry angles (any scalar or vector in MATLAB syntax):\ne.g. -30 0 30 or -180:10:180\n0 is PA') ...
    sprintf('Couch angles (any scalar or vector in MATLAB syntax):\ne.g. -30 0 30 or -180:10:180\n0 for coplanar beam') ...
    };
defAns = {...
    'SS', ...
    '0.1', ...
    '0.6', ...
    '0.5 0.5 0.3', ...
    '4', ...
    '-90 90', ...
    '0 0'};

%% Prompt input dialog
if isempty(dosecalcSetup)
    answer = inputdlg(prompt, 'Dose calculation parameters', 1, defAns);
else
    answer = inputdlg(prompt, 'Dose calculation parameters', 1, inpAns);
end

%% TODO if hit cancel

%% Convert answer to struct
dosecalcSetup.treat_mode = answer{1};
dosecalcSetup.det_margin = str2double(answer{2});
dosecalcSetup.spread_parameters = str2double(answer{3});
dosecalcSetup.spacing_parameters = sscanf(answer{4}, '%f%f%f');
dosecalcSetup.aperture_size = str2double(answer{5});
dosecalcSetup.beam_angles = eval(['[' answer{6} ']']);
dosecalcSetup.couch_angles = eval(['[' answer{7} ']']);

%% Check input data integrity
dosecalcSetup.treat_mode = strtrim(dosecalcSetup.treat_mode);
if ~strcmpi(dosecalcSetup.treat_mode, 'SS') && ~strcmpi(dosecalcSetup.treat_mode, 'DET')
    msgbox('Unknow treatment mode specified');
end

if isempty(dosecalcSetup.det_margin) || dosecalcSetup.det_margin < 0 || dosecalcSetup.det_margin > 10
    msgbox('DET margin out of range');
end

if isempty(dosecalcSetup.spread_parameters) || dosecalcSetup.spread_parameters < 0 || dosecalcSetup.spread_parameters > 10
    msgbox('Beam width out of range');
end

if isempty(dosecalcSetup.spacing_parameters) || any(dosecalcSetup.spacing_parameters < 0) || any(dosecalcSetup.spacing_parameters > 10)
    msgbox('Spot spacing out of range');
end

if isempty(dosecalcSetup.aperture_size) || dosecalcSetup.aperture_size < 0 || dosecalcSetup.aperture_size > 100
    msgbox('Aperture out of range');
end

if isempty(dosecalcSetup.beam_angles) || any(~isnumeric(dosecalcSetup.beam_angles))
    msgbox('Gantry angle out of range');
end

%% Choose output directory
if isempty(patient_dir)
    patient_dir = uifile('getdir', 'Select the patient directory');
end

try
    load(fullfile(patient_dir, 'matlab_files', 'Geometry.mat'));
catch
    msgbox('Error loading Geometry file');
    return;
end

%% prepare variables for input file
geometry_rhomw_filename     = fullfile(patient_dir, 'geometry_files', 'rhomw.bin');
geometry_Smw_filename       = fullfile(patient_dir, 'geometry_files', 'Smw.bin');
geometry_tumormask_filename = fullfile(patient_dir, 'geometry_files', 'target_mask.bin');
geometry_Fmw2_filename      = fullfile(patient_dir, 'geometry_files', 'Fmw2.bin');
batch_dose = fullfile(patient_dir, 'batch_dose.bin');
beam_spec = fullfile(patient_dir, 'beam_spec_DAVE.txt');

% size of the CT data grid
[M,N,Q] = size(Geometry.rhomw);

% voxel size
dx = double(Geometry.voxel_size(1));
dy = double(Geometry.voxel_size(2));
dz = double(Geometry.voxel_size(3));

% grid start location
x0 = double(Geometry.start(1));
y0 = double(Geometry.start(2));
z0 = double(Geometry.start(3));

% tumor zstart, zend. NOT used any more
tumor_zstart = 0;
tumor_zend = 1;

read_beamspec = 0;
write_beamspec = 1;
use_scale2D = 1;
Mxp = 512;
Nyp = 512;
SAD = 100;
num_energies = 576;
pencil_beam_filenames_file = 'pencil_beam_files/pencil_beam_filenames.txt';
pencil_beam_path = 'pencil_beam_files';
    
sig_xp = dosecalcSetup.spread_parameters;
sig_yp = dosecalcSetup.spread_parameters;
xSpace = dosecalcSetup.spacing_parameters(1);
ySpace = dosecalcSetup.spacing_parameters(2);
zSpace = dosecalcSetup.spacing_parameters(3);
del_xp = dosecalcSetup.aperture_size;
del_yp = dosecalcSetup.aperture_size;
num_angles = numel(dosecalcSetup.beam_angles);
phi = dosecalcSetup.beam_angles;
couch_phi = dosecalcSetup.couch_angles;

%% write input file
fprintf('Generating input file.\n');

fid = fopen(fullfile(patient_dir, 'dosecalc_input.txt'), 'w');
fprintf(fid,'#######################################################################\n');
fprintf(fid,['## Proton dose calculator input file for: ' patient_dir '\n']);
fprintf(fid,'#######################################################################\n\n');
fprintf(fid,'# Text entries not preceeded by comments must not be changed!!\n\n');
fprintf(fid,'# GEOMETRY DATA\n\n');
fprintf(fid,'# CT density grid file\n');
fprintf(fid,'rhomw_grid_file\n');
fprintf(fid,'%s\n\n',geometry_rhomw_filename);
fprintf(fid,'# Stopping power ratio grid file\n');
fprintf(fid,'Smw_grid_file\n');
fprintf(fid,'%s\n\n',geometry_Smw_filename);
fprintf(fid,'# Tumor mask file\n');
fprintf(fid,'mask_file\n');
fprintf(fid,'%s\n\n',geometry_tumormask_filename);
fprintf(fid,'# Grid containing 2D scaling F^2 factors\n');
fprintf(fid,'Fmw2_grid_file\n');
fprintf(fid,'%s\n\n',geometry_Fmw2_filename);
fprintf(fid,'# Output dose file\n');
fprintf(fid,'batch_dose_file\n');
fprintf(fid,'%s\n\n',batch_dose);
fprintf(fid,'# Beam specification file\n');
fprintf(fid,'beam_spec_file\n');
fprintf(fid,'%s\n\n',beam_spec);
fprintf(fid,'# Data grid dimensions (Xdim,Ydim,Zdim)\n');
fprintf(fid,'grid_dimensions\n');
fprintf(fid,'%d %d %d\n\n',M,N,Q);
fprintf(fid,'# Voxel dimensions (dx,dy,dz) [cm]\n');
fprintf(fid,'voxel_dimensions\n');
fprintf(fid,'%f %f %f\n\n',dx,dy,dz);
fprintf(fid,'# Grid start coordinates (shifted) (x0,y0,z0) [cm]\n');
fprintf(fid,'start_coordinates\n');
fprintf(fid,'%f %f %f\n\n',x0,y0,z0);
fprintf(fid,'# Tumor Z coordinates (tumor_zstart, tumor_zend) [cm]\n');
fprintf(fid,'tumor_coordinates\n');
fprintf(fid,'%f %f\n\n',tumor_zstart,tumor_zend);
fprintf(fid,'# USER DEFINED BEAM DATA\n\n');
fprintf(fid,'# Calculation mode\n');
fprintf(fid,'treat_mode\n');
fprintf(fid,'%s\n\n', dosecalcSetup.treat_mode);
fprintf(fid,'# DET margin\n');
fprintf(fid,'det_margin\n');
fprintf(fid,'%f\n\n', dosecalcSetup.det_margin);
fprintf(fid,'# Load beam specification file [1 for yes, 0 for no]\n');
fprintf(fid,'read_beam_spec_file\n');
fprintf(fid,'%d\n\n',read_beamspec);
fprintf(fid,'# Write beam specification file [1 for yes, 0 for no]\n');
fprintf(fid,'write_beam_spec_file\n');
fprintf(fid,'%d\n\n',write_beamspec);
fprintf(fid,'# Use 2D scaling method [1 for yes, 0 for no]\n');
fprintf(fid,'use_2D_scaling\n');
fprintf(fid,'%d\n\n',use_scale2D); 
fprintf(fid,'# Gaussian spread parameters (sig_xp,sig_yp) [cm]\n');
fprintf(fid,'spread_parameters\n');
fprintf(fid,'%f %f\n\n',sig_xp,sig_yp);
fprintf(fid,'# Gaussian spacing parameters (xSpace,ySpace,zSpace) [cm]\n');
fprintf(fid,'spacing_parameters\n');
fprintf(fid,'%f %f %f\n\n',xSpace,ySpace,zSpace);
fprintf(fid,'# Rectangular aperture size in beams eye view (del_xp,del_yp) [cm]\n');
fprintf(fid,'aperture_size\n');
fprintf(fid,'%f %f\n\n',del_xp,del_yp);
fprintf(fid,'# Number of beam angles\n');
fprintf(fid,'num_beam_angles\n');
fprintf(fid,'%d\n\n',num_angles);
fprintf(fid,'# Beam angles for calculation [degrees]\n');
fprintf(fid,'gantry_angles\n');
fprintf(fid,'%f\n',phi);
fprintf(fid,'\n');
fprintf(fid,'# Couch angles for calculation [degrees]\n');
fprintf(fid,'couch_angles\n');
fprintf(fid,'%f\n',couch_phi);
fprintf(fid,'\n');
fprintf(fid,'# Number of beamlets filling each beam angle aperture (Mxp,Nyp)\n');
fprintf(fid,'number_beamlets\n');
fprintf(fid,'%d %d\n\n',Mxp,Nyp);
fprintf(fid,'# Source to axis distance (SAD) [cm]\n');
fprintf(fid,'SAD\n');
fprintf(fid,'%f\n\n',SAD);
fprintf(fid,'# PENCIL BEAM DATA\n\n');
fprintf(fid,'# Number of pencil beam energies\n');
fprintf(fid,'number_energies\n');
fprintf(fid,'%d\n\n',num_energies);
fprintf(fid,'# Pencil beam filenames file\n');
fprintf(fid,'pencil_beam_filenames\n');
fprintf(fid,'%s\n\n',pencil_beam_filenames_file);
fprintf(fid,'# Location of pencil beam files\n');
fprintf(fid,'pencil_beam_path\n');
fprintf(fid,'%s\n',pencil_beam_path);
fclose(fid);

% save data to matlab structure
save(fullfile(patient_dir, 'matlab_files', 'dosecalcSetup.mat'), 'dosecalcSetup');