/*
* プロジェクトのプロパティ
* [XNA Game Studio]タブ
* Use HiDef to access the complete API
*/
using System;
using Microsoft.Xna.Framework;
using Microsoft.Xna.Framework.Graphics;
using Microsoft.Xna.Framework.Input;
namespace AstroSim3
{
class Game1 : Game
{
GraphicsDeviceManager graphics;
SpriteBatch sprite;
SpriteFont font;
Effect effect;
VertexBuffer vertexBuffer;
Texture2D texture;
KeyboardState kStateOld = new KeyboardState();
// fps
int fpsSec = -1;
int fpsDraw = 0;
int fpsCount = 0;
const double solarYear = 365.24219; // 太陽年(day) 365d5h48m45s
const double anomalisticYear = 365.259643; // 近点年(day) 365d6h13m53.1552s
const double earthRotationPerSec = (2 * Math.PI) / 86400; // 地球が1秒間に回転する角度
const double e = 0.01671022; // 離心率(Orbital eccentricity)
readonly double K = Math.Sqrt((1 + e) / (1 - e)); // ケプラー方程式の定数
const double epsilon = 1.0e-14;
DateTime dt = DateTime.UtcNow;
int tz = 0;
readonly string[] tzName = { "UTC", "JST-9" };
int cursor = 0;
readonly int[] cursorPos = { 0, 5, 8, 11, 14, 15, 17 };
readonly int[] cursorLen = { 4, 2, 2, 2, 1, 1, 3 };
int repeatFrame = 0;
double MJD; // 修正ユリウス日
double MJD_date;
double MJD_time;
double T; // 2000/1/1 12:00(UT)からのユリウス世紀(36525日)
double obliquity; // 黄道傾斜角
double eclipticLon; // 黄径
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");
texture = Content.Load<Texture2D>("earthmap1k");
effect = Content.Load<Effect>("Effect1");
effect.Parameters["WorldMap"].SetValue(texture);
effect.Parameters["revision"].SetValue(
(float)Math.Sin((0.533 / 2 + (35 * 60 + 8) / 3600.0) * Math.PI / 180));
VertexPositionTexture[] vertices = new VertexPositionTexture[4];
vertices[0] = new VertexPositionTexture(new Vector3(-1, 1, 0), new Vector2(0, 0));
vertices[1] = new VertexPositionTexture(new Vector3(1, 1, 0), new Vector2(1, 0));
vertices[2] = new VertexPositionTexture(new Vector3(-1, -1, 0), new Vector2(0, 1));
vertices[3] = new VertexPositionTexture(new Vector3(1, -1, 0), new Vector2(1, 1));
vertexBuffer = new VertexBuffer(GraphicsDevice,
typeof(VertexPositionTexture), 4, BufferUsage.WriteOnly);
vertexBuffer.SetData(vertices);
CalcParameter();
base.LoadContent();
}
void CalcParameter()
{
// 修正ユリウス日
int y = dt.Year;
int m = dt.Month;
int d = dt.Day;
if (m < 3)
{
y--;
m += 12;
}
MJD_date = (int)(365.25 * y) + (y / 400) - (y / 100) + (int)(30.59 * (m - 2)) + d - 678912;
MJD_time = (dt.Hour * 60 + dt.Minute) / 1440.0;
MJD = MJD_date + MJD_time;
// 黄道傾斜角
T = (MJD - 51544.5) / 36525;
obliquity = (84381.406 - 46.836769 * T - 0.00059 * T * T + 0.001813 * T * T * T) / 3600;
// 平均近点角（概算）近日点から次の近日点までの角度
double Ma = ModAngle((2 * Math.PI) * ((MJD / anomalisticYear - 0.1242853) % 1.0));
double E;
double Ta; // 真近点角
KeplersEquation(Ma, out E, out Ta);
// 春分点(vernal equinox)の真近点角
double MJDv = ((int)(MJD / solarYear - 0.3399541) + 0.3399541) * solarYear;
double Mv = ModAngle((2 * Math.PI) * ((MJDv / anomalisticYear - 0.1242853) % 1.0));
double Tv;
KeplersEquation(Mv, out E, out Tv);
// 黄径（概算）春分点から次の春分点までの角度
eclipticLon = ModAngle(Ta - Tv);
// 楕円効果と傾斜効果
int ellipseEffect = (int)Math.Round(Ma / earthRotationPerSec - Ta / earthRotationPerSec);
int obliquityEffect = CalcObliquityEffect(obliquity, eclipticLon);
int equationOfTime = ellipseEffect + obliquityEffect; // 均時差
float transit = (43200 - equationOfTime) / 86400.0f; // 南中時
// 天球上の太陽軌道の高さと半径
double solarDecl = Math.Sin(eclipticLon) * obliquity * Math.PI / 180; // 太陽の赤緯
double solarAlt = Math.Sin(solarDecl);
double solarRad = Math.Cos(solarDecl);
effect.Parameters["MJD_time"].SetValue((float)MJD_time);
effect.Parameters["solarDecl"].SetValue((float)solarDecl);
effect.Parameters["solarAlt"].SetValue((float)solarAlt);
effect.Parameters["solarRad"].SetValue((float)solarRad);
effect.Parameters["transit"].SetValue(transit);
}
double ModAngle(double angle)
{
while (angle <= -Math.PI) angle += (2 * Math.PI);
while (Math.PI < angle) angle -= (2 * Math.PI);
return angle;
}
// 傾斜効果の計算
int CalcObliquityEffect(double obliquity, double eclipticLon)
{
// 地球を基準とした太陽の公転
double x = Math.Cos(eclipticLon);
double r = Math.Sin(eclipticLon);
double y = Math.Cos(obliquity * Math.PI / 180) * r;
double celestialEquator = Math.Atan2(y, x); // 天の赤道上の角度
return (int)Math.Round(eclipticLon / earthRotationPerSec - celestialEquator / earthRotationPerSec);
}
// 漸化式によりケプラー方程式を解く
// 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;
}
protected override void Update(GameTime gameTime)
{
int delta = 0;
KeyboardState kState = Keyboard.GetState();
if (kState.IsKeyDown(Keys.Escape)) Exit();
if (kState.IsKeyDown(Keys.Left) && kStateOld.IsKeyUp(Keys.Left))
{
cursor = (cursor + 6) % 7;
}
if (kState.IsKeyDown(Keys.Right) && kStateOld.IsKeyUp(Keys.Right))
{
cursor = (cursor + 1) % 7;
}
if (kState.IsKeyDown(Keys.Up)) delta = 1;
if (kState.IsKeyDown(Keys.Down)) delta = -1;
kStateOld = kState;
if (delta == 0)
{
repeatFrame = 0;
}
else
{
if (--repeatFrame <= 0)
{
switch (cursor)
{
case 0: dt = dt.AddYears(delta); break;
case 1: dt = dt.AddMonths(delta); break;
case 2: dt = dt.AddDays(delta); break;
case 3: dt = dt.AddHours(delta); break;
case 4: dt = dt.AddMinutes(delta * 10); break;
case 5: dt = dt.AddMinutes(delta); break;
case 6: tz = (tz + 1) % 2; break;
}
if (cursor != 6)
{
CalcParameter();
}
repeatFrame = (repeatFrame == 0) ? 6 : 30;
}
}
base.Update(gameTime);
}
protected override void Draw(GameTime gameTime)
{
GraphicsDevice.Clear(Color.CornflowerBlue);
GraphicsDevice.SamplerStates[0] = SamplerState.LinearClamp;
foreach (EffectPass pass in effect.CurrentTechnique.Passes)
{
pass.Apply();
GraphicsDevice.SetVertexBuffer(vertexBuffer);
GraphicsDevice.DrawPrimitives(PrimitiveType.TriangleStrip, 0, 2);
}
// fps
fpsDraw++;
if (gameTime.TotalGameTime.Seconds != fpsSec)
{
fpsCount = fpsDraw;
fpsDraw = 0;
fpsSec = gameTime.TotalGameTime.Seconds;
}
sprite.Begin();
string text = new string(' ', cursorPos[cursor]) + new string('_', cursorLen[cursor]);
sprite.DrawString(font, text, new Vector2(0, 0), Color.Red);
text = string.Format("{0} {1}", dt.AddHours(tz * 9).ToString("yyyy/MM/dd HH:mm"), tzName[tz]);
sprite.DrawString(font, text, new Vector2(0, 0), Color.Red);
text = string.Format("fps={0} MJD={1:f4}", fpsCount, MJD);
sprite.DrawString(font, text, new Vector2(0, 20), Color.Red);
text = string.Format("eLon={0:f6} T={1:f7} obliq={2:f3}", eclipticLon, T, obliquity);
sprite.DrawString(font, text, new Vector2(0, 40), Color.Red);
sprite.End();
base.Draw(gameTime);
}
}
}
texture WorldMap;
float MJD_time; // 修正ユリウス日の時刻
float solarDecl; // 太陽の赤緯
float solarAlt; // 天球上の太陽軌道の高さ
float solarRad; // 天球上の太陽軌道の半径
float transit; // 南中時
float revision; // 太陽の視角(0.533deg)と大気差(35m8s)による昼の長さの補正
sampler TextureSampler = sampler_state
{
texture = <WorldMap>;
mipfilter = linear;
minfilter = linear;
magfilter = linear;
};
struct VertexShaderInput
{
float4 Position : POSITION0;
float2 TexCoord : TEXCOORD0;
};
struct VertexShaderOutput
{
float4 Position : POSITION0;
float2 TexCoord : TEXCOORD0;
};
VertexShaderOutput VertexShaderFunction(VertexShaderInput input)
{
VertexShaderOutput output;
output.Position = input.Position;
output.TexCoord = input.TexCoord;
return output;
}
// input.TexCoord.x 経度 Left=0.0 - Right=1.0
// input.TexCoord.y 緯度 Top=0.0 - Bottom=1.0
float4 PixelShaderFunction(VertexShaderOutput input) : COLOR0
{
float lat = radians(90 - input.TexCoord.y * 180); // 緯度 Top=90 - Bottom=-90
float lon = -180 + input.TexCoord.x * 360; // 経度 Left=-180 - Right=180
// 黄緯線
if (abs(lat - solarDecl) < 0.0025) return float4(1, 1, 0, 1);
// 赤道
if (abs(lat) < 0.0025) return float4(1, 0, 0, 1);
// ローカル時刻
float localtime = MJD_time + lon / 360;
localtime = frac(localtime);
if (localtime < 0.001) return float4(1, 1, 1, 1); // 00:00
// テクスチャサンプリング
float4 output;
output = tex2D(TextureSampler, input.TexCoord);
// 天球上の太陽軌道と地平面の交点＝日出・日没
float ra = sin(lat) * revision;
float rx = cos(lat) * revision;
float x = ((solarAlt + ra) * -tan(lat) - rx) / solarRad;
if (x <= -1) // 白夜
{
return output;
}
if (1 <= x) // 極夜
{
output.rgb *= 0.5;
return output;
}
float halfDaytime = degrees(acos(x)) / 360;
// 日出・日没時刻（南中時－半分の昼の長さ）
float rising = transit - halfDaytime;
float setting = transit + halfDaytime;
if (localtime <= rising || setting <= localtime)
{
output.rgb *= 0.5;
}
return output;
}
technique Technique1
{
pass Pass1
{
VertexShader = compile vs_3_0 VertexShaderFunction();
PixelShader = compile ps_3_0 PixelShaderFunction();
}
}
