using System;
using System.Collections.Generic;
using Microsoft.Xna.Framework;
using Microsoft.Xna.Framework.Graphics;
using Microsoft.Xna.Framework.Input;
namespace AstroSim2
{
class Game1 : Game
{
GraphicsDeviceManager graphics;
SpriteBatch sprite;
SpriteFont font;
BasicEffect effect;
VertexBuffer vbLine;
IndexBuffer ibLine;
List<VertexPositionColor> lineVertices = new List<VertexPositionColor>();
List<short> lineIndices = new List<short>();
VertexBuffer vbTriangle;
IndexBuffer ibTriangle;
List<VertexPositionColor> triangleVertices = new List<VertexPositionColor>();
List<short> triangleIndices = new List<short>();
Vector3[] ellipticalOrbit = new Vector3[365 + 1]; // 楕円軌道
int[] ellipseEffect = new int[365 + 1]; // 楕円効果
int[] obliquityEffect = new int[365 + 1]; // 傾斜効果
const int meanSolarDay = 24 * 60 * 60; // 平均太陽日(s)
const double siderealYear = 365.25636; // 恒星年(day) 365d6h9m9.504s
const double solarYear = 365.24219; // 太陽年(day) 365d5h48m45s
const double anomalisticYear = 365.259643; // 近点年(day) 31558433.1552s 365d6h13m53.1552s
const int siderealYearSec = (int)(siderealYear * meanSolarDay); // 恒星年(s)
const int solarYearSec = (int)(solarYear * meanSolarDay); // 太陽年(s)
const int anomalisticYearSec = (int)(anomalisticYear * meanSolarDay); // 近点年(s)
const double au = 149597870700; // 天文単位(m)
const double a = au * 1.00000011; // 軌道長半径(Semimajor axis)(m)
const double e = 0.01671022; // 離心率(Orbital eccentricity)
readonly double K = Math.Sqrt((1 + e) / (1 - e)); // ケプラー方程式の定数
const double siderealDay = 86164.0905; // 恒星日(s) 23h56m4.0905s
const double earthRotationPerSec = (2 * Math.PI) / siderealDay; // 1秒間に地球が回転する角度
const double obliquity = 84381.406; // 黄道傾斜角(degsec) 23deg26m21.406s
const float Sr = 1392000 * 1000 / 2; // 太陽の半径(m)
const double epsilon = 1.0e-14;
Vector3 camPos = new Vector3(0, (float)(au * -2), 0);
float camLat = 0;
float camLon = 90;
public Game1()
{
graphics = new GraphicsDeviceManager(this);
graphics.PreferredBackBufferWidth = 1280;
graphics.PreferredBackBufferHeight = 720;
Content.RootDirectory = "Content";
IsMouseVisible = true;
}
protected override void LoadContent()
{
sprite = new SpriteBatch(GraphicsDevice);
font = Content.Load<SpriteFont>("SpriteFont1");
effect = new BasicEffect(GraphicsDevice);
effect.VertexColorEnabled = true;
effect.Projection = Matrix.CreatePerspectiveFieldOfView(MathHelper.ToRadians(45),
GraphicsDevice.Viewport.AspectRatio, (float)(au * 0.001), (float)(au * 5));
CalcEllipticalOrbit(); // 楕円軌道の計算
CalcEllipseEffect(); // 楕円効果の計算
CalcObliquityEffect(); // 傾斜効果の計算
// 各種ラインの生成
GenerateLines();
vbLine = new VertexBuffer(GraphicsDevice,
typeof(VertexPositionColor), lineVertices.Count, BufferUsage.WriteOnly);
vbLine.SetData(lineVertices.ToArray());
ibLine = new IndexBuffer(GraphicsDevice,
typeof(short), lineIndices.Count, BufferUsage.WriteOnly);
ibLine.SetData(lineIndices.ToArray());
// 太陽の生成
GenerateSun((float)(a * e), 0, 0, Color.Yellow);
vbTriangle = new VertexBuffer(GraphicsDevice,
typeof(VertexPositionColor), triangleVertices.Count, BufferUsage.WriteOnly);
vbTriangle.SetData(triangleVertices.ToArray());
ibTriangle = new IndexBuffer(GraphicsDevice,
typeof(short), triangleIndices.Count, BufferUsage.WriteOnly);
ibTriangle.SetData(triangleIndices.ToArray());
base.LoadContent();
}
// 漸化式によりケプラー方程式を解く
// M 平均近点角(mean anomaly)
// E 離心近点角(Eccentric anomaly)
// T 真近点角(true anomaly)
void KeplersEquation(double M, out double E, out double T)
{
double E0 = M; // 初項
for (int i = 0; ; )
{
i++;
E = M + e * Math.Sin(E0);
if ((E0 - epsilon < E) && (E < E0 + epsilon))
{
break;
}
if (10 <= i)
{
Console.WriteLine(string.Format("計算打ち切り M={0} E={1}", M, E));
break;
}
E0 = E;
}
T = Math.Atan(K * Math.Tan(E / 2)) * 2;
if (T < 0) T += (2 * Math.PI);
}
// 楕円軌道の計算
void CalcEllipticalOrbit()
{
int day = 0;
for (int t = 0; t < siderealYearSec; t += meanSolarDay) // 恒星年で近似
{
double E, T;
KeplersEquation((2 * Math.PI) * t / siderealYearSec, out E, out T);
double r = a * (1 - e * Math.Cos(E));
float x = (float)(Math.Cos(T) * r + a * e); // 焦点との距離を加算
float y = (float)(Math.Sin(T) * r);
ellipticalOrbit[day++] = new Vector3(x, y, 0);
}
}
// 楕円効果の計算
void CalcEllipseEffect()
{
int day = 0;
for (int t = 0; t < anomalisticYearSec; t++) // 近点年
{
double E, T;
KeplersEquation((2 * Math.PI) * t / anomalisticYearSec, out E, out T);
// 地球の自転
double earthRotation = (2 * Math.PI) * ((t / siderealDay) % 1.0);
if (T <= earthRotation)
{
ellipseEffect[day] = t - meanSolarDay * day;
day++;
t += 86164; // 1回転する間は追い越さないので飛ばす
}
else
{
// 角度の差だけ秒を進める。小数点以下は切捨て。桁落ちによる誤差増大に注意
t += (int)((T / earthRotationPerSec) - (earthRotation / earthRotationPerSec));
}
}
}
// 傾斜効果の計算
void CalcObliquityEffect()
{
double radObliq = obliquity * Math.PI / (180 * 60 * 60);
double cosObliq = Math.Cos(radObliq);
int day = 0;
for (int t = 0; t < solarYearSec; t++) // 太陽年
{
// 地球を基準とした太陽の公転
double rad = (2 * Math.PI) * t / solarYearSec;
double x = Math.Cos(rad);
double r = Math.Sin(rad);
double y = cosObliq * r;
double celestialEquator = Math.Atan2(y, x); // 天の赤道上の角度
if (celestialEquator < 0) celestialEquator += (2 * Math.PI);
// 地球の自転
double earthRotation = (2 * Math.PI) * ((t / siderealDay) % 1.0);
if (celestialEquator <= earthRotation)
{
obliquityEffect[day] = t - meanSolarDay * day;
day++;
t += 86164; // 1回転する間は追い越さないので飛ばす
}
else
{
t += (int)((celestialEquator / earthRotationPerSec)
- (earthRotation / earthRotationPerSec));
}
}
}
// 各種ラインの生成
void GenerateLines()
{
double aspectRatio = Math.Sqrt(1 - e * e); // 楕円軌道の長短比
double b = a * aspectRatio; // 軌道短半径(m)
// 長径と短径
lineVertices.Add(new VertexPositionColor(new Vector3((float)-a, 0, 0), Color.White));
lineVertices.Add(new VertexPositionColor(new Vector3((float)a, 0, 0), Color.White));
lineVertices.Add(new VertexPositionColor(new Vector3(0, (float)-b, 0), Color.White));
lineVertices.Add(new VertexPositionColor(new Vector3(0, (float)b, 0), Color.White));
for (int i = 0; i < 4; i++)
{
lineIndices.Add((short)i);
}
// 楕円軌道
int vertexNum = lineVertices.Count;
for (int i = 0; i < 365; i++)
{
lineVertices.Add(new VertexPositionColor(ellipticalOrbit[i], Color.White));
lineIndices.Add((short)(vertexNum + i));
lineIndices.Add((short)(vertexNum + (i + 1) % 365));
}
// 楕円効果
vertexNum = lineVertices.Count;
for (int i = 0; i < 365; i++)
{
Vector3 pos = ellipticalOrbit[i];
pos.Z = (float)(Sr * ellipseEffect[i] / 30); // 30秒進むと太陽半個分上昇
lineVertices.Add(new VertexPositionColor(pos, Color.Cyan));
lineIndices.Add((short)(vertexNum + i));
lineIndices.Add((short)(vertexNum + (i + 1) % 365));
}
// 傾斜効果
vertexNum = lineVertices.Count;
for (int i = 0; i < 365; i++)
{
Vector3 pos = ellipticalOrbit[i];
int j = (290 + i) % 365; // 春分点3/20 近日点1/4 と仮定
pos.Z = (float)(Sr * obliquityEffect[j] / 30);
lineVertices.Add(new VertexPositionColor(pos, Color.Yellow));
lineIndices.Add((short)(vertexNum + i));
lineIndices.Add((short)(vertexNum + (i + 1) % 365));
}
// 均時差
vertexNum = lineVertices.Count;
for (int i = 0; i < 365; i++)
{
Vector3 pos = ellipticalOrbit[i];
int j = (290 + i) % 365;
pos.Z = (float)(Sr * (ellipseEffect[i] + obliquityEffect[j]) / 30);
lineVertices.Add(new VertexPositionColor(pos, Color.Red));
lineIndices.Add((short)(vertexNum + i));
lineIndices.Add((short)(vertexNum + (i + 1) % 365));
}
}
void GenerateSun(float x, float y, float z, Color c)
{
int i = triangleVertices.Count;
triangleVertices.Add(new VertexPositionColor(new Vector3(x, y, z + Sr), c));
triangleVertices.Add(new VertexPositionColor(new Vector3(x + Sr, y, z), c));
triangleVertices.Add(new VertexPositionColor(new Vector3(x, y + Sr, z), c));
triangleVertices.Add(new VertexPositionColor(new Vector3(x - Sr, y, z), c));
triangleVertices.Add(new VertexPositionColor(new Vector3(x, y - Sr, z), c));
triangleVertices.Add(new VertexPositionColor(new Vector3(x, y, z - Sr), c));
GenerateTriangle(i, 0, 2, 1);
GenerateTriangle(i, 0, 3, 2);
GenerateTriangle(i, 0, 4, 3);
GenerateTriangle(i, 0, 1, 4);
GenerateTriangle(i, 1, 2, 5);
GenerateTriangle(i, 2, 3, 5);
GenerateTriangle(i, 3, 4, 5);
GenerateTriangle(i, 4, 1, 5);
}
void GenerateTriangle(int vertexNum, int v1, int v2, int v3)
{
triangleIndices.Add((short)(vertexNum + v1));
triangleIndices.Add((short)(vertexNum + v2));
triangleIndices.Add((short)(vertexNum + v3));
}
protected override void Update(GameTime gameTime)
{
KeyboardState kState = Keyboard.GetState();
if (kState.IsKeyDown(Keys.Escape)) Exit();
if (kState.IsKeyDown(Keys.W)) Move(0, 0);
if (kState.IsKeyDown(Keys.S)) Move(180, 0);
if (kState.IsKeyDown(Keys.A)) Move(0, 90);
if (kState.IsKeyDown(Keys.D)) Move(0, -90);
if (kState.IsKeyDown(Keys.Up)) camLat = Math.Min(camLat + 0.5f, 89.9f);
if (kState.IsKeyDown(Keys.Down)) camLat = Math.Max(camLat - 0.5f, -89.9f);
if (kState.IsKeyDown(Keys.Left)) camLon = (camLon + 1) % 360;
if (kState.IsKeyDown(Keys.Right)) camLon = (camLon + 359) % 360;
if (kState.IsKeyDown(Keys.PageUp)) Move(90, 0);
if (kState.IsKeyDown(Keys.PageDown)) Move(-90, 0);
base.Update(gameTime);
}
void Move(float lat, float lon)
{
float rad = MathHelper.ToRadians(camLat + lat);
float z = (float)(Math.Sin(rad) * au * 0.0025);
float r = (float)(Math.Cos(rad) * au * 0.0025);
if (lon == 0) camPos.Z += z;
rad = MathHelper.ToRadians(camLon + lon);
float x = (float)Math.Cos(rad) * r;
float y = (float)Math.Sin(rad) * r;
camPos.X += x;
camPos.Y += y;
}
protected override void Draw(GameTime gameTime)
{
GraphicsDevice.Clear(Color.CornflowerBlue);
GraphicsDevice.DepthStencilState = DepthStencilState.Default;
//GraphicsDevice.BlendState = BlendState.AlphaBlend;
// カメラ
float rad = MathHelper.ToRadians(camLat);
float z = (float)(Math.Sin(rad) * au);
float r = (float)(Math.Cos(rad) * au);
rad = MathHelper.ToRadians(camLon);
float x = (float)Math.Cos(rad) * r;
float y = (float)Math.Sin(rad) * r;
effect.View = Matrix.CreateLookAt(camPos, camPos + new Vector3(x, y, z), Vector3.UnitZ);
foreach (EffectPass pass in effect.CurrentTechnique.Passes)
{
pass.Apply();
GraphicsDevice.SetVertexBuffer(vbLine);
GraphicsDevice.Indices = ibLine;
GraphicsDevice.DrawIndexedPrimitives(PrimitiveType.LineList,
0, 0, vbLine.VertexCount, 0, ibLine.IndexCount / 2);
GraphicsDevice.SetVertexBuffer(vbTriangle);
GraphicsDevice.Indices = ibTriangle;
GraphicsDevice.DrawIndexedPrimitives(PrimitiveType.TriangleList,
0, 0, vbTriangle.VertexCount, 0, ibTriangle.IndexCount / 3);
}
sprite.Begin();
string text = string.Format("x={0:f2} y={1:f2} z={2:f2}",
camPos.X / au, camPos.Y / au, camPos.Z / au);
sprite.DrawString(font, text, new Vector2(0, 0), Color.White);
text = string.Format("lat={0:f0} lon={1:f0}", camLat, camLon);
sprite.DrawString(font, text, new Vector2(0, 20), Color.White);
sprite.End();
base.Draw(gameTime);
}
}
}
