Rendering the sky like Stellarium with JPARSEC

Rendering the sky like Stellarium with JPARSEC

There is an old project called Stellarium for Java. It is a port of the old 0.8 version of Stellarium, which is a very popular planetarium due to its high quality, realistic renderings. This project uses the JOGL library to offer a similar experience to the original using Java, but unfortunately it was later abandoned.

I have recently analized the original code to look how this port renders the sky, in particular how it computes the atmosphere brightness and the sky colors considering the position of the Sun or the Moon. I have taken just this code and connected it to the rendering process used in JPARSEC, so that similar renderings can be done (pixel by pixel, so slower) without using the GPU. For this task I have updated JPARSEC with a few improvements (availables from the JPARSEC repository at bitbucket) to improve the horizon rendering with terrain, adding two more more textures.

The idea is to render the sky with JPARSEC using a given texture for the horizon and a transparent image, and reuse the SkyRendering object to compute, pixel by pixel, the illumination and color of the background atmosphere using the formulae from Stellarium. With this you only have to draw the sky image on top of the atmosphere image for the final image. So here are the great results…




 

I leave to you the task of analyzing the code to try to understand the underlying maths and physics, I prefer not to describe it since I'm not an expert, I have just borrowed the code. Here is the rather long Java class that can be used for such renderings in case you want to play with it.

StellariumAtmosphere.java
import static java.lang.StrictMath.PI;
import static java.lang.StrictMath.abs;
import static java.lang.StrictMath.acos;
import static java.lang.StrictMath.asin;
import static java.lang.StrictMath.cos;
import static java.lang.StrictMath.exp;
import static java.lang.StrictMath.log;
import static java.lang.StrictMath.log10;
import static java.lang.StrictMath.pow;
import static java.lang.StrictMath.sin;
import static java.lang.StrictMath.tan;
 
import java.awt.Graphics2D;
import java.awt.image.BufferedImage;
 
import jparsec.astronomy.TelescopeElement;
import jparsec.astronomy.CoordinateSystem.COORDINATE_SYSTEM;
import jparsec.ephem.Ephem;
import jparsec.ephem.EphemerisElement;
import jparsec.ephem.Functions;
import jparsec.ephem.Target.TARGET;
import jparsec.ephem.planets.EphemElement;
import jparsec.graph.chartRendering.AWTGraphics;
import jparsec.graph.chartRendering.Graphics;
import jparsec.graph.chartRendering.PlanetRenderElement;
import jparsec.graph.chartRendering.SkyRenderElement;
import jparsec.graph.chartRendering.Projection.PROJECTION;
import jparsec.graph.chartRendering.RenderPlanet;
import jparsec.graph.chartRendering.SkyRenderElement.HORIZON_TEXTURE;
import jparsec.graph.chartRendering.SkyRenderElement.LEYEND_POSITION;
import jparsec.graph.chartRendering.SkyRenderElement.MILKY_WAY_TEXTURE;
import jparsec.graph.chartRendering.SkyRenderElement.REALISTIC_STARS;
import jparsec.graph.chartRendering.frame.SkyRendering;
import jparsec.io.ConsoleReport;
import jparsec.io.image.Picture;
import jparsec.math.Constant;
import jparsec.math.FastMath;
import jparsec.observer.City;
import jparsec.observer.LocationElement;
import jparsec.observer.ObserverElement;
import jparsec.time.AstroDate;
import jparsec.time.TimeElement;
 
public class StellariumAtmosphere {
 
	public static SkyRendering render;
 
	// Returns an image showing the sky for a given fixed set of properties: observer, appearance, and so on
	public static BufferedImage getStarsImage() throws Exception {        	
    	// Parameters for the rendering: geographical location and elevation, time zone,
    	// camera pointing direction and inclination, sensor size in pixels, field of view, ...
    	String locName = "Madrid";
    	double cameraAcimut = 90 * Constant.DEG_TO_RAD; // 0 = North
    	double cameraElevation = 50 * Constant.DEG_TO_RAD; // 0 = horizon, 90 = zenith
    	double cameraInclination = 0 * Constant.DEG_TO_RAD;
    	int cameraPixelX = 1200, cameraPixelY = 800;
    	double fov = 120 * Constant.DEG_TO_RAD;
    	TimeElement.SCALE ts = TimeElement.SCALE.LOCAL_TIME;
    	AstroDate astro = new AstroDate(2018, 2, 19, 10, 0, 0);
    	TimeElement date = new TimeElement(astro, ts);
 
    	EphemerisElement eph = new EphemerisElement(TARGET.SUN,
    		EphemerisElement.COORDINATES_TYPE.APPARENT, EphemerisElement.EQUINOX_OF_DATE,
    		EphemerisElement.TOPOCENTRIC, EphemerisElement.REDUCTION_METHOD.IAU_2006,
    		EphemerisElement.FRAME.DYNAMICAL_EQUINOX_J2000,
    		EphemerisElement.ALGORITHM.MOSHIER);
    	eph.optimizeForSpeed();
    	ObserverElement observer = ObserverElement.parseCity(City.findCity(locName));
 
    	SkyRenderElement sky = getSky(cameraAcimut, cameraElevation, cameraPixelX, cameraPixelY, fov, cameraInclination);
 
	sky.drawMilkyWayContoursWithTextures = MILKY_WAY_TEXTURE.OPTICAL;
	EphemElement ephemSun = Ephem.getEphemeris(date, observer, eph, true);
	if (ephemSun.elevation > -5 * Constant.DEG_TO_RAD) sky.drawMilkyWayContoursWithTextures = MILKY_WAY_TEXTURE.NO_TEXTURE;
	sky.drawMilkyWayContours = sky.drawMilkyWayContoursWithTextures == MILKY_WAY_TEXTURE.NO_TEXTURE ? false : true;
	int yMargin = 0;
 
	// These lines are used to change the vertical position of the horizon, to create panoramas
//		sky.height = sky.width/2+150;
//		yMargin = sky.height/2-100;
 
    	render = new SkyRendering(date, observer, eph, sky, locName, yMargin);  
    	// JPARSEC supports transparency using these two lines
    	AWTGraphics.setBufferedImageType(BufferedImage.TYPE_INT_ARGB);
    	BufferedImage out = render.createTransparentBufferedImage();
    	return out;
    }
 
