Upload code file.

Author: phcoding

README.md

@@ -0,0 +1,4 @@
+# 晶格系统可视化
+1. 代码文件中 `main_surface2.m` 和 `main_wigner.m` 为可执行脚本，其它皆为使用到的函数。两文件中，前者为生成三种晶系对应正格子和倒格子的代码，后者为生成三种晶系对应的WS原胞的代码。绘制某种晶格的正格子或者倒格子需要将对应部分代码取消注释。整体代码稍显繁杂，没有做重构工作，具体细节自行研究。
+
+2. gif文件夹下为上述3D图动态旋转动画，可用im2gif函数生成。

coor_trans.m

@@ -0,0 +1,8 @@
+function [ cc ] = coor_trans( c, trans_mat )
+%COOR_TRANS Get Cartesian coordinate from another coordinate by
+%transformation matrix.
+%   Detailed explanation goes here
+
+cc = c*trans_mat;
+end
+

create_atom.m

@@ -0,0 +1,19 @@
+function [ X, Y, Z, C ] = create_atom( coordinate, r)
+%CREATE_ATOM Create one atom and return surface matrix data.
+%   Detailed explanation goes here
+
+range = 100;
+% phi =linspace(0,2*pi,range);
+% thet = linspace(0,pi,range);
+% [phi,thet] = meshgrid(phi, thet);
+% X = r*sin(thet).*cos(phi);
+% Y = r*sin(thet).*sin(phi);
+% Z = r*cos(thet);
+% C = sin(X+Y+Z);
+[X,Y,Z] = sphere(range);
+C = sin((X+Y+Z)/3);
+X = r*X + coordinate(1);
+Y = r*Y + coordinate(2);
+Z = r*Z + coordinate(3);
+end
+

create_bone.m

@@ -0,0 +1,29 @@
+function [ X, Y, Z, C ] = create_bone( cstart, cend, radius, varargin )
+%CREATE_BONE Create one bone between two point on a special coordinate
+%described by transformation matrix and return surface matrix data on
+%Cartesian  coordinates.
+%   Detailed explanation goes here
+
+if length(varargin) == 1
+    trans_mat = varargin{1};
+else
+    trans_mat = [1,0,0;0,1,0;0,0,1];
+    disp('Warning: trans_mat param error ! defult to using cartesian coordinates !');
+end
+cstart_c = cstart*trans_mat;
+cend_c = cend*trans_mat;
+vector_c = cend_c - cstart_c;
+l = norm(vector_c);
+t = acos(vector_c(3)/l);
+if norm(vector_c(1:2)) ~= 0
+    if vector_c(2) >= 0
+        p = acos(vector_c(1)/norm(vector_c(1:2)));
+    else
+        p = -acos(vector_c(1)/norm(vector_c(1:2)));
+    end
+else
+    p = 0;
+end
+[X, Y, Z, C] = create_bones(cstart_c, t, p, l, radius);
+end
+

create_bones.m

@@ -0,0 +1,25 @@
+function [ X, Y, Z, C ] = create_bones( coordinate, t, p, l, r)
+%CREATE_BONES Create one bone whitch starts from a point and described by
+%Ball coordinates with radius of r.
+%   Detailed explanation goes here
+
+range = 20;
+phi = linspace(0,2*pi, range);
+h = linspace(0,l,range);
+[phi, h] = meshgrid(phi, h);
+X = r*cos(phi);
+Y = r*sin(phi);
+Z = h;
+
+X_ = cos(t)*X+sin(t)*Z;
+Z = -sin(t)*X+cos(t)*Z;
+
+X = cos(p)*X_-sin(p)*Y;
+Y = sin(p)*X_+cos(p)*Y;
+
+X = X+coordinate(1);
+Y = Y+coordinate(2);
+Z = Z+coordinate(3);
+C = phi;
+end
+

create_crystal_system.m

@@ -0,0 +1,12 @@
+function [ crystal_system ] = create_crystal_system( trans_mat, limits )
+%CREATE_CRYSTAL_SYSTEM Summary of this function goes here
+%   Detailed explanation goes here
+
+[x,y,z] = meshgrid(limits(1,1):limits(1,2), limits(2,1):limits(2,2), ...
+    limits(3,1):limits(3,2));
+x = reshape(x, [numel(x),1]);
+y = reshape(y, [numel(y),1]);
+z = reshape(z, [numel(z),1]);
+crystal_system = [x,y,z]*trans_mat;
+end
+

create_parallelogram.m

@@ -0,0 +1,7 @@
+function [ X,Y,Z,C ] = create_parallelogram( points, trans_mat )
+%CREATE_PARALLELOGRAM Summary of this function goes here
+%   Detailed explanation goes here
+[X,Y,Z,C] = create_surface([points(1,:);points(2,:);...
+    points(2,:)+points(3,:)-points(1,:);points(3,:)], trans_mat);
+end
+

create_surface.m

