%% 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/19
clearvars -except reuse;
clc;
digits(32);
format long;

%% user interaction
simid = 1; %user-defined number to differentiate between batch runs
istart = 1; %how many tips to generate?
iend = 100;

%% end of user interaction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% set up deterministic stream of random numbers outputting seed
prefix = ['BatchMilling.SimID.' num2str(simid) '.' num2str(istart) '.' num2str(iend)];
diary( [prefix '.txt'] );
s = RandStream('mt19937ar','Seed','shuffle'); %rand draws from global stream
RandStream.setGlobalStream(s);
s.get
diary('off');

%% 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();

XYZ = load_template_lattice('D:\sc_01.xyz');

%% batch-processing generate pillars
for pillarid=istart:1:iend
    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);


    %% rotate template lattice
    Rthis = ROTATIONS(pillarid).R; %sm.UserR;
    RR(1:3,:) = mtimes(Rthis,XYZ(1:3,:));

    FINALGRID = mill_tip_with_vacuum_capsule( RR );
    
    visualize_result( FINALGRID );
    
    write_meshgen_conformant_result( FINALGRID, [prefix '.PillarID.' num2str(pillarid) '.txt']);
end
    
    


% %binary version - faulty
% header = ['BINARY ' num2str(uint32(natoms)) ' 0 0\n'];
% fprintf(fileID,header);
% ftell(fileID)
% for i=1:1:natoms
%     fwrite(fileID,[FINALGRID(1,i),FINALGRID(2,i),FINALGRID(3,i)],'float32');
%     %ftell(fileID)
%     fwrite(fileID,FINALGRID(4,i),'short');
%     if mod(i,10000) == 0
%         i
%     end
%     %ftell(fileID)
% end
% fseek(fileID,0,'cof');
% ftell(fileID)
% fwrite(fileID,FINALGRID(1,:),'float32',(2*4+1*2+0*4),'n');
% ftell(fileID)
% fseek(fileID,length(header)-1+1*4,'bof');
% ftell(fileID)
% fwrite(fileID,FINALGRID(2,:),'float32',(1*4+1*2+1*4),'n');
% ftell(fileID)
% fseek(fileID,length(header)-1+2*4,'bof');
% ftell(fileID)
% fwrite(fileID,FINALGRID(3,:),'float32',(2*4+1*2+2*4),'n');
% fclose(fileID);
[mynodefile_name ' output written, happy evaporating!']


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




cylinderfront = AABB(2,1);
cylinderrear = AABB(2,3);
cylinderleft = AABB(1,1);
cylinderright = AABB(1,2);
cylinderbottom = AABB(3,1);
cylindertop = cylinderbottom + tip.H;
tiptop = cylindertop + tip.RB;
%vectorized preselection
kick = RXYZ(3,:) < cylinderbottom | RXYZ(3,:) > tiptop | ...
    RXYZ(1,:) < cylinderleft | RXYZ(1,:) > cylinderright | ...
    RXYZ(2,:) < cylinderfront | RXYZ(2,:) > cylinderrear;
%hlfsph = kick(1,:) == 0 & XYZ(3,:) >= cylindertop;
%cyl = kick(1,:) == 0 & XYZ(3,:) < cylindertop;
natoms-sum(kick)
['Candidate ions in tip ' num2str(natoms-sum(kick))]
clearvars sXYZ;

N = 0;
Nvac = 0;
Nbase = 0;
Ntip = 0;
for i=1:natoms
    %not even in vacuum capsule irrelevant
    if varkick(1,i) == 1
        continue;
    end
%     %case separation most likely first
%     %below AABB
%     if XYZ(3,i) < cylinderbottom
%         continue;
%     end
%     %above AABB
%     if XYZ(3,i) > tiptop
%         continue;
%     end
%     %too left AABB
%     if XYZ(1,i) < cylinderleft
%         continue;
%     end
%     %too right AABB
%     if XYZ(1,i) > cylinderright
%         continue;
%     end
%     %too front AABB
%     if XYZ(2,i) < cylinderfront
%         continue;
%     end
%     %too rear
%     if XYZ(2,i) > cylinderrear
%         continue;
%     end
    % not continued? so more detailed check
    % in cylinder part?
    if RXYZ(3,i) < cylindertop
        %inside infinitesimal circular slab of truncated cone?
        Rz = tip.RB - ((tip.RB-tip.RT)/tip.H * (RXYZ(3,i)-cylinderbottom));
        %cylinder axis is z-axis, most likely in circle
        if ((RXYZ(1,i)-0)^2+(RXYZ(2,i)-0)^2) <= Rz^2 
            N = N + 1;
            sXYZ(1:3,N) = RXYZ(1:3,i);
        end %else in AABB but not in cylinder
        N
    else
        %inside halfsphere?
        Rhlfsph = tip.RT;
        if ((RXYZ(1,i)-0)^2+(RXYZ(2,i)-0)^2+(RXYZ(3,i)-cylindertop)^2) <= Rhlfsph^2
            N = N + 1;
            sXYZ(1:3,N) = RXYZ(1:3,i);
        end %else in AABB but not above halfsphere dome
        N
    end
