核数据检索——中子宏观反应截面

MT项选择很多，有专门的列表，是编号。你可以参考ENDF数据库的“能谱绘图”，点开后MT选项后面有个help下拉菜单，里面是详细的列表，根据需要选吧。数据格式如下所示，我示意的是Th-232，MT取1，N-TOT

CENDL21

ENDFB6

9.023200+4 2.300400+2 0 0 0 09040 3 1

0.000000+0 0.000000+0 0 0 24379040 3 1

208 5437 2 9040 3 1

1.000000-5 3.896260+2 1.000000-4 1.320680+2 1.200000-4 1.216890+29040 3 1

1.400000-4 1.136240+2 1.600000-4 1.071180+2 1.800000-4 1.017310+29040 3 1

2.000000-4 9.717550+1 2.200000-4 9.325540+1 2.400000-4 8.983570+19040 3 1

2.530000-4 8.783380+1 2.600000-4 8.681930+1 2.800000-4 8.413140+19040 3 1

=====================================================================

JEF22

9.023200+4 2.300400+2 0 0 0 09040 3 1

0.000000+0 0.000000+0 0 0 23859040 3 1

277 5385 2 0 09040 3 1

1.000000-5 3.875200+2 1.000000-4 1.306300+2 1.200000-4 1.202800+29040 3 1

1.400000-4 1.122300+2 1.600000-4 1.057500+2 1.800000-4 1.003700+29040 3 1

2.000000-4 9.583000+1 2.200000-4 9.192000+1 2.400000-4 8.851000+19040 3 1

======================================================================

JENDL

9.02320+ 4 2.30045+ 2 0 0 0 09040 3 1

0.00000+ 0 0.00000+ 0 0 0 31519040 3 1

3 5 17 2151 59040 3 1

1.00000- 5 4.72560+ 0 2.53000- 2 9.39500- 2 3.50000+ 3 2.52594- 49040 3 1

3.50000+ 3-1.70000- 1 4.00000+ 3-2.10000- 1 5.00000+ 3-2.50000- 19040 3 1

首先，需要说明两点。

第一，本案例使用的是单月潮汐观测数据，处理方法则是基于长期观测资料的调和分析来进行处理。中期观测资料的分析需要分别求主要分潮、随从分潮，短期观测资料分析则还需要计算不同观测序列的权重，但是核心算法与长期观测资料分析是一致的，都是建立矛盾方程，然后使用最小二乘法建立法方程，求出法方程系数，再求出矩阵X、Y。

第二，本文主要介绍方法步骤，所用代码大多为关键步骤实现，仅供参考。如需完整代码，请关注博主的另一篇资源。

下面开始介绍本案例的处理，从理论上讲，首先，我们需要选取分潮，确立我们所要分析的天文分潮，本案例用8个主要分潮——M2、S2、N2、K2、K1、O1、P1、Q1，四个半日潮，四个全日潮（其实去看潮汐相关研究的文献就会发现，基本都是以这八个分潮为主的）。相邻数据时间间隔为一小时，以所有数据的中间数对应时刻作为时间原点，然后对观测记录数据进行排序。

从实际出发，在matlab中，首先清理空间、准备环境（这是一个良好习惯），然后需要导入数据，对分潮进行排序。

%%导入数据

clear allclose allclc

data=importdata('C:\Users\STAR\Desktop\TideData_01.txt')

%%分潮排序（八个）

M2=1

S2=2

N2=3

K2=4

K1=5

O1=6

P1=7

Q1=8

登录后复制

常量的准备，输入杜德森数，以及天文元素随时间的变化速度。

%%杜德森数（八个分潮各有七个杜德森数）

miu{M2}=[2,0,0,0,0,0,0]

miu{S2}=[2,2,-2,0,0,0,0]

miu{N2}=[2,-1,0,1,0,0,0]

miu{K2}=[2,2,0,0,0,0,0]

miu{K1}=[1,1,0,0,0,0,1]

miu{O1}=[1,-1,0,0,0,0,-1]

miu{P1}=[1,1,-2,0,0,0,-1]

miu{Q1}=[1,-2,0,1,0,0,-1]

%%天文元素随时间的变化速度

rateOfChange=[14.49205211,0.54901653,0.04106864,0.00464183,0.00220641,0.00000196]

登录后复制

然后需要计算时间原点。个人习惯将代码模块化，本处使用自编函数TimeCalculation。

year=2003

month=3

day=1

hour=0

[year,month,day,hour]=TimeCalculation(2003,3,1,0,360)

登录后复制

计算分潮角速度

for i=M2:Q1

sigma(i)=AngularVelocity(miu{i},rateOfChange)

end

sigma=deg2rad(sigma) %注意，这里涉及到一个角度转弧度的 *** 作

登录后复制

然后是计算各天文元素

[tao,ss,hhh,pp,NNN,ppp]=AstronomicalElements(year,month,day,hour)

tao=deg2rad(tao)

ss=deg2rad(ss)

hhh=deg2rad(hhh)

pp=deg2rad(pp)

NNN=deg2rad(NNN)

ppp=deg2rad(ppp)

astronomicalElements=[tao,ss,hhh,pp,NNN,ppp]

登录后复制

计算分潮初始位相，同样使用自编函数

for i=M2:Q1

v0(i)=InitialPhase(miu{i},astronomicalElements)

end

v0(2) = 6.2832%由于S2分潮的特殊性，直接赋值

登录后复制

计算交点因子及订正角

deltaMiu4{M2}=[0,0,0,2,2]

deltaMiu5{M2}=[-2,-1,0,0,1]

rho{M2}=[0.0005,-0.0373,1,0.0006,0.0002]

deltaMiu4{K2}=[0,0,0,0]

deltaMiu5{K2}=[-1,0,1,2]

rho{K2}=[-0.0128,1,0.2980,0.0324]

deltaMiu4{K1}=[-2,0,0,0,0,0]

deltaMiu5{K1}=[-1,-2,-1,0,1,2]

rho{K1}=[0.0002,0.0001,-0.0198,1,0.1356,-0.0029]

deltaMiu4{P1}=[0,0,0,2,2]

deltaMiu5{P1}=[-2,-1,0,0,1]

rho{P1}=[0.0008,-0.0112,1,-0.0015,-0.0003]

deltaMiu4{O1}=[0,0,0,2,2,2]

deltaMiu5{O1}=[-2,-1,0,-1,0,1]

rho{O1}=[-0.0058,0.1885,1,0.0002,-0.0064,-0.0010]

for i=[M2,K2,K1,P1,O1]

[f(i),u(i)]=IntersectionFactorAndCorrectionAngle(deltaMiu4{i},deltaMiu5{i},rho{i},pp,NNN)

end

f(Q1)=f(O1)

u(Q1)=u(O1)

f(N2)=f(M2)

u(N2)=u(M2)

f(S2)=1

u(S2)=0

登录后复制

建立法方程并求x和y

deltaT=1

N=data(end,1)

NN=(N-1)/2

A(0+1,0+1)=N

for i=M2:Q1

A(0+1,i+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2))

A(i+1,0+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2))

A(i+1,i+1)=(N+(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2

B(i,i)=(N-(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2

end

for i=M2:Q1

for j=M2:Q1

if ~(i==j)

A(i+1,j+1)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))+(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2

B(i,j)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))-(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2

end

end

end

F1(0+1)=sum(data(:,2))

for i=M2:Q1

F1(i+1)=sum(data(:,2).*cos((data(:,1)-361)*sigma(i)*deltaT))

F2(i)=sum(data(:,2).*sin((data(:,1)-361)*sigma(i)*deltaT))

end

X=F1/A(:,1:end)

Y=F2/B(:,1:end)

登录后复制

计算准调和振幅R和位相theta

R=(X(2:end).^2+Y.^2).^0.5

for i=M2:Q1

theta(i)=CalculatedPhase(R(i),X(i+1),Y(i))

end

登录后复制

计算分潮的调和常数

H=R./f

g=theta+v0+u

for i=M2:Q1

g(i)=rad2pi(g(i))

end

登录后复制

潮汐预报、计算自报余差

S0=X(0+1)

for i=1:721

h(i)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g))

end

r=data(1:end,2)'-h(1:end)

delta=sum(r.^2)^0.5/data(end,1)

登录后复制

计算预报潮位

starttime=(31-1+30+1)*24+1

endtime=starttime+31*24+1

forecastTime=1:31*24+1+1

for i=starttime:endtime

forecastTide(i-starttime+1)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g))

end

登录后复制

绘制预报图

figure(2)

plot(forecastTime,forecastTide)

title("Forecast(May 1 to June 1)")

xlabel("Serial number")

ylabel("Stage")

legend("Forecast",'Location','Best')

登录后复制

写在文末，潮汐体现的是海洋动力及水文要素的变化规律和控制机制，从古至今，关于潮汐的研究从未停止，远有沈括著作《梦溪笔谈》，近有牛顿平衡潮理论之广泛应用，希望本文能为有需要的人提供一点思路及方法。

---