3
* This file Copyright (C) 2004 Robert Spearman
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
// This is an ad hoc meteor model
21
// Could use a simple ablation physics model in the future
24
NOTE: Here the radiant is always along the ecliptic at the apex of the Earth's way.
25
In reality, individual meteor streams have varying velocity vectors and therefore radiants
26
which are generally not at the apex of the Earth's way, such as the Perseids shower.
29
// Improved realism and efficiency 2004-12
34
Meteor::Meteor(Projector *proj, Navigator* nav, ToneReproductor* eye, double v)
36
// velocity = 11+(double)rand()/((double)RAND_MAX+1)*v; // abs range 11-72 km/s
41
// determine meteor model view matrix (want z in dir of travel of earth, z=0 at center of earth)
42
// meteor life is so short, no need to recalculate
43
double equ_rotation; // rotation needed to align with path of earth
44
Vec3d sun_dir = nav->helio_to_earth_equ( Vec3d(0,0,0) );
46
Mat4d tmat = Mat4d::xrotation(-23.45f*M_PI/180.f); // ecliptical tilt
47
sun_dir.transfo4d(tmat); // convert to ecliptical coordinates
49
equ_rotation = acos( sun_dir.dot( Vec3d(1,0,0) ) );
50
if( sun_dir[1] < 0 ) equ_rotation = 2*M_PI - equ_rotation;
52
equ_rotation -= M_PI_2;
54
mmat = Mat4d::xrotation(23.45f*M_PI/180.f) * Mat4d::zrotation(equ_rotation) * Mat4d::yrotation(M_PI_2);
57
// select random trajectory using polar coordinates in XY plane, centered on observer
58
xydistance = (double)rand()/((double)RAND_MAX+1)*(VISIBLE_RADIUS);
59
double angle = (double)rand()/((double)RAND_MAX+1)*2*M_PI;
61
// find observer position in meteor coordinate system
62
obs = nav->local_to_earth_equ(Vec3d(0,0,EARTH_RADIUS));
63
obs.transfo4d(mmat.transpose());
65
// set meteor start x,y
66
pos_internal[0] = pos_train[0] = position[0] = xydistance*cos(angle) +obs[0];
67
pos_internal[1] = pos_train[1] = position[1] = xydistance*sin(angle) +obs[1];
69
// determine life of meteor (start and end z value based on atmosphere burn altitudes)
71
// D is distance from center of earth
72
double D = sqrt( position[0]*position[0] + position[1]*position[1] );
74
if( D > EARTH_RADIUS+HIGH_ALTITUDE ) {
80
start_h = sqrt( pow(EARTH_RADIUS+HIGH_ALTITUDE,2) - D*D);
82
// determine end of burn point, and nearest point to observer for distance mag calculation
83
// mag should be max at nearest point still burning
84
if( D > EARTH_RADIUS+LOW_ALTITUDE ) {
85
end_h = -start_h; // earth grazing
86
min_dist = xydistance;
88
end_h = sqrt( pow(EARTH_RADIUS+LOW_ALTITUDE,2) - D*D);
89
min_dist = sqrt( xydistance*xydistance + pow( end_h - obs[2], 2) );
92
if(min_dist > VISIBLE_RADIUS ) {
93
// on average, not visible (although if were zoomed ...)
99
// limit lifetime to 0.5-3.0 sec
100
double tmp_h = start_h - velocity * (0.5 + (double)rand()/((double)RAND_MAX+1) * 2.5);
101
if( tmp_h > end_h ) {
106
pos_train[2] = position[2] = start_h;
108
// printf("New meteor: %f %f s:%f e:%f v:%f\n", position[0], position[1], start_h, end_h, velocity);
113
// Determine drawing color given magnitude and eye
114
// (won't be visible during daylight)
116
// *** color varies somewhat based on velocity, plus atmosphere reddening
118
// determine intensity
119
float Mag1 = (double)rand()/((double)RAND_MAX+1)*6.75f - 3;
120
float Mag2 = (double)rand()/((double)RAND_MAX+1)*6.75f - 3;
121
float Mag = (Mag1 + Mag2)/2.0f;
123
mag = (5. + Mag) / 256.0;
124
if (mag>250) mag = mag - 256;
126
float term1 = expf(-0.92103f*(mag + 12.12331f)) * 108064.73f;
131
// Compute the equivalent star luminance for a 5 arc min circle and convert it
132
// in function of the eye adaptation
133
rmag = eye->adapt_luminance(term1);
134
rmag = rmag/powf(proj->get_fov(),0.85f)*50.f;
136
// if size of star is too small (blink) we put its size to 1.2 --> no more blink
137
// And we compensate the difference of brighteness with cmag
139
cmag=rmag*rmag/1.44f;
142
mag = cmag; // assumes white
144
// most visible meteors are under about 180km distant
145
// scale max mag down if outside this range
147
if(min_dist!=0) scale = 180*180/(min_dist*min_dist);
148
if( scale < 1 ) mag *= scale;
156
// returns true if alive
157
bool Meteor::update(int delta_time)
160
if(!alive) return(0);
162
if( position[2] < end_h ) {
163
// burning has stopped so magnitude fades out
164
// assume linear fade out
166
mag -= max_mag * (double)delta_time/500.0f;
167
if( mag < 0 ) alive=0; // no longer visible
171
// *** would need time direction multiplier to allow reverse time replay
172
position[2] = position[2] - velocity/1000.0f*(double)delta_time;
174
// train doesn't extend beyond start of burn
175
if( position[2] + velocity*0.5f > start_h ) {
176
pos_train[2] = start_h ;
178
pos_train[2] -= velocity*(double)delta_time/1000.0f;
181
//printf("meteor position: %f delta_t %d\n", position[2], delta_time);
183
// determine visual magnitude based on distance to observer
184
double dist = sqrt( xydistance*xydistance + pow( position[2]-obs[2], 2) );
186
if( dist == 0 ) dist = .01; // just to be cautious (meteor hits observer!)
188
dist_multiplier = min_dist*min_dist / (dist*dist);
194
// returns true if visible
195
bool Meteor::draw(Projector *proj, Navigator* nav)
198
if(!alive) return(0);
202
Vec3d spos = position;
203
Vec3d epos = pos_train;
206
spos.transfo4d(mmat);
207
epos.transfo4d(mmat);
209
// convert to local and correct for earth radius [since equ and local coordinates in stellarium use same 0 point!]
210
spos = nav->earth_equ_to_local( spos );
211
epos = nav->earth_equ_to_local( epos );
212
spos[2] -= EARTH_RADIUS;
213
epos[2] -= EARTH_RADIUS;
215
int t1 = proj->project_local_check(spos/1216, start); // 1216 is to scale down under 1 for desktop version
216
int t2 = proj->project_local_check(epos/1216, end);
218
// don't draw if not visible (but may come into view)
219
if( t1 + t2 == 0 ) return 1;
221
// printf("[%f %f %f] (%d, %d) (%d, %d)\n", position[0], position[1], position[2], (int)start[0], (int)start[1], (int)end[0], (int)end[1]);
224
// connect this point with last drawn point
226
double tmag = mag*dist_multiplier;
228
// compute an intermediate point so can curve slightly along projection distortions
230
Vec3d posi = pos_internal;
231
posi[2] = position[2] + (pos_train[2] - position[2])/2;
232
posi.transfo4d(mmat);
233
posi = nav->earth_equ_to_local( posi );
234
posi[2] -= EARTH_RADIUS;
235
proj->project_local(posi/1216, intpos);
237
// draw dark to light
238
glBegin(GL_LINE_STRIP);
240
glVertex3f(end[0],end[1],0);
241
glColor4f(1,1,1,tmag/2);
242
glVertex3f(intpos[0],intpos[1],0);
243
glColor4f(1,1,1,tmag);
244
glVertex3f(start[0],start[1],0);
247
//glPointSize(1); // glpoint is not portable...
249
glVertex3f(start[0],start[1],0);
254
// TEMP - show radiant
255
Vec3d radiant = Vec3d(0,0,0.5f);
256
radiant.transfo4d(mmat);
257
if( projection->project_earth_equ(radiant, start) ) {
260
glVertex3f(start[0]-10,start[1],0);
261
glVertex3f(start[0]+10,start[1],0);
265
glVertex3f(start[0],start[1]-10,0);
266
glVertex3f(start[0],start[1]+10,0);
276
bool Meteor::is_alive(void)
added
removed
src/meteor.cpp