    public static BufferedImage getSkyBrightness() throws Exception {
	EphemElement ephemSun = render.getRenderSkyObject().calcPlanet(TARGET.SUN, true, true);
	EphemElement ephemMoon = render.getRenderSkyObject().calcPlanet(TARGET.Moon, true, false);
 
	ConsoleReport.fullEphemReportToConsole(ephemSun);
 
	ToneReproductor eye = new ToneReproductor();
	StellariumAtmosphere atmosphere = new StellariumAtmosphere();
	float atmBr = 0.0001f, turbidity = 1f, atmI = 1000f;
	if (ephemSun.elevation < 0) {
		double elev = Math.min(18, Math.abs(ephemSun.elevation * Constant.RAD_TO_DEG));
		atmI = (float) (1000 + 5000 * elev / 15.0);
		//turbidity = (float) (2 - 1 * Math.abs(ephemSun.elevation * Constant.RAD_TO_DEG) / 15.0);
	}
        eye.setWorldAdaptationLuminance(3.75f + atmBr * 40000.f);
        atmosphere.atmIntensity = 255 * atmBr * atmI;			
 
        int year = render.getTimeObject().astroDate.getYear(), month = render.getTimeObject().astroDate.getMonth();
        double latitude = render.getObserverObject().getLatitudeDeg(), altitude = render.getObserverObject().getHeight(), 
        		temperature = render.getObserverObject().getTemperature(), relativeHumidity = render.getObserverObject().getHumidity();
        double moonPhase = Math.PI - ephemMoon.elongation;
 
        LocationElement sunLoc = new LocationElement(ephemSun.azimuth, ephemSun.elevation, 1);
        LocationElement moonLoc = new LocationElement(ephemMoon.azimuth, ephemMoon.elevation, 1);
        LocationElement zenLoc = new LocationElement(0, Constant.PI_OVER_TWO, 1);
	double moonZenD = LocationElement.getAngularDistance(moonLoc, zenLoc);
	double sunZenD = LocationElement.getAngularDistance(sunLoc, zenLoc);
 
	atmosphere.skyLight.setParams((float) sunZenD, turbidity);
	atmosphere.skyBright.setLoc(Math.toRadians(latitude), altitude, temperature, relativeHumidity);
	atmosphere.skyBright.setSunMoon(Math.cos(moonZenD), Math.cos(sunZenD));
 
	atmosphere.skyBright.setDate(year, month, moonPhase);
 
	int w = render.getRenderSkyObject().render.width;
	int h = render.getRenderSkyObject().render.height;
	BufferedImage out = new BufferedImage(w, h, BufferedImage.TYPE_INT_RGB);
	int black = (new java.awt.Color(0, 0, 0)).getRGB();
	for (int x = 0; x < w; ++x) {
	for (int y = 0; y < h; ++y) {
		LocationElement skyLoc = render.getRenderSkyObject().getSkyLocation(x, y);
		if (skyLoc.getLatitude() <= 0) {
			out.setRGB(x, y, black);
			continue;
		}
		double sunD = LocationElement.getApproximateAngularDistance(skyLoc, sunLoc); 
		double moonD = LocationElement.getApproximateAngularDistance(skyLoc, moonLoc); 
		double zenD = LocationElement.getApproximateAngularDistance(skyLoc, zenLoc); 
 
		double cosDistSun = FastMath.cos(sunD), cosDistZen = FastMath.cos(zenD);
		float color[] = atmosphere.skyLight.get_xyY_valuev(cosDistSun, 1.0 / cosDistZen);
		color[2] = (float) atmosphere.skyBright.getLuminance(FastMath.cos(moonD), 
			cosDistSun, cosDistZen);
 
		eye.xyYToRGB(color);
 
		int r = (int) (atmosphere.atmIntensity * color[0]);
		int g = (int) (atmosphere.atmIntensity * color[1]);
		int b = (int) (atmosphere.atmIntensity * color[2]);
		if (r > 255) r = 255;
		if (g > 255) g = 255;
		if (b > 255) b = 255;
		if (r < 0) r = 0;
		if (g < 0) g = 0;
		if (b < 0) b = 0;
		int rgb = (new java.awt.Color(r, g, b)).getRGB();
		out.setRGB(x, y, rgb);
	}
	}
	return out;
    }
 
