%?%?Create?matrix?will?all?fluctuating?velocity?components?for?each?snapshot?in?a?column
%?%
%?%?Do?POD?analysis
%?n?速度場(chǎng)的總的節(jié)點(diǎn)數(shù)
%?m?時(shí)間文件的個(gè)數(shù)，生成矩陣的維數(shù)
%?Uall?將所有文件的速度U和速度Ｖ重構(gòu)得到的POD分析的矩陣
%?R?為二階相關(guān)函數(shù)構(gòu)成的矩陣
%?eValue?為矩陣R特征值
%?eVec?為矩陣R特征向量
%?menergy為各特征值所對(duì)應(yīng)的能量百分比
%?phi為得到的基函數(shù)
clear;clc;

kx=92;
ky=99;
kz=22;
lx1=11;
lx2=82;
lz1=4;
lz2=19;
lx=lx2-lx1+1;
lz=lz2-lz1+1;
n=lx*ky;
dr=0.6872;
x=-35.5357:dr:-35.5357+dr*（kx-1）;
y=-31.3456:dr:-31.3456+dr*（ky-1）;
z=0.6096:dr:0.6096+dr*（kz-1）;
dt=0.002;
t0=fix（lx*0.75）;
lt=200;
%?file1=t0+1;
%?file2=t0+1+lt;
file1=1;
file2=t0+lt+t0;
lf=file2-file1+1;

m=importdata（‘mean?velocity_flipxy.dat‘）;
m1=m.data（:?4）;
m2=reshape（m1?kx?ky?kz）;
m3=m2（lx1:lx2?:?:）;

U=zeros（?n?lf）;
for?iz=4
????mu?=?reshape（m3（:?:?iz）?[]?1）;
????for?i=file1:file2
????????fname??=?sprintf（‘d%04d.dat‘?i）;
????????a??????=?importdata（fname）;
????????temp???=?a.data（:4）;
????????temp???=?reshape（temp?kx?ky?kz）;
????????iu?????=?-temp（:?:?iz）;
????????iu?????=?flipud（iu）;
????????iu?????=?fliplr（iu）;
????????iu?????=?iu（lx1:lx2?:）;
%?????????iu?????=?iu‘;
????????iu?????=?reshape（iu?[]?1）;
????????
????????U（:?i-file1+1）?=?iu-mu;
????????disp（fname）;
????end
????
????R=U‘*U;??%?solve:?eV?is?eigenvectors?D?is?eigenvalues?in?diagonal?matrix
????[eVD]=eig（R）;?%?sort?eigenvalues?in?ascending?order-I?is?sorted?index?vector
????disp（‘eigenvalues?have?been?done.‘）;
????[LI]=sort（diag（D））;
????for?i=1:length（D）
????????eValue（length（D）+1-i）=L（i）;???%?Eigenvalues?sorted?in?descending?order
????????eVec（:length（D）+1-i）=eV（:I（i））;??%?Eigenvectors?sorted?in?the?same?order
????end;
????for??k=1:length（eValue）
????????if?eValue（k）<=0
????????????eValue（k）=0;??????????????%?last?eigenvalue?should?be?zero
????????end
????end
????menergy=eValue/sum（eValue）;???????%?relative?energy?associated?with?mode?m
????
????%計(jì)算能量百分比
????for?j=1:lf
????????cumulant（j）=sum（menergy（1:j））;
????end
????fname=sprintf（‘pod_energy_z%02d.dat‘?iz）;
????fenergy=fopen（fname?‘w‘）;
????fprintf（fenergy‘VARIABLES?=?“mode“?“Percent“?“Cumulant?Percent?“\r\n‘）;
????for?i=1:length（eValue）
????????fprintf（fenergy‘%d\t%g\t%g\r\n‘?i?menergy（i）?cumulant（i））;
????end
????fclose（fenergy）;
????disp（‘energy?percent?have?been?done.‘）;
????
????%計(jì)算歸一化的基函數(shù)
????for?i=1:lf
????????tmp=U*eVec（:i）;
????????phi（:i）=tmp/norm（tmp）;
????end;
????disp（‘normalizations?have?been?done.‘）;
????
????%?輸出基函數(shù)的等值線圖
????fname=sprintf（‘pod_phi_z%02d.dat‘?iz）;
????fphi=fopen（fname‘w‘）;
????fprintf（fphi‘VARIABLES?=?“x“?“y“?“phi?u“\r\n‘）;
????for?k=1:50
????????fprintf（fphi‘\r\nZone?T=“mode?%d“i=%d?j=%d\r\n‘?k?lx?ky）;
????????for?iy=1:ky
????????????for?ix=1:lx
????????????????fprintf（fp