3
* Copyright (C) 2003 Fabien Chereau
5
* This program is free software; you can redistribute it and/or
6
* modify it under the terms of the GNU General Public License
7
* as published by the Free Software Foundation; either version 2
8
* of the License, or (at your option) any later version.
10
* This program is distributed in the hope that it will be useful,
11
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
* GNU General Public License for more details.
15
* You should have received a copy of the GNU General Public License
16
* along with this program; if not, write to the Free Software
17
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20
// Class which compute and display the daylight sky color using openGL
21
// the sky is computed with the Skylight class.
23
// TODO : Adaptative resolution for optimization
25
#include "stellarium.h"
26
#include "stel_utility.h"
27
#include "stellastro.h"
28
#include "stel_atmosphere.h"
30
Atmosphere::Atmosphere() : sky_resolution(48), tab_sky(NULL), world_adaptation_luminance(0.f), atm_intensity(0)
32
// Create the vector array used to store the sky color on the full field of view
33
tab_sky = new Vec3f*[sky_resolution+1];
34
for (int k=0; k<sky_resolution+1 ;k++)
36
tab_sky[k] = new Vec3f[sky_resolution+1];
41
Atmosphere::~Atmosphere()
43
for (int k=0; k<sky_resolution+1 ;k++)
45
if (tab_sky[k]) delete [] tab_sky[k];
47
if (tab_sky) delete [] tab_sky;
50
void Atmosphere::compute_color(double JD, Vec3d sunPos, Vec3d moonPos, float moon_phase,
51
ToneReproductor * eye, Projector* prj,
52
float latitude, float altitude, float temperature, float relative_humidity)
55
float min_mw_lum = 0.13;
57
// no need to calculate if not visible
58
if(!fader.getInterstate())
61
world_adaptation_luminance = 3.75f;
62
milkyway_adaptation_luminance = min_mw_lum; // brighter than without atm, since no drawing addition of atm brightness
67
atm_intensity = fader.getInterstate();
74
// these are for radii
75
double sun_angular_size = atan(696000./AU/sunPos.length());
76
double moon_angular_size = atan(1738./AU/moonPos.length());
78
double touch_angle = sun_angular_size + moon_angular_size;
79
double dark_angle = moon_angular_size - sun_angular_size;
84
// determine luminance falloff during solar eclipses
85
double separation_angle = acos( sunPos.dot( moonPos )); // angle between them
87
// printf("touch at %f\tnow at %f (%f)\n", touch_angle, separation_angle, separation_angle/touch_angle);
89
// bright stars should be visible at total eclipse
90
// TODO: correct for atmospheric diffusion
91
// TODO: use better coverage function (non-linear)
92
// because of above issues, this algorithm darkens more quickly than reality
93
if( separation_angle < touch_angle)
99
float asun = sun_angular_size*sun_angular_size;
100
min = (asun - moon_angular_size*moon_angular_size)/asun; // minimum proportion of sun uncovered
103
else min = 0.004; // so bright stars show up at total eclipse
105
if(separation_angle < dark_angle) atm_intensity = min;
106
else atm_intensity *= min + (1.-min)*(separation_angle-dark_angle)/(touch_angle-dark_angle);
108
// printf("atm int %f (min %f)\n", atm_intensity, min);
112
sun_pos[0] = sunPos[0];
113
sun_pos[1] = sunPos[1];
114
sun_pos[2] = sunPos[2];
117
moon_pos[0] = moonPos[0];
118
moon_pos[1] = moonPos[1];
119
moon_pos[2] = moonPos[2];
121
sky.set_paramsv(sun_pos, 5.f);
123
skyb.set_loc(latitude * M_PI/180., altitude, temperature, relative_humidity);
124
skyb.set_sun_moon(moon_pos[2], sun_pos[2]);
126
// Calculate the date from the julian day.
130
skyb.set_date(date.years, date.months, moon_phase);
132
float stepX = (float)prj->getViewportWidth() / sky_resolution;
133
float stepY = (float)prj->getViewportHeight() / sky_resolution;
134
float viewport_left = (float)prj->getViewportPosX();
135
float viewport_bottom = (float)prj->getViewportPosY();
137
Vec3d point(1., 0., 0.);
139
// Variables used to compute the average sky luminance
141
unsigned int nb_lum = 0;
143
// Compute the sky color for every point above the ground
144
for (int x=0; x<=sky_resolution; ++x)
146
for(int y=0; y<=sky_resolution; ++y)
148
prj->unproject_local((double)viewport_left+x*stepX, (double)viewport_bottom+y*stepY,point);
153
point[2] = -point[2];
154
// The sky below the ground is the symetric of the one above :
155
// it looks nice and gives proper values for brightness estimation
158
b2.pos[0] = point[0];
159
b2.pos[1] = point[1];
160
b2.pos[2] = point[2];
162
// Use the Skylight model for the color
163
sky.get_xyY_valuev(b2);
165
// Use the Skybright.cpp 's models for brightness which gives better results.
166
b2.color[2] = skyb.get_luminance(moon_pos[0]*b2.pos[0]+moon_pos[1]*b2.pos[1]+
167
moon_pos[2]*b2.pos[2], sun_pos[0]*b2.pos[0]+sun_pos[1]*b2.pos[1]+
168
sun_pos[2]*b2.pos[2], b2.pos[2]);
171
sum_lum+=b2.color[2];
173
eye->xyY_to_RGB(b2.color);
174
tab_sky[x][y].set(atm_intensity*b2.color[0],atm_intensity*b2.color[1],atm_intensity*b2.color[2]);
178
world_adaptation_luminance = 3.75f + 3.5*sum_lum/nb_lum*atm_intensity;
179
milkyway_adaptation_luminance = min_mw_lum*(1-atm_intensity) + 30*sum_lum/nb_lum*atm_intensity;
187
// Draw the atmosphere using the precalc values stored in tab_sky
188
void Atmosphere::draw(Projector* prj, int delta_time)
190
if(fader.getInterstate())
192
// printf("Atm int: %f\n", atm_intensity);
193
glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_COLOR);
195
float stepX = (float)prj->getViewportWidth() / sky_resolution;
196
float stepY = (float)prj->getViewportHeight() / sky_resolution;
197
float viewport_left = (float)prj->getViewportPosX();
198
float view_bottom = (float)prj->getViewportPosY();
200
glDisable(GL_TEXTURE_2D);
202
prj->set_orthographic_projection(); // set 2D coordinate
203
for (int y2=0; y2<sky_resolution; ++y2)
205
glBegin(GL_QUAD_STRIP);
206
for(int x2=0; x2<sky_resolution+1; ++x2)
208
glColor3fv(tab_sky[x2][y2]);
209
glVertex2i((int)(viewport_left+x2*stepX),(int)(view_bottom+y2*stepY));
210
glColor3fv(tab_sky[x2][y2+1]);
211
glVertex2i((int)(viewport_left+x2*stepX),(int)(view_bottom+(y2+1)*stepY));
215
prj->reset_perspective_projection();