    // Returns the sky render object containing all properties to define the rendering appearance
    public static SkyRenderElement getSky(double acimut, double elevation, int width, int height, double fov, double incl) {
    	PlanetRenderElement planet = new PlanetRenderElement(false, true, true, false);
	SkyRenderElement sky = new SkyRenderElement(COORDINATE_SYSTEM.HORIZONTAL,
		PROJECTION.STEREOGRAPHICAL, acimut, elevation, width, height, planet, 
		TelescopeElement.OBJECTIVE_50mm_f1_4);
	sky.telescope.ocular.focalLength = TelescopeElement.getOcularFocalLengthForCertainField(fov, sky.telescope);
	sky.drawObjectsLimitingMagnitude = 3f;        
	sky.drawStarsLimitingMagnitude = 5.0f;
	sky.drawPlanetsMoonSun = true;
	sky.drawSkyCorrectingLocalHorizon = true;        
	sky.drawSkyBelowHorizon = false;        
	sky.drawFastLabels = SkyRenderElement.SUPERIMPOSED_LABELS.AVOID_SUPERIMPOSING_ACCURATE;        
	sky.drawFastLabelsInWideFields = false;        
	sky.drawConstellationLimits = false;
	sky.drawLeyend = LEYEND_POSITION.TOP;
	sky.drawExternalGrid = false;
	sky.drawCoordinateGrid = false;
	sky.drawDeepSkyObjectsAllMessierAndCaldwell = false;
	sky.drawStarsRealistic = REALISTIC_STARS.SPIKED;
	sky.drawMeteorShowers = false;
	sky.drawNebulaeContours = false;
	sky.drawComets = sky.drawAsteroids = false;
	sky.drawConstellationNames = false;
	sky.poleAngle = (float) incl;
	sky.drawStarsLabelsLimitingMagnitude = 1.1f;
	sky.drawPlanetsMoonSun = true;
	sky.drawPlanetsLabels = true;
	sky.drawIcons = false;
	sky.drawSpaceProbes = false;
	sky.drawDeepSkyObjects = false;
	sky.drawConstellationNames = true;
	//sky.drawConstellationNamesType = CONSTELLATION_NAME.SPANISH;
	sky.drawConstellationNamesFont = Graphics.FONT.SANS_SERIF_ITALIC_22;
	sky.drawPlanetsNamesFont = Graphics.FONT.SANS_SERIF_BOLD_16;
	sky.drawStarsNamesFont = Graphics.FONT.DIALOG_PLAIN_12;
	sky.drawDeepSkyObjectsTextures = false;
	sky.drawCoordinateGridEclipticLabels = false;
	sky.planetRender.textures = false;
	sky.drawMilkyWayContoursWithTextures = MILKY_WAY_TEXTURE.OPTICAL;
	sky.drawMilkyWayContours = sky.drawMilkyWayContoursWithTextures == MILKY_WAY_TEXTURE.NO_TEXTURE ? false : true;
	sky.drawHorizonTexture = HORIZON_TEXTURE.VELETA_30m;
	sky.drawLeyend = LEYEND_POSITION.NO_LEYEND;
	sky.planetRender.highQuality = true;
	RenderPlanet.MAXIMUM_TEXTURE_QUALITY_FACTOR = 4f;
 
	sky.setColorMode(SkyRenderElement.COLOR_MODE.BLACK_BACKGROUND);
	sky.background = Functions.getColor(0, 0, 0, 0);
	return sky;
    }
 
    private Skylight skyLight = new Skylight();
 
    private SkyBright skyBright = new SkyBright();
 
    double worldAdaptationLuminance;
 
    double milkywayAdaptationLuminance;
 
    private float atmIntensity;
 
    public static void main(String args[]) {
    	try {
        	//Translate.setDefaultLanguage(LANGUAGE.SPANISH);
 
    		BufferedImage sky = getStarsImage();
    		Picture pic = new Picture(getSkyBrightness());
    		Graphics2D g = pic.getImage().createGraphics();
    		AWTGraphics.enableAntialiasing(g);
    		g.drawImage(sky, 0, 0, null);
    		g.dispose();
    		pic.show("");
    		pic.write("/home/alonso/test.png");
 
    	} catch (Exception exc) {
    		exc.printStackTrace();
    	}
    }
}
 
class ToneReproductor {
 
