Matlab计算最外层(球面)

记号:x 代表宇宙膨胀系数μ
____________________________________________________________________________
函数文件1 aqm.m   (最外层的半径a , 球面qm
function a=aqm(x,H,M)
global H M
a=27.8592*(M*(1/x^3-3/x+2)/H^2)^(1/3)
____________________________________________________________________________
函数文件2    (文件名为zvx1.m ,是速度红移的一次项)
function zv1=zvx1(x,H,M)
global H M
zv1=0.092928258*(H*M*(1/x^3-3/x+2))^(1/3)
____________________________________________________________________________
函数文件3 zvx2.m   (文件名为zvx2.m ,是速度红移的二次项,即相对论修正项)
function zv2=zvx2(x,H,M)
global H M
zv2=0.5*
0.092928258*(H*M*(1/x^3-3/x+2))^(1/3)^2
____________________________________________________________________________

函数文件4 zvx12   (文件名为zvx2.m ,是速度红移的一次项与二次项之和)
function zv12=zvx12(x,H,M)
global H M
zv12=zvx1(x)+zvx2(x)
____________________________________________________________________________
函数文件5 zex.m   (斥力红移)
function ze=zex(x,H,M)
global H M
k=0.008635661;
ze=0.5*k*(H*M)^(2/3)*(1/x^3-1)/((1/x^3-3/x+2)^(1/3)-k*(H*M)^(2/3)*(0.5/x^3+1))
____________________________________________________________________________
函数文件6   zvex.m   (速度红移与斥力红移之和)
function zve=zvex(x,H,M)
global H M
zve=zvx12(x)+zex(x)
____________________________________________________________________________
函数文件7 zqm.m   (球面红移的非线性方程,即:  速度红移 + 斥力红移 - 2.3 = 0
function x=zqm(u,H,M)
global H M
x=((0.092828258*(H*M*(1/u^3-3/u+2))^(1/3))+0.5*(0.092828258*
(H*M*(1/u^3-3/u+2))^(1/3))^2+(0.5*0.008636*(H*M)^(2/3)*(1/u^3-1))/
(1/u^3-3/u+2)^(1/3)/(1-0.008636*(H*M)^(2/3)*(0.5/u^3+1)/(1/u^3-3/u+2)^(1/3)))-2.30
____________________________________________________________________________
函数文件8 p.m   (用字母p代表宇宙密度ρ
function p=p(a,M)
global M
p=81.2568*M/a^3
____________________________________________________________________________
函数文件9 v0xa.m   (用字母V代表希腊字母Λv0Λ0
function v0=v0xa(x,a)
v0=22.70933/x^3/a^3
______________________________________________________________________________
函数文件10 vxa   (用字母V代表宇宙常数Λ
function V=vxa(x,a)
V=22.7093347*(1/x^3-1)/a^3
______________________________________________________________________________
M_文件:  qm.m
global H M     
 % 设置 HM 为全局变量
H=input('H=');    % 键盘输入H
M=input('M=');
u=input('u=');    
% u=0.4 μ的近似值,这里用u代替μ
x=fzero('zqm',u)   % 解非线性方程,求宇宙膨胀系数μ
x=input('x=');      % 键盘输入x
a=aqm(x)        % 求宇宙半径
zvex(x)          % 求速度红移和斥力红移之和
v0xa(x,a)         % 求斥力常数Λ0
vxa(x,a)          % 求宇宙常数Λ
p(a,M)           % 求宇宙密度ρ
_____________________________________________________________________________