Real sun position


photo

Recommended Posts

Hello,

How should I code the rotation of the sun (LightWorld), in order to simulate its real position in the sky under a given latitude/longitude and date/time ? Do you have a real code maybe ? (btw, it would be a great addition to have this by default when a LightWorld is child of a geodetic node)

Thanks

Link to post

Hello,

I'm afraid we have no any code examples or demos that shows how to move LightWorld object accorrding to real Sun position.
However it is possible.

You can use geodetic pivot (this and this could be helpful) to compare real coordinates and decart coordnates.

How to submit a good bug report
---
FTP server for test scenes and user uploads:

Link to post

Hello,

Thanks to corrected code posted here https://stackoverflow.com/questions/8708048/position-of-the-sun-given-time-of-day-latitude-and-longitude . I ported it in script. Here it is for future reference, and maybe get included in a future version of the SDK :-)

void sunPosition(double & el, double& az, double lat, double lon, int year=2017, int month=12, int day=31, double hour=23, int min=59,int sec=59) {
    double pi = 3.14159265;
    double twopi = 2 * pi;
    double deg2rad = pi / 180.;

    // Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
    int month_days[13] = ( 0,31,28,31,30,31,30,31,31,30,31,30 );
    for (int i=0; i<month; i++)
        day += month_days[i];
    int leapdays = (year % 4) == 0 && ((year % 400) == 0 || (year % 100) != 0) && day >= 60 && !(month==2 && day==60);
    if (leapdays>0) day++;

    // Get Julian date - 2400000
    hour += min / 60. + sec / 3600.; // hour plus fraction
    double delta = year - 1949.;
    double leap = int(delta / 4.); // former leapyears
    double jd = 32916.5 + delta * 365 + leap + day + hour / 24.;

    // The input to the Atronomer's almanach is the difference between
    // the Julian date and JD 2451545.0 (noon, 1 January 2000)
    double time = jd - 51545.;

    // Ecliptic coordinates

    // Mean longitude
    double mnlong = 280.460 + .9856474 * time;
    mnlong = mnlong % 360;
    if (mnlong<0) mnlong+= 360;

    // Mean anomaly
    double mnanom = 357.528 + .9856003 * time;
    mnanom = mnanom %360;
    if (mnanom<0) mnanom+= 360;
    mnanom *=deg2rad;

    // Ecliptic longitude and obliquity of ecliptic
    double eclong = mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom);
    eclong = eclong % 360;
    if (eclong<0) eclong+= 360;
    double oblqec = 23.439 - 0.0000004 * time;
    eclong *=deg2rad;
    oblqec *=deg2rad;

    // Celestial coordinates
    // Right ascension and declination
    double num = cos(oblqec) * sin(eclong);
    double den = cos(eclong);
    double ra = atan(num / den);
    if (den<0) ra+= pi;
    if (den >= 0 && num < 0) ra+= twopi;
    double dec = asin(sin(oblqec) * sin(eclong));

    // Local coordinates
    // Greenwich mean sidereal time
    double gmst = 6.697375 + .0657098242 * time + hour;
    gmst = gmst % 24;
    if (gmst<0) gmst+= 24.;

    // Local mean sidereal time
    double lmst = gmst + lon / 15.;
    lmst = lmst %24.;
    if (lmst<0) lmst += 24.;
    lmst = lmst * 15. * deg2rad;

    // Hour angle
    double ha = lmst - ra;
    if (ha<-pi) ha+= twopi;
    if (ha>pi) ha-= twopi;

    // Latitude to radians
    lat = lat * deg2rad;

    // Azimuth and elevation
    el = asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha));
    az = asin(-cos(dec) * sin(ha) / cos(el));

    // For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
    int cosAzPos = (0 <= sin(dec) - sin(el) * sin(lat));
    int sinAzNeg = (sin(az) < 0);
    if (cosAzPos && sinAzNeg) az +=twopi;
    if (!cosAzPos) az=  pi - az;

    // return el and az
    el = el / deg2rad;
    az = az / deg2rad;
    lat = lat / deg2rad;
}


#include <core/unigine.h>

Node sun = engine.editor.getNodeByName("sun");
if (sun!=NULL) {
	double lat = 43.15;
	double lon = 5.1;
	double el, az;
	sunPosition(el,az, lat, lon, 2017,6,21, 12,00,00);
	log.message("el=%d , az=%d\n", el,az);
	sun.setRotation( quat(90,270,270) * quat(az,0,0) * quat(0,90,0) * quat(el,0,0) * quat(90,0,0) );
}

The code is quite costly, so you certainly don't want to run it every frame, but probably once every minutes or so for a "real time" sun movement.

:)

  • Like 3
Link to post