%% testing virtip
% a functionalityfunctionality to generate virtual tip for setting up 
% a TapSIM simulation for a single-crystalline pillar of arbitrary crystal structure
% Markus Kühbach, 2017/12/17-19
clearvars -except XYZ;
%clear;
clc;
digits(32);
format long;

%% user interaction

simid = 4; %user-defined number to differentiate between batch runs
%simid 1 http://www.crystallography.net/cod/ Alu fcc
%simid 2 http://www.crystallography.net/cod/9008506.html Magnesium hcp
%simid 3 http://www.crystallography.net/cod/9008523.html Zirconium hcp
%simid 4 http://www.crystallography.net/cod/9008558.html Tungsten bcc
istart = 0; %how many tips to generate?
iend = 0;
prefix = ['BatchMill.SimID.' num2str(simid) '.' num2str(istart) '.' num2str(iend) ]; %istart iend

use_existent_rotations = 0;
use_existent_rotations_plus_scatter = 1;
use_identity_rotations = 0;

%% end of user interaction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if use_existent_rotations == 1 | use_existent_rotations_plus_scatter == 1
    load(['H:/Paper/Paper15_MachineLearningAPTTapSim/bb_simulation/' ...
        'VirtualTapSimTips/BatchMill.SimID.' num2str(simid) '/' ...
        'BatchMill.SimID.' num2str(simid) '.1.100.Meta.mat'],'ROTATIONS'); 
%% set up deterministic stream of random numbers outputting seed
else
    diary( [prefix '.PRNG.txt'] );
    s = RandStream('mt19937ar','Seed','shuffle'); %rand draws from global stream
    RandStream.setGlobalStream(s);
    s.get
    diary('off');
end

%% set up for all pillars fixed tip and vacuum capsule geometry
global vac;
global tip;
global base;
global lat;

set_pillar_geometry();
get_aabb();
shift_aabb_origin();
specify_crystal_lattice();
clearvars lat;

if simid == 1
    XYZ = load_template_lattice('D:\fccal_03.xyz'); %'D:\sc_01.xyz');
elseif simid == 2
    XYZ = load_template_lattice('D:\mghcp_01.xyz'); %'D:\sc_01.xyz');
elseif simid == 3
    XYZ = load_template_lattice('D:\zrhcp_01.xyz');
elseif simid == 4
    XYZ = load_template_lattice('D:\wbcc_01.xyz');
else
    disp('No simulation ID found');
end


%% batch-processing generate pillars
for pillarid=istart:1:iend
    
    %manual overwrite unrotated pattern...
    if use_identity_rotations == 1
        pillarid = 0;
    end
    
    if use_existent_rotations ~= 1 & use_existent_rotations_plus_scatter == 0
        ROTATIONS(pillarid).R = sample_random_so3(); %random rotor SO3
        %ROTATIONS(pillarid).R = eye(3,3); %no rotation
        %ROTATIONS(pillarid).R = rotation_axis(45.0, 3);
    end
    
    if use_existent_rotations_plus_scatter == 1
        %generate fixed rotation quaternion rotor about random axis
        ROTATIONS_PLUS_SCATTER(pillarid).R = mtimes(sample_random_rotor(0.5),ROTATIONS(pillarid).R);
    end

    %% rotate template lattice into temporary copy
    if use_identity_rotations == 1
        Rthis = eye(3,3);
        
        Rthis = mtimes(sample_random_rotor(0.5),eye(3,3));
        
    elseif use_identity_rotations == 1
         Rthis = ROTATIONS(pillarid).R;
    elseif use_existent_rotations_plus_scatter == 1
        Rthis = ROTATIONS_PLUS_SCATTER(pillarid).R;
    end
    RR(1:3,:) = mtimes(Rthis, XYZ(1:3,:));
    ['PillarID ' num2str(pillarid) ' template lattice rotated']

    %% build a tip and vacuum capsule
    [FINALGRID, ThinningThreshold, GuardThickness, nov, niv] = mill_tip_with_vacuum_capsule2( RR );
    
    MESHSTATS(pillarid).novac = sum(FINALGRID(4,:)==0);
    MESHSTATS(pillarid).nivac = sum(FINALGRID(4,:)==1);
    MESHSTATS(pillarid).nbot = sum(FINALGRID(4,:)==2);
    MESHSTATS(pillarid).ntip = sum(FINALGRID(4,:)>=10);
    MESHSTATS(pillarid).total = MESHSTATS(pillarid).novac + MESHSTATS(pillarid).nivac + ...
        MESHSTATS(pillarid).nbot + MESHSTATS(pillarid).ntip;
        
    visualize_result( FINALGRID, [prefix '.PillarID.' num2str(pillarid)] );

    %FINALGRID(4,:) == 1 are not accepted by TAPSIM as it internally
    %refers to the macroscopic domain wall parts, therefore reset 1s to 0s
    FINALGRID(4,FINALGRID(4,:)==1) = 0;
    sum(FINALGRID(4,:)==0)
    
    write_meshgen_conformant_result( FINALGRID, [prefix '.PillarID.' num2str(pillarid) '.txt']);
    clearvars Rthis RR FINALGRID ans niv nov 
end

clearvars XYZ RR FINALGRID ans niv nov;
save([prefix '.MetaScatter' num2str(pillarid) '.mat']);




%Rthis = ROTATIONS(1).R;
%a = [1.0; 2.0; 3.0];
%mtimes(Rthis,a)
%mtimes(eye(3,3),a)