惑星や星の回転力と重力をシミュレートできる「太陽系」シミュレーターに興味があります。
私たちの太陽系をシミュレートし、さまざまな速度でシミュレートできるようにしたいと思います (つまり、地球や他の惑星が太陽の周りを日、年などで回転するのを観察します)。惑星を追加したり、惑星の質量を変更したりして、システムにどのように影響するかを確認したいと思います.
この種のシミュレーターを作成するための正しい方向に私を向けるリソースを誰かが持っていますか?
この目的のために設計された既存の物理エンジンはありますか?
それはここにあるすべてであり、一般的に、Jean Meeus が書いたものすべてです。
ニュートンの万有引力の法則とケプラーの惑星運動の法則を知り、理解する必要があります。この 2 つは単純で、高校で勉強したことがなければ、聞いたことがあると思います。最後に、シミュレーターを可能な限り正確にしたい場合は、 n-Body 問題に慣れておく必要があります。
シンプルに始めるべきです。オブジェクトとその周りを回転Sunするオブジェクトを作成してみてください。Earthこれにより、非常に堅実なスタートが得られ、そこから拡張するのはかなり簡単です。惑星オブジェクトは次のようになります。
Class Planet {
float x;
float y;
float z; // If you want to work in 3D
double velocity;
int mass;
}
それを覚えておいてF = MA、残りはただの退屈な数学です:P
これは、一般的な N 体の問題に関する優れたチュートリアルです。
http://www.artcompsci.org/#msa
Ruby を使用して作成されていますが、他の言語などにマッピングするのは非常に簡単です。一般的な統合アプローチのいくつかをカバーしています。Forward-Euler、Leapfrog、Hermite。
実装する必要があるのは、適切な微分方程式 (ケプラーの法則) とルンゲクッタを使用することだけです。（少なくともこれは私にとってはうまくいきましたが、おそらくもっと良い方法があります）
そのようなシミュレーターはオンラインでたくさんあります。
これは、500 行の C コードで実装された単純なものの 1 つです。(モーションアルゴリズムははるかに少ない) http://astro.berkeley.edu/~dperley/programs/ssms.html .
これも確認してください:
http://en.wikipedia.org/wiki/Kepler_problem
http://en.wikipedia.org/wiki/Two-body_problem
http://en.wikipedia.org/wiki/N-body_problem
物理学では、これはN 体問題として知られています。惑星が 3 つ以上ある系では、これを手で解くことができないことで有名です。幸いなことに、コンピューターを使って非常に簡単に近似解を求めることができます。
このコードをゼロから作成するための優れた論文は、ここにあります。
ただし、ここで警告の言葉が重要だと思います。期待した結果が得られない場合があります。方法を確認したい場合:
問題はこれです。
ええ、現代の天文学者は、土星の質量が太陽の周りの地球の軌道をどのように変化させるかに関心を持っています. しかし、これは非常にマイナーな効果です。太陽の周りの惑星の経路をプロットする場合、太陽系に他の惑星があるかどうかはほとんど問題になりません。太陽はとても大きいので、他のすべての重力をかき消してしまいます。これに対する唯一の例外は次のとおりです。
明確にするために、はい、惑星間のいくつかの相互作用を計算できます。しかし、現実的な太陽系を作成する場合、これらの相互作用は肉眼では重要ではありません。
でも試してみて、見つけてください！
C++ で記述され、OpenGL を使用する n-body モデリング ツールキットであるnModを確認してください。かなり人口の多い太陽系モデルが付属しており、簡単に変更できるはずです。また、彼は n-body シミュレーション全般に関する非常に優れた wiki を持っています。これを作成したのと同じ人物がMoodyという新しいプログラムも作成していますが、それほど進んでいないようです。
さらに、数個以上のオブジェクトで n 体シミュレーションを行う場合は、高速多極子法(高速多極子アルゴリズムとも呼ばれます) を検討する必要があります。計算の数を O(N^2) から O(N) に減らして、シミュレーションを本当に高速化できます。この記事の著者によると、これは 20 世紀で最も成功したアルゴリズムのトップ 10の 1 つでもあります。
惑星物理学をシミュレートするアルゴリズム。
これは、私の Android アプリでの Keppler パーツの実装です。ソース全体をダウンロードできるように、主要な部分は私の Web サイトにあります: http://www.barrythomas.co.uk/keppler.html
これは、軌道の「次の」位置に惑星を描くための私の方法です。追跡しようとしている惑星と同じ周期を持つ円で、一度に 1 度ずつ、円を 1 周するようなステップを考えてみてください。このメソッドの外では、回転の度数を含む dTime と呼ばれるグローバル double をステップ カウンターとして使用します。
メソッドに渡される主要なパラメーターは、dEccentricty、dScalar (軌道がすべてディスプレイに収まるようにするスケーリング係数)、dYear (地球年での軌道の持続時間)、および近日点が適切な場所になるように軌道を方向付けるためのものです。文字盤には、いわばdLongPeri - 近日点の経度。
ドロープラネット:
public void drawPlanet (double dEccentricity, double dScalar, double dYear, Canvas canvas, Paint paint,
String sName, Bitmap bmp, double dLongPeri)
{
double dE, dr, dv, dSatX, dSatY, dSatXCorrected, dSatYCorrected;
float fX, fY;
int iSunXOffset = getWidth() / 2;
int iSunYOffset = getHeight() / 2;
// get the value of E from the angle travelled in this 'tick'
dE = getE (dTime * (1 / dYear), dEccentricity);
// get r: the length of 'radius' vector
dr = getRfromE (dE, dEccentricity, dScalar);
// calculate v - the true anomaly
dv = 2 * Math.atan (
Math.sqrt((1 + dEccentricity) / (1 - dEccentricity))
*
Math.tan(dE / 2)
);
// get X and Y coords based on the origin
dSatX = dr / Math.sin(Math.PI / 2) * Math.sin(dv);
dSatY = Math.sin((Math.PI / 2) - dv) * (dSatX / Math.sin(dv));
// now correct for Longitude of Perihelion for this planet
dSatXCorrected = dSatX * (float)Math.cos (Math.toRadians(dLongPeri)) -
dSatY * (float)Math.sin(Math.toRadians(dLongPeri));
dSatYCorrected = dSatX * (float)Math.sin (Math.toRadians(dLongPeri)) +
dSatY * (float)Math.cos(Math.toRadians(dLongPeri));
// offset the origin to nearer the centre of the display
fX = (float)dSatXCorrected + (float)iSunXOffset;
fY = (float)dSatYCorrected + (float)iSunYOffset;
if (bDrawOrbits)
{
// draw the path of the orbit travelled
paint.setColor(Color.WHITE);
paint.setStyle(Paint.Style.STROKE);
paint.setAntiAlias(true);
// get the size of the rect which encloses the elliptical orbit
dE = getE (0.0, dEccentricity);
dr = getRfromE (dE, dEccentricity, dScalar);
rectOval.bottom = (float)dr;
dE = getE (180.0, dEccentricity);
dr = getRfromE (dE, dEccentricity, dScalar);
rectOval.top = (float)(0 - dr);
// calculate minor axis from major axis and eccentricity
// http://www.1728.org/ellipse.htm
double dMajor = rectOval.bottom - rectOval.top;
double dMinor = Math.sqrt(1 - (dEccentricity * dEccentricity)) * dMajor;
rectOval.left = 0 - (float)(dMinor / 2);
rectOval.right = (float)(dMinor / 2);
rectOval.left += (float)iSunXOffset;
rectOval.right += (float)iSunXOffset;
rectOval.top += (float)iSunYOffset;
rectOval.bottom += (float)iSunYOffset;
// now correct for Longitude of Perihelion for this orbit's path
canvas.save();
canvas.rotate((float)dLongPeri, (float)iSunXOffset, (float)iSunYOffset);
canvas.drawOval(rectOval, paint);
canvas.restore();
}
int iBitmapHeight = bmp.getHeight();
canvas.drawBitmap(bmp, fX - (iBitmapHeight / 2), fY - (iBitmapHeight / 2), null);
// draw planet label
myPaint.setColor(Color.WHITE);
paint.setTextSize(30);
canvas.drawText(sName, fX+20, fY-20, paint);
}
上記のメソッドは、さらに 2 つのメソッドを呼び出して、E (平均異常値) と r (惑星が見つかった端のベクトルの長さ) の値を提供します。
getE:
public double getE (double dTime, double dEccentricity)
{
// we are passed the degree count in degrees (duh)
// and the eccentricity value
// the method returns E
double dM1, dD, dE0, dE = 0; // return value E = the mean anomaly
double dM; // local value of M in radians
dM = Math.toRadians (dTime);
int iSign = 1;
if (dM > 0) iSign = 1; else iSign = -1;
dM = Math.abs(dM) / (2 * Math.PI); // Meeus, p 206, line 110
dM = (dM - (long)dM) * (2 * Math.PI) * iSign; // line 120
if (dM < 0)
dM = dM + (2 * Math.PI); // line 130
iSign = 1;
if (dM > Math.PI) iSign = -1; // line 150
if (dM > Math.PI) dM = 2 * Math.PI - dM; // line 160
dE0 = Math.PI / 2; // line 170
dD = Math.PI / 4; // line 170
for (int i = 0; i < 33; i++) // line 180
{
dM1 = dE0 - dEccentricity * Math.sin(dE0); // line 190
dE0 = dE0 + dD * Math.signum((float)(dM - dM1));
dD = dD / 2;
}
dE = dE0 * iSign;
return dE;
}
getRfromE:
public double getRfromE (double dE, double dEccentricty, double dScalar)
{
return Math.min(getWidth(), getHeight()) / 2 * dScalar * (1 - (dEccentricty * Math.cos(dE)));
}
非常に難しく、物理学の深い知識が必要なように見えますが、実際には非常に簡単です。2 つの式とベクトルの基本的な理解が必要です。
質量が m1 と m2 で、それらの間の距離が d の、planet1 と Planet2 の間の引力 (または重力): Fg = G*m1*m2/d^2; Fg = m*a。G は定数です。加速度 "a" が小さすぎず、大きすぎず、およそ "0.01" または "0.1" になるようにランダムな値を代入して求めます。
その瞬間に現在の惑星に作用している全ベクトル力がある場合、瞬間加速度 a=(総力)/(現在の惑星の質量) を求めることができます。そして、現在の加速度と現在の速度と現在の位置があれば、新しい速度と新しい位置を見つけることができます
実際に見たい場合は、次の超簡単なアルゴリズム (疑似コード) を使用できます。
int n; // # of planets
Vector2D planetPosition[n];
Vector2D planetVelocity[n]; // initially set by (0, 0)
double planetMass[n];
while (true){
for (int i = 0; i < n; i++){
Vector2D totalForce = (0, 0); // acting on planet i
for (int j = 0; j < n; j++){
if (j == i)
continue; // force between some planet and itself is 0
Fg = G * planetMass[i] * planetMass[j] / distance(i, j) ^ 2;
// Fg is a scalar value representing magnitude of force acting
// between planet[i] and planet[j]
// vectorFg is a vector form of force Fg
// (planetPosition[j] - planetPosition[i]) is a vector value
// (planetPosition[j]-planetPosition[i])/(planetPosition[j]-plantetPosition[i]).magnitude() is a
// unit vector with direction from planet[i] to planet[j]
vectorFg = Fg * (planetPosition[j] - planetPosition[i]) /
(planetPosition[j] - planetPosition[i]).magnitude();
totalForce += vectorFg;
}
Vector2D acceleration = totalForce / planetMass[i];
planetVelocity[i] += acceleration;
}
// it is important to separate two for's, if you want to know why ask in the comments
for (int i = 0; i < n; i++)
planetPosition[i] += planetVelocity[i];
sleep 17 ms;
draw planets;
}
物理学をシミュレートしている場合は、Box2Dを強くお勧めします。
これは優れた物理シミュレータであり、物理シミュレーションを使用して、必要なボイラープレートの量を実際に削減します。
Bate、Muller、および White によるAstrodynamics の基礎は、学部生の航空宇宙エンジニアのために母校で読む必要があります。これは、地球軌道の物体の軌道力学をカバーする傾向があります...しかし、それはおそらく、理解を開始するために必要な物理学と数学のレベルです.
「Jean Meeusが書いたすべてのもの」に対する@Stefano Boriniの提案に対して+1。
Dear Friend here is the graphics code that simulate solar system
Kindly refer through it
/*Arpana*/
#include<stdio.h>
#include<graphics.h>
#include<conio.h>
#include<math.h>
#include<dos.h>
void main()
{
int i=0,j=260,k=30,l=150,m=90;
int n=230,o=10,p=280,q=220;
float pi=3.1424,a,b,c,d,e,f,g,h,z;
int gd=DETECT,gm;
initgraph(&gd,&gm,"c:\tc\bgi");
outtextxy(0,10,"SOLAR SYSTEM-Appu");
outtextxy(500,10,"press any key...");
circle(320,240,20); /* sun */
setfillstyle(1,4);
floodfill(320,240,15);
outtextxy(310,237,"sun");
circle(260,240,8);
setfillstyle(1,2);
floodfill(258,240,15);
floodfill(262,240,15);
outtextxy(240,220,"mercury");
circle(320,300,12);
setfillstyle(1,1);
floodfill(320,298,15);
floodfill(320,302,15);
outtextxy(335,300,"venus");
circle(320,160,10);
setfillstyle(1,5);
floodfill(320,161,15);
floodfill(320,159,15);
outtextxy(332,150, "earth");
circle(453,300,11);
setfillstyle(1,6);
floodfill(445,300,15);
floodfill(448,309,15);
outtextxy(458,280,"mars");
circle(520,240,14);
setfillstyle(1,7);
floodfill(519,240,15);
floodfill(521,240,15);
outtextxy(500,257,"jupiter");
circle(169,122,12);
setfillstyle(1,12);
floodfill(159,125,15);
floodfill(175,125,15);
outtextxy(130,137,"saturn");
circle(320,420,9);
setfillstyle(1,13);
floodfill(320,417,15);
floodfill(320,423,15);
outtextxy(310,400,"urenus");
circle(40,240,9);
setfillstyle(1,10);
floodfill(38,240,15);
floodfill(42,240,15);
outtextxy(25,220,"neptune");
circle(150,420,7);
setfillstyle(1,14);
floodfill(150,419,15);
floodfill(149,422,15);
outtextxy(120,430,"pluto");
getch();
while(!kbhit()) /*animation*/
{
a=(pi/180)*i;
b=(pi/180)*j;
c=(pi/180)*k;
d=(pi/180)*l;
e=(pi/180)*m;
f=(pi/180)*n;
g=(pi/180)*o;
h=(pi/180)*p;
z=(pi/180)*q;
cleardevice();
circle(320,240,20);
setfillstyle(1,4);
floodfill(320,240,15);
outtextxy(310,237,"sun");
circle(320+60*sin(a),240-35*cos(a),8);
setfillstyle(1,2);
pieslice(320+60*sin(a),240-35*cos(a),0,360,8);
circle(320+100*sin(b),240-60*cos(b),12);
setfillstyle(1,1);
pieslice(320+100*sin(b),240-60*cos(b),0,360,12);
circle(320+130*sin(c),240-80*cos(c),10);
setfillstyle(1,5);
pieslice(320+130*sin(c),240-80*cos(c),0,360,10);
circle(320+170*sin(d),240-100*cos(d),11);
setfillstyle(1,6);
pieslice(320+170*sin(d),240-100*cos(d),0,360,11);
circle(320+200*sin(e),240-130*cos(e),14);
setfillstyle(1,7);
pieslice(320+200*sin(e),240-130*cos(e),0,360,14);
circle(320+230*sin(f),240-155*cos(f),12);
setfillstyle(1,12);
pieslice(320+230*sin(f),240-155*cos(f),0,360,12);
circle(320+260*sin(g),240-180*cos(g),9);
setfillstyle(1,13);
pieslice(320+260*sin(g),240-180*cos(g),0,360,9);
circle(320+280*sin(h),240-200*cos(h),9);
setfillstyle(1,10);
pieslice(320+280*sin(h),240-200*cos(h),0,360,9);
circle(320+300*sin(z),240-220*cos(z),7);
setfillstyle(1,14);
pieslice(320+300*sin(z),240-220*cos(z),0,360,7);
delay(20);
i++;
j++;
k++;
l++;
m++;
n++;
o++;
p++;
q+=2;
}
getch();
}