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即Λ₀）
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       % 设置 H、M 为全局变量
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)         % 求斥力常数Λ₀
vxa(x,a)          % 求宇宙常数Λ
p(a,M)           % 求宇宙密度ρ
_____________________________________________________________________________