gtFillCrystFields.m

function parameters = gtFillCrystFields(parameters, phase_ndx)
% parameters = gtFillCrystFields(parameters, phase_ndx)
% fill in the cryst(phase_ndx)
% handle all crystallography stuff in one easy function...
%
% concept - we need certain information about each phase, and we have to get
% this from somewhere.  This could be from a .dat file created by XOP, a
% .cif file from Diamond or somewhere else, or could be supplied from the
% head of the user - we should not be fussy.
%
% The user should not have to copy anything - as long as the files are on
% disk, we should move them to the right place
% For the moment at least, keep the concept of archiving a dat and cif
% file for every dataset in id11/graintracking/
%
% Information to deal with:
%
% We already have:
% energy, max and min two theta values
%
% We need to get:
% lattice parameters
% spacegroup
% {hkl} or {hkil} reflections
%
% Optional, but nice, are:
% calculated intensities
% symmetry operators (although we don't use this info yet...)
%
% We generate ourselves:
% crystal system and lattice type we can get from spacegroup number
% signed hkls we generate by applying the (hard coded) symmetry operators
% d spaces from hkls and lattice parameters
% theta from dspacing and energy
%

disp([' ~~~~~~~~~~~~~ Get crystallography information for phase ' ...
    num2str(phase_ndx) ' ~~~~~~~~~~~~~ '])

% now fill in this information for this phase
p = make_parameters(2);
phase_ok = false;
while (~phase_ok)
    % is this part of a series of measurements?  In this case we can simply
    % copy from another parameters file
    disp(['If you have already measured this material with same geometry and conditions, ' ...
        'you can copy from another parameters file.']);
    disp(['If the experimental geometry is different, '...
        'it may be safer to load from a data file prepared for the specific conditions.']);

    test = inputwdefault('Do you want to load a parameters file? [y/n]', 'n');
    if (ismember(lower(test), {'y', 'yes'}))
        %%%% get complete information from another parameters file %%%%%
        parameters = sfGetInfoFromParameters(parameters, phase_ndx);
		
		% Update reflections for actual photon energy
		% !!! intensity missing, but it is not yet used in the processing 
		parameters.cryst = sfUpdateReflections(parameters.cryst, parameters.acq.energy);
		
        % this is a "new" experiment, set up the names and read from the datafile
        fprintf('Please modify the description for this phase (%d)\n', phase_ndx);
        parameters = sfModifyMaterialDescription(parameters, phase_ndx);
    else
        if (~isempty(parameters.xop(phase_ndx).filename))
            xop_filename = fullfile(parameters.xop(phase_ndx).xop_dir, ...
                parameters.xop(phase_ndx).filename);
            question = sprintf(['Do you want to reset crystallographic information' ...
                ' previously loaded from file "%s", for phase: %02d? [y/n]'], ...
                xop_filename, phase_ndx);
            if (~ismember(inputwdefault(question, 'y'), {'y', 'yes'}))
                parameters.xop(phase_ndx)   = p.xop;
                parameters.cryst(phase_ndx) = p.cryst;
            end
        end

        % this is a "new" experiment, set up the names and read from the datafile
        fprintf('Please give a description for this phase (%d)\n', phase_ndx);
        parameters = sfModifyMaterialDescription(parameters, phase_ndx);

        % create the xop directory in which to save the files used
        [xop_dir, ~] = gtCrystXopInitializeDirectory(parameters.acq.name, ...
            parameters.cryst(phase_ndx).name);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Get crystal system and lattice from the spacegroup number
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Get crystallographic data
        disp('You can either load crystallographic info (lattice pars, symmetry) from a file or enter it manually.');
        query = 'Do you want to load it from a .cif file? [y/n]';
        test = inputwdefault(query, 'n');
        if (ismember(lower(test), {'y', 'yes'}))
            parameters = sfGetCrystInfoFromDataFiles(parameters, phase_ndx);
        else
            parameters = sfGetCrystInfoFromUser(parameters, phase_ndx);
        end

        [lattice_system, crystal_system, herm_symbol] = gtReadSpaceGroup(parameters.cryst(phase_ndx).spacegroup);
        parameters.cryst(phase_ndx).lattice_system = lattice_system;
        parameters.cryst(phase_ndx).crystal_system = crystal_system;
        parameters.cryst(phase_ndx).hermann_mauguin = herm_symbol;
        clear lattice_system crystal_system herm_symbol;

        %%%% get information from files or user %%%%%
        % Get reflection data

        query = 'Do you want to load reflection data from .csv or .dat files? (recommended) [y/n]';
        test = inputwdefault(query, 'n');
        if (ismember(lower(test), {'y', 'yes'}))
            parameters = sfGetReflectionInfoFromDataFiles(parameters, phase_ndx, xop_dir);
            % expand the input data
            parameters = sfExpandData(parameters, phase_ndx);
            parameters = sfGetSymmCrystInfo(parameters, phase_ndx);
        else
            query = 'Do you want to generate reflection data using FABLE? [y/n]';
            test = inputwdefault(query, 'y');
            if (ismember(lower(test), {'y', 'yes'}))
                parameters = sfGetReflectionInfoFromFable(parameters, phase_ndx);
                parameters = sfGetSymmCrystInfoFromFable(parameters, phase_ndx);
            else
                parameters = sfGetReflectionInfoFromUser(parameters, phase_ndx);
                % expand the input data
                parameters = sfExpandData(parameters, phase_ndx);
                parameters = sfGetSymmCrystInfo(parameters, phase_ndx);
            end
        end

    end % not loading parameters file
    
    % Check things look okay
    query = 'Do you want to look at the contents of cryst for this phase? [y/n]';
    test = inputwdefault(query, 'y');
    if (ismember(lower(test), {'y', 'yes'}))
        parameters = sfCheckWithUser(parameters, phase_ndx);
    end

    % checking parameters for this phase - to be checked
    tmp_pars.cryst = parameters.cryst(phase_ndx);
    tmp_pars.xop   = parameters.xop(phase_ndx);

    tmp_pars = gtCheckParameters(tmp_pars, 'cryst', 'verbose', true);
    tmp_pars = gtCheckParameters(tmp_pars, 'xop', 'verbose', true);
    
    parameters.cryst(phase_ndx) = tmp_pars.cryst;
    parameters.xop(phase_ndx)   = tmp_pars.xop;
    clear tmp_pars
    
    % does this look ok?
    query = 'Is this phase OK? [y/n] ';
    test = input(query, 's');
    if (ismember(lower(test), {'y', 'yes'}))
        phase_ok = true;
        disp('');
        
        disp('Saving parameters.mat in the analysis directory...')
        parameters_name = fullfile(parameters.acq.dir,'parameters.mat');
        save(parameters_name, 'parameters')
        disp([' ~~~~~~~~~~~~~ Finished getting crystallography information for phase ' ...
            num2str(phase_ndx) ' ~~~~~~~~~~~~~ '])
    end
end % end while loop

end % end of the main function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function parameters = sfCheckWithUser(parameters, phase_id)
% show the user the cryst field so they can decide if it looks ok
    disp(parameters.cryst(phase_id));
    
end % end of sfCheckWithUser

function parameters = sfGetInfoFromParameters(parameters, phase_id)
% in this case, will fill cryst directly

    % Get the master parameters
    [masterfile, masterpath] = uigetfile({'*.mat', 'MAT-files (*.mat)'}, ...
        'Select a parameters file');
    if ~(strcmp(masterfile, 'parameters.mat'))
        warning('gtCrystFillFiels:sfGetInfoFromParameters:strange_file', ...
            ['File should be named ''parameters.mat'', but got ''' masterfile '''!']);
    end
    master = load(fullfile(masterpath, masterfile));
    if ~isfield(master, 'parameters')
        gtError('gtFillCrystFields:wrong_file', ...
            'Given file doesn''t DCT parameters!');
    elseif ~isfield(master.parameters, 'cryst')
        gtError('gtFillCrystFields:wrong_file', ...
            'Given parameters file doesn''t contains any crystallographic info!');
    end
    % If there is more than one phase
    nPhases = length(master.parameters.cryst);
    if (nPhases > 1)
        disp('Which phase to use?');
        for m = 1:length(master.parameters.cryst)
            sprintf(' (%d) %s\n', m, master.parameters.cryst(m).name);
        end
        mm = input(['[1' sprintf('/%d', 2:nPhases) ']']);
    else
        mm = 1;
    end
    % Add these values to our parameters field
    parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id),master.parameters.cryst(mm),true,true,true);
    % keep track of where we got the data from

    % create the xop directory in which to save the files used
    [xop_dir, ~] = gtCrystXopInitializeDirectory(parameters.acq.name, parameters.cryst(phase_id).name);
    parameters.xop(phase_id).xop_dir = xop_dir;
    parameters.xop(phase_id).filename = 'master_parameters.mat';
    % copy the parameters file that we used as the reference
    copyfile(fullfile(masterpath, masterfile), fullfile(xop_dir, 'master_parameters.mat'))
    % return to the main loop to allow user to check manually

end % end of sfGetInfoFromParameters

function parameters = sfModifyMaterialDescription(parameters, phase_id)
% get this stuff via gtModfiyStructure
% Usertype = 1 for material, name, composition

    list = build_list_v2();
    header = sprintf('Manual input of crystallographic parameters for phase %d:', phase_id);
    parameters.cryst(phase_id) = gtModifyStructure(parameters.cryst(phase_id), list.cryst, 1, header);

end % end of sfModifyMaterialDescription

function parameters = sfGetReflectionInfoFromDataFiles(parameters, phase_id, xop_dir)
% get a data file to try
    
    disp('Please select a data file from which to load powder diffraction data');
    % pick a suitable datafile
    [datafile, datapath] = uigetfile('*.csv; *.dat', ...
        'Select a powder diffraction data file', fullfile('/data','id11','archive','file_xop'));
    % try to read from it
    try
        xop = gtLoadReflectionsFromFile(fullfile(datapath, datafile));
        parameters.xop(phase_id) = gtAddMatFile(parameters.xop(phase_id), xop, true, false, false);
        clear xop;
        % add this data to parameters.cryst
    catch mexc
        gtPrintException(mexc, 'Did not manage to read from the selected file!');
        disp('Keep trying - possibly with a different datafile or manual input.');
        return;
    end

    try
        % copy the file that we used
        copyfile(fullfile(datapath, datafile), fullfile(xop_dir, datafile))
        parameters.xop(phase_id).xop_dir = xop_dir;
        parameters.xop(phase_id).filename = datafile;
        disp('Read some powder diffraction data');
    catch mexc
        gtPrintException(mexc, ...
            'Did not manage to copy the file to the shared directory! This is not a big issue...');
    end

end % end of sfGetReflectionInfoFromDataFiles

function parameters = sfGetCrystInfoFromDataFiles(parameters, phase_id)
% need to get spacegroup and lattice parameters from a .cif file
% this goes into parameters.cryst

    [ciffile, cifpath] = uigetfile('*.cif', ...
        'Select a suitable crystallography file', fullfile('/data','id11','archive','file_cif'));
    cif = gtCrystCifInitialize(fullfile(cifpath, ciffile));
    parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id), cif, true, true, false);
    clear cif;

    % copy the file that we used to the xop_dir
    unix(sprintf('cp %s %s',fullfile(cifpath, ciffile), fullfile(parameters.acq.dir, ciffile)));

end % end of sfGetCrystInfoFromDataFiles

function parameters = sfGetReflectionInfoFromFable(parameters, phase_id)
% need to get spacegroup and lattice parameters
% this goes into parameters.cryst
    disp('Using fable calculation of theoretical reflections')
    disp('The only missing information is the intensity')
    % reflection list
    list = gtCrystCalculateReflections(parameters.cryst(phase_id), parameters.detgeo, parameters.acq.energy, 'reflections_');
    parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id), list, true, false, false, false, false);

    clear list
end% end of sfGetReflectionInfoFromFable

function parameters = sfGetSymmCrystInfoFromFable(parameters, phase_id)
    % symmetry operators from fable calculation
    [symm_unique, symm_all] = gtCrystCalculateSymmetryOperators(phase_id, ...
        parameters.cryst(phase_id).spacegroup, ...
        parameters.cryst(phase_id).crystal_system);

    parameters.cryst(phase_id).symm_unique = symm_unique;
    parameters.cryst(phase_id).symm_all = symm_all;
end

function parameters = sfGetReflectionInfoFromUser(parameters, phase_id)
% get info interactively

    % We need to get:
    % {hkl} or {hkil} reflections
    % lattice parameters
    % spacegroup
    % get the list for the fields we are interested in
    % in parameters.cryst

    % no datafile to copy in this case

    % in parameters .xop
    hkl = input('Manually input the hkl families as an nx3 matrix: [h1 k1 l1; h2 k2 l2; etc]');
    parameters.xop(phase_id).hkl = hkl;
    disp('Filling parameters.xop.int and mult with ones');
    parameters.xop(phase_id).int = ones(size(hkl));
    parameters.xop(phase_id).mult = ones(size(hkl));

end % end of sfGetReflectionInfoFromUser

function parameters = sfGetCrystInfoFromUser(parameters, phase_id)
% need to get spacegroup and lattice parameters from a .cif file
% Usertype = 2 for lattice parameters and spacegroup

    list = build_list_v2();
    header = sprintf('Manual input of crystallographic parameters for phase %d', phase_id);
    parameters.cryst(phase_id) = gtModifyStructure(parameters.cryst(phase_id), list.cryst, 2, header);

end % end of sfGetCrystInfoFromUser

function parameters = sfExpandData(parameters, phase_id)
% file parameters.cryst from reflection data in parameters.xop and symmetry
% and lattice parameters in parameters.cryst

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % getting reflection on the detector, basing on theta max and min
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % after xop
    % this cleans the xop output into distinct hkl families
    % and recalculates dspacings and thetas according to parameters
    % discards anything not on the detector
    % this info goes into the parameters.cryst field
    detector = gtDetectorTwotheta(parameters, phase_id);
    parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id), detector, true, true, false);
    clear detector;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using symmetry operators getting all the possible reflections
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % this expands the hkl families in parameters.cryst
    hklsymmetry = gtSymmetricReflections(parameters.cryst(phase_id), parameters.acq.energy);
    parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id), hklsymmetry, true, true, false);
    clear hklsymmetry;
end % end of sfExpandData

function parameters = sfGetSymmCrystInfo(parameters, phase_id)
    % Calculating symmetry operators
    parameters.cryst(phase_id).symm = gtCrystGetSymmetryOperators( ...
        parameters.cryst(phase_id).crystal_system, ...
        parameters.cryst(phase_id).spacegroup );
end

function cryst = sfUpdateReflections(cryst, energy)
% Updates the reflection data according to a new photon energy 
    % Update the d-spacing info in cryst according to new lattice parameters.
    % Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
	% cryst.dspacing   = gtCrystDSpacing(cryst.hkl, Bmat);
	% cryst.dspacingsp = gtCrystDSpacing(cryst.hklsp, Bmat);

	% Update the theta info according to new photon energy
    lambda = gtConvEnergyToWavelength(energy);
	cryst.theta   = asind(0.5*lambda./cryst.dspacing);
	cryst.thetasp = asind(0.5*lambda./cryst.dspacingsp);
	
	% Update the intensity info: MISSING !
	cryst.int   = [];
	cryst.intsp = [];
end