1 /*
2 OpenUniverse 1.0
3 Copyright (C) 2000 Raul Alonso <amil@las.es>
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
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.
14
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., 675 Mass Ave, Cambridge, MA 02139, USA.
18 */
19
20 #include "ou.h"
21 #include <string.h>
22
23
gettime(void)24 float gettime(void)
25 {
26 static float told = 0.0f;
27 float tnew, ris;
28
29 tnew = glutGet((GLenum) GLUT_ELAPSED_TIME);
30
31 ris = tnew - told;
32
33 told = tnew;
34
35 return ris / 1000.0;
36 }
37
38 /* Rotates point around the axis, 'angle' radians,
39 borrowed from Mesa */
Rotation(double angle,double axis[3],double point[3])40 void Rotation(double angle, double axis[3], double point[3])
41 {
42
43 double mag, s, c;
44 double x, y, z, xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;
45
46 s = sin(angle);
47 c = cos(angle);
48
49 mag = MODULE(axis);
50
51 if (mag == 0.0)
52 return;
53
54 x = axis[X] / mag;
55 y = axis[Y] / mag;
56 z = axis[Z] / mag;
57
58 xx = x * x;
59 yy = y * y;
60 zz = z * z;
61 xy = x * y;
62 yz = y * z;
63 zx = z * x;
64 xs = x * s;
65 ys = y * s;
66 zs = z * s;
67 one_c = 1.0 - c;
68
69 x = point[X] * ((one_c * xx) + c);
70 y = point[X] * ((one_c * xy) - zs);
71 z = point[X] * ((one_c * zx) + ys);
72
73 x += point[Y] * ((one_c * xy) + zs);
74 y += point[Y] * ((one_c * yy) + c);
75 z += point[Y] * ((one_c * yz) - xs);
76
77 x += point[Z] * ((one_c * zx) - ys);
78 y += point[Z] * ((one_c * yz) + xs);
79 z += point[Z] * ((one_c * zz) + c);
80
81 point[X] = x;
82 point[Y] = y;
83 point[Z] = z;
84 }
85
86
MakeQuat(double * q,double angle,double ax,double ay,double az)87 void MakeQuat(double *q, double angle, double ax, double ay, double az)
88 {
89 double d;
90 double s;
91
92 angle /= 2.0;
93 d = DISTANCE(ax, ay, az);
94 ax /= d;
95 ay /= d;
96 az /= d;
97
98 s = -sin(angle);
99
100 q[0] = cos(angle);
101 q[1] = s * ax;
102 q[2] = s * ay;
103 q[3] = s * az;
104 }
105
106
Quat2Matrix(double * q,float m[16])107 void Quat2Matrix(double *q, float m[16])
108 {
109 float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
110
111 x2 = q[1] + q[1];
112 y2 = q[2] + q[2];
113 z2 = q[3] + q[3];
114 xx = q[1] * x2;
115 xy = q[1] * y2;
116 xz = q[1] * z2;
117 yy = q[2] * y2;
118 yz = q[2] * z2;
119 zz = q[3] * z2;
120 wx = q[0] * x2;
121 wy = q[0] * y2;
122 wz = q[0] * z2;
123
124 m[0] = 1.0 - (yy + zz);
125 m[1] = xy - wz;
126 m[2] = xz + wy;
127 m[3] = 0.0;
128
129 m[4] = xy + wz;
130 m[5] = 1.0 - (xx + zz);
131 m[6] = yz - wx;
132 m[7] = 0.0;
133
134 m[8] = xz - wy;
135 m[9] = yz + wx;
136 m[10] = 1.0 - (xx + yy);
137 m[11] = 0.0;
138
139 m[12] = 0;
140 m[13] = 0;
141 m[14] = 0;
142 m[15] = 1;
143 }
144
145
Mat2Quat(float m[4][4],double * q)146 void Mat2Quat(float m[4][4], double *q)
147 {
148 float tr, s, q2[4];
149 int i, j, k;
150 int nxt[3] = { 1, 2, 0 };
151
152 tr = m[0][0] + m[1][1] + m[2][2];
153
154 if (tr > 0.0) {
155 s = sqrt(tr + 1.0);
156 q[0] = s / 2.0;
157 s = 0.5 / s;
158 q[1] = (m[1][2] - m[2][1]) * s;
159 q[2] = (m[2][0] - m[0][2]) * s;
160 q[3] = (m[0][1] - m[1][0]) * s;
161 } else {
162 i = 0;
163 if (m[1][1] > m[0][0])
164 i = 1;
165 if (m[2][2] > m[i][i])
166 i = 2;
167 j = nxt[i];
168 k = nxt[j];
169
170 s = sqrt((m[i][i] - (m[j][j] + m[k][k])) + 1.0);
171
172 q2[i] = s * 0.5;
173
174 if (s != 0.0)
175 s = 0.5 / s;
176
177 q2[3] = (m[j][k] - m[k][j]) * s;
178 q2[j] = (m[i][j] + m[j][i]) * s;
179 q2[k] = (m[i][k] + m[k][i]) * s;
180
181 q[1] = q2[0];
182 q[2] = q2[1];
183 q[3] = q2[2];
184 q[0] = q2[3];
185 }
186 }
187
Euler2Quat(double roll,double pitch,double yaw,double * q)188 void Euler2Quat(double roll, double pitch, double yaw, double *q)
189 {
190 float cr, cp, cy, sr, sp, sy, cpcy, spsy;
191
192 cr = cos(roll / 2);
193 cp = cos(pitch / 2);
194 cy = cos(yaw / 2);
195
196 sr = sin(roll / 2);
197 sp = sin(pitch / 2);
198 sy = sin(yaw / 2);
199
200 cpcy = cp * cy;
201 spsy = sp * sy;
202
203 q[0] = cr * cpcy + sr * spsy;
204 q[1] = sr * cpcy - cr * spsy;
205 q[2] = cr * sp * cy + sr * cp * sy;
206 q[3] = cr * cp * sy - sr * sp * cy;
207
208 }
209
210
211
212 #define SWAP_ROWS(a, b) { double *_tmp = a; (a)=(b); (b)=_tmp; }
213 #define MAT(m,r,c) m[(c)*4+(r)]
214
215 /*
216 * Compute inverse of 4x4 transformation matrix.
217 * Code contributed by Jacques Leroy jle@star.be
218 * Return TRUE for success, FALSE for failure (singular matrix)
219 *
220 * Aug '99: Adapted for openuniverse from the Mesa 3.1 source code
221 */
invert_matrix_general(double mat[16],double inv[16])222 static int invert_matrix_general(double mat[16], double inv[16])
223 {
224 double *m = mat;
225 double *out = inv;
226 double wtmp[4][8];
227 double m0, m1, m2, m3, s;
228 double *r0, *r1, *r2, *r3;
229
230 r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
231
232 r0[0] = MAT(m, 0, 0), r0[1] = MAT(m, 0, 1),
233 r0[2] = MAT(m, 0, 2), r0[3] = MAT(m, 0, 3),
234 r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
235 r1[0] = MAT(m, 1, 0), r1[1] = MAT(m, 1, 1),
236 r1[2] = MAT(m, 1, 2), r1[3] = MAT(m, 1, 3),
237 r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
238 r2[0] = MAT(m, 2, 0), r2[1] = MAT(m, 2, 1),
239 r2[2] = MAT(m, 2, 2), r2[3] = MAT(m, 2, 3),
240 r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
241 r3[0] = MAT(m, 3, 0), r3[1] = MAT(m, 3, 1),
242 r3[2] = MAT(m, 3, 2), r3[3] = MAT(m, 3, 3),
243 r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
244
245 /* choose pivot - or die */
246 if (fabs(r3[0]) > fabs(r2[0]))
247 SWAP_ROWS(r3, r2);
248 if (fabs(r2[0]) > fabs(r1[0]))
249 SWAP_ROWS(r2, r1);
250 if (fabs(r1[0]) > fabs(r0[0]))
251 SWAP_ROWS(r1, r0);
252 if (0.0 == r0[0])
253 return 0;
254
255 /* eliminate first variable */
256 m1 = r1[0] / r0[0];
257 m2 = r2[0] / r0[0];
258 m3 = r3[0] / r0[0];
259 s = r0[1];
260 r1[1] -= m1 * s;
261 r2[1] -= m2 * s;
262 r3[1] -= m3 * s;
263 s = r0[2];
264 r1[2] -= m1 * s;
265 r2[2] -= m2 * s;
266 r3[2] -= m3 * s;
267 s = r0[3];
268 r1[3] -= m1 * s;
269 r2[3] -= m2 * s;
270 r3[3] -= m3 * s;
271 s = r0[4];
272 if (s != 0.0) {
273 r1[4] -= m1 * s;
274 r2[4] -= m2 * s;
275 r3[4] -= m3 * s;
276 }
277 s = r0[5];
278 if (s != 0.0) {
279 r1[5] -= m1 * s;
280 r2[5] -= m2 * s;
281 r3[5] -= m3 * s;
282 }
283 s = r0[6];
284 if (s != 0.0) {
285 r1[6] -= m1 * s;
286 r2[6] -= m2 * s;
287 r3[6] -= m3 * s;
288 }
289 s = r0[7];
290 if (s != 0.0) {
291 r1[7] -= m1 * s;
292 r2[7] -= m2 * s;
293 r3[7] -= m3 * s;
294 }
295
296 /* choose pivot - or die */
297 if (fabs(r3[1]) > fabs(r2[1]))
298 SWAP_ROWS(r3, r2);
299 if (fabs(r2[1]) > fabs(r1[1]))
300 SWAP_ROWS(r2, r1);
301 if (0.0 == r1[1])
302 return 0;
303
304 /* eliminate second variable */
305 m2 = r2[1] / r1[1];
306 m3 = r3[1] / r1[1];
307 r2[2] -= m2 * r1[2];
308 r3[2] -= m3 * r1[2];
309 r2[3] -= m2 * r1[3];
310 r3[3] -= m3 * r1[3];
311 s = r1[4];
312 if (0.0 != s) {
313 r2[4] -= m2 * s;
314 r3[4] -= m3 * s;
315 }
316 s = r1[5];
317 if (0.0 != s) {
318 r2[5] -= m2 * s;
319 r3[5] -= m3 * s;
320 }
321 s = r1[6];
322 if (0.0 != s) {
323 r2[6] -= m2 * s;
324 r3[6] -= m3 * s;
325 }
326 s = r1[7];
327 if (0.0 != s) {
328 r2[7] -= m2 * s;
329 r3[7] -= m3 * s;
330 }
331
332 /* choose pivot - or die */
333 if (fabs(r3[2]) > fabs(r2[2]))
334 SWAP_ROWS(r3, r2);
335 if (0.0 == r2[2])
336 return 0;
337
338 /* eliminate third variable */
339 m3 = r3[2] / r2[2];
340 r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
341 r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6], r3[7] -= m3 * r2[7];
342
343 /* last check */
344 if (0.0 == r3[3])
345 return 0;
346
347 s = 1.0 / r3[3]; /* now back substitute row 3 */
348 r3[4] *= s;
349 r3[5] *= s;
350 r3[6] *= s;
351 r3[7] *= s;
352
353 m2 = r2[3]; /* now back substitute row 2 */
354 s = 1.0 / r2[2];
355 r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
356 r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
357 m1 = r1[3];
358 r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
359 r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
360 m0 = r0[3];
361 r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
362 r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
363
364 m1 = r1[2]; /* now back substitute row 1 */
365 s = 1.0 / r1[1];
366 r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
367 r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
368 m0 = r0[2];
369 r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
370 r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
371
372 m0 = r0[1]; /* now back substitute row 0 */
373 s = 1.0 / r0[0];
374 r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
375 r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
376
377 MAT(out, 0, 0) = r0[4];
378 MAT(out, 0, 1) = r0[5], MAT(out, 0, 2) = r0[6];
379 MAT(out, 0, 3) = r0[7], MAT(out, 1, 0) = r1[4];
380 MAT(out, 1, 1) = r1[5], MAT(out, 1, 2) = r1[6];
381 MAT(out, 1, 3) = r1[7], MAT(out, 2, 0) = r2[4];
382 MAT(out, 2, 1) = r2[5], MAT(out, 2, 2) = r2[6];
383 MAT(out, 2, 3) = r2[7], MAT(out, 3, 0) = r3[4];
384 MAT(out, 3, 1) = r3[5], MAT(out, 3, 2) = r3[6];
385 MAT(out, 3, 3) = r3[7];
386
387 return 1;
388 }
389
390 #undef SWAP_ROWS
391
392
393
394 /* Based upon the gluLookAt function from the Mesa 3.1 source tree.
395 The very only differences with the original function are:
396 -Removed the call to glTranlate
397 -Rotation matrix is inverted
398
399 This function is used to rotate spaceships so they orientation
400 matches their dir and up vectors.
401
402 */
MyLookAt(double dir[3],double up[3])403 void MyLookAt(double dir[3], double up[3])
404 {
405 GLdouble m[16], inv[16];
406 GLdouble x[3], y[3], z[3];
407 GLdouble mag;
408
409 /* Make rotation matrix */
410
411 /* Z vector */
412 z[0] = dir[0];
413 z[1] = dir[1];
414 z[2] = dir[2];
415 mag = sqrt(z[0] * z[0] + z[1] * z[1] + z[2] * z[2]);
416 if (mag) { /* mpichler, 19950515 */
417 z[0] /= mag;
418 z[1] /= mag;
419 z[2] /= mag;
420 }
421
422 /* Y vector */
423 y[0] = up[0];
424 y[1] = up[1];
425 y[2] = up[2];
426
427 /* X vector = Y cross Z */
428 x[0] = y[1] * z[2] - y[2] * z[1];
429 x[1] = -y[0] * z[2] + y[2] * z[0];
430 x[2] = y[0] * z[1] - y[1] * z[0];
431
432 /* Recompute Y = Z cross X */
433 y[0] = z[1] * x[2] - z[2] * x[1];
434 y[1] = -z[0] * x[2] + z[2] * x[0];
435 y[2] = z[0] * x[1] - z[1] * x[0];
436
437 /* mpichler, 19950515 */
438 /* cross product gives area of parallelogram, which is < 1.0 for
439 * non-perpendicular unit-length vectors; so normalize x, y here
440 */
441
442 mag = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
443 if (mag) {
444 x[0] /= mag;
445 x[1] /= mag;
446 x[2] /= mag;
447 }
448
449 mag = sqrt(y[0] * y[0] + y[1] * y[1] + y[2] * y[2]);
450 if (mag) {
451 y[0] /= mag;
452 y[1] /= mag;
453 y[2] /= mag;
454 }
455 #define M(row,col) m[col*4+row]
456 M(0, 0) = x[0];
457 M(0, 1) = x[1];
458 M(0, 2) = x[2];
459 M(0, 3) = 0.0;
460 M(1, 0) = y[0];
461 M(1, 1) = y[1];
462 M(1, 2) = y[2];
463 M(1, 3) = 0.0;
464 M(2, 0) = z[0];
465 M(2, 1) = z[1];
466 M(2, 2) = z[2];
467 M(2, 3) = 0.0;
468 M(3, 0) = 0.0;
469 M(3, 1) = 0.0;
470 M(3, 2) = 0.0;
471 M(3, 3) = 1.0;
472 #undef M
473
474 if (invert_matrix_general(m, inv))
475 glMultMatrixd(inv);
476
477 }
478
479
SmoothTrans(double cam[3],double eye[3],double up[3],double tocam[3],double toeye[3],double toup[3],double mark)480 int SmoothTrans(double cam[3], double eye[3], double up[3],
481 double tocam[3], double toeye[3], double toup[3],
482 double mark)
483 {
484 double v[3], d, trans_speed;
485
486
487
488 d = 0.0;
489 trans_speed = 20.0;
490
491 if (!paused)
492 switch (cmode) {
493 default:
494 case 0:
495 case 1:
496 break;
497 case 2:
498 case 3:
499 ADDVECTORS(cam, cam, planets[currtarget].vel);
500 break;
501 }
502
503 SUBVECTORS(v, tocam, cam);
504 d += MODULE(v);
505 DIVVECTOR(v, v, trans_speed);
506 ADDVECTORS(cam, cam, v);
507
508 SUBVECTORS(v, toeye, eye);
509 d += MODULE(v);
510 DIVVECTOR(v, v, trans_speed);
511 ADDVECTORS(eye, eye, v);
512
513 SUBVECTORS(v, toup, up);
514 d += MODULE(v);
515 DIVVECTOR(v, v, trans_speed);
516 ADDVECTORS(up, up, v);
517
518 if (d <= mark)
519 return 0;
520 else
521 return 1;
522 }
523
texture_LOD(GLubyte * image,int * width,int * height,int components)524 GLubyte *texture_LOD(GLubyte * image, int *width, int *height,
525 int components)
526 {
527 int i, j, k, nw, nh;
528
529 switch (LOD) {
530 case LOW:
531 nh = *height / 4;
532 nw = *width / 4;
533 for (j = 0; j < nh; j++)
534 for (i = 0; i < nw; i++)
535 for (k = 0; k < components; k++)
536 image[j * nw * components + i * components + k] =
537 image[j * 4 * (*width) * components +
538 i * 4 * components + k];
539 *height /= 4;
540 *width /= 4;
541 break;
542 case MEDIUM:
543 nh = *height / 2;
544 nw = *width / 2;
545 for (j = 0; j < nh; j++)
546 for (i = 0; i < nw; i++)
547 for (k = 0; k < components; k++)
548 image[j * nw * components + i * components + k] =
549 image[j * 2 * (*width) * components +
550 i * 2 * components + k];
551 *height /= 2;
552 *width /= 2;
553 break;
554 case HIGH:
555 default:
556 break;
557 }
558 return image;
559 }
560
log(char * message)561 void log(char *message)
562 {
563 if (logfile) {
564 fprintf(logfile, "%s", message);
565 fflush(logfile);
566 }
567 }
568
shutdown(int status)569 void shutdown(int status)
570 {
571 fclose(logfile);
572 exit(status);
573 }
574
575