Wrote matlab code to solve Parker Wind

Here is a matlab routine to solve for the parker wind given lambda x is in units of planet radius, and y is velocity in units of sound speed.

Also notified Titos about typo in arxiv preprint that made initial attempts to solve the equation nonsensical.

Also note that the original version of this post assumed that

when in fact

So here is the correct version of the matlab routine 

lambda=11.5;
N=81;
decs=4;
xx=logspace(0,decs,N);
yy=zeros(N,1);
syms y positive;
for i=1:N
  yy(i)=vpasolve(y - log(y) == -3 - 4*log(lambda/2)+4*log(xx(i))+2*lambda/xx(i),y,(xx(i) > lambda/2) * 100+.001);
end
fid=fopen(['parker_wind_',num2str(lambda),'.dat'],'w')
fprintf(fid, '%i\n',N);
fprintf(fid,'radius, density, velocity\n')
for i=1:N
    fprintf(fid,'%e %e %e\n',xx(i),1d0/xx(i)^2/sqrt(yy(i)/yy(1)), sqrt(yy(i)));
end
fclose(fid);

Attachments (1)

Download all attachments as: .zip

Comments

No comments.