地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その５）

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その５）

---

        public double tjdCalc(Dtime a)                  //Tjdの計算
        {
            int jm = a.Month;
            int jy = a.Year;
            if (a.Month < 3)
            {
                jy = a.Year - 1;
                jm += 12;
            }
            double jd = (int)(365.25 * jy) + (int)(jy / 400) - (int)(jy / 100)
                + (int)(30.59 * (jm - 2)) + a.Day + 1721088.5
                + (double)a.Hour / 24 + (double)a.Minute / 1440 + (double)a.Second / 86400;
            double tjd = jd - 2440000.5;
            return tjd;
        }
        private void orbital_sphere()                   //平面地図に軌道を表示する
        {
            keido = keido - 1;
            xl = mx;
            if (mb2 == true)
            {
                mx = System.Windows.Forms.Cursor.Position.X;
                if (mb == false)
                {
                    xl = mx;
                    mb = true;
                }
            }
            int dx = xl - mx;
            if (Math.Abs(dx) > 360)
                dx = 0;
            keido = keido + dx;
            if (keido >= 360)
                keido -= 360;
            if (keido < 0)
                keido += 360;
            int rkeido = keido - 90;
            if (rkeido > 180)
                rkeido -= 360;
            string rkeidoDirection = "E";
            if (rkeido < 0)
            {
                rkeidoDirection = "W";
                rkeido *= -1;
            }
            this.label2.Text = rkeido.ToString() + rkeidoDirection;

            DateTime pDateJ = DateTime.Now;             //現地日時表示
            this.label1.Text = pDateJ.ToString();
            
            DateTime crnt1 = DateTime.UtcNow;           //現在日時をTJDへ
            Dtime crnt2;
            crnt2.Year = crnt1.Year;
            crnt2.Month = crnt1.Month;
            crnt2.Day = crnt1.Day;
            crnt2.Hour = crnt1.Hour;
            crnt2.Minute = crnt1.Minute;
            crnt2.Second = crnt1.Second;
            double tjdCurrent = tjdCalc(crnt2);
            
            Bitmap canvas = new Bitmap(bmp3);           //軌道計算
            Graphics g = Graphics.FromImage(canvas);
            Pen p = new Pen(Color.FromArgb(128, Color.Red), 1);
            Font fnt = new Font("MS UI Gothic", 9);
            int iff = 120;
            int[] xll = new int[iff + 1];
            int[] yl = new int[iff + 1];
            double[] hsa = new double[iff + 1];
            for (int i = 0; i < iff + 1; i++)
            {
                                                        //グリニッジ平均恒星時計算
                double tg = 24 * (0.671262 + 1.0027379094 * (tjdCurrent + (double)(iff - i) / 60 / 24));
                tg = tg % 24;
                double tl = tg + keido / 15;
                if (tl > 24)
                    tl -= 24;
                                                        //真近点離角計算
                double fai = ((tjdCurrent + (double)(iff - i) / 60 / 24 - t1.tjdep)
                    * t1.mm * 2 * Math.PI + t1.mk / 180 * Math.PI) % (Math.PI * 2);
                double eec = fai + t1.ec * Math.Sin(fai);
                double ek, de;
                do
                {
                    ek = fai + t1.ec * Math.Sin(eec);
                    de = Math.Abs(eec - ek);
                    eec = ek;
                }
                while (de > 0.000000000001);
                double tanvk = Math.Sqrt(1 - t1.ec * t1.ec) * Math.Sin(ek) / (Math.Cos(ek) - t1.ec);
                double vk = Math.Atan(tanvk);
                double vkj = Math.Cos(ek) - t1.ec;
                if (vkj < 0)
                    vk = vk + Math.PI;
                                                        //軌道位置計算
                double tt = 24 * 3600 / t1.mm;          //公転周期
                double ar = (6.67384E-11 * 5.972E+24 * tt * tt / 4 / Math.PI / Math.PI);
                ar = Math.Pow(ar, (double)1 / (double)3);       //軌道長半径
                double rk = ar * (1 - t1.ec * Math.Cos(ek));
                double xk = rk * Math.Cos(vk);
                double yk = rk * Math.Sin(vk);
                                                        //地心赤道直交座標変換
                double xs = xk * t1.sx1 - yk * t1.sx2;
                double ys = xk * t1.sy1 - yk * t1.sy2;
                double zs = xk * t1.sz1 + yk * t1.sz2;
                                                        //G系地心直交座標変換
                tg = tg / 24 * 2 * Math.PI;
                double xg = xs * Math.Cos(tg) + ys * Math.Sin(tg);
                double yg = ys * Math.Cos(tg) - xs * Math.Sin(tg);
                double zg = zs;
                                                        //G系地心直交座標から経度・緯度及び楕円体高変換
                double skeido = Math.Atan(yg / xg);
                if (xg < 0)
                    skeido = skeido + Math.PI;
                string LongitudeDirection = "E";
                double skeidox = skeido;
                if (skeido > Math.PI)
                {
                    skeido = 2 * Math.PI - skeido;
                    skeidox = -skeido;
                    LongitudeDirection = "W";
                }
                if (skeido < 0)
                {
                    skeidox = skeido;
                    skeido *= -1;
                    LongitudeDirection = "W";
                }
                double pg = Math.Sqrt(xg * xg + yg * yg);
                double sido0 = Math.Atan(zg / pg);
                double ds, sido, n;
                double ag = 6377397.155;
                double e2 = 0.006674372230614;
                do
                {
                    n = ag / Math.Sqrt(1 - e2 * e2 * Math.Sin(sido0) * Math.Sin(sido0));
                    sido = Math.Atan(zg / (pg - e2 * e2 * n * Math.Cos(sido0)));
                    ds = sido - sido0;
                    sido0 = sido;
                }
                while (ds > 0.000000000001);
                string LatitudeDirection = "N";
                double sidoy = sido;
                if (sido < 0)
                {
                    LatitudeDirection = "S";
                    sido *= -1;
                }
                double hs = pg / Math.Cos(sido);
                skeido = skeido / Math.PI * 180;
                sido = sido / Math.PI * 180;
                                                        //表示する
                xll[i] = (int)((double)canvas.Width / 2 * (1 + skeidox / Math.PI));
                yl[i] = (int)((double)canvas.Height / 2 * (1 - 2 * sidoy / Math.PI));
                hsa[i] = hs / n;
                hs -= n;
                hs /= 1000;
                int rl = 4;
                if (i == iff)
                {
                    p = new Pen(Color.FromArgb(255, Color.AntiqueWhite), 2);
                    rl = 6;
                    int cxl = xll[i];
                    int cyl = yl[i];
                    if (cxl > h + h - 180)
                        cxl = h + h - 180;
                    g.DrawString("LAT=" + sido.ToString("###") + LatitudeDirection
                        + ",LON=" + skeido.ToString("###") + LongitudeDirection
                        + ",ALT=" + hs.ToString("###KM"), fnt, Brushes.Coral, cxl, cyl + 5);
                }
                g.DrawEllipse(p, xll[i] - rl / 2, yl[i] - rl / 2, rl, rl);
            }
            fnt.Dispose();
            g.Dispose();
            this.pictureBox7.Image = canvas;
            Bitmap sphereBmp = new Bitmap(h, h);
            sphereBmp = sphere(xll, yl, hsa, iff);
            this.pictureBox9.Image = sphereBmp;
        }