   public ToneReproductor() {
        lda = 50.f;
        lwa = 40000.f;
        maxDL = 100.f;
        gamma = 2.3f;
 
        // Update alpha_da and beta_da values
        double log10Lwa = log10(lwa);
        alphaWa = 0.4f * log10Lwa + 1.519f;
        betaWa = -0.4f * log10Lwa * log10Lwa + 0.218f * log10Lwa + 6.1642f;
 
        setDisplayAdaptationLuminance(lda);
        setWorldAdaptationLuminance(lwa);
    }
 
    /**
     * Set the eye adaptation luminance for the display and precompute what can be
     * Usual luminance range is 1-100 cd/m^2 for a CRT screen
     *
     * @param _Lda
     */
    void setDisplayAdaptationLuminance(double _Lda) {
        lda = _Lda;
 
        // Update alpha_da and beta_da values
        double log10Lda = log10(lda);
        alphaDa = 0.4f * log10Lda + 1.519f;
        betaDa = -0.4f * log10Lda * log10Lda + 0.218f * log10Lda + 6.1642f;
 
        // Update terms
        alphaWaOverAlphaDa = alphaWa / alphaDa;
        term2 = pow(10.f, (betaWa - betaDa) / alphaDa) / (PI * 0.0001f);
    }
 
    /**
     * Set the eye adaptation luminance for the world and precompute what can be
     *
     * @param _Lwa
     */
    public void setWorldAdaptationLuminance(double _Lwa) {
        lwa = _Lwa;
 
        // Update alpha_da and beta_da values
        double log10Lwa = log10(lwa);
        alphaWa = 0.4f * log10Lwa + 1.519f;
        betaWa = -0.4f * log10Lwa * log10Lwa + 0.218f * log10Lwa + 6.1642f;
 
        // Update terms
        alphaWaOverAlphaDa = alphaWa / alphaDa;
        term2 = pow(10.f, (betaWa - betaDa) / alphaDa) / (PI * 0.0001f);
 
    }
 
    /**
     * Convert from xyY color system to RGB according to the adaptation
     * The Y component is in cd/m^2
     */
    public void xyYToRGB(float[] color) {
        // TODO: Fred the parameter should an SColor object
        // 1. Hue conversion
        float log10Y = (float) log10(color[2]);
        // if log10Y>0.6, photopic vision only (with the cones, colors are seen)
        // else scotopic vision if log10Y<-2 (with the rods, no colors, everything blue),
        // else mesopic vision (with rods and cones, transition state)
        if (log10Y < 0.6) {
            // Compute s, ratio between scotopic and photopic vision
            float s = 0.f;
            if (log10Y > -2.f) {
                float op = (log10Y + 2.f) / 2.6f;
                s = 3.f * op * op - 2 * op * op * op;
            }
 
            // Do the blue shift for scotopic vision simulation (night vision) [3]
            // The "night blue" is x,y(0.25, 0.25)
            color[0] = (1.f - s) * 0.25f + s * color[0];// Add scotopic + photopic components
            color[1] = (1.f - s) * 0.25f + s * color[1];// Add scotopic + photopic components
 
            // Take into account the scotopic luminance approximated by V [3] [4]
            double V = color[2] * (1.33f * (1.f + color[1] / color[0] + color[0] * (1.f - color[0] - color[1])) - 1.68f);
            color[2] = (float) (0.4468f * (1.f - s) * V + s * color[2]);
        }
 
        // 2. Adapt the luminance value and scale it to fit in the RGB range [2]
        color[2] = (float) pow(adaptLuminance(color[2]) / maxDL, 1.d / gamma);
 
        // Convert from xyY to XZY
        double X = color[0] * color[2] / color[1];
        double Y = color[2];
        double Z = (1.f - color[0] - color[1]) * color[2] / color[1];
 
        // Use a XYZ to Adobe RGB (1998) matrix which uses a D65 reference white
        color[0] = (float) (2.04148f * X - 0.564977f * Y - 0.344713f * Z);
        color[1] = (float) (-0.969258f * X + 1.87599f * Y + 0.0415557f * Z);
        color[2] = (float) (0.0134455f * X - 0.118373f * Y + 1.01527f * Z);
    }
 
    /**
     * Set the maximum display luminance : default value = 100 cd/m^2
     * This value is used to scale the RGB range
     */
    void set_max_display_luminance(float _maxdL) {
        maxDL = _maxdL;
    }
 
    /**
     * Set the display gamma : default value = 2.3
     */
    void setDisplayGamma(float _gamma) {
        gamma = _gamma;
    }
 
    /**
     * Return adapted luminance from world to display
     */
    public double adaptLuminance(double worldLluminance) {
        return pow(worldLluminance * PI * 0.0001d, alphaWaOverAlphaDa) * term2;
    }
 
    private double lda;// Display luminance adaptation (in cd/m^2)
 
    double lwa;// World   luminance adaptation (in cd/m^2)
 
    double maxDL;// Display maximum luminance (in cd/m^2)
 
    double gamma;// Screen gamma value
 
    // Precomputed variables
    double alphaDa;
 
    double betaDa;
 
    double alphaWa;
 
    double betaWa;
 
