PCOX-自定义函数wfnread分析
发完上期的推文，过冷水中觉得“庖丁解牛”的典故用好牵强，明明“愚公移山”更加符合我们的故事，具体内容见愚公移山。
本期我们来看一看wfnread函数
[wfn,pot,rho,wavedata,XYZ]=wfnread(fileName, h)还记得我们是在讲readoutput的过程中插入的wfnread。为何要在
[data, wfn, D_valid, W_valid] = readOutput(outName, datName)中自定义一个函数来处理datName文件而不是用处理outName文件的方法？这是因为两个文件类型不一样，读取数据的方法差别很大，如果都放在一起，容易造成混乱。我们来看代码：
myFile = fopen('*.dat','r', 'l');
prefix = fread(myFile, 1, 'uint32');
if (prefix > 100)
myFile = fopen('*.dat','r','b');
prefix = fread(myFile, 1, 'uint32');
End这里有两个函数fopen、fread读者有可能不太熟悉。
fileID = fopen(filename,permission,machinefmt,encodingIn)：
A = fread(fileID,sizeA,precision,skip,machinefmt)Fopen：打开一个二进制文件，filename：文件名；permission：表示文件操作模式，是只读呢？可写呢？等等，machinefmt：二进制字节顺序；encodingIn：字节编码模式
Fread：读取二进制文件，fileID：打开的变量名；sizeA：读取内容的规模；skp:表示跳过几个元素。A:表示读取的具体内容。
两个函数的参数设置很复杂，一两句讲不清。而过冷水又不是那种解牛而不全的人，预知详情当然得且看下回分解。这里我们只需说清当下参数的含义即可。
大小端判断完毕就开始文件的读写，我们要读写什么内容呢？
prefix = fread(myFile, 1, 'uint32');
date = fread(myFile, prefix, 'uchar');
date = (char(date'))
date =
[11-JUL-2005 16:30:28 -0500]它是文件的生成日期，并没有什么*用。
prefix = fread(myFile, 1, 'uint32');
prefix = fread(myFile, 1, 'uint32');
wfn.step = fread(myFile, 1, 'double');
wfn.radius = fread(myFile, 1, 'double');
wfn.ndim = fread(myFile, 1, 'uint32');
wfn.num_wedges = fread(myFile, 1, 'uint32');
wfn.trans=fread(myFile, 1, 'uint32')%这个数据是我手动添加的；
wfn.spin=fread(myFile, 1, 'uint32')%这个数据也是手动添加；
fread(myFile, 2, 'uint32'); % i_1's含义是不是读取两行数据？
prefix = fread(myFile, 1, 'uint32');
prefix = fread(myFile, 1, 'uint32');
wfn.shift = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');这段代码和之前的readoutput函数中的：
inMain = load('pTemp/main.dat');
data.timeDiag = inMain(:,1);
data.FermiLevel = inMain(:,2);
data.timeHart = inMain(:,3);有相似性，需要注意的点有:
prefix = fread(myFile, 1, 'uint32');这段代码的的生成值prefix并没有在后续代码调试中使用，是不是有点智障设置？你猜这段代码想干什么？过冷水起初以为是很重要的一个量，后来明这段代码只是起到调节指针的作用。二进制文件写入也不完全写入都是我们想要的数据，空格或者有些数据是我们不需要写入的，就用该方法跳过。这段代码用的就很技巧性，同时也对生成文件的格式有严格要求，文件格式不一致，多一个空行或者多一个数据，会引起很大的麻烦。
prefix = fread(myFile, 1, 'uint32');
fread(myFile, prefix/4, 'uint32'); % i_1's
prefix = fread(myFile, 1, 'uint32');
prefix = fread(myFile, 1, 'uint32');
fread(myFile, prefix/8, 'double'); % one's
prefix = fread(myFile, 1, 'uint32');
%read in (one,kk=1,3)
prefix = fread(myFile, 1, 'uint32');
fread(myFile, prefix/8, 'double'); % one's
prefix = fread(myFile, 1, 'uint32');
%read in (zero,kk=1,9*3)
prefix = fread(myFile, 1, 'uint32');
fread(myFile, prefix/8, 'double'); % zero's
prefix = fread(myFile, 1, 'uint32');
%read in (i_1)
prefix = fread(myFile, 1, 'uint32');
fread(myFile, prefix/4, 'uint32'); % i_1
prefix = fread(myFile, 1, 'uint32');
%read in (wfn%kx(kk), wfn%ky(kk), wfn%kz(kk), kk=1, wfn%nwedge)
prefix = fread(myFile, 1, 'uint32');
if (prefix/4 ~= wfn.num_wedges*3);
out = 'ERROR: reading wfn.k* from wfn.dat'
end这么一大段代码都是不断在重复读文件调节指针，最后做了判断。What you want do？
一顿神乎其技的操作，让人不知道文件指针跑到什么地方去了，为了确定文件指针在我们想要的位置，做了一个判断指针位置的操作。正常人会这么写代码吗？为什么不锁定关键语句？写到这才发现这段代码写的这么怪能够完成特定任务，那是因为代码参照*.dat的生成代码来写的。生成代码怎么将数据转换成二进制文件，那matlab就按照格式来转化成可看数据。不了解数据生成代码前提下，自然好多操作有点费解。
XYZ = fread(myFile, [3,wfn.num_wedges], 'uint32');
XYZ = transpose(XYZ);
prefix = fread(myFile, 1, 'uint32');这代代码的含义是读取三维坐标
prefix = fread(myFile, 1, 'uint32');
wfn.num_states = fread(myFile, prefix/4, 'uint32');
prefix = fread(myFile, 1, 'uint32');
wfn.first_state = 1;
if (wfn.num_states > 35)
global tmp_states first_tmp;
tmp_states = wfn.num_states;
uiwait(wfn_warning);
wfn.first_state = first_tmp;
wfn.num_states = tmp_states;
clear global tmp_states first_tmp
end这段代码是对状态数进行判断如果wfn.num_states > 35会进行重新赋值，这段赋值操作用到自定义函数 uiwait(wfn_warning)，因为当前不影响我们的运行，所以就不讲该自定义函数，过冷水担心不断深入展开，到时读者都不知道转到什么地方去了.
prefix = fread(myFile, 1, 'uint32');
out = fread(myFile, prefix/4, 'uint32'); % i_1's
prefix = fread(myFile, 1, 'uint32');
prefix = fread(myFile, 1, 'uint32');
wfn.en_tmp = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');
%read in (occ_tmp(kk),kk=1,nstate)
prefix = fread(myFile, 1, 'uint32');
wfn.occ_tmp = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');
%read in (pot%vnew(kk,1),kk=1,grid%nwedge)
prefix = fread(myFile, 1, 'uint32');
pot = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');
%read in (elec_st%rho(kk,1),kk=1,grid%nwedge)
prefix = fread(myFile, 1, 'uint32');
rho = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');这段代码还是读取数据。我们有新的数据pot、rho涉及到的内容和前段内容一致，就不在此讲解了。记性好的读者应该还能记得我们这段自定义函数的目的是要获取：wfn, pot, rho, wavedata, XYZ，现在唯独没有 wavedata,数据，所以现在我们开始完成最后一步提取 wavedata:
wavedata = zeros(wfn.num_wedges, wfn.num_states);
for (state = 1:wfn.num_states);
prefix = fread(myFile, 1, 'uint32');
wavedata(:,state) = fread(myFile, prefix/8, 'double');
prefix = fread(myFile, 1, 'uint32');
end先自定义wavedata变量，它的格式根据前面的参数确定然后就是读取每一个state下的数据。
这样我们这段自定义函数的功能就完成了。
腾讯云开发者