---

次の記事へ続く。

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その４）

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その４）

---

        public void marble()                            //昼夜混合画像変換
        {
            //
            //
            //System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
            //sw.Start();
            //
            //
            DateTime pDate = DateTime.Now;              //現在日時と季節判定
            double sDay = pDate.DayOfYear;
            double pSeason = (sDay / 365.25 - 0.72) * 2 * Math.PI;
            if (pSeason < 0)
                pSeason = pSeason + 2 * Math.PI;
            double aSeason = Math.Tan((90 - 23.4 * Math.Sin(pSeason)) / 180 * Math.PI);
            int sHur = pDate.Hour;
            int pUTC = 9;
            if (sHur < pUTC)
                sHur = sHur - pUTC;
            else
                sHur = sHur + (24 - pUTC);
            double sTime = (double)(sHur * 3600 + pDate.Minute * 60 + pDate.Second) / 24 / 3600 * 360;

            Bitmap canvas3 = new Bitmap(w, h);          //画像セット
            
            int[] yDay = new int[w - 1];                //明暗境界線計算
            for (int j = 0; j < w - 1; j++)
            {
                double yDayD = Math.Atan(aSeason * Math.Cos(((double)j
                    / (double)w * 360 + sTime - 180) / 180 * Math.PI))
                    / Math.PI * h + h2;
                yDay[j] = (int)yDayD;
            }
            Color c, c1, c2;                            //昼夜計算
            byte cR, cG, cB, cA;
            int jkk, ikk;
            int[,] DorN = new int[w, h];
            int Dsum;
            int Eve = w / 48;
            int Eve2 = (Eve + 1) * (Eve + 1);
            for (int j = 0; j < w - 1; j++)
            {
                for (int i = 0; i < h - 1; i++)
                {
                    if (pSeason < Math.PI)
                    {
                        if (i > yDay[j])
                        {
                            DorN[j, i] = 1;
                        }
                        else
                        {
                            DorN[j, i] = 0;
                        }
                    }
                    else
                    {
                        if (i > yDay[j])
                        {
                            DorN[j, i] = 0;
                        }
                        else
                        {
                            DorN[j, i] = 1;
                        }
                    }
                }
            }
            for (int j = 0; j < w - 1; j++)             //グラデーション生成
            {
                for (int i = 0; i < h - 1; i++)
                {
                    Dsum = 0;
                    for (int jk = 0; jk < Eve; jk++)
                    {
                        for (int ik = 0; ik < Eve; ik++)
                        {
                            jkk = j + jk + jk + jk + jk - Eve - Eve;
                            if (jkk < 0)
                                jkk += w;
                            else if (jkk > w - 1)
                                jkk -= w;
                            ikk = i + ik + ik + ik + ik - Eve - Eve;
                            if (ikk < 0)
                            {
                                ikk = -ikk;
                                jkk += w2;
                            }
                            else if (ikk > h - 1)
                            {
                                ikk = 2 * h - 2 - ikk;
                                jkk += w2;
                            }
                            if (jkk > w - 1)
                            {
                                jkk -= w;
                            }
                            Dsum += DorN[jkk, ikk];
                        }
                    }
                    c1 = bmp1.GetPixel(j, i);
                    c2 = bmp2.GetPixel(j, i);
                    double a2 = (double)Dsum / Eve2;
                    double a1 = 1 - a2;
                    cR = (byte)((double)c2.R * a1 + (double)c1.R * a2);
                    cG = (byte)((double)c2.G * a1 + (double)c1.G * a2);
                    cB = (byte)((double)c2.B * a1 + (double)c1.B * a2);
                    cA = (byte)((double)c2.A * a1 + (double)c1.A * a2);
                    c = Color.FromArgb(cA, cR, cG, cB);
                    bmp3.SetPixel(j, i, c);
                }
            }
            //
            //
            //sw.Stop();
            //MessageBox.Show(sw.Elapsed.ToString());
            //
            //
        }

