+#include <cmath>
#include "beamgate.h"
#include "layout.h"
#include "signal.h"
using namespace std;
using namespace Msp;
+namespace {
+
+double optical_thickness(double zenith_cosine)
+{
+ double r = 635; // Radius of Earth divided by atmospheric thickness at zenith
+ double rc = r*zenith_cosine;
+ return sqrt(2*r+1+rc*rc)-rc;
+}
+
+GL::Color spectrum_to_rgb(const double *spectrum)
+{
+ // Data from http://cvrl.ioo.ucl.ac.uk/cmfs.htm
+ static double color_matching[] =
+ {
+ 3.769647E-03, 4.146161E-04, 1.847260E-02,
+ 9.382967E-03, 1.059646E-03, 4.609784E-02,
+ 2.214302E-02, 2.452194E-03, 1.096090E-01,
+ 4.742986E-02, 4.971717E-03, 2.369246E-01,
+ 8.953803E-02, 9.079860E-03, 4.508369E-01,
+ 1.446214E-01, 1.429377E-02, 7.378822E-01,
+ 2.035729E-01, 2.027369E-02, 1.051821E+00,
+ 2.488523E-01, 2.612106E-02, 1.305008E+00,
+ 2.918246E-01, 3.319038E-02, 1.552826E+00,
+ 3.227087E-01, 4.157940E-02, 1.748280E+00,
+ 3.482554E-01, 5.033657E-02, 1.917479E+00,
+ 3.418483E-01, 5.743393E-02, 1.918437E+00,
+ 3.224637E-01, 6.472352E-02, 1.848545E+00,
+ 2.826646E-01, 7.238339E-02, 1.664439E+00,
+ 2.485254E-01, 8.514816E-02, 1.522157E+00,
+ 2.219781E-01, 1.060145E-01, 1.428440E+00,
+ 1.806905E-01, 1.298957E-01, 1.250610E+00,
+ 1.291920E-01, 1.535066E-01, 9.991789E-01,
+ 8.182895E-02, 1.788048E-01, 7.552379E-01,
+ 4.600865E-02, 2.064828E-01, 5.617313E-01,
+ 2.083981E-02, 2.379160E-01, 4.099313E-01,
+ 7.097731E-03, 2.850680E-01, 3.105939E-01,
+ 2.461588E-03, 3.483536E-01, 2.376753E-01,
+ 3.649178E-03, 4.277595E-01, 1.720018E-01,
+ 1.556989E-02, 5.204972E-01, 1.176796E-01,
+ 4.315171E-02, 6.206256E-01, 8.283548E-02,
+ 7.962917E-02, 7.180890E-01, 5.650407E-02,
+ 1.268468E-01, 7.946448E-01, 3.751912E-02,
+ 1.818026E-01, 8.575799E-01, 2.438164E-02,
+ 2.405015E-01, 9.071347E-01, 1.566174E-02,
+ 3.098117E-01, 9.544675E-01, 9.846470E-03,
+ 3.804244E-01, 9.814106E-01, 6.131421E-03,
+ 4.494206E-01, 9.890228E-01, 3.790291E-03,
+ 5.280233E-01, 9.994608E-01, 2.327186E-03,
+ 6.133784E-01, 9.967737E-01, 1.432128E-03,
+ 7.016774E-01, 9.902549E-01, 8.822531E-04,
+ 7.967750E-01, 9.732611E-01, 5.452416E-04,
+ 8.853376E-01, 9.424569E-01, 3.386739E-04,
+ 9.638388E-01, 8.963613E-01, 2.117772E-04,
+ 1.051011E+00, 8.587203E-01, 1.335031E-04,
+ 1.109767E+00, 8.115868E-01, 8.494468E-05,
+ 1.143620E+00, 7.544785E-01, 5.460706E-05,
+ 1.151033E+00, 6.918553E-01, 3.549661E-05,
+ 1.134757E+00, 6.270066E-01, 2.334738E-05,
+ 1.083928E+00, 5.583746E-01, 1.554631E-05,
+ 1.007344E+00, 4.895950E-01, 1.048387E-05,
+ 9.142877E-01, 4.229897E-01, 0.000000E+00,
+ 8.135565E-01, 3.609245E-01, 0.000000E+00,
+ 6.924717E-01, 2.980865E-01, 0.000000E+00,
+ 5.755410E-01, 2.416902E-01, 0.000000E+00,
+ 4.731224E-01, 1.943124E-01, 0.000000E+00,
+ 3.844986E-01, 1.547397E-01, 0.000000E+00,
+ 2.997374E-01, 1.193120E-01, 0.000000E+00,
+ 2.277792E-01, 8.979594E-02, 0.000000E+00,
+ 1.707914E-01, 6.671045E-02, 0.000000E+00,
+ 1.263808E-01, 4.899699E-02, 0.000000E+00,
+ 9.224597E-02, 3.559982E-02, 0.000000E+00,
+ 6.639960E-02, 2.554223E-02, 0.000000E+00,
+ 4.710606E-02, 1.807939E-02, 0.000000E+00,
+ 3.292138E-02, 1.261573E-02, 0.000000E+00,
+ 2.262306E-02, 8.661284E-03, 0.000000E+00,
+ 1.575417E-02, 6.027677E-03, 0.000000E+00,
+ 1.096778E-02, 4.195941E-03, 0.000000E+00,
+ 7.608750E-03, 2.910864E-03, 0.000000E+00,
+ 5.214608E-03, 1.995557E-03, 0.000000E+00,
+ 3.569452E-03, 1.367022E-03, 0.000000E+00,
+ 2.464821E-03, 9.447269E-04, 0.000000E+00,
+ 1.703876E-03, 6.537050E-04, 0.000000E+00,
+ 1.186238E-03, 4.555970E-04, 0.000000E+00,
+ 8.269535E-04, 3.179738E-04, 0.000000E+00,
+ 5.758303E-04, 2.217445E-04, 0.000000E+00,
+ 4.058303E-04, 1.565566E-04, 0.000000E+00,
+ 2.856577E-04, 1.103928E-04, 0.000000E+00,
+ 2.021853E-04, 7.827442E-05, 0.000000E+00,
+ 1.438270E-04, 5.578862E-05, 0.000000E+00,
+ 1.024685E-04, 3.981884E-05, 0.000000E+00,
+ 7.347551E-05, 2.860175E-05, 0.000000E+00,
+ 5.259870E-05, 2.051259E-05, 0.000000E+00,
+ 3.806114E-05, 1.487243E-05, 0.000000E+00,
+ 2.758222E-05, 1.080001E-05, 0.000000E+00,
+ 2.004122E-05, 7.863920E-06, 0.000000E+00,
+ 1.458792E-05, 5.736935E-06, 0.000000E+00,
+ 1.068141E-05, 4.211597E-06, 0.000000E+00,
+ 7.857521E-06, 3.106561E-06, 0.000000E+00,
+ 5.768284E-06, 2.286786E-06, 0.000000E+00,
+ 4.259166E-06, 1.693147E-06, 0.000000E+00,
+ 3.167765E-06, 1.262556E-06, 0.000000E+00,
+ 2.358723E-06, 9.422514E-07, 0.000000E+00,
+ 1.762465E-06, 7.053860E-07, 0.000000E+00,
+ };
+
+ // Compute CIE XYZ tristimulus values
+ double x = 0;
+ double y = 0;
+ double z = 0;
+ for(unsigned i=0; i<89; ++i)
+ {
+ x += spectrum[i]*color_matching[i*3];
+ y += spectrum[i]*color_matching[i*3+1];
+ z += spectrum[i]*color_matching[i*3+2];
+ }
+
+ // Compute linear RGB color
+ return GL::Color(3.2406*x-1.5372*y-0.4986*z, -0.9689*x+1.8758*y+0.0415*z, 0.0557*x-0.2040*y+1.0570*z);
+}
+
+}
+
namespace R2C2 {
Layout3D::Layout3D(Layout &l):
{
// South, 15° from zenith
sun.set_position(GL::Vector4(0, -0.259, 0.966, 0));
- sun.set_diffuse(GL::Color(0.7));
+ sun.set_diffuse(GL::Color(0.0));
+ sun.set_specular(GL::Color(0.0));
lighting.set_ambient(GL::Color(0.2));
lighting.attach(0, sun);
Vector diff = Vector(sun.get_position())-sun_dir;
if(diff.norm()>0.0025f)
sun.set_position(GL::Vector4(sun_dir, 0.0f));
- // TODO replace these very unscientific guesses with correct formulas
- lighting.set_ambient(GL::Color(pow(0.2f*(sun_dir.z+1.0f), 2.0f)));
- sun.set_diffuse(GL::Color(0.7f*max(sun_dir.z, 0.0f)));
+
+ double T = 5777;
+ double h = 6.62606957e-34; // Planck constant
+ double k = 1.3806488e-23; // Boltzmann constant
+ double c = 299792458; // Speed of light in vacuum
+ double wavelengths[89];
+ for(unsigned i=0; i<89; ++i)
+ wavelengths[i] = (390+i*5)/1e9;
+
+ // Compute the spectrum of sunlight in space
+ double incoming_spectrum[89];
+ for(unsigned i=0; i<89; ++i)
+ incoming_spectrum[i] = (2*h*pow(c, 2))/pow(wavelengths[i], 5)/(exp(h*c/(wavelengths[i]*k*T))-1);
+
+ // Apply rayleigh scattering to get the spectrum at the surface
+ double sun_thickness = optical_thickness(sun_dir.z);
+ double transmitted_spectrum[89];
+ for(unsigned i=0; i<89; ++i)
+ {
+ double f = 7.9e-27/pow(wavelengths[i], 4);
+ transmitted_spectrum[i] = incoming_spectrum[i]*exp(-f*sun_thickness);
+ }
+
+ GL::Color sunlight_color = spectrum_to_rgb(transmitted_spectrum)*2e-15;
+ /* When the sun is setting, the amount of direct illlumination is
+ proportional to the fraction of sun still over the horizon. This formula
+ assumes a square sun, but the difference is hardly noticeable. */
+ double cutoff = sun_dir.z>0.005 ? 1 : sun_dir.z>-0.005 ? (sun_dir.z+0.005)*100 : 0;
+ sun.set_diffuse(sunlight_color*cutoff);
+
+ double skylight_spectrum[89] = {};
+ double total_weight = 0;
+ for(int x=-2; x<=2; ++x)
+ for(int y=-2; y<=2; ++y)
+ for(int z=0; z<=2; ++z)
+ {
+ if(abs(x)!=2 && abs(y)!=2 && abs(z)!=2)
+ continue;
+
+ double l_sq = x*x+y*y+z*z;
+ double l = sqrt(l_sq);
+ double view_thickness = optical_thickness(z/l);
+ double w = z/l/l_sq;
+ for(unsigned i=0; i<89; ++i)
+ {
+ double f = 7.9e-27/pow(wavelengths[i], 4);
+ skylight_spectrum[i] += w*incoming_spectrum[i]*(exp(-f*sun_thickness)-exp(-f*view_thickness))/(1-sun_thickness/view_thickness);
+ }
+ total_weight += w;
+ }
+
+ GL::Color skylight_color = spectrum_to_rgb(skylight_spectrum)*(2e-15/total_weight);
+ lighting.set_ambient(skylight_color*0.5);
}
void Layout3D::object_added(Object &o)
projects / r2c2.git / commitdiff