end


%% inspect TAPSim example first bcc1105**.txt

idtypes = unique(bcc1105(:,4));
which = bcc1105(:,4) == 0; %ID zero is vacuum 
type0 = bcc1105(which,1:3);
which = bcc1105(:,4) == 2;
type2 = bcc1105(which,1:3);
which = bcc1105(:,4) == 10; %IDs >= are ions to evaporate
type10 = bcc1105(which,1:3);
figure('Position',[100 100 500 500]);
hold on
plot3(type2(:,1),type2(:,2),type2(:,3),'.','MarkerSize',2,'Color',[0 1 0]); %conductive base plate with explicit z == 0.0
hold on
plot3(type0(:,1),type0(:,2),type0(:,3),'.','MarkerSize',2,'Color',[1 0 0]); %vacuum highly resolved protection... layer
hold on
plot3(type10(:,1),type10(:,2),type10(:,3),'.','MarkerSize',2,'Color',[0 0 1]); %sample


%% inspect TAPSim example binary bcc1105**.bin
clearvars file_name fileID prefixchar floatpart shortpart np xyz id;
file_name = 'Y:\TAPSIM\tapsim_v1.0b_r3225\Test\bcc_110_5.00e-10.bin';
%see meshgen writeTapSim_Nodefile
%ASCII header has 18 characters, i.e. skip first 18*1Byte then read
%C structs float32*3 4B,short16 2B
%for the first after the char+newline header we read always the first four byte skip two fours plus one short
%immediately thereafter read again
fileID = fopen(file_name);
prefixchar = 18;
fseek(fileID, (prefixchar*1), 'bof');
floatpart = fread(fileID, [1, inf], '1*float32', (2*4+1*2+0*4), 'n'); 
np = length(floatpart(1,:));
xyz(:,1) = single(reshape(floatpart(1,:),[np 1]));
clearvars floatpart;
fseek(fileID,(prefixchar*1+1*4), 'bof');
floatpart = fread(fileID, [1, inf], '1*float32', (1*4+1*2+1*4), 'n'); %one trailing float32 + then a short surplus a leading float32
np = length(floatpart(1,:));
xyz(:,2) = single(reshape(floatpart(1,:),[np 1]));
clearvars floatpart;
fseek(fileID,(prefixchar*1+2*4), 'bof');
floatpart = fread(fileID, [1, inf], '1*float32', (0*4+1*2+2*4), 'n'); %one trailing float32 + then a short surplus a leading float32
np = length(floatpart(1,:));
xyz(:,3) = single(reshape(floatpart(1,:),[np 1]));
clearvars floatpart;
fseek(fileID,(prefixchar*1+3*4), 'bof');
shortpart = fread(fileID, [1, inf], '1*int16', (3*4), 'n'); %one trailing float32 + then a short surplus a leading float32
np = length(shortpart(1,:));
id(:,1) = uint16(reshape(shortpart(1,:),[np 1]));
clearvars shortpart;
clearvars bcc1105bin;
bcc1105bin(:,1:3) = xyz;
bcc1105bin(:,4) = id;
clearvars idtypes which type0 type2 type10

idtypes = unique(bcc1105bin(:,4));
which = bcc1105bin(:,4) == 0; %ID zero is vacuum 
type0 = bcc1105bin(which,1:3);
which = bcc1105bin(:,4) == 2;
type2 = bcc1105bin(which,1:3);
which = bcc1105bin(:,4) == 10; %IDs >= are ions to evaporate
type10 = bcc1105bin(which,1:3);

which = bcc1105bin(:,3) < 1.0e-10;
typebase = bcc1105bin(which,1:4);
hold on
plot3(typebase(:,1),typebase(:,2),typebase(:,3),'.','MarkerSize',4,'Color',[1 0 0]); 
min(typebase(1432,3))