---

次の記事へ続く。

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その３）

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その３）

---

         public Image resizeImage(Image jpg, int resizeHeight)      //画像サイズ変換
        {
            int resizeWidth = (int)(jpg.Width * ((double)resizeHeight / (double)jpg.Height));
            Bitmap resizeBmp = new Bitmap(resizeWidth, resizeHeight);
            Graphics g = Graphics.FromImage(resizeBmp);
            g.InterpolationMode = System.Drawing.Drawing2D.InterpolationMode.HighQualityBicubic;
            g.DrawImage(jpg, 0, 0, resizeWidth, resizeHeight);
            g.Dispose();
            return resizeBmp;
        }
        public int[] matY(int imageHeight)              //平面→球面投影座標変換Y配列
        {
            int[] ymat = new int[imageHeight - 1];

            for (int i = 0; i < imageHeight - 1; i++)
            {
                ymat[i] = (int)((double)imageHeight / Math.PI * Math.Acos(1 - 2 * i / (double)imageHeight));
            }
            return ymat;
        }
        public int[,] matX(int[] ymat, int iH)          //平面→球面投影座標変換X配列
        {
            int[,] xmat = new int[iH - 1, iH - 1];
            double f;
            for (int i = 0; i < iH - 1; i++)
            {
                for (int j = 0; j < iH - 1; j++)
                {
                    xmat[j, i] = iH + 1;
                    double ch = (double)(Math.Sqrt((i - iH / 2) * (i - iH / 2) + (j - iH / 2) * (j - iH / 2)));
                    if (ch < iH / 2)
                    {
                        f = (1 - 2 * (double)j / (double)iH) / Math.Sin((double)ymat[i] * Math.PI / (double)iH);
                        if (f > 0.99)
                            f = 0.99;
                        if (f < -1)
                            f = -1;
                        xmat[j, i] = (int)((double)iH / Math.PI * Math.Acos(f));
                    }
                }
            }
            return xmat;
        }

---

 次の記事へ続く。

 

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その２）

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その２）