    double alphaWaOverAlphaDa;
 
    double term2;
 
}
 
class Skylight {
    public class SkylightStruct {
        double zenith_angle;// zenith_angle : angular distance to the zenith in radian
 
        double dist_sun;// dist_sun     : angular distance to the sun in radian
 
        double color[] = new double[3];// 3 component color, can be RGB or CIE color system
    }
 
    public static class SkylightStruct2 {
        public double pos[] = new double[3];// Vector to the position (vertical = pos[2])
 
        public float color[] = new float[3];// 3 component color, can be RGB or CIE color system
    }
 
    public Skylight() {
    }
 
    void setParams(float _sun_zenith_angle, float _turbidity) {
        // Set the two Main variables
        thetas = _sun_zenith_angle;
        T = _turbidity;
 
        // Precomputation of the distribution coefficients and zenith luminances/color
        compute_zenith_luminance();
        compute_zenith_color();
        compute_luminance_distribution_coefs();
        compute_color_distribution_coefs();
 
        // Precompute everything possible to increase the get_CIE_value() function speed
        double cos_thetas = cos(thetas);
        termX = (float) (zenithColorX / ((1.f + aX * exp(bX)) * (1.f + cX * exp(dX * thetas) + eX * cos_thetas * cos_thetas)));
        termY = (float) (zenith_color_y / ((1.f + Ay * exp(By)) * (1.f + Cy * exp(Dy * thetas) + Ey * cos_thetas * cos_thetas)));
        term_Y = (float) (zenith_luminance / ((1.f + AY * exp(BY)) * (1.f + CY * exp(DY * thetas) + EY * cos_thetas * cos_thetas)));
 
    }
 
    public void setParamsV(float[] _sun_pos, float _turbidity) {
        // Store sun position
        sunPos[0] = _sun_pos[0];
        sunPos[1] = _sun_pos[1];
        sunPos[2] = _sun_pos[2];
 
        // Set the two Main variables
        thetas = (float) (PI / 2 - asin(sunPos[2]));
        T = _turbidity;
 
        // Precomputation of the distribution coefficients and zenith luminances/color
        compute_zenith_luminance();
        compute_zenith_color();
        compute_luminance_distribution_coefs();
        compute_color_distribution_coefs();
 
        // Precompute everything possible to increase the get_CIE_value() function speed
        double cos_thetas = sunPos[2];
        termX = (float) (zenithColorX / ((1.f + aX * exp(bX)) * (1.f + cX * exp(dX * thetas) + eX * cos_thetas * cos_thetas)));
        termY = (float) (zenith_color_y / ((1.f + Ay * exp(By)) * (1.f + Cy * exp(Dy * thetas) + Ey * cos_thetas * cos_thetas)));
        term_Y = (float) (zenith_luminance / ((1.f + AY * exp(BY)) * (1.f + CY * exp(DY * thetas) + EY * cos_thetas * cos_thetas)));
    }
 
    /**
     * Compute CIE luminance for zenith in cd/m^2
     */
    void compute_zenith_luminance() {
        zenith_luminance = (float) (1000.f * ((4.0453f * T - 4.9710f) * tan((0.4444f - T / 120.f) * (PI - 2.f * thetas)) -
                0.2155f * T + 2.4192f));
        if (zenith_luminance <= 0.f) zenith_luminance = 0.00000000001f;
    }
 
    /**
     * Compute CIE x and y color components
     */
    void compute_zenith_color() {
        thetas2 = thetas * thetas;
        thetas3 = thetas2 * thetas;
        T2 = T * T;
 
        zenithColorX = (0.00166f * thetas3 - 0.00375f * thetas2 + 0.00209f * thetas) * T2 +
                (-0.02903f * thetas3 + 0.06377f * thetas2 - 0.03202f * thetas + 0.00394f) * T +
                (0.11693f * thetas3 - 0.21196f * thetas2 + 0.06052f * thetas + 0.25886f);
 
        zenith_color_y = (0.00275f * thetas3 - 0.00610f * thetas2 + 0.00317f * thetas) * T2 +
                (-0.04214f * thetas3 + 0.08970f * thetas2 - 0.04153f * thetas + 0.00516f) * T +
                (0.15346f * thetas3 - 0.26756f * thetas2 + 0.06670f * thetas + 0.26688f);
 
    }
 
    /**
     * Compute the luminance distribution coefficients
     */
    void compute_luminance_distribution_coefs() {
        AY = 0.1787f * T - 1.4630f;
        BY = -0.3554f * T + 0.4275f;
        CY = -0.0227f * T + 5.3251f;
        DY = 0.1206f * T - 2.5771f;
        EY = -0.0670f * T + 0.3703f;
    }
 
    /**
     * Compute the color distribution coefficients
     */
    void compute_color_distribution_coefs() {
        aX = -0.0193f * T - 0.2592f;
        bX = -0.0665f * T + 0.0008f;
        cX = -0.0004f * T + 0.2125f;
        dX = -0.0641f * T - 0.8989f;
        eX = -0.0033f * T + 0.0452f;
 
        Ay = -0.0167f * T - 0.2608f;
        By = -0.0950f * T + 0.0092f;
        Cy = -0.0079f * T + 0.2102f;
        Dy = -0.0441f * T - 1.6537f;
        Ey = -0.0109f * T + 0.0529f;
    }
 
