参考文献“New views of solar wind with the Lambert W function" by Steven R.Cranmer,精读SectionIV.THE parker solar wind promblem部分，利用Lambert函数来求解太阳风等温模型Equations 13-16，以图表显示太阳风解,由于采用了近似函数，所以部分解计算值为负，无法显示。

源代码：

clear;

K=1.38*10^(-23)

T=[0.1*10^6,0.2*10^6,0.5*10^6,1*10^6,2*10^6,5*10^6]

m=1.6726231*10^(-27)

RT=K*T/m

Vc=sqrt(RT)

a=Vc

G=6.67*10^(-11)

M=1.989*10^30

sunR=6.96*10^8

e=exp(1)

for m=1:1:6

   rc(m)=G*M/(2*Vc(m)^2)

   for i=1:1:200

       r=sunR*i

       D(m)=(r/rc(m))^(-4)*exp(4*(1-rc(m)/r)-1)

       if r<=rc

          u(m,i)=sqrt(-a(m)^2*(-1+sqrt(2-2*e*D(m))))/1000

       else

          u(m,i)=sqrt(-a(m)^2*(-1-sqrt(2-2*e*D(m))))/1000

       end

    end

end

figure(1);

l=1:1:200;

%画图

plot(l,u(1,:),'g',l,u(2,:),'b',l,u(3,:),'k',l,u(4,:),'r',l,u(5,:),'y',l,u(6,:),'c');

set(gca,'xscale','log');

set(gca,'yscale','log');

legend('T=0.1','T=0.2','T=0.5','T=1','T=2','T=5');

grid;

title('parker太阳风解');

xlabel('r/R');ylabel('u');

hold on;

%标记临界点

plot(rc(1)/sunR,Vc(1)/1000,'go',rc(2)/sunR,Vc(2)/1000,'bo',rc(3)/sunR,Vc(3)/1000,'ko',rc(4)/sunR,Vc(4)/1000,'ro',rc(5)/sunR,Vc(5)/1000,'yo',rc(6)/sunR,Vc(6)/1000,'co');

hold off;