---

        private void timer3_Tick(object sender, EventArgs e)    //昼夜再描画タイマー
        {
            marble();
        }
        private void timer4_Tick(object sender, EventArgs e)    //地球回転再描画タイマー
        {
            orbital_sphere();
        }
        private Bitmap makeBmp_desktop(string filename)         //ビットマップ読み込み＆サイズ変換
        {
            Bitmap bmp = new Bitmap(h + h, h);
            System.Drawing.Image jpg;
            try
            {
                jpg = System.Drawing.Bitmap.FromFile(filename);
                jpg = resizeImage(jpg, h);
                bmp = new Bitmap(jpg);
                bmpGetTry = true;
            }
            catch (System.IO.FileNotFoundException)
            {
                string eMsg = "File not found";
                this.label2.Text = eMsg;
                bmpGetTry = false;
            }
            return bmp;
        }
        private Bitmap makeBmp(string url, string filename)         //ビットマップ読み込み＆サイズ変換
        {
            Bitmap bmp = new Bitmap(h + h, h);
            WebClient wc = new WebClient();
            System.Drawing.Image jpg;
            Stream st2 = wc.OpenRead(url);
            jpg = new Bitmap(st2);
            jpg.Save(filename);
            jpg = resizeImage(jpg, h);
            bmp = new Bitmap(jpg);
            st2.Close();
            return bmp;
        }
        private double telDate()                        //TLE形式シリアル日時計算
        {
            DateTime pDate = DateTime.UtcNow;
            double cuday2 = pDate.DayOfYear;
            cuday2 = cuday2 + (double)(pDate.Hour * 3600 + pDate.Minute * 60 + pDate.Second) / 24 / 3600;
            double cu = (pDate.Year - 2000) * 1000 + cuday2;
            return cu;
        }
        private TLEdata ISStle(string url)              //TLEデータ読み込み
        {
            double cu = telDate();                      //現在日時設定
            WebClient wc = new WebClient();             //webからTLE読み込み・抽出
            TLEdata t;
            t = new TLEdata();
            try
            {
                Stream st = wc.OpenRead(url);
                StreamReader sr = new StreamReader(st, Encoding.GetEncoding("Shift_JIS"));
                Boolean readstatus = false;
                while (readstatus == false)
                {
                    string cText0 = sr.ReadLine();
                    if (cText0 == null)
                        readstatus = true;
                    if (cText0.Length == 7)
                    {
                        string wordPick = cText0.Substring(4, 3);
                        if (wordPick == "ISS")
                        {
                            string cText1 = sr.ReadLine();
                            string cText2 = sr.ReadLine();
                            double ep = double.Parse(cText1.Substring(22, 14));
                            double incc = double.Parse(cText2.Substring(12, 8));
                            double ra = double.Parse(cText2.Substring(21, 8));
                            t.ec = double.Parse(cText2.Substring(30, 7)) / 10000000;
                            double ap = double.Parse(cText2.Substring(38, 8));
                            t.mk = double.Parse(cText2.Substring(47, 8));
                            t.mm = double.Parse(cText2.Substring(56, 11));
                            if (ep > cu - 1)
                                readstatus = true;
                            incc = incc / 180 * Math.PI;
                            ra = ra / 180 * Math.PI;
                            ap = ap / 180 * Math.PI;
                                                        //JD translation for Epoch
                            Dtime epc;
                            int[] mday = { 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
                            epc.Year = (int)(ep / 1000) + 2000;
                            double epday2 = ep - (epc.Year - 2000) * 1000;
                            int rep = epc.Year % 4;
                            if (rep == 0)
                                mday[2] = 29;
                            epc.Month = 0;
                            epc.Day = 0;
                            int i = 0;
                            while (epday2 > 0)
                            {
                                epc.Month++;
                                epc.Day = (int)epday2;
                                epday2 = epday2 - mday[++i];
                            }
                            epday2 = epday2 + mday[i];
                            epc.Hour = (int)((epday2 - (int)epday2) * 24);
                            epc.Minute = (int)((epday2 - (int)epday2 - (double)epc.Hour / 24) * 24 * 60);
                            epc.Second = (int)((epday2 - (int)epday2 - (double)epc.Hour / 24
                                - (double)epc.Minute / 24 / 60) * 24 * 60 * 60);
                            t.tjdep = tjdCalc(epc);
                                                        //地心赤道直交座標変換用係数
                            t.sx1 = Math.Cos(ra) * Math.Cos(ap) - Math.Sin(ra) * Math.Cos(incc) * Math.Sin(ap);
                            t.sx2 = Math.Cos(ra) * Math.Sin(ap) + Math.Sin(ra) * Math.Cos(incc) * Math.Cos(ap);
                            t.sy1 = Math.Sin(ra) * Math.Cos(ap) + Math.Cos(ra) * Math.Cos(incc) * Math.Sin(ap);
                            t.sy2 = Math.Sin(ra) * Math.Sin(ap) - Math.Cos(ra) * Math.Cos(incc) * Math.Cos(ap);
                            t.sz1 = Math.Sin(incc) * Math.Sin(ap);
                            t.sz2 = Math.Sin(incc) * Math.Cos(ap);
                        }
                    }
                }
                sr.Close();
                st.Close();
            }
            catch (System.IO.FileNotFoundException)
            {
                string eMsg = "File not found";
                this.label2.Text = eMsg;
            }
            return t;
        }

---

次の記事へ続く。

 

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その１）

地球儀上でISS(国際宇宙ステーション)の位置を表示するプログラム（その１）

ソースの備忘録です。

---

using System;                                           //My sphere with ISS orbital 2014/09/05
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using System.Net;
using System.IO;
using System.Threading;

namespace WindowsFormsApplication1
{
    public partial class Form1 : Form
    {
        private int xl, mx;                             //マウス移動量
        private Boolean mb = false;                     //マウス状態
        private Boolean mb2 = false;                    //マウス状態
        private Boolean bmpGetTry = true;               //画像有無ステータス
        private Point mMousePoint;                      //マウス座標
        private double expandSize = 0.9;                //球体拡大率
        private int keido = 135;                        //経度初期値(日本)
        private const int h = 500;                      //画像縦サイズ
        private const int h2 = h / 2;
        private const int w = 2 * h;
        private const int w2 = h;
        private Bitmap bmp1 = new Bitmap(h + h, h);     //昼画像
        private Bitmap bmp2 = new Bitmap(h + h, h);     //夜画像
        private Bitmap bmp3 = new Bitmap(h + h, h);     //昼夜ミックス
        private int[] y = new int[h];                   //平面→球面変換配列Y
        private int[,] x = new int[h, h];               //平面→球面変換配列X
        public struct Dtime                             //Datetime同等構造体
        {
            public int Year;
            public int Month;
            public int Day;
            public int Hour;
            public int Minute;
            public int Second;
        }
        public class TLEdata                            //2行フォーマット
        {
            public double ec;                           //離心率(Eccentricity)
            public double mk;                           //平均近点角(Mean anomaly)
            public double mm;                           //平均運動(Mean motion)
            public double sx1;
            public double sx2;
            public double sy1;
            public double sy2;
            public double sz1;
            public double sz2;
            public double tjdep;                        //元期(Epoch)
        }
        private TLEdata t1;                             //ISSのTLEデータ派生クラス

        public Form1()
        {
            InitializeComponent();
            this.MouseDown += new MouseEventHandler(Form1_MouseDown);
            this.MouseMove += new MouseEventHandler(Form1_MouseMove);
            orbital_main();
        }
        private void orbital_main()
        {
                                                        //TLE読み込み(下記NASA URLから1日以内前のTLEを抽出する)
            string tle = "http://spaceflight.nasa.gov/realdata/sightings/SSapplications/Post/JavaSSOP/orbit/ISS/SVPOST.html";
            t1 = ISStle(tle);

            string loadurlname;
            string loadfilename = @"earth_noon.jpg";    //昼画像読み込み(from Desktop)
            bmp1 = makeBmp_desktop(loadfilename);
            if (bmpGetTry == false)
            {
                loadurlname                            //(from NASA;Desktopになかった場合)
                    = "http://eoimages.gsfc.nasa.gov/images/imagerecords/73000/73801/world.topo.bathy.200409.3x5400x2700.jpg";
                bmp1 = makeBmp(loadurlname, loadfilename);
            }
            loadfilename = @"earth_night.jpg";          //夜画像読み込み(from Desktop)
            bmp2 = makeBmp_desktop(loadfilename);
            if (bmpGetTry == false)
            {
                loadurlname                            //(from NASA;Desktopになかった場合)
                    = "http://eoimages.gsfc.nasa.gov/images/imagerecords/55000/55167/earth_lights_lrg.jpg";
                bmp2 = makeBmp(loadurlname, loadfilename);
            }
            y = matY(h);                                //平面→球面投影配列準備
            x = matX(y, h);
            marble();
            //タイマーセット
            timer3.Interval = 600000;                    //昼夜画像書き換え
            timer3.Start();
            timer4.Interval = 100;                      //地球回転
            timer4.Start();
        }

---

次の記事に続く。