    /**
     * Compute the sky color at the given position in the CIE color system and store it in p.color
     * p.color[0] is CIE x color component
     * p.color[1] is CIE y color component
     * p.color[2] is CIE Y color component (luminance)
     */
    void get_xyY_value(SkylightStruct p) {
        double cos_dist_sun = cos(p.dist_sun);
        double one_over_cos_zenith_angle = 1.f / cos(p.zenith_angle);
        p.color[0] = termX * (1.f + aX * exp(bX * one_over_cos_zenith_angle)) * (1.f + cX * exp(dX * p.dist_sun) +
                eX * cos_dist_sun * cos_dist_sun);
        p.color[1] = termY * (1.f + Ay * exp(By * one_over_cos_zenith_angle)) * (1.f + Cy * exp(Dy * p.dist_sun) +
                Ey * cos_dist_sun * cos_dist_sun);
        p.color[2] = term_Y * (1.f + AY * exp(BY * one_over_cos_zenith_angle)) * (1.f + CY * exp(DY * p.dist_sun) +
                EY * cos_dist_sun * cos_dist_sun);
    }
 
    /**
     * Compute the sky color at the given position in the CIE color system and store it in p.color
     * p.color[0] is CIE x color component
     * p.color[1] is CIE y color component
     * p.color[2] is CIE Y color component (luminance)
     */
    public void get_xyY_valuev(SkylightStruct2 p) {
        //	if (p.pos[2]<0.)
        //	{
        //		p.color[0] = 0.25;
        //		p.color[1] = 0.25;
        //		p.color[2] = 0;
        //		return;
        //	}
 
        double cosDistSun = sunPos[0] * (p.pos[0]) + sunPos[1] * (p.pos[1]) + sunPos[2] * (p.pos[2]) - 0.0000001f;
        double oneOverCosZenithAngle = 1.f / p.pos[2];
        float distSun = (float) acos(cosDistSun);
 
        p.color[0] = (float) (termX * (1.f + aX * exp(bX * oneOverCosZenithAngle)) * (1.f + cX * exp(dX * distSun) +
                eX * cosDistSun * cosDistSun));
        p.color[1] = (float) (termY * (1.f + Ay * exp(By * oneOverCosZenithAngle)) * (1.f + Cy * exp(Dy * distSun) +
                Ey * cosDistSun * cosDistSun));
        p.color[2] = (float) (term_Y * (1.f + AY * exp(BY * oneOverCosZenithAngle)) * (1.f + CY * exp(DY * distSun) +
                EY * cosDistSun * cosDistSun));
 
        if (p.color[2] < 0 || p.color[0] < 0 || p.color[1] < 0) {
            p.color[0] = 0.25f;
            p.color[1] = 0.25f;
            p.color[2] = 0;
        }
    }
 
    /**
     * Compute the sky color at the given position in the CIE color system and store it in p.color
     * p.color[0] is CIE x color component
     * p.color[1] is CIE y color component
     * p.color[2] is CIE Y color component (luminance)
     */
    public float[] get_xyY_valuev(double cosDistSun, double oneOverCosZenithAngle) {
        //	if (p.pos[2]<0.)
        //	{
        //		p.color[0] = 0.25;
        //		p.color[1] = 0.25;
        //		p.color[2] = 0;
        //		return;
        //	}
 
        float distSun = (float) FastMath.acos(cosDistSun);
 
        float color[] = new float[3];
        color[0] = (float) (termX * (1.f + aX * exp(bX * oneOverCosZenithAngle)) * (1.f + cX * exp(dX * distSun) +
                eX * cosDistSun * cosDistSun));
        color[1] = (float) (termY * (1.f + Ay * exp(By * oneOverCosZenithAngle)) * (1.f + Cy * exp(Dy * distSun) +
                Ey * cosDistSun * cosDistSun));
        color[2] = (float) (term_Y * (1.f + AY * exp(BY * oneOverCosZenithAngle)) * (1.f + CY * exp(DY * distSun) +
                EY * cosDistSun * cosDistSun));
 
        if (color[2] < 0 || color[0] < 0 || color[1] < 0) {
            color[0] = 0.25f;
            color[1] = 0.25f;
            color[2] = 0;
        }
        return color;
    }
 
    /**
     * Return the current zenith color in xyY color system
     */
    void get_zenith_color(double[] v) {
        v[0] = zenithColorX;
        v[1] = zenith_color_y;
        v[2] = zenith_luminance;
    }
 
    private float thetas;// angular distance between the zenith and the sun in radian
 
    float T;// Turbidity : i.e. sky "clarity"
 
    //  1 : pure air
    //  2 : exceptionnally clear
    //  4 : clear
    //  8 : light haze
    // 25 : haze
    // 64 : thin fog
 
