mirror of
https://github.com/curioustorvald/Terrarum.git
synced 2026-03-07 12:21:52 +09:00
partially working sky model
This commit is contained in:
minjaesong
@@ -0,0 +1,34 @@
|
|||||||
|
fun sky(direction, location, date, time, conditions) -> spectral radiance
|
||||||
|
|
||||||
|
T<turbidity> = (t_m + t_h) / t_m; t_m = molecular atmos; t_h = haze atmos; See Fig 3
|
||||||
|
|
||||||
|
// CIE model for parametric model for clear skies
|
||||||
|
Y_C = Y_z * (0.91 + 10 * e pow (-3 * gamma) + 0.45 * cos2(gamma)) * (1 - e pow (-0.32/cos(theta)))
|
||||||
|
/((0.91 + 10 * e pow (-3 * theta) + 0.45 * cos2(theta_s)) * (1 - e pow (-0.32)));
|
||||||
|
Y_z = luminance at the zenith (see Fig 4)
|
||||||
|
|
||||||
|
// for overcast skies
|
||||||
|
Y_OC = Y_z * (1 + 2 * cos(theta)) / 3
|
||||||
|
|
||||||
|
// Perez model
|
||||||
|
fun F(theta, gamma) = (1 + A * e pow (B / cos(theta))) * (1 + C * e pow (D * gamma) + E * cos2(gamma));
|
||||||
|
A, B, C, D, E = distributional coefficients;
|
||||||
|
gamma, theta = angles as in Fig 4
|
||||||
|
|
||||||
|
// luminance Y for sky in any viewing direction
|
||||||
|
Y_P = Y_z * F(theta, gamma) / F(0, theta_s)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
x_z and y_z are the sun's spectral radiance, given in the Appendix
|
||||||
|
x<Chromacity> = x_z * F(theta, gamma) / F(0, theta_s),
|
||||||
|
y<Chromacity> = y_z * F(theta, gamma) / F(0, theta_s)
|
||||||
|
|
||||||
|
|
||||||
|
// Aerial perspective model
|
||||||
|
// attenuation
|
||||||
|
// TODO
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
J. Preetham, A & Shirley, Peter & E. Smits, Brian. (1999). A Practical Analytic Model for Daylight. Proceedings of ACM SIGGRAPH. 99. 91-100. 10.1145/311535.311545.
|
||||||
196
src/net/torvald/parametricsky/Application.kt
Normal file
196
src/net/torvald/parametricsky/Application.kt
Normal file
@@ -0,0 +1,196 @@
|
|||||||
|
package net.torvald.parametricsky
|
||||||
|
|
||||||
|
import com.badlogic.gdx.Game
|
||||||
|
import com.badlogic.gdx.Gdx
|
||||||
|
import com.badlogic.gdx.Screen
|
||||||
|
import com.badlogic.gdx.backends.lwjgl.LwjglApplication
|
||||||
|
import com.badlogic.gdx.backends.lwjgl.LwjglApplicationConfiguration
|
||||||
|
import com.badlogic.gdx.graphics.Pixmap
|
||||||
|
import com.badlogic.gdx.graphics.Texture
|
||||||
|
import com.badlogic.gdx.graphics.g2d.SpriteBatch
|
||||||
|
import com.badlogic.gdx.math.Affine2
|
||||||
|
import net.torvald.colourutil.CIEYXY
|
||||||
|
import net.torvald.colourutil.CIEXYZUtil.toXYZ
|
||||||
|
import net.torvald.colourutil.CIEXYZUtil.toColorRaw
|
||||||
|
import net.torvald.colourutil.CIEXYZUtil.toColor
|
||||||
|
import net.torvald.terrarum.inUse
|
||||||
|
import java.awt.BorderLayout
|
||||||
|
import java.awt.Dimension
|
||||||
|
import javax.swing.*
|
||||||
|
import kotlin.math.pow
|
||||||
|
|
||||||
|
|
||||||
|
const val WIDTH = 720
|
||||||
|
const val HEIGHT = 720
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Created by minjaesong on 2018-08-01.
|
||||||
|
*/
|
||||||
|
class Application : Game() {
|
||||||
|
|
||||||
|
/* Variables:
|
||||||
|
* 1. Canvas Y (theta)
|
||||||
|
* 2. Gamma (180deg - solar_azimuth; Canvas X)
|
||||||
|
* 3. Solar angle (theta_s)
|
||||||
|
* 4. Turbidity
|
||||||
|
*
|
||||||
|
* Sampling rate:
|
||||||
|
* theta in 70 downTo 0 step 10 (8 entries) // canvas
|
||||||
|
* gamma in -90..90 step 12 (16 entries) // canvas
|
||||||
|
* theta_s in 0..70 step 10 (8 entries) // time of the day
|
||||||
|
* turbidity in {1.5, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64} (12 entries) // weather of the day
|
||||||
|
*/
|
||||||
|
|
||||||
|
|
||||||
|
private lateinit var oneScreen: Pixmap
|
||||||
|
private lateinit var batch: SpriteBatch
|
||||||
|
|
||||||
|
private lateinit var testTex: Texture
|
||||||
|
|
||||||
|
var turbidity = 7.0
|
||||||
|
var thetaOfSun = 0.0
|
||||||
|
|
||||||
|
override fun getScreen(): Screen {
|
||||||
|
return super.getScreen()
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun setScreen(screen: Screen?) {
|
||||||
|
super.setScreen(screen)
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun render() {
|
||||||
|
Gdx.graphics.setTitle("Daylight Model — F: ${Gdx.graphics.framesPerSecond}")
|
||||||
|
|
||||||
|
genTexLoop(turbidity, thetaOfSun)
|
||||||
|
|
||||||
|
|
||||||
|
val tex = Texture(oneScreen)
|
||||||
|
tex.setFilter(Texture.TextureFilter.Linear, Texture.TextureFilter.Linear)
|
||||||
|
|
||||||
|
batch.inUse {
|
||||||
|
batch.draw(tex, 0f, 0f, WIDTH.toFloat(), HEIGHT.toFloat())
|
||||||
|
}
|
||||||
|
|
||||||
|
tex.dispose()
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun pause() {
|
||||||
|
super.pause()
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun resume() {
|
||||||
|
super.resume()
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun resize(width: Int, height: Int) {
|
||||||
|
super.resize(width, height)
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun dispose() {
|
||||||
|
oneScreen.dispose()
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Generated texture is as if you took the panorama picture of sky: up 70deg to horizon, east-south-west;
|
||||||
|
* with sun not moving (sun is at exact south, sun's height is adjustable)
|
||||||
|
*/
|
||||||
|
private fun genTexLoop(T: Double, theta_s: Double) {
|
||||||
|
// loop thru gamma and theta
|
||||||
|
for (y in 0..90) { // theta
|
||||||
|
for (x in 0..90) { // gamma
|
||||||
|
val theta = Math.toRadians(y.toDouble()) // of observer
|
||||||
|
val gamma = Math.toRadians(90 - x.toDouble()) // of observer
|
||||||
|
|
||||||
|
val Y_z = Model.getAbsoluteZenithLuminance(T, theta_s)
|
||||||
|
val x_z = Model.getZenithChromaX(T, theta_s)
|
||||||
|
val y_z = Model.getZenithChromaY(T, theta_s)
|
||||||
|
|
||||||
|
val Y_p = Y_z * Model.getFforLuma(theta, gamma, T) / Model.getFforLuma(0.0, theta_s, T)
|
||||||
|
val x_p = (x_z * Model.getFforChromaX(theta, gamma, T) / Model.getFforChromaX(0.0, theta_s, T)).coerceIn(0.0, 1.0)
|
||||||
|
val y_p = (y_z * Model.getFforChromaY(theta, gamma, T) / Model.getFforChromaY(0.0, theta_s, T)).coerceIn(0.0, 1.0)
|
||||||
|
|
||||||
|
val normalisedY = Y_p.toFloat().pow(0.5f).div(10f)
|
||||||
|
|
||||||
|
//println("$Y_p -> $normalisedY, $x_p, $y_p")
|
||||||
|
|
||||||
|
val rgbColour = CIEYXY(normalisedY, x_p.toFloat(), y_p.toFloat()).toXYZ().toColor()
|
||||||
|
//val rgbColour = CIEYXY(normalisedY, 0.3128f, 0.3290f).toXYZ().toColorRaw()
|
||||||
|
|
||||||
|
oneScreen.setColor(rgbColour)
|
||||||
|
oneScreen.drawPixel(x, y)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// end loop
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun create() {
|
||||||
|
batch = SpriteBatch()
|
||||||
|
testTex = Texture(Gdx.files.internal("assets/test_texture.tga"))
|
||||||
|
|
||||||
|
oneScreen = Pixmap(90, 90, Pixmap.Format.RGBA8888)
|
||||||
|
|
||||||
|
|
||||||
|
ApplicationController(this)
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
class ApplicationController(app: Application) : JFrame() {
|
||||||
|
|
||||||
|
val mainPanel = JPanel()
|
||||||
|
|
||||||
|
val turbidityControl = JSlider(2, 25, 7)
|
||||||
|
val theta_sControl = JSlider(0, 85, 0)
|
||||||
|
|
||||||
|
val turbidityValueDisp = JLabel()
|
||||||
|
val theta_sValueDisp = JLabel()
|
||||||
|
|
||||||
|
init {
|
||||||
|
val turbidityPanel = JPanel()
|
||||||
|
val theta_sPanel = JPanel()
|
||||||
|
|
||||||
|
turbidityPanel.add(JLabel("Turbidity"))
|
||||||
|
turbidityPanel.add(turbidityControl)
|
||||||
|
turbidityPanel.add(turbidityValueDisp)
|
||||||
|
|
||||||
|
turbidityValueDisp.text = turbidityControl.value.toString()
|
||||||
|
theta_sValueDisp.text = theta_sControl.value.toString()
|
||||||
|
|
||||||
|
theta_sPanel.add(JLabel("Theta_s"))
|
||||||
|
theta_sPanel.add(theta_sControl)
|
||||||
|
theta_sPanel.add(theta_sValueDisp)
|
||||||
|
|
||||||
|
mainPanel.add(turbidityPanel)
|
||||||
|
mainPanel.add(theta_sPanel)
|
||||||
|
|
||||||
|
this.isVisible = true
|
||||||
|
this.defaultCloseOperation = WindowConstants.EXIT_ON_CLOSE
|
||||||
|
this.size = Dimension(300, 400)
|
||||||
|
|
||||||
|
this.add(mainPanel)
|
||||||
|
|
||||||
|
|
||||||
|
turbidityControl.addChangeListener {
|
||||||
|
turbidityValueDisp.text = turbidityControl.value.toString()
|
||||||
|
app.turbidity = turbidityControl.value.toDouble()
|
||||||
|
}
|
||||||
|
|
||||||
|
theta_sControl.addChangeListener {
|
||||||
|
theta_sValueDisp.text = theta_sControl.value.toString()
|
||||||
|
app.thetaOfSun = Math.toRadians(theta_sControl.value.toDouble())
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
fun main(args: Array<String>) {
|
||||||
|
val config = LwjglApplicationConfiguration()
|
||||||
|
config.width = WIDTH
|
||||||
|
config.height = HEIGHT
|
||||||
|
config.foregroundFPS = 0
|
||||||
|
|
||||||
|
LwjglApplication(Application(), config)
|
||||||
|
}
|
||||||
390
src/net/torvald/parametricsky/Model.kt
Normal file
390
src/net/torvald/parametricsky/Model.kt
Normal file
@@ -0,0 +1,390 @@
|
|||||||
|
package net.torvald.parametricsky
|
||||||
|
|
||||||
|
import net.torvald.dataclass.Matrix
|
||||||
|
import kotlin.math.pow
|
||||||
|
|
||||||
|
/**
|
||||||
|
* This model is about the backdrop; to paint the background using realistic colour of the sky.
|
||||||
|
* This model does not provide the global light in itself; It's up to you.
|
||||||
|
*
|
||||||
|
* J. Preetham, A & Shirley, Peter & E. Smits, Brian. (1999). *A Practical Analytic Model for Daylight*. Proceedings of ACM SIGGRAPH. 1999. 91-100. 10.1145/311535.311545.
|
||||||
|
*
|
||||||
|
* This implementation is NOT FOR realtime calculation, especially in video games.
|
||||||
|
* You are advised to precalculate necessary colours and put the results as a colour map.
|
||||||
|
*
|
||||||
|
* Created by minjaesong on 2018-08-01.
|
||||||
|
*/
|
||||||
|
|
||||||
|
object Model {
|
||||||
|
|
||||||
|
/** Skylight Distribution Coefficients and Zenith Values */
|
||||||
|
data class DistributionCoeff(
|
||||||
|
var A: Double, var B: Double, var C: Double, var D: Double, var E: Double
|
||||||
|
)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param theta_s solar angle from zenith in radians. 0 means at zenith; 0.5pi means at horizon, >0.5pi indicates below horizon
|
||||||
|
* @param phi_s solar azimuth in radians. 0: East; 0.5pi: South; pi: West
|
||||||
|
*/
|
||||||
|
data class SolarPosition(var theta_s: Double, var phi_s: Double)
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
|
fun getFforLuma(theta: Double, gamma: Double, T: Double) = _getFbyPerez(theta, gamma, getLuminanceDistributionFun(T))
|
||||||
|
fun getFforChromaX(theta: Double, gamma: Double, T: Double) = _getFbyPerez(theta, gamma, getChromaXDistributionFun(T))
|
||||||
|
fun getFforChromaY(theta: Double, gamma: Double, T: Double) = _getFbyPerez(theta, gamma, getChromaYDistributionFun(T))
|
||||||
|
|
||||||
|
/*private*/ fun _getFbyPerez(theta: Double, gamma: Double, dc: DistributionCoeff): Double {
|
||||||
|
val A = dc.A; val B = dc.B; val C = dc.C; val D = dc.D; val E = dc.E
|
||||||
|
val e = Math.E
|
||||||
|
return (1.0 + A * e.pow(B / cos(theta))) *
|
||||||
|
(1.0 + C * e.pow(D * gamma) + E * cos2(gamma))
|
||||||
|
}
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
*/
|
||||||
|
fun getLuminanceDistributionFun(T: Double): DistributionCoeff {
|
||||||
|
val mat = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( 0.1787, -1.4630),
|
||||||
|
doubleArrayOf(-0.3554, 0.4275),
|
||||||
|
doubleArrayOf(-0.0227, 5.3251),
|
||||||
|
doubleArrayOf( 0.1206, -2.5771),
|
||||||
|
doubleArrayOf(-0.0670, 0.3703)
|
||||||
|
))
|
||||||
|
val mat2 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( T ),
|
||||||
|
doubleArrayOf(1.0)
|
||||||
|
))
|
||||||
|
val mat3 = mat * mat2
|
||||||
|
return DistributionCoeff(
|
||||||
|
mat3.data[0][0],
|
||||||
|
mat3.data[1][0],
|
||||||
|
mat3.data[2][0],
|
||||||
|
mat3.data[3][0],
|
||||||
|
mat3.data[4][0]
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
*/
|
||||||
|
fun getChromaXDistributionFun(T: Double): DistributionCoeff {
|
||||||
|
val mat = Matrix(arrayOf(
|
||||||
|
doubleArrayOf(-0.0193, -0.2592),
|
||||||
|
doubleArrayOf(-0.0665, 0.0008),
|
||||||
|
doubleArrayOf(-0.0004, 0.2125),
|
||||||
|
doubleArrayOf(-0.0641, -0.8989),
|
||||||
|
doubleArrayOf(-0.0033, 0.0452)
|
||||||
|
))
|
||||||
|
val mat2 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( T ),
|
||||||
|
doubleArrayOf(1.0)
|
||||||
|
))
|
||||||
|
val mat3 = mat * mat2
|
||||||
|
return DistributionCoeff(
|
||||||
|
mat3.data[0][0],
|
||||||
|
mat3.data[1][0],
|
||||||
|
mat3.data[2][0],
|
||||||
|
mat3.data[3][0],
|
||||||
|
mat3.data[4][0]
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
*/
|
||||||
|
fun getChromaYDistributionFun(T: Double): DistributionCoeff {
|
||||||
|
val mat = Matrix(arrayOf(
|
||||||
|
doubleArrayOf(-0.0167, -0.2608),
|
||||||
|
doubleArrayOf(-0.0950, 0.0092),
|
||||||
|
doubleArrayOf(-0.0079, 0.2102),
|
||||||
|
doubleArrayOf(-0.0441, -1.6537),
|
||||||
|
doubleArrayOf(-0.0109, 0.0529)
|
||||||
|
))
|
||||||
|
val mat2 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( T ),
|
||||||
|
doubleArrayOf(1.0)
|
||||||
|
))
|
||||||
|
val mat3 = mat * mat2
|
||||||
|
return DistributionCoeff(
|
||||||
|
mat3[0][0],
|
||||||
|
mat3[1][0],
|
||||||
|
mat3[2][0],
|
||||||
|
mat3[3][0],
|
||||||
|
mat3[4][0]
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
* @param theta_s angle from zenith in radians
|
||||||
|
* @return Luminance in candela per metre squared
|
||||||
|
*/
|
||||||
|
fun getAbsoluteZenithLuminance(T: Double, theta_s: Double): Double {
|
||||||
|
return (4.0453 * T - 4.9710) * tan((4.0 / 9.0 - T / 120.0) * (PI - 2 * theta_s)) - 0.2155 * T + 2.4192
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
* @param theta_s angle from zenith in radians
|
||||||
|
* @return X value of the chroma in CIEYxy colour space
|
||||||
|
*/
|
||||||
|
fun getZenithChromaX(T: Double, theta_s: Double): Double {
|
||||||
|
val mat1 = Matrix(arrayOf(doubleArrayOf(T.pow(2), T, 1.0)))
|
||||||
|
val mat2 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( 0.00166, -0.00375, 0.00209, 0.0),
|
||||||
|
doubleArrayOf(-0.02903, 0.06377, -0.03202, 0.00394),
|
||||||
|
doubleArrayOf( 0.11693, -0.21196, 0.06052, 0.25886)
|
||||||
|
))
|
||||||
|
val mat3 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf(theta_s.pow(3)),
|
||||||
|
doubleArrayOf(theta_s.pow(2)),
|
||||||
|
doubleArrayOf(theta_s),
|
||||||
|
doubleArrayOf(1.0)
|
||||||
|
))
|
||||||
|
|
||||||
|
return (mat1 * mat2 * mat3)[0][0]
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param T turbidity
|
||||||
|
* @param theta_s angle from zenith in radians
|
||||||
|
* @return Y value of the chroma in CIEYxy colour space
|
||||||
|
*/
|
||||||
|
fun getZenithChromaY(T: Double, theta_s: Double): Double {
|
||||||
|
val mat1 = Matrix(arrayOf(doubleArrayOf(T.pow(2), T, 1.0)))
|
||||||
|
val mat2 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf( 0.00275, -0.00610, 0.00317, 0.0),
|
||||||
|
doubleArrayOf(-0.04214, 0.08970, -0.04153, 0.00516),
|
||||||
|
doubleArrayOf( 0.15346, -0.26756, 0.06670, 0.26688)
|
||||||
|
))
|
||||||
|
val mat3 = Matrix(arrayOf(
|
||||||
|
doubleArrayOf(theta_s.pow(3)),
|
||||||
|
doubleArrayOf(theta_s.pow(2)),
|
||||||
|
doubleArrayOf(theta_s),
|
||||||
|
doubleArrayOf(1.0)
|
||||||
|
))
|
||||||
|
|
||||||
|
return (mat1 * mat2 * mat3)[0][0]
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/** refractive index of air */
|
||||||
|
private val RYL_n = 1.0003
|
||||||
|
/** number of molecules per unit volume */
|
||||||
|
private val RYL_N = 2.545 * (10 pow 25)
|
||||||
|
/** depolarisation factor */
|
||||||
|
private val RYL_p_n = 0.035
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Rayleigh angular scattering
|
||||||
|
*/
|
||||||
|
fun getBetaMofTheta(theta: Double, lambda: Double): Double {
|
||||||
|
return (PI_SQR * (RYL_n.pow(2.0) - 1) pow 2.0) / (2 * RYL_N * lambda.pow(4.0)) *
|
||||||
|
((6.0 + 3.0 * RYL_p_n) / (6.0 - 7.0 * RYL_p_n)) * (1.0 + cos2(theta))
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Rayleigh total scattering
|
||||||
|
*/
|
||||||
|
fun getBetaM(lambda: Double): Double {
|
||||||
|
return (8 * PI_SQR * (RYL_n.pow(2.0) - 1) pow 2.0) / (3 * RYL_N * lambda.pow(4.0)) *
|
||||||
|
((6.0 + 3.0 * RYL_p_n) / (6.0 - 7.0 * RYL_p_n))
|
||||||
|
}
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
|
/** Table 1 */
|
||||||
|
private val mieEtaLUT = arrayOf(
|
||||||
|
doubleArrayOf(4.192, 4.193, 4.177, 4.147, 4.072),
|
||||||
|
doubleArrayOf(3.311, 3.319, 3.329, 3.335, 3.339),
|
||||||
|
doubleArrayOf(2.860, 2.868, 2.878, 2.883, 2.888),
|
||||||
|
doubleArrayOf(2.518, 2.527, 2.536, 2.542, 2.547),
|
||||||
|
doubleArrayOf(1.122, 1.129, 1.138, 1.142, 1.147),
|
||||||
|
doubleArrayOf(0.3324, 0.3373, 0.3433, 0.3467, 0.3502),
|
||||||
|
doubleArrayOf(0.1644, 0.1682, 0.1730, 0.1757, 0.1785),
|
||||||
|
doubleArrayOf(0.1239, 0.1275, 0.1320, 0.1346, 0.1373),
|
||||||
|
doubleArrayOf(0.08734, 0.09111, 0.09591, 0.09871, 0.10167),
|
||||||
|
doubleArrayOf(0.08242, 0.08652, 0.09179, 0.09488, 0.09816),
|
||||||
|
doubleArrayOf(0.08313, 0.08767, 0.09352, 0.09697, 0.10065),
|
||||||
|
doubleArrayOf(0.09701, 0.1024, 0.1095, 0.1137, 0.1182),
|
||||||
|
doubleArrayOf(0.1307, 0.1368, 0.1447, 0.1495, 0.1566)
|
||||||
|
)
|
||||||
|
private val mieEtaLUTThetas = doubleArrayOf(1.0, 4.0, 7.0, 10.0, 30.0, 60.0, 80.0, 90.0, 110.0, 120.0, 130.0, 150.0, 180.0)
|
||||||
|
private val mieEtaLUTLambdas = doubleArrayOf(400.0, 450.0, 550.0, 650.0, 850.0)
|
||||||
|
private val solarSpectralQuantities = arrayOf(
|
||||||
|
/* K */doubleArrayOf(0.650393,0.653435,0.656387,0.657828,0.660644,0.662016,0.663365,0.665996,0.667276,0.668532,0.669765,0.670974,0.67216,0.673323,0.674462,0.675578,0.67667,0.677739,0.678784,0.678781,0.679802,0.6808,0.681775,0.681771,0.682722,0.683649,0.683646,0.68455,0.684546,0.685426,0.686282,0.686279,0.687112,0.687108,0.687917,0.687913,0.688699,0.688695,0.688691,0.689453,0.689449),
|
||||||
|
/* S_0 */doubleArrayOf(63.4,65.8,94.8,104.8,105.9,96.8,113.9,125.6,125.5,121.3,121.3,113.5,113.1,110.8,106.5,108.8,105.3,104.4,100.0,96.0,95.1,89.1,90.5,90.3,88.4,84.0,85.1,81.9,82.6,84.9,81.3,71.9,74.3,76.4,63.3,71.7,77.0,65.2,47.7,68.6,65.0),
|
||||||
|
/* S_1 */doubleArrayOf(38.5,35.0,43.4,46.3,43.9,37.1,36.7,35.9,32.6,27.9,24.3,20.1,16.2,13.2,8.6,6.1,4.2,1.9,0.0,-1.6,-3.5,-3.5,-5.8,-7.2,-8.6,-9.5,-10.9,-10.7,-12.0,-14.0,-13.6,-12.0,-13.3,-12.9,-10.6,-11.6,-12.2,-10.2,-7.8,-11.2,-10.4),
|
||||||
|
/* S_2 */doubleArrayOf(3.0,1.2,-1.1,-0.5,-0.7,-1.2,-2.6,-2.9,-2.8,-2.6,-2.6,-1.8,-1.5,-1.3,-1.2,-1.0,-0.5,-0.3,0.0,0.2,0.5,2.1,3.2,4.1,4.7,5.1,6.7,7.3,8.6,9.8,10.2,8.3,9.6,8.5,7.0,7.6,8.0,6.7,5.2,7.4,6.8),
|
||||||
|
/* Sun */doubleArrayOf(1655.90,1623.37,2112.75,2588.82,2582.91,2423.23,2676.05,2965.83,3054.54,3005.75,3066.37,2883.04,2871.21,2782.50,2710.06,2723.36,2636.13,2550.38,2506.02,2531.16,2535.59,2513.42,2463.15,2417.32,2368.53,2321.21,2282.77,2233.98,2197.02,2152.67,2109.79,2072.83,2024.04,1987.08,1942.72,1907.24,1862.89,1825.92,0.0,0.0,0.0),
|
||||||
|
/* k_o */doubleArrayOf(0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.003,0.006,0.009,0.014,0.021,0.030,0.040,0.048,0.063,0.075,0.085,0.103,0.120,0.120,0.115,0.125,0.120,0.105,0.090,0.079,0.067,0.057,0.048,0.036,0.028,0.023,0.018,0.014,0.011,0.010,0.009,0.007,0.004,0.000),
|
||||||
|
/* k_wa */doubleArrayOf(0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.01600,0.02400,0.01250,1.00000,0.87000,0.06100,0.00100,0.00001,0.00001,0.00060),
|
||||||
|
/* K_g */doubleArrayOf(0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,3.00,0.21,0.00)
|
||||||
|
)
|
||||||
|
private val mieKLUT = solarSpectralQuantities[0]
|
||||||
|
private val mieKLUTLambdas = doubleArrayOf( // also doubles as lambda for solar spectral quantities
|
||||||
|
380.0, 390.0, 400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, 480.0, 490.0, 500.0, 510.0, 520.0, 530.0, 540.0, 550.0, 560.0, 570.0, 580.0, 590.0, 600.0, 610.0, 620.0, 630.0, 640.0, 650.0, 660.0, 670.0, 680.0, 690.0, 700.0, 710.0, 720.0, 730.0, 740.0, 750.0, 760.0, 770.0, 780.0
|
||||||
|
)
|
||||||
|
/**
|
||||||
|
* e.g.
|
||||||
|
*
|
||||||
|
* 0 2 4 5 7 , find 3
|
||||||
|
*
|
||||||
|
* will return (1, 2), which corresponds value (2, 4) of which input value 3 is in between.
|
||||||
|
*/
|
||||||
|
private fun binarySearchInterval(value: Double, array: DoubleArray): Pair<Int, Int> {
|
||||||
|
var low: Int = 0
|
||||||
|
var high: Int = array.size - 1
|
||||||
|
|
||||||
|
while (low <= high) {
|
||||||
|
val mid = (low + high).ushr(1)
|
||||||
|
val midVal = array[mid]
|
||||||
|
|
||||||
|
if (value < midVal)
|
||||||
|
high = mid - 1
|
||||||
|
else if (value > midVal)
|
||||||
|
low = mid + 1
|
||||||
|
else
|
||||||
|
return Pair(mid, mid)
|
||||||
|
}
|
||||||
|
|
||||||
|
val first = Math.max(high, 0)
|
||||||
|
val second = Math.min(low, array.size - 1)
|
||||||
|
return Pair(first, second)
|
||||||
|
}
|
||||||
|
/** Just a table lookup with linear interpolation */
|
||||||
|
private fun getMieEta(theta: Double, lambda: Double): Double {
|
||||||
|
val lambdaIndex = binarySearchInterval(lambda, mieEtaLUTLambdas)
|
||||||
|
val lambdaInterval = Pair(mieEtaLUTLambdas[lambdaIndex.first], mieEtaLUTLambdas[lambdaIndex.second])
|
||||||
|
val thetaIndex = binarySearchInterval(theta, mieEtaLUTThetas)
|
||||||
|
val thetaInterval = Pair(mieEtaLUTThetas[thetaIndex.first], mieEtaLUTThetas[thetaIndex.second])
|
||||||
|
|
||||||
|
val lambdaStep = (lambda - mieEtaLUTLambdas[0]) / (lambdaInterval.second - lambdaInterval.first)
|
||||||
|
val thetaStep = (theta - mieEtaLUTThetas[0]) / (thetaInterval.second - thetaInterval.first)
|
||||||
|
|
||||||
|
return interpolateBilinear(lambdaStep, thetaStep,
|
||||||
|
mieEtaLUT[thetaIndex.first][lambdaIndex.first],
|
||||||
|
mieEtaLUT[thetaIndex.first][lambdaIndex.second],
|
||||||
|
mieEtaLUT[thetaIndex.second][lambdaIndex.first],
|
||||||
|
mieEtaLUT[thetaIndex.second][lambdaIndex.second]
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun getMieK(lambda: Double): Double {
|
||||||
|
val lambdaIndex = binarySearchInterval(lambda, mieKLUTLambdas)
|
||||||
|
val lambdaInterval = Pair(mieKLUTLambdas[lambdaIndex.first], mieKLUTLambdas[lambdaIndex.second])
|
||||||
|
|
||||||
|
val lambdaStep = (lambda - mieKLUTLambdas[0]) / (lambdaInterval.second - lambdaInterval.first)
|
||||||
|
|
||||||
|
return interpolateLinear(lambdaStep, mieKLUT[lambdaIndex.first], mieKLUT[lambdaIndex.second])
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Mie angular scattering
|
||||||
|
* @param T turbidity
|
||||||
|
* @param theta angle from zenith in radians
|
||||||
|
* @param lambda monochromatic light, valid between [400,850] nm
|
||||||
|
*/
|
||||||
|
fun getBetaPofTheta(T: Double, theta: Double, lambda: Double): Double {
|
||||||
|
val c = (0.6544 * T - 0.6510) * (10.0 pow -16.0)
|
||||||
|
return 0.434 * c * (TWOPI / lambda).pow(2.0) * 0.5 * getMieEta(theta, lambda)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Mie total scattering
|
||||||
|
* @param T turbidity
|
||||||
|
* @param lambda monochromatic light, valid between [380,780] nm
|
||||||
|
*/
|
||||||
|
fun getBetaP(T: Double, lambda: Double): Double {
|
||||||
|
val c = (0.6544 * T - 0.6510) * (10.0 pow -16.0)
|
||||||
|
return 0.434 * c * (TWOPI / lambda).pow(2.0) * getMieK(lambda)
|
||||||
|
}
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Get a solar time (t)
|
||||||
|
*
|
||||||
|
* @param t_s standard time in decimal hours
|
||||||
|
* @param J Julian date [1..365]
|
||||||
|
* @param SM standard meridian for the time zone in radians
|
||||||
|
* @param L obresver's longitude in radians
|
||||||
|
*/
|
||||||
|
fun t_solarTime(t_s: Double, J: Double, SM: Double = 0.0, L: Double = 0.0): Double {
|
||||||
|
return t_s + 0.170 * sin(FOURPI * (J - 80) / 373.0) - 0.129 * sin(TWOPI * (J - 8) / 355.0) + (12 * (SM - L) / PI)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Approximated solar declination
|
||||||
|
*
|
||||||
|
* @param J Julian date [1..365]
|
||||||
|
*/
|
||||||
|
fun delta_solarDeclination(J: Double): Double {
|
||||||
|
return 0.4093 * sin(TWOPI * (J - 81) / 368.0)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @param l obresver's latitude in radians
|
||||||
|
* @param delta solar declination in radians
|
||||||
|
* @param t solar time in decimal hours
|
||||||
|
*/
|
||||||
|
fun getSolarPosition(l: Double, delta: Double, t: Double): SolarPosition {
|
||||||
|
return SolarPosition(
|
||||||
|
HALFPI - arcsin(sin(l) * sin(delta) - cos(l) * cos(delta) * cos(PI * t / 12)),
|
||||||
|
arctan((-cos(delta) * sin(PI * t / 12)) / (cos(l) * sin(delta) - sin(l) * cos(delta) * cos(PI * t / 12)))
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
|
private val E = Math.E
|
||||||
|
private val HALFPI = 0.5 * Math.PI
|
||||||
|
private val PI = Math.PI
|
||||||
|
private val TWOPI = 2.0 * Math.PI
|
||||||
|
private val FOURPI = 4.0 * Math.PI
|
||||||
|
private val PI_SQR = PI * PI
|
||||||
|
private fun sin(a: Double) = Math.sin(a)
|
||||||
|
private fun cos(a: Double) = Math.cos(a)
|
||||||
|
private fun tan(a: Double) = Math.tan(a)
|
||||||
|
private fun arcsin(a: Double) = Math.asin(a)
|
||||||
|
private fun arctan(a: Double) = Math.atan(a)
|
||||||
|
private fun cos2(a: Double) = 0.5 * (1.0 + cos(2 * a))
|
||||||
|
private infix fun Double.pow(other: Double) = Math.pow(this, other)
|
||||||
|
private infix fun Int.pow(exp: Int): Int {
|
||||||
|
var exp = exp
|
||||||
|
var base = this
|
||||||
|
var result = 1
|
||||||
|
while (true) {
|
||||||
|
if (exp and 1 != 0)
|
||||||
|
result *= base
|
||||||
|
exp = exp shr 1
|
||||||
|
if (exp.inv() != 0)
|
||||||
|
break
|
||||||
|
base *= base
|
||||||
|
}
|
||||||
|
|
||||||
|
return result
|
||||||
|
}
|
||||||
|
private fun interpolateLinear(step: Double, start: Double, end: Double): Double {
|
||||||
|
if (start == end) {
|
||||||
|
return start
|
||||||
|
}
|
||||||
|
if (step <= 0f) {
|
||||||
|
return start
|
||||||
|
}
|
||||||
|
return if (step >= 1f) {
|
||||||
|
end
|
||||||
|
}
|
||||||
|
else (1f - step) * start + step * end
|
||||||
|
}
|
||||||
|
/** X and Y starts at top left; X goes down, Y goes right as value increases positively */
|
||||||
|
private fun interpolateBilinear(stepX: Double, stepY: Double, topLeft: Double, topRight: Double, bottomLeft: Double, bottomRight: Double): Double {
|
||||||
|
val val1 = interpolateLinear(stepX, topLeft, topRight)
|
||||||
|
val val2 = interpolateLinear(stepX, bottomLeft, bottomRight)
|
||||||
|
return interpolateLinear(stepY, val1, val2)
|
||||||
|
}
|
||||||
|
}
|
||||||
Reference in New Issue
Block a user