Ps)+0.0008sin(2S-3h+Ps)+0.0115sin(2S-2h)+0.0037sin(2S- 3.2 自定义函数部分代码

2P)-0.002sin(2S-2N)+0.0009sin(3S-2h-P) (式7)

function deltag=gravity_earthtide(a,L,phi)

β=-0.0048sin(P-N)-0.0008sin(2h-P-N)+0.003sin(S-%输入参量a为MATLAB时间，L为台站经度，东经为正，m phi为台站纬度，北纬为正，输出变量为修正值，单位μ2h+N)+0.0895sin(S-N)+0.001sin(2S-2h+P-N)+0.0049sin(2S-Gal。

P-N)+0.0006sin(3S-2h-N) (式8)

timevec=datevec(a)；%转化成向量形式，以便读取参数Sinδ=sinεsinλcosβ+cosεsinβ (式9)

m m m y=timevec(1)；CosδcosH=cosβcosλcosθ+sinθ(cosεcosβ

m m m=timevec(2)；sinλ-sinεsinβ) (式10)

m

m m d=timevec(3)；

CosZ =sinφsinδ+cosφcosδcosH (式11)m t=timevec(4)； %y，m，d分别为年、月、日，t为北

θ=(t-8)×15°+h+L-180° (式12)

京时间的小时数

式中λ为月亮的黄经，β为月亮的黄纬，δ为月亮的

m m …

赤纬，H为月亮的时角，θ为地方恒星时，φ为台站地心纬… %中间部分为理论公式，与计算步骤相同度，L为台站的经度，东经为正，t为北京时间的小时数。Deltag=deltag*1.16 %最后结果乘以潮汐因子

计算 、cosZs 4 计算实例 =1+0.0168cos(h-Ps)+0.0003cos(2h-2Ps) (式13)以成都(L=104.0°，φ=30.6°)2008年1月1日为λ=h+0.0335sin(h-Ps)+0.0004sin(2h-2Ps) s 例，计算结果如下：β=0 (式14)

s CosZ =sinφsinεsinλ+cosφ(cosλcosθ+sinθs s s cosεsinλ) (式15)

s 3 MATLAB程序设计3.1 m文件

编制m文件earthtide_day1_day2，实现两个时间段之

表1 成都2008年1月1日重力固体潮理论值

间整点时刻的修正值的自动计算，代码如下：

以下为与已知理论值的对比。

day1=input('输入计算开始时间：')；%输入时间为时间向量，一行六列

day2=input('输入计算截止时间：')；L=input('输入经度：')；phi=input('输入纬度：')；

Day1num=datenum(day1)；%将时间向量转化为数字表示的MATLAB时间

day2num=datenum(day2)；x1=[2008 8 8 20 0 0]；x2=[2008 8 8 21 0 0]；

delta=datenum(x2)-datenum(x1)； noh=round((day2num-day1num)/delta)；

daynum=linspace(day1num,day2num,noh+1)；%将两时刻之间按一小时间隔段均分

daynum=daynum'；n=length(daynum)；表2 2008.11.01～2008.11.10成都地区重力固体潮理论值对比amend=[]；可以看出，计算结果与已知理论值在时间和数值上均for i=1:n 较吻合，证明了计算结果的可靠性。amend=[amend;gravity_earthtide(daynum(i),L,phi)] 5 结束语；%amend为修正值，MATLAB的语言结构决定了它在科学计算方面的优势，

end % gravity_earthtide为自定义函数相比其他语言，程序易于编写且更加简洁精练，用于重力dayvec=datevec(daynum)；固体潮理论值的计算可以实现一月或一年的输出，并可实gravitytide=[dayvec amend]；现结果的可视化，数据以Excel文件格式成表备用，方便查name1=input('输入保存数据的xls表格名称:')；找和计算。Xlswrite(name1,gravitytide)；%输出形式为Excel表格参考文献：(略)plot(daynum,amend) %绘制理论曲线图

中国西部科技 2009年1月(上旬)第08卷第01期第162期

总34