    // Computed variables depending on the 2 above
 
    float zenith_luminance;// Y color component of the CIE color at zenith (luminance)
 
    float zenithColorX;// x color component of the CIE color at zenith
 
    float zenith_color_y;// y color component of the CIE color at zenith
 
    double eye_lum_conversion;// luminance conversion for an eye adapted to screen luminance (around 40 cd/m^2)
 
    double AY, BY, CY, DY, EY;// Distribution coefficients for the luminance distribution function
 
    float aX, bX, cX, dX, eX;// Distribution coefficients for x distribution function
 
    double Ay, By, Cy, Dy, Ey;// Distribution coefficients for y distribution function
 
    float termX;// Precomputed term for x calculation
 
    float termY;// Precomputed term for y calculation
 
    float term_Y;// Precomputed term for luminance calculation
 
    float sunPos[] = new float[3];
 
    static float thetas2;
 
    static float thetas3;
 
    static float T2;
}
 
class SkyBright {
    public SkyBright() {
        setDate(2003, 8, 0);
        setLoc(Constant.PI_OVER_FOUR, 1000.d, 25.d, 40.d);
        setSunMoon(0.5, 0.5);
    }
 
    /**
     * @param year
     * @param month     1=Jan, 12=Dec
     * @param moonPhase in radian 0=Full Moon, PI/2=First Quadrant/Last Quadran, PI=No Moon
     */
    public void setDate(int year, int month, double moonPhase) {
        magMoon = -12.73d + 1.4896903d * abs(moonPhase) + 0.04310727d * pow(moonPhase, 4.d);
 
        RA = (month - 3.d) * 0.52359878d;
 
        // Term for dark sky brightness computation
        bNightTerm = 1.0e-13 + 0.3e-13 * cos(0.56636d * (year - 1992.d));
    }
 
 
    public void setLoc(double latitude, double altitude, double temperature, double relativeHumidity) {
        double signLatitude = (latitude >= 0.d) ? 2.d - 1.d : 0 - 1.d;
 
        // extinction Coefficient for V band
        double KR = 0.1066d * exp(-altitude / 8200.d);
        double KA = 0.1d * exp(-altitude / 1500.d) * pow(1.d - 0.32d / log(relativeHumidity / 100.d), 1.33d) *
                (1.d + 0.33d * signLatitude * sin(RA));
        double KO = 0.031d * (3.d + 0.4d * (latitude * cos(RA) - cos(3.d * latitude))) / 3.d;
        double KW = 0.031d * 0.94d * (relativeHumidity / 100.d) * exp(temperature / 15.d) * exp(-altitude / 8200.d);
        K = KR + KA + KO + KW;
    }
 
    /**
     * Set the moon and sun zenith angular distance (cosin given) and precompute what can be
     *
     * @param cosDistMoonZenith
     * @param cosDistSunZenith
     */
    public void setSunMoon(double cosDistMoonZenith, double cosDistSunZenith) {
        // Air mass for Moon
        if (cosDistMoonZenith < 0) airMassMoon = 40.f;
        else airMassMoon = 1.f / (cosDistMoonZenith + 0.025f * exp(-11.f * cosDistMoonZenith));
 
        // Air mass for Sun
        if (cosDistSunZenith < 0) airMassSun = 40;
        else airMassSun = 1.f / (cosDistSunZenith + 0.025f * exp(-11.f * cosDistSunZenith));
 
        bMoonTerm1 = pow(10.f, -0.4 * (magMoon + 54.32f));
 
        C3 = pow(10.f, -0.4f * K * airMassMoon);// Term for moon brightness computation
 
        bTwilightTerm = -6.724f + 22.918312f * (Constant.PI_OVER_TWO - acos(cosDistSunZenith));
 
        C4 = pow(10.f, -0.4f * K * airMassSun);// Term for sky brightness computation
    }
 