@@ -0,0 +1,11 @@
+function [ X, Y, Z, C ] = create_surface( points, trans_mat )
+%CREATE_SURFACE Summary of this function goes here
+%   Detailed explanation goes here
+
+points = points*trans_mat;
+X = [points(:,1)',points(1,1)];
+Y = [points(:,2)',points(1,2)];
+Z = [points(:,3)',points(1,3)];
+C = sqrt(pow2(Z)+pow2(Y)+pow2(X));
+end
+

frames.m

@@ -0,0 +1,37 @@
+function [ fs ] = frames( varargin )
+%FRAMES 将输入资源转换为frame数据组
+%   fs = frames(...)
+%   输入可以为图片文件路径、FIG文件路径、图片数据和frame影片数据
+len = length(varargin);
+fs = struct('cdata',num2cell(zeros(1,len,'uint8')),'colormap',num2cell(zeros(1,len,'double')));
+for i=1:len
+    e = varargin{1,i};
+    switch(class(e))
+        case 'char'
+            filename = strsplit(e,'.');
+            fmt = upper(filename{end});
+            try
+                switch(fmt)
+                    case {'PNG','JPG','JPEG','BMP'}
+                        fs(1,i) = im2frame(imread(e));
+                    case {'FIG'}
+                        open(e);
+                        fs(1,i) = getframe;
+                        close gcf;
+                    otherwise
+                        warning(['filetype ' fmt ' not surpported !']);
+                end
+            catch
+                warning(['resource "' e '" load failed ![INDEX ' num2str(i) ']']);
+            end
+        case 'struct'
+            fs(1,i) = e;
+        case 'uint8'
+            try
+                fs(1,i) = im2frame(e);
+            catch
+                warning(['imdata load failed !' '[INDEX ' num2str(i) ']']);
+            end
+    end
+end
+end

get_rl.m

@@ -0,0 +1,14 @@
+function [ rl_tm ] = get_rl( pl_tm )
+%GET_RL Get Reciprocal Lattice from Positive crystal lattice
+%transformation matrix.
+%   RL_TM = GET_RL( PL_TM )
+%   argsin:
+%       pl_tm       Positive crystal lattice transformation matrix.
+%   argsout:
+%       rl_tm       Reciprocal Lattice transformation matrix.
+
+rl_tm(1,:) = 2*pi*cross(pl_tm(2,:), pl_tm(3,:))/dot(pl_tm(1,:), cross(pl_tm(2,:), pl_tm(3,:)));
+rl_tm(2,:) = 2*pi*cross(pl_tm(3,:), pl_tm(1,:))/dot(pl_tm(2,:), cross(pl_tm(3,:), pl_tm(1,:)));
+rl_tm(3,:) = 2*pi*cross(pl_tm(1,:), pl_tm(2,:))/dot(pl_tm(3,:), cross(pl_tm(1,:), pl_tm(2,:)));
+end
+

get_surface_symmetry_point.m

@@ -0,0 +1,14 @@
+function [ s ] = get_surface_symmetry_point( o, surface )
+%GET_SURFACE_SYMMETRY_POINT Summary of this function goes here
+%   Detailed explanation goes here
+
+a = surface(1,:);
+b = surface(2,:);
+v_ao = o-a;
+v_ab = a-b;
+l_ap = dot(v_ao,v_ab)/norm(v_ab);
+v_ap = l_ap*v_ab/norm(v_ab);
+v_as = v_ao-2*v_ap;
+s = v_as+a;
+end
+

get_surface_symmetry_points.m

@@ -0,0 +1,12 @@
+function [ sps ] = get_surface_symmetry_points( points, surface )
+%GET_SURFACE_SYMMETRY_POINTS Summary of this function goes here
+%   Detailed explanation goes here
+
+count = 0;
+sps = zeros(size(points));
+for o = points'
+    count  = count + 1;
+    sps(count,:) = get_surface_symmetry_point(o', surface);
+end
+end
+

get_symmetry_point.m

@@ -0,0 +1,15 @@
+function [ s ] = get_symmetry_point( o, line )
+%GET_SYMMETRY_POINT Summary of this function goes here
+%   Detailed explanation goes here
+
+ a = line(1,:);
+ b = line(2,:);
+ v_ab = b-a;
+ v_ao = o-a;
+ l_ap = dot(v_ao,v_ab)/norm(v_ab);
+ v_ap = l_ap*v_ab/norm(v_ab);
+ v_po = v_ao-v_ap;
+ v_as = v_ao-2*v_po;
+ s = v_as + a;
+end
+

get_symmetry_points.m

@@ -0,0 +1,12 @@
+function [ sps ] = get_symmetry_points( points, line )
+%GET_SYMMETRY_POINTS Summary of this function goes here
+%   Detailed explanation goes here
+
+count = 0;
+sps = zeros(size(points));
+for o = points'
+    count  = count + 1;
+    sps(count,:) = get_symmetry_point(o', line);
+end
+end
+

gif/bcc_pl.gif

@@ -0,0 +1,48 @@
+function [  ] = im2gif( varargin )
+%IM2GIF 将图片数据集保存为GIF动图
+%   IM2GIF(frames, filename, delay, size, n)
+%   frames  影片数据组，格式为 1xn，每一个元素为{cdata,colormap}结构体
+%   filename图片保存文件名
+%   delay   图片合成间隔，单位秒（s）
+%   size    图片尺寸size=[wight,height]
+%   n       图片保存的颜色位数
+vlen = length(varargin);
+if vlen < 2
+    error('number of arguments error !');
+else
+    frames = varargin{1,1};
+    filename = varargin{1,2};
+    % 默认参数配置
+    delay = .5;
+    size = false;
+    n = 32;
+    if vlen >2
+        delay = varargin{1,3};
+    end
+    if vlen > 3
+        size = varargin{1,4};
+    end
+    if vlen > 4
+        n = varargin{1,5};
+    end
+end
+
+flen = length(frames);
+if flen < 2
+    warning('No dynamic photo created: frames too less !');
+end
+for i=1:flen
+    im = frame2im(frames(1,i));
+    % 改变图片尺寸
+    if size~=false
+        im = imresize(im, fliplr(size));
+    end
+    % 转换数据格式
+    [I,map] = rgb2ind(im, n);
+    if i==1
+        imwrite(I,map,filename, 'gif', 'Loopcount', inf, 'DelayTime', delay);
+    else
+        imwrite(I,map,filename, 'gif', 'WriteMode', 'append', 'DelayTime', delay);
+    end
+end
+end