开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > GIS编程 > 查看源码
sunlight.c
资源名称：glut23-src.zip　[点击查看]
上传用户：xk288cn
上传日期：2007-05-28
资源大小：4876k
文件大小：13k
源码类别：

GIS编程

开发平台：

Visual C++

sunlight.c：源码内容
  1. /* sunlight.c
  2.  * 
  3.  * To compile:
  4.  * cc -o sunlight sunlight.c -lglut -lGLU -lGL -lXmu -lXext -lX11 -lm
  5.  *
  6.  * Usage: sunlight
  7.  *
  8.  * Requires: globe.raw
  9.  *
  10.  * This program demonstrates the sun lighting the earth, displayed
  11.  * either as a globe or a flat map.  The time of day and day of year
  12.  * can be adjusted interactively, and the globe can be viewed from
  13.  * any angle.
  14.  *
  15.  * The lighting is done using OpenGL's lighting.  The flat map is
  16.  * lighted by using the normals of globe.  This is a rather unique
  17.  * use of normals, making a flat surface lit as if it were round.
  18.  *
  19.  * The left mouse adjusts the globe orientation, the time of day, or
  20.  * the time of year.  The right mouse activates a menu.
  21.  *
  22.  * Chris Schoeneman - 9/6/98
  23.  * crs23@bigfoot.com
  24.  * http://www.bigfoot.com/~crs23/
  25.  *        
  26.  */
  28. #include <stdlib.h>
  29. #include <stdio.h>
  30. #include <math.h>
  31. #include <GL/glut.h>
  33. #ifndef M_PI
  34. #define M_PI 3.14159265358979323846
  35. #endif
  37. /* menu constants */
  38. #define MENU_SHOW_GLOBE 1
  39. #define MENU_SHOW_MAP 2
  40. #define MENU_ADJUST_GLOBE 3
  41. #define MENU_ADJUST_DAY 4
  42. #define MENU_ADJUST_TIME 5
  43. #define MENU_QUIT 6
  45. static int mapMode = 0;
  46. static int adjustMode = 0;
  47. static GLdouble aspectRatio = 1.0;
  48. static double timeOfDay = 0.0;
  49. static double timeOfYear = 0.0;
  51. /*
  52.  * sun direction functions
  53.  */
  55. static const double radPerDeg = M_PI / 180.0;
  56. static const double radPerHour = M_PI / 12.0;
  57. static const double siderealHoursPerHour = 1.002737908;
  58. static const double epoch = 2415020.0;
  60. static double getGreenwichSideral(double julianDay)
  61. {
  62.   double jdMidnight;
  63.   double dayFraction;
  64.   double T; /* time (fractions of a century) */
  65.   double greenwichMidnight;
  67.   /* true position requires sidereal time of midnight at prime meridian.
  68.      get midnight of given julian day (midnight has decimal of .5). */
  69.   jdMidnight = floor(julianDay);
  70.   if (julianDay - jdMidnight >= 0.5) jdMidnight += 0.5;
  71.   else jdMidnight -= 0.5;
  73.   T = (jdMidnight - epoch) / 36525.0;
  75.   /* get fraction of a day */
  76.   dayFraction = (julianDay - jdMidnight) * 24.0;
  78.   /* get Greenwich midnight sidereal time (in hours) */
  79.   greenwichMidnight =
  80. fmod((0.00002581 * T + 2400.051262) * T + 6.6460656, 24.0);
  82.   /* return Greenwich sidereal time */
  83.   return radPerHour * (greenwichMidnight + dayFraction * siderealHoursPerHour);
  84. }
  86. static void getTruePosition(double julianDay,
  87. float latitude, float longitude,
  88. double sx, double sy, double sz,
  89. float pos[3])
  90. {
  91.   double tx, ty, tz;
  92.   double localSidereal;
  94.   /* convert to radians */
  95.   latitude *= radPerDeg;
  96.   longitude *= radPerDeg;
  98.   /* get local sidereal time */
  99.   localSidereal = getGreenwichSideral(julianDay) - longitude;
  101.   /* rotate around polar axis (y-axis) by local sidereal time */
  102.   tx = sx * cos(localSidereal) - sz * sin(localSidereal);
  103.   ty = sy;
  104.   tz = sz * cos(localSidereal) + sx * sin(localSidereal);
  106.   /* rotate by latitude to local position */
  107.   pos[0] = (float)tx;
  108.   pos[1] = (float)(ty * cos(latitude) - tz * sin(latitude));
  109.   pos[2] = (float)(tz * cos(latitude) + ty * sin(latitude));
  110. }
  112. static void getSunDirection(double julianDay, float latitude,
  113. float longitude, float pos[3])
  114. {
  115.   double C;
  116.   double T; /* time (fractions of a century) */
  117.   double sx, sy, sz; /* sun direction */
  118.   double obliquity; /* earth's tilt */
  119.   double meanAnomaly;
  120.   double geometricMeanLongitude;
  121.   double trueLongitude;
  123.   T = (julianDay - epoch) / 36525.0;
  125.   obliquity = radPerDeg *
  126. (23.452294 + (-0.0130125 + (-0.00000164 +
  127. 0.000000503 * T) * T) * T);
  129.   meanAnomaly = radPerDeg * (358.47583 +
  130. ((0.0000033 * T + 0.000150) * T + 35999.04975) * T);
  132.   geometricMeanLongitude = radPerDeg *
  133. ((0.0003025 * T + 36000.76892) * T + 279.69668);
  135.   C = radPerDeg *
  136. (sin(meanAnomaly) * (1.919460 - (0.004789 + 0.000014 * T) * T) +
  137. sin(2.0 * meanAnomaly) * (0.020094 - 0.0001 * T) +
  138. sin(3.0 * meanAnomaly) * 0.000293);
  140.   trueLongitude = geometricMeanLongitude + C;
  142.   /* position of sun if earth didn't rotate: */
  143.   sx = sin(trueLongitude) * cos(obliquity);
  144.   sy = sin(trueLongitude) * sin(obliquity);
  145.   sz = cos(trueLongitude);
  147.   /* get true position */
  148.   getTruePosition(julianDay, latitude, longitude, sx, sy, sz, pos);
  149. }
  151. /*
  152.  * globe/map calculation
  153.  */
  155. #define LAT_GRID 20
  156. #define LON_GRID 40
  157. static GLfloat sphere[LAT_GRID + 1][LON_GRID + 1][3];
  158. static GLfloat map[LAT_GRID + 1][LON_GRID + 1][2];
  159. static GLfloat uv[LAT_GRID + 1][LON_GRID + 1][2];
  161. static void initSphere(void)
  162. {
  163.   int lat, lon;
  165.   for (lat = 0; lat <= LAT_GRID; lat++) {
  166.     const float y = (float)-cos(M_PI * (double)lat / (double)LAT_GRID);
  167.     const float r = (float)sin(M_PI * (double)lat / (double)LAT_GRID);
  168.     for (lon = 0; lon <= LON_GRID; lon++) {
  169.       sphere[lat][lon][0] = r * (float)sin(-M_PI + 2.0 * M_PI *
  170. (double)lon / (double)LON_GRID);
  171.       sphere[lat][lon][1] = y;
  172.       sphere[lat][lon][2] = r * (float)cos(-M_PI + 2.0 * M_PI *
  173. (double)lon / (double)LON_GRID);
  174.     }
  175.   }
  176. }
  178. static void initMap(void)
  179. {
  180.   int lat, lon;
  182.   for (lat = 0; lat <= LAT_GRID; lat++) {
  183.     const float y = 1.0f * ((float)lat / (float)LAT_GRID - 0.5f);
  184.     for (lon = 0; lon <= LON_GRID; lon++) {
  185.       map[lat][lon][0] = 2.0f * ((float)lon / (float)LON_GRID - 0.5f);
  186.       map[lat][lon][1] = y;
  187.     }
  188.   }
  189. }
  191. static int isPowerOfTwo(int x)
  192. {
  193.   if (x <= 0) return 0;
  194.   while ((x & 1) == 0) x >>= 1;
  195.   return x == 1;
  196. }
  198. static void initTexture(const char* filename)
  199. {
  200.   FILE* imageFile;
  201.   int lat, lon;
  202.   unsigned char buffer[4];
  203.   unsigned char* image;
  204.   int dx, dy;
  206.   /* initialize texture coordinates */
  207.   for (lat = 0; lat <= LAT_GRID; lat++) {
  208.     const float y = (float)lat / (float)LAT_GRID;
  209.     for (lon = 0; lon <= LON_GRID; lon++) {
  210.       uv[lat][lon][0] = (float)lon / (float)LON_GRID;
  211.       uv[lat][lon][1] = y;
  212.     }
  213.   }
  215.   /* open image file */
  216.   imageFile = fopen(filename, "rb");
  217.   if (!imageFile) {
  218.     fprintf(stderr, "cannot open image file %s for readingn", filename);
  219.     exit(1);
  220.   }
  222.   /* read image size */
  223.   fread(buffer, 1, 4, imageFile);
  224.   dx = (int)((buffer[0] << 8) | buffer[1]);
  225.   dy = (int)((buffer[2] << 8) | buffer[3]);
  226.   if (!isPowerOfTwo(dx) || !isPowerOfTwo(dy)) {
  227.     fprintf(stderr, "image %s has an illegal size: %dx%dn", filename, dx, dy);
  228.     exit(1);
  229.   }
  231.   /* read image */
  232.   image = (unsigned char*)malloc(3 * dx * dy * sizeof(unsigned char));
  233.   fread(image, 3 * dx, dy, imageFile);
  235.   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  236.   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  237.   glTexImage2D(GL_TEXTURE_2D, 0, 3, dx, dy, 0,
  238. GL_RGB, GL_UNSIGNED_BYTE, image);
  240.   /* tidy up */
  241.   free(image);
  242.   fclose(imageFile);
  243. }
  245. /*
  246.  * misc OpenGL stuff
  247.  */
  249. static const GLfloat sunColor[] = { 1.0f, 1.0f, 1.0f, 1.0f };
  250. static GLfloat sunDirection[4] = { 0.0f, 0.0f, 1.0f, 0.0f };
  251. static const GLfloat surfaceColor[] = { 0.8f, 0.8f, 0.8f, 1.0f };
  253. static void initProjection(void)
  254. {
  255.   glMatrixMode(GL_PROJECTION);
  256.   glLoadIdentity();
  258.   if (mapMode) {
  259.     /* orthographic for map */
  260.     if (aspectRatio < 2.0)
  261.       glOrtho(-1.0, 1.0, -1.0 / aspectRatio, 1.0 / aspectRatio, 2.0, 4.0);
  262.     else
  263.       glOrtho(-0.5 * aspectRatio, 0.5 * aspectRatio, -0.5, 0.5, 2.0, 4.0);
  264.   }
  266.   else {
  267.     /* prespective for globe */
  268.     gluPerspective(40.0, aspectRatio, 1.0, 5.0);
  269.   }
  271.   glMatrixMode(GL_MODELVIEW);
  272. }
  274. static void initSunlight(void)
  275. {
  276.   getSunDirection(epoch + timeOfDay + 365.0 * timeOfYear,
  277. 0.0f, 0.0f, sunDirection);
  279.   glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
  280.   glLightfv(GL_LIGHT0, GL_DIFFUSE, sunColor);
  281.   glLightfv(GL_LIGHT0, GL_POSITION, sunDirection);
  282.   glEnable(GL_LIGHT0);
  284.   glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, surfaceColor);
  285. }
  287. /*
  288.  * glut callbacks
  289.  */
  291. static GLfloat angle = 0; /* in degrees */
  292. static GLfloat angle2 = 0; /* in degrees */
  293. static int moving = 0, startx, starty;
  295. static void reshapeCB(GLsizei w, GLsizei h)
  296. {
  297.   aspectRatio = (GLdouble)w / (GLdouble)h;
  298.   glViewport(0, 0, w, h);
  299.   initProjection();
  300.   glutPostRedisplay();
  301. }
  303. static void redrawCB(void)
  304. {
  305.   glClear(GL_COLOR_BUFFER_BIT);
  307.   glPushMatrix();
  309.     /* Perform scene rotations based on user mouse input. */
  310.     if (!mapMode) {
  311.       glRotatef(angle2, 1.0f, 0.0f, 0.0f);
  312.       glRotatef(angle, 0.0f, 1.0f, 0.0f);
  313.     }
  315.     glLightfv(GL_LIGHT0, GL_POSITION, sunDirection);
  317.     if (mapMode) {
  318.       /* use the normals for a sphere on a flat surface to get
  319.        * the lighting as if we unwrapped a lighted sphere. */
  321.       int lat, lon;
  322.       for (lat = 0; lat < LAT_GRID; lat++) {
  323. glBegin(GL_TRIANGLE_STRIP);
  325.   for (lon = 0; lon <= LON_GRID; lon++) {
  326.     glTexCoord2fv(uv[lat + 1][lon]);
  327.     glNormal3fv(sphere[lat + 1][lon]);
  328.     glVertex2fv(map[lat + 1][lon]);
  330.     glTexCoord2fv(uv[lat][lon]);
  331.     glNormal3fv(sphere[lat][lon]);
  332.     glVertex2fv(map[lat][lon]);
  333.   }
  335. glEnd();
  336.       }
  337.     }
  339.     else {
  340.       /* draw a sphere */
  341.       int lat, lon;
  342.       for (lat = 0; lat < LAT_GRID; lat++) {
  343. glBegin(GL_TRIANGLE_STRIP);
  345.   for (lon = 0; lon <= LON_GRID; lon++) {
  346.     glTexCoord2fv(uv[lat + 1][lon]);
  347.     glNormal3fv(sphere[lat + 1][lon]);
  348.     glVertex3fv(sphere[lat + 1][lon]);
  350.     glTexCoord2fv(uv[lat][lon]);
  351.     glNormal3fv(sphere[lat][lon]);
  352.     glVertex3fv(sphere[lat][lon]);
  353.   }
  355. glEnd();
  356.       }
  357.     }
  359.   glPopMatrix();
  361.   glutSwapBuffers();
  362. }
  364. static void mouseCB(int button, int state, int x, int y)
  365. {
  366.   if (button == GLUT_LEFT_BUTTON) {
  367.     if (state == GLUT_DOWN) {
  368.       moving = 1;
  369.       startx = x;
  370.       starty = y;
  371.     }
  372.     else if (state == GLUT_UP) {
  373.       moving = 0;
  374.     }
  375.   }
  376. }
  378. static void motionCB(int x, int y)
  379. {
  380.   if (moving) {
  381.     switch (adjustMode) {
  382.       case 0:
  383. /* move globe */
  384. if (!mapMode) {
  385.   angle = angle + (x - startx);
  386.   angle2 = angle2 + (y - starty);
  387.   glutPostRedisplay();
  388. }
  389. break;
  391.       case 1:
  392. /* change day of year */
  393. timeOfYear = timeOfYear + (y - starty) / 365.0;
  394. while (timeOfYear < 0) timeOfYear += 1.0;
  395. while (timeOfYear >= 1.0) timeOfYear -= 1.0;
  396. getSunDirection(epoch + timeOfDay + 365.0 * timeOfYear,
  397. 0.0f, 0.0f, sunDirection);
  398. glutPostRedisplay();
  399. break;
  401.       case 2:
  402. /* change time of day (by 4 minute increments (24 hrs / 360)) */
  403. timeOfDay = timeOfDay + (y - starty) / 360.0;
  404. while (timeOfDay < 0) timeOfDay += 1.0;
  405. while (timeOfDay >= 1.0) timeOfDay -= 1.0;
  406. getSunDirection(epoch + timeOfDay + 365.0 * timeOfYear,
  407. 0.0f, 0.0f, sunDirection);
  408. glutPostRedisplay();
  409. break;
  410.     }
  412.     startx = x;
  413.     starty = y;
  414.   }
  415. }
  417. static void keyCB(unsigned char c, int x, int y)
  418. {
  419.   if (c == 27) {
  420.     exit(0);
  421.   }
  422.   glutPostRedisplay();
  423. }
  426. /*
  427.  * menu handling
  428.  */
  430. static void handleMenu(int value)
  431. {
  432.   switch (value) {
  433.     default:
  434.       break;
  436.     case MENU_SHOW_GLOBE:
  437.       mapMode = 0;
  438.       initProjection();
  439.       glutPostRedisplay();
  440.       glutChangeToMenuEntry(2, "Show Globe *", MENU_SHOW_GLOBE);
  441.       glutChangeToMenuEntry(3, "Show Map", MENU_SHOW_MAP);
  442.       break;
  444.     case MENU_SHOW_MAP:
  445.       mapMode = 1;
  446.       initProjection();
  447.       glutPostRedisplay();
  448.       glutChangeToMenuEntry(2, "Show Globe", MENU_SHOW_GLOBE);
  449.       glutChangeToMenuEntry(3, "Show Map *", MENU_SHOW_MAP);
  450.       break;
  452.     case MENU_ADJUST_GLOBE:
  453.       adjustMode = 0;
  454.       glutChangeToMenuEntry(4, "Adjust Globe *", MENU_ADJUST_GLOBE);
  455.       glutChangeToMenuEntry(5, "Adjust Day", MENU_ADJUST_DAY);
  456.       glutChangeToMenuEntry(6, "Adjust Time", MENU_ADJUST_TIME);
  457.       break;
  459.     case MENU_ADJUST_DAY:
  460.       adjustMode = 1;
  461.       glutChangeToMenuEntry(4, "Adjust Globe", MENU_ADJUST_GLOBE);
  462.       glutChangeToMenuEntry(5, "Adjust Day *", MENU_ADJUST_DAY);
  463.       glutChangeToMenuEntry(6, "Adjust Time", MENU_ADJUST_TIME);
  464.       break;
  466.     case MENU_ADJUST_TIME:
  467.       adjustMode = 2;
  468.       glutChangeToMenuEntry(4, "Adjust Globe", MENU_ADJUST_GLOBE);
  469.       glutChangeToMenuEntry(5, "Adjust Day", MENU_ADJUST_DAY);
  470.       glutChangeToMenuEntry(6, "Adjust Time *", MENU_ADJUST_TIME);
  471.       break;
  473.     case MENU_QUIT:
  474.       exit(0);
  475.   }
  476. }
  478. /*
  479.  * main
  480.  */
  482. int main(int argc, char** argv)
  483. {
  484.   /*
  485.    * initialize GLUT and open a window
  486.    */
  487.   glutInit(&argc, argv);
  488.   glutInitWindowSize(512, 512);
  489.   glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
  490.   glutCreateWindow("Sunlight");
  491.   glutDisplayFunc(redrawCB);
  492.   glutReshapeFunc(reshapeCB);
  493.   glutMouseFunc(mouseCB);
  494.   glutMotionFunc(motionCB);
  495.   glutKeyboardFunc(keyCB);
  497.   /*
  498.    * make the menu
  499.    */
  500.   glutCreateMenu(handleMenu);
  501.   glutAddMenuEntry("SUNLIGHT", 0);
  502.   glutAddMenuEntry("Show Globe *", MENU_SHOW_GLOBE);
  503.   glutAddMenuEntry("Show Map", MENU_SHOW_MAP);
  504.   glutAddMenuEntry("Adjust Globe *", MENU_ADJUST_GLOBE);
  505.   glutAddMenuEntry("Adjust Day", MENU_ADJUST_DAY);
  506.   glutAddMenuEntry("Adjust Time", MENU_ADJUST_TIME);
  507.   glutAddMenuEntry("Quit", MENU_QUIT);
  508.   glutAttachMenu(GLUT_RIGHT_BUTTON);
  510.   /*
  511.    * initialize GL
  512.    */
  513.   initProjection();
  514.   gluLookAt(0.0, 0.0, 3.0,
  515.     0.0, 0.0, 0.0,
  516.     0.0, 1.0, 0.0);
  517.   initSunlight();
  518.   glEnable(GL_LIGHTING);
  519.   glEnable(GL_CULL_FACE);
  520.   glEnable(GL_TEXTURE_2D);
  522.   /*
  523.    * initialize data structures
  524.    */
  525.   initSphere();
  526.   initMap();
  527.   initTexture("globe.raw");
  529.   glutMainLoop();
  530.   return 0;             /* ANSI C requires main to return int. */
  531. }