    /**
     * Compute the luminance at the given position
     *
     * @param cosDistMoon cos(angular distance between moon and the position)
     * @param cosDistSun  cos(angular distance between sun  and the position)
     * @param cosDistZenithcos(angulardistancebetweenzenithandtheposition)
     *
     * @return
     */
    public double getLuminance(double cosDistMoon, double cosDistSun, double cosDistZenith) {
        // catch rounding errors here or end up with white flashes in some cases
        if (cosDistMoon < -1.d) cosDistMoon = -1.d;
        if (cosDistMoon > 1.d) cosDistMoon = 1.d;
        if (cosDistSun < -1.d) cosDistSun = -1.d;
        if (cosDistSun > 1.d) cosDistSun = 1.d;
        if (cosDistZenith < -1.d) cosDistZenith = -1.d;
        if (cosDistZenith > 1.d) cosDistZenith = 1.d;
 
        double distMoon = FastMath.acos(cosDistMoon);
        double distSun = FastMath.acos(cosDistSun);
 
        // Air mass
        double X = 1.d / (cosDistZenith + 0.025f * FastMath.exp(-11.d * cosDistZenith));
        double bKX = pow(10.d, -0.4f * K * X);
 
        // Dark night sky brightness
        bNight = 0.4f + 0.6f / FastMath.sqrt(0.04f + 0.96f * cosDistZenith * cosDistZenith);
        bNight *= bNightTerm * bKX;
 
        // Moonlight brightness
        double FM = 18886.28 / (distMoon * distMoon + 0.0007f) + pow(10.d, 6.15f - (distMoon + 0.001) * 1.43239f);
        FM += 229086.77f * (1.06f + cosDistMoon * cosDistMoon);
        bMoon = bMoonTerm1 * (1.d - bKX) * (FM * C3 + 440000.d * (1.d - C3));
 
        //Twilight brightness
        bTwilight = pow(10.d, bTwilightTerm + 0.063661977f * FastMath.acos(cosDistZenith) / K) *
                (1.7453293f / distSun) * (1.d - bKX);
 
        // Daylight brightness
        double FS = 18886.28f / (distSun * distSun + 0.0007f) + pow(10.d, 6.15f - (distSun + 0.001) * 1.43239f);
        FS += 229086.77f * (1.06f + cosDistSun * cosDistSun);
        bDaylight = 9.289663e-12 * (1.d - bKX) * (FS * C4 + 440000.d * (1.d - C4));
 
        // 27/08/2003 : Decide increase moonlight for more halo effect...
        bMoon *= 2.;
 
        // Total sky brightness
        bTotal = bDaylight > bTwilight ? bNight + bTwilight + bMoon : bNight + bDaylight + bMoon;
 
        return (bTotal < 0.d) ? 0.d : bTotal * 900900.9f * PI * 1e-4 * 3239389 * 2;
        //5;	// In cd/m^2 : the 32393895 is empirical term because the
        // lambert -> cd/m^2 formula seems to be wrong...
    }
 
    /*
250 REM  Visual limiting magnitude
260 BL=B(3)/1.11E-15 : REM in nanolamberts*/
 
    // Airmass for each component
    //cos_dist_zenith =cos(dist_zenith);
    //double gaz_mass = 1.f / ( cos_dist_zenith + 0.0286f *exp(-10.5f * cos_dist_zenith) );
    //double aerosol_mass = 1.f / ( cos_dist_zenith + 0.0123f *exp(-24.5f * cos_dist_zenith) );
    //double ozone_mass = 1.f /sqrt( 0.0062421903f - cos_dist_zenith * cos_dist_zenith / 1.0062814f );
    // Total extinction for V band
    //double DM = KR*gaz_mass + KA*aerosol_mass + KO*ozone_mass + KW*gaz_mass;
 
    /*
	// Visual limiting magnitude
	if (BL>1500.0)
	{
		C1 = 4.466825e-9;
		C2 = 1.258925e-6;
	}
	else
	{
		C1 = 1.584893e-10;
		C2 = 0.012589254;
	}
 
	double TH = C1*Math.pow(1.f+Math.sqrt(C2*BL),2.f); // in foot-candles
	double MN = -16.57-2.5*Math.log10(TH)-DM+5.0*Math.log10(SN); // Visual Limiting Magnitude
	*/
 
    /**
     * Air mass for the Moon
     */
    private double airMassMoon;
 
    /**
     * Air mass for the Sun
     */
    double airMassSun;
 
    /**
     * Total brightness
     */
    double bTotal;
 
    /**
     * Dark night brightness
     */
    double bNight;
 
    /**
     * Twilight brightness
     */
    double bTwilight;
 
    /**
     * Daylight sky brightness
     */
    double bDaylight;
 
    /**
     * Moon brightness
     */
    double bMoon;
 
    /**
     * Moon magnitude
     */
    double magMoon;
 
    /**
     * Something related with date
     */
    double RA;
 
    /**
     * Useful coef...
     */
    double K;
 
    /**
     * Term for moon brightness computation
     */
    double C3;
 
    /**
     * Term for sky brightness computation
     */
    double C4;
 
    /**
     * Snellen Ratio (20/20=1.0, good 20/10=2.0)
     */
    double SN = 1;
 
    // Optimisation variables
    double bNightTerm;
 
    double bMoonTerm1;
 
    double bTwilightTerm;
}

The method getSkyBrightness contains some variables that can be adapted in case you want to exagerate the brightness (atmBr, atmI) or the red glow close to the horizon (turbidity). I have set them in a fast way after a few tests to try to use acceptable values depending on the Sun elevation.

These kind of images are already being used for showing the sky during the day at the web page of the OAN ephemerides server.

Discussion

Enter your comment
   __    ____   ____   ___    __  ___
  / /   / __/  /  _/  / _ )  /  |/  /
 / /__ _\ \   _/ /   / _  | / /|_/ / 
/____//___/  /___/  /____/ /_/  /_/
 
 
blog/stellarium_like_renderings.txt · created: 2018/02/19 16:10 (Last modified 2018/11/21 11:19) by Tomás Alonso Albi
 
Recent changes RSS feed Creative Commons License Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki