15 double optical_thickness(double zenith_cosine)
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;
22 GL::Color spectrum_to_rgb(const double *spectrum)
24 // Data from http://cvrl.ioo.ucl.ac.uk/cmfs.htm
25 static double color_matching[] =
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,
118 // Compute CIE XYZ tristimulus values
122 for(unsigned i=0; i<89; ++i)
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];
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);
137 Layout3D::Layout3D(Layout &l):
139 catalogue(layout.get_catalogue())
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));
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));
153 const set<Object *> &lobjs = layout.get_all<Object>();
154 for(set<Object *>::iterator i=lobjs.begin(); i!=lobjs.end(); ++i)
158 Layout3D::~Layout3D()
160 while(!utilities.empty())
161 delete *utilities.begin();
162 while(!objects.empty())
163 delete objects.begin()->second;
166 void Layout3D::add(Object3D &o)
168 insert_unique(objects, &o.get_object(), &o);
171 Object3D &Layout3D::get_3d(Object &o) const
173 return *get_item(objects, &o);
176 void Layout3D::remove(Object3D &o)
178 objects.erase(&o.get_object());
181 void Layout3D::add(Utility3D &u)
183 utilities.insert(&u);
186 void Layout3D::remove(Utility3D &u)
191 void Layout3D::tick()
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));
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;
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);
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)
223 double f = 7.9e-27/pow(wavelengths[i], 4);
224 transmitted_spectrum[i] = incoming_spectrum[i]*exp(-f*sun_thickness);
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);
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)
240 if(abs(x)!=2 && abs(y)!=2 && abs(z)!=2)
243 double l_sq = x*x+y*y+z*z;
244 double l = sqrt(l_sq);
245 double view_thickness = optical_thickness(z/l);
247 for(unsigned i=0; i<89; ++i)
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);
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);
260 void Layout3D::object_added(Object &o)
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);
274 void Layout3D::object_removed(Object &o)
276 ObjectMap::iterator i = objects.find(&o);