source/3d/layout.cpp

   1 #include <cmath>
   2 #include "beamgate.h"
   3 #include "layout.h"
   4 #include "signal.h"
   5 #include "terrain.h"
   6 #include "track.h"
   7 #include "utility.h"
   8 #include "vehicle.h"
   9
  10 using namespace std;
  11 using namespace Msp;
  12
  13 namespace {
  14
  15 double optical_thickness(double zenith_cosine)
  16 {
  17         double r = 635;  // Radius of Earth divided by atmospheric thickness at zenith
  18         double rc = r*zenith_cosine;
  19         return sqrt(2*r+1+rc*rc)-rc;
  20 }
  21
  22 GL::Color spectrum_to_rgb(const double *spectrum)
  23 {
  24         // Data from http://cvrl.ioo.ucl.ac.uk/cmfs.htm
  25         static double color_matching[] =
  26         {
  27                 3.769647E-03, 4.146161E-04, 1.847260E-02,
  28                 9.382967E-03, 1.059646E-03, 4.609784E-02,
  29                 2.214302E-02, 2.452194E-03, 1.096090E-01,
  30                 4.742986E-02, 4.971717E-03, 2.369246E-01,
  31                 8.953803E-02, 9.079860E-03, 4.508369E-01,
  32                 1.446214E-01, 1.429377E-02, 7.378822E-01,
  33                 2.035729E-01, 2.027369E-02, 1.051821E+00,
  34                 2.488523E-01, 2.612106E-02, 1.305008E+00,
  35                 2.918246E-01, 3.319038E-02, 1.552826E+00,
  36                 3.227087E-01, 4.157940E-02, 1.748280E+00,
  37                 3.482554E-01, 5.033657E-02, 1.917479E+00,
  38                 3.418483E-01, 5.743393E-02, 1.918437E+00,
  39                 3.224637E-01, 6.472352E-02, 1.848545E+00,
  40                 2.826646E-01, 7.238339E-02, 1.664439E+00,
  41                 2.485254E-01, 8.514816E-02, 1.522157E+00,
  42                 2.219781E-01, 1.060145E-01, 1.428440E+00,
  43                 1.806905E-01, 1.298957E-01, 1.250610E+00,
  44                 1.291920E-01, 1.535066E-01, 9.991789E-01,
  45                 8.182895E-02, 1.788048E-01, 7.552379E-01,
  46                 4.600865E-02, 2.064828E-01, 5.617313E-01,
  47                 2.083981E-02, 2.379160E-01, 4.099313E-01,
  48                 7.097731E-03, 2.850680E-01, 3.105939E-01,
  49                 2.461588E-03, 3.483536E-01, 2.376753E-01,
  50                 3.649178E-03, 4.277595E-01, 1.720018E-01,
  51                 1.556989E-02, 5.204972E-01, 1.176796E-01,
  52                 4.315171E-02, 6.206256E-01, 8.283548E-02,
  53                 7.962917E-02, 7.180890E-01, 5.650407E-02,
  54                 1.268468E-01, 7.946448E-01, 3.751912E-02,
  55                 1.818026E-01, 8.575799E-01, 2.438164E-02,
  56                 2.405015E-01, 9.071347E-01, 1.566174E-02,
  57                 3.098117E-01, 9.544675E-01, 9.846470E-03,
  58                 3.804244E-01, 9.814106E-01, 6.131421E-03,
  59                 4.494206E-01, 9.890228E-01, 3.790291E-03,
  60                 5.280233E-01, 9.994608E-01, 2.327186E-03,
  61                 6.133784E-01, 9.967737E-01, 1.432128E-03,
  62                 7.016774E-01, 9.902549E-01, 8.822531E-04,
  63                 7.967750E-01, 9.732611E-01, 5.452416E-04,
  64                 8.853376E-01, 9.424569E-01, 3.386739E-04,
  65                 9.638388E-01, 8.963613E-01, 2.117772E-04,
  66                 1.051011E+00, 8.587203E-01, 1.335031E-04,
  67                 1.109767E+00, 8.115868E-01, 8.494468E-05,
  68                 1.143620E+00, 7.544785E-01, 5.460706E-05,
  69                 1.151033E+00, 6.918553E-01, 3.549661E-05,
  70                 1.134757E+00, 6.270066E-01, 2.334738E-05,
  71                 1.083928E+00, 5.583746E-01, 1.554631E-05,
  72                 1.007344E+00, 4.895950E-01, 1.048387E-05,
  73                 9.142877E-01, 4.229897E-01, 0.000000E+00,
  74                 8.135565E-01, 3.609245E-01, 0.000000E+00,
  75                 6.924717E-01, 2.980865E-01, 0.000000E+00,
  76                 5.755410E-01, 2.416902E-01, 0.000000E+00,
  77                 4.731224E-01, 1.943124E-01, 0.000000E+00,
  78                 3.844986E-01, 1.547397E-01, 0.000000E+00,
  79                 2.997374E-01, 1.193120E-01, 0.000000E+00,
  80                 2.277792E-01, 8.979594E-02, 0.000000E+00,
  81                 1.707914E-01, 6.671045E-02, 0.000000E+00,
  82                 1.263808E-01, 4.899699E-02, 0.000000E+00,
  83                 9.224597E-02, 3.559982E-02, 0.000000E+00,
  84                 6.639960E-02, 2.554223E-02, 0.000000E+00,
  85                 4.710606E-02, 1.807939E-02, 0.000000E+00,
  86                 3.292138E-02, 1.261573E-02, 0.000000E+00,
  87                 2.262306E-02, 8.661284E-03, 0.000000E+00,
  88                 1.575417E-02, 6.027677E-03, 0.000000E+00,
  89                 1.096778E-02, 4.195941E-03, 0.000000E+00,
  90                 7.608750E-03, 2.910864E-03, 0.000000E+00,
  91                 5.214608E-03, 1.995557E-03, 0.000000E+00,
  92                 3.569452E-03, 1.367022E-03, 0.000000E+00,
  93                 2.464821E-03, 9.447269E-04, 0.000000E+00,
  94                 1.703876E-03, 6.537050E-04, 0.000000E+00,
  95                 1.186238E-03, 4.555970E-04, 0.000000E+00,
  96                 8.269535E-04, 3.179738E-04, 0.000000E+00,
  97                 5.758303E-04, 2.217445E-04, 0.000000E+00,
  98                 4.058303E-04, 1.565566E-04, 0.000000E+00,
  99                 2.856577E-04, 1.103928E-04, 0.000000E+00,
 100                 2.021853E-04, 7.827442E-05, 0.000000E+00,
 101                 1.438270E-04, 5.578862E-05, 0.000000E+00,
 102                 1.024685E-04, 3.981884E-05, 0.000000E+00,
 103                 7.347551E-05, 2.860175E-05, 0.000000E+00,
 104                 5.259870E-05, 2.051259E-05, 0.000000E+00,
 105                 3.806114E-05, 1.487243E-05, 0.000000E+00,
 106                 2.758222E-05, 1.080001E-05, 0.000000E+00,
 107                 2.004122E-05, 7.863920E-06, 0.000000E+00,
 108                 1.458792E-05, 5.736935E-06, 0.000000E+00,
 109                 1.068141E-05, 4.211597E-06, 0.000000E+00,
 110                 7.857521E-06, 3.106561E-06, 0.000000E+00,
 111                 5.768284E-06, 2.286786E-06, 0.000000E+00,
 112                 4.259166E-06, 1.693147E-06, 0.000000E+00,
 113                 3.167765E-06, 1.262556E-06, 0.000000E+00,
 114                 2.358723E-06, 9.422514E-07, 0.000000E+00,
 115                 1.762465E-06, 7.053860E-07, 0.000000E+00,
 116         };
 117
 118         // Compute CIE XYZ tristimulus values
 119         double x = 0;
 120         double y = 0;
 121         double z = 0;
 122         for(unsigned i=0; i<89; ++i)
 123         {
 124                 x += spectrum[i]*color_matching[i*3];
 125                 y += spectrum[i]*color_matching[i*3+1];
 126                 z += spectrum[i]*color_matching[i*3+2];
 127         }
 128
 129         // Compute linear RGB color
 130         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);
 131 }
 132
 133 }
 134
 135 namespace R2C2 {
 136
 137 Layout3D::Layout3D(Layout &l):
 138         layout(l),
 139         catalogue(layout.get_catalogue())
 140 {
 141         // South, 15° from zenith
 142         sun.set_position(GL::Vector4(0, -0.259, 0.966, 0));
 143         sun.set_diffuse(GL::Color(0.0));
 144         sun.set_specular(GL::Color(0.0));
 145         lighting.set_ambient(GL::Color(0.2));
 146         lighting.attach(0, sun);
 147         lighting.set_sky_color(GL::Color(0.2));
 148         lighting.set_horizon_angle(Geometry::Angle<float>::from_degrees(5));
 149
 150         layout.signal_object_added.connect(sigc::mem_fun(this, &Layout3D::object_added));
 151         layout.signal_object_removed.connect(sigc::mem_fun(this, &Layout3D::object_removed));
 152
 153         const set<Object *> &lobjs = layout.get_all<Object>();
 154         for(set<Object *>::iterator i=lobjs.begin(); i!=lobjs.end(); ++i)
 155                 object_added(**i);
 156 }
 157
 158 Layout3D::~Layout3D()
 159 {
 160         while(!utilities.empty())
 161                 delete *utilities.begin();
 162         while(!objects.empty())
 163                 delete objects.begin()->second;
 164 }
 165
 166 void Layout3D::add(Object3D &o)
 167 {
 168         insert_unique(objects, &o.get_object(), &o);
 169 }
 170
 171 Object3D &Layout3D::get_3d(Object &o) const
 172 {
 173         return *get_item(objects, &o);
 174 }
 175
 176 void Layout3D::remove(Object3D &o)
 177 {
 178         objects.erase(&o.get_object());
 179 }
 180
 181 void Layout3D::add(Utility3D &u)
 182 {
 183         utilities.insert(&u);
 184 }
 185
 186 void Layout3D::remove(Utility3D &u)
 187 {
 188         utilities.erase(&u);
 189 }
 190
 191 void Layout3D::tick()
 192 {
 193         Time::TimeDelta t = layout.get_clock().get_current_time();
 194         Angle time_of_day = Angle::from_turns(t/Time::day);
 195         Angle latitude = Angle::from_degrees(52.5);
 196         Angle axial_tilt = Angle::from_degrees(23.45);
 197         Transform trans = Transform::rotation(Angle::quarter_turn()-latitude, Vector(-1, 0, 0))*
 198                 Transform::rotation(Angle::half_turn()-time_of_day, Vector(0, 0, 1))*
 199                 Transform::rotation(axial_tilt, Vector(-1, 0, 0));
 200         Vector sun_dir = trans.transform_linear(Vector(0, -1, 0));
 201         Vector diff = sun.get_position().slice<3>(0)-sun_dir;
 202         if(diff.norm()>0.0025f)
 203                 sun.set_position(compose(sun_dir, 0.0f));
 204
 205         double T = 5777;
 206         double h = 6.62606957e-34;  // Planck constant
 207         double k = 1.3806488e-23;   // Boltzmann constant
 208         double c = 299792458;       // Speed of light in vacuum
 209         double wavelengths[89];
 210         for(unsigned i=0; i<89; ++i)
 211                 wavelengths[i] = (390+i*5)/1e9;
 212
 213         // Compute the spectrum of sunlight in space
 214         double incoming_spectrum[89];
 215         for(unsigned i=0; i<89; ++i)
 216                 incoming_spectrum[i] = (2*h*pow(c, 2))/pow(wavelengths[i], 5)/(exp(h*c/(wavelengths[i]*k*T))-1);
 217
 218         // Apply rayleigh scattering to get the spectrum at the surface
 219         double sun_thickness = optical_thickness(sun_dir.z);
 220         double transmitted_spectrum[89];
 221         for(unsigned i=0; i<89; ++i)
 222         {
 223                 double f = 7.9e-27/pow(wavelengths[i], 4);
 224                 transmitted_spectrum[i] = incoming_spectrum[i]*exp(-f*sun_thickness);
 225         }
 226
 227         GL::Color sunlight_color = spectrum_to_rgb(transmitted_spectrum)*2e-15;
 228         /* When the sun is setting, the amount of direct illlumination is
 229         proportional to the fraction of sun still over the horizon.  This formula
 230         assumes a square sun, but the difference is hardly noticeable. */
 231         double cutoff = sun_dir.z>0.005 ? 1 : sun_dir.z>-0.005 ? (sun_dir.z+0.005)*100 : 0;
 232         sun.set_diffuse(sunlight_color*cutoff);
 233
 234         double skylight_spectrum[89] = {};
 235         double total_weight = 0;
 236         for(int x=-2; x<=2; ++x)
 237                 for(int y=-2; y<=2; ++y)
 238                         for(int z=0; z<=2; ++z)
 239                         {
 240                                 if(abs(x)!=2 && abs(y)!=2 && abs(z)!=2)
 241                                         continue;
 242
 243                                 double l_sq = x*x+y*y+z*z;
 244                                 double l = sqrt(l_sq);
 245                                 double view_thickness = optical_thickness(z/l);
 246                                 double w = z/l/l_sq;
 247                                 for(unsigned i=0; i<89; ++i)
 248                                 {
 249                                         double f = 7.9e-27/pow(wavelengths[i], 4);
 250                                         skylight_spectrum[i] += w*incoming_spectrum[i]*(exp(-f*sun_thickness)-exp(-f*view_thickness))/(1-sun_thickness/view_thickness);
 251                                 }
 252                                 total_weight += w;
 253                         }
 254
 255         GL::Color skylight_color = spectrum_to_rgb(skylight_spectrum)*(2e-15/total_weight);
 256         lighting.set_sky_color(skylight_color*0.65);
 257         lighting.set_ambient(skylight_color*0.35);
 258 }
 259
 260 void Layout3D::object_added(Object &o)
 261 {
 262         if(Track *t = dynamic_cast<Track *>(&o))
 263                 new Track3D(*this, *t);
 264         else if(Signal *s = dynamic_cast<Signal *>(&o))
 265                 new Signal3D(*this, *s);
 266         else if(Vehicle *v = dynamic_cast<Vehicle *>(&o))
 267                 new Vehicle3D(*this, *v);
 268         else if(BeamGate *g = dynamic_cast<BeamGate *>(&o))
 269                 new BeamGate3D(*this, *g);
 270         else if(Terrain *r = dynamic_cast<Terrain *>(&o))
 271                 new Terrain3D(*this, *r);
 272 }
 273
 274 void Layout3D::object_removed(Object &o)
 275 {
 276         ObjectMap::iterator i = objects.find(&o);
 277         if(i!=objects.end())
 278                 delete i->second;
 279 }
 280
 281 } // namespace R2C2