figure('Position',[100 100 500 500]);
hold on
plot3(type2(:,1),type2(:,2),type2(:,3),'.','MarkerSize',4,'Color',[1 0 0]); %conductive base plate with explicit z == 0.0
hold on
plot3(type0(:,1),type0(:,2),type0(:,3),'.','MarkerSize',2,'Color',[0 1 0]); %vacuum highly resolved protection... layer
view(-0,-90)
view(90,0)
ylim([-15e-9 15e-9])
zlim([+50e-9 +80e-9])
hold on
view(0,90)
xlim([-15e-9 15e-9])
ylim([-15e-9 15e-9])
pbaspect([1 1 1])
hold on
these = type10(:,3) < 1.0e-8;
plot3(type10(these,1),type10(these,2),type10(these,3),'.','MarkerSize',2,'Color',[0 0 1]); %sample
min(type10(:,1))*1.0e9
max(type10(:,1))*1.0e9
min(type10(:,2))*1.0e9
max(type10(:,2))*1.0e9
min(type10(:,3))*1.0e9
max(type10(:,3))*1.0e9
BCCDIFF = bcc1105 - bcc1105bin;
max(BCCDIFF(:,1))
max(BCCDIFF(:,2))
max(BCCDIFF(:,3))
max(BCCDIFF(:,4))
% the example in fact shows the same files and thus how to also operate
% with the text file
min(type10(:,1))*1.0e9
max(type10(:,1))*1.0e9
min(type10(:,2))*1.0e9
max(type10(:,2))*1.0e9
min(type10(:,3))*1.0e9
max(type10(:,3))*1.0e9
%% write user-defined new node file currently ASCII style



% what we learn from this example?
% the input file provides the supporting points for the numerical fine
% resolved potential computation mesh, we can access it, have unstood its
% structure we can modify it to our will and needs
% to begin with we choose the build a tip with a vacuum and base plate of
% the same dimensions and enclose inside a truncated cone to cap
% so limitations currently:
% single ion type only
% fixed maximum cylinder radius and tip height

% view(0,0)
% xlim([-15e-9 15e-9])
% ylim([-2.0e-9 2.0e-9])
% what about the base plate?




% xlim([-2 2])
% ylim([-2 2])
% zlim([-2 2])


h = 10;
k = 0;
l = 0;
normhkl = h*ar + k*br + l*cr;
normhkl = normhkl ./ norm(normhkl);
normhkl
%rotate plane about origin simply by rotating plane normal vector normhkl
phi = 0.0/180.0*pi;
R = [cos(phi), -1*sin(phi), 0; sin(phi), cos(phi), 0; 0, 0, 1];
mtimes(R,normhkl)

% heuristic approach to require not so many excess ions if in practice a
% sphere centered about the center of gravity of the tip suffices
% shift grid appropriately...










% generate random rotation quaternion in SO3
phi = 0.0/180.0*pi;
R = [cos(phi), -1*sin(phi), 0; sin(phi), cos(phi), 0; 0, 0, 1];

% quaternion to matrix
% ### accept for now default


% get axis aligned bounding box to assure having enough ions to erode
AABB(1:3,1) = [min(RXYZ(1,:));min(RXYZ(2,:));min(RXYZ(3,:))];
AABB(1:3,2) = [max(RXYZ(1,:));max(RXYZ(2,:));max(RXYZ(3,:))];
AABB





%simple cubic test
sc(1:3,1) = [0;0;0];
sc(1:3,2) = [1;0;0];
sc(1:3,3) = [0;1;0];
sc(1:3,4) = [0;0;1];
sc(1:3,5) = [1;1;0];
sc(1:3,6) = [0;1;1];
sc(1:3,7) = [1;0;1];
sc(1:3,8) = [1;1;1];

phi = 45.0/180.0*pi;
R = [cos(phi), -1*sin(phi), 0; sin(phi), cos(phi), 0; 0, 0, 1];

for e=1:8
    psc(1:3,e) = mtimes(R,sc(1:3,e));
end

figure('Position',[100 100 500 500]);
grid('on')
box('on')
for e=1:8
    hold on
    plot3(sc(1,:),sc(2,:),sc(3,:),'.','MarkerSize',30,'Color',[0 0 1]);
    hold on
    plot3(psc(1,:),psc(2,:),psc(3,:),'.','MarkerSize',30,'Color',[1 0 0]);
end
pbaspect([ 1 1 1])







