Arma 3 Night Sky Planets

/*
    Place the file in a mission folder and execute it for example for player's init field:
        0 = [] execvm "fn_drawPlanets.sqf";

    Calculations based on https://aa.quae.nl/en/reken/hemelpositie.html
*/
bis_fnc_drawPanets_timeRefresh = 0;
draw3DPlanets = {

    if (sunormoon < 1) then {
        //--- Get world coordinates
        _lat = -getnumber (configfile >> "CfgWorlds" >> worldname >> "latitude");   //latitude = -35.097;
        _lon = getnumber (configfile >> "CfgWorlds" >> worldname >> "longitude");   //longitude = 16.482;
        _timeZone = -2;
        _time = daytime;

        if (time > bis_fnc_drawPanets_timeRefresh) then {
            //--- orbital params at 01/01/2000
            _planets = [
                //Name      a       e       i       AP      LAN     M0      Color       Size
                ["Earth",   1.00000,    0.01671,    0.000,      288.064,    174.873,    357.529,    [0,0,0],    0.0],
                ["Mercury", 0.38710,    0.20563,    7.005,      29.125,     48.331,     174.795,    [255,243,239],  0.2],
                ["Venus",   0.72333,    0.00677,    3.395,      54.884,     76.680,     50.416,     [253,254,255],  0.25],
                ["Mars",    1.52368,    0.09340,    1.850,      286.502,    49.558,     19.373,     [245,151,118],  0.25],
                ["Jupiter", 5.20260,    0.04849,    1.303,      273.687,    100.464,    20.020,     [249,250,232],  0.25]
                //--- Not visible without a telescope
                //["Saturn",    9.55491,    0.05551,    2.489,      339.391,    113.666,    317.021,    [253,254,209],  0.2],
                //["Uranus",    19.21845,   0.04630,    0.773,      98.999,     74.006,     141.050,    [188,253,250],  0.2],
                //["Neptune",   30.11039,   0.00899,    1.770,      276.340,    131.784,    256.225,    [190,226,242],  0.2]
            ];
            draw3DPlanets_data = [];

            //--- Get number of days since 01/01/2000
            _date = date;//[2004,1,1,0,0];
            _year = _date select 0;
            _day = ((datetonumber _date) * 365);
            //_time = 1;
            _diffYear = _year - 2000;
            _days = abs((_diffYear * 365) + floor(_diffYear / 4) + _day); // ToDo: Correct leap year calculation

            _hXE = 0;
            _hYE = 0;
            _hZE = 0;
            _ME = 0;
            {
                _name = _x select 0;
                _a = _x select 1;
                _e = _x select 2;
                _i = _x select 3;
                _AP = _x select 4;
                _LAN = _x select 5;
                _M0 = _x select 6;
                _color = _x select 7;
                _size = _x select 8;

                //--- Get mean anomaly
                _n = 0.9856076686 / (_a^(3/2));
                _M = (_M0 + _n * _days) % 360;
                _TA = _M + 2 * _e * sin(_M) + 1.25 * _e^2 * sin(2 * _M);

                //--- Get eccentric anomaly
                _EA = _M;
                _f = 0;
                _ratio = 0;
                for "_i" from 0 to 10 do {
                    _ff = 1 - _e * cos(_f);
                    _f = _EA - _e * sin(_EA) - _M;
                    _EA = _EA - _ratio;
                    _ratio = _ratio + _f / _ff;
                };

                //--- Get true anomaly
                _TA = 2 * atan(sqrt((1 + _e) / (1 - _e)) * tan(_EA / 2));

                //--- Get distance from the sun
                _r = (_a * (1 - _e^2)) / (1 + _e * cos(_TA));

                //--- Get heliocentric elliptical coordinates
                _hX = _r * ((cos(_LAN) * cos(_AP + _TA)) - sin(_LAN) * cos(_i) * sin(_AP + _TA));
                _hY = _r * ((sin(_LAN) * cos(_AP + _TA)) + cos(_LAN) * cos(_i) * sin(_AP + _TA));
                _hZ = _r * sin(_i) * sin(_AP + _TA);
                if (_foreachindex == 0) then {
                    _hXE = _hX;
                    _hYE = _hY;
                    _hZE = _hZ;
                    _ME = _M;
                } else {
                    //--- Get rectangular geocentric ecliptical coordinates
                    _posX = _hX - _hXE;
                    _posY = _hY - _hYE;
                    _posZ = _hZ - _hZE;

                    //--- Get geocentric ecliptical longitude and latitude
                    _lambda = _posY atan2 _posX;
                    _beta = asin(_posZ / sqrt(_posX^2 + _posY^2 + _posZ^2));

                    //--- The right ascension and declination
                    _epsilon = 23.4397;
                    _alpha = (sin(_lambda) * cos(_epsilon) - tan(_beta) * sin(_epsilon)) atan2 (cos(_lambda));
                    _delta = asin(sin(_beta) * cos(_epsilon) + cos(_beta) * sin(_epsilon) * sin(_lambda));

                    //--- The sidereal time
                    //_theta = (_ME + 102.937 + 15 * (_time + _timeZone) + _lon) % 360;

                    //--- Hour angle
                    //_H = _theta - _alpha;

                    //--- The height and azimuth
                    //_dirH = (sin(_H)) atan2 (cos(_H) * sin(_lat) - tan(_delta) * cos(_lat));
                    //_dirV = asin(sin(_lat) * sin(_delta) + cos(_lat) * cos(_delta) * cos(_H));

                    draw3DPlanets_data pushback [_name,_ME,_alpha,_delta,_color,_size];
                };
            } foreach _planets;
            bis_fnc_drawPanets_timeRefresh = time + 1;
        };

        {

            _name = _x select 0;
            _ME = _x select 1;
            _alpha = _x select 2;
            _delta = _x select 3;
            _color = _x select 4;
            _size = _x select 5;
            _color = [
                (_color select 0) / 255,
                (_color select 1) / 255,
                (_color select 2) / 255,
                1
            ];

            if (cameraview == "external") then {_name = "";};

            _theta = (_ME + 102.937 + 15 * (_time + _timeZone) + _lon) % 360;
            _H = _theta - _alpha;
            _dirH = (sin(_H)) atan2 (cos(_H) * sin(_lat) - tan(_delta) * cos(_lat));
            _dirV = asin(sin(_lat) * sin(_delta) + cos(_lat) * cos(_delta) * cos(_H));

            _dirH = 282 - _dirH; //--- Magic constant

            if (_dirV > 5) then {
                _dis = viewdistance;
                _pos = [
                    (position player select 0) + (cos _dirH * _dis),
                    (position player select 1) + (sin _dirH * _dis),
                    (position player select 2) + (tan _dirV * _dis)
                ];
                drawicon3d ["\a3\Data_f\light_flare_ca.paa",_color,_pos,_size,_size,0,_name,0];
            };
        } foreach draw3DPlanets_data;
    };
};

if (isnil "draw3DPlanetsHandler") then {
    draw3DPlanetsHandler = addmissioneventhandler ["draw3D",{[] call draw3DPlanets}];
};