here the C codes:
/* Include files */
#include "SimpleSH.h"
#include "sh.h"
/* OpenGL display mode: points, triangles and fill */
int g_ViewMode;
/* Current indexes for the SH plotting */
int currentL, currentM;
/* Window position */
int iXWinPos, iYWinPos;
/* Rotation angles based on the mouse movements */
int iYAngle, iYBegin;
int iXAngle, iXBegin;
/* Window width and height */
int iDisplayWidth, iDisplayHeight;
/* Mouse movement status */
int iMovingStatus;
/* Animation flag */
int iAnimatingStatus;
/* Number of triangles for the current SH displayed */
int iResolution;
/* To scale up/down the SH plotting */
double scaleLevel;
/* program's main entry point */
int main(int argc, char** argv)
{
/* Default values */
/* Initialize l and m indexes */
currentL = currentM = 0;
/* Default resolution is 800 by 600*/
iDisplayWidth = 800;
iDisplayHeight = 600;
/* Rotation angles */
iYBegin = iYAngle = 0;
iXBegin = iXAngle = 0;
/* Mouse not moving by default */
iMovingStatus = 0;
/* No animation by default */
iAnimatingStatus = 0;
/* Initial window position */
iXWinPos = iYWinPos = 10;
/* Display wireframe model by default */
g_ViewMode = GL_LINE;
/* Default model resolution */
iResolution = 64;
/* Default scale resolution */
scaleLevel = 1.0;
/* Show credits */
showCredits();
/* Parse the command-line arguments */
parseArgs(argc, argv);
/* Glut initialization */
glutInit(&argc, argv);
/* Double buffering and RGB display */
glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGB);
/* Window size */
glutInitWindowSize(iDisplayWidth, iDisplayHeight);
/* Window initial position */
glutInitWindowPosition(iXWinPos, iYWinPos);
/* Window creations */
glutCreateWindow("SimpleSH");
/* Initialize OpenGL settings */
initOpenGL();
/* Assign GLUT call-back functions */
glutDisplayFunc(displayScene);
glutReshapeFunc(reshapeWindow);
glutMouseFunc(mouseInput);
glutMotionFunc(mouseMotion);
glutKeyboardFunc(keyboard);
/* Enter GLUT loop */
glutMainLoop();
}
/* Keyboard function */
void keyboard(unsigned char key, int x, int y)
{
switch (key)
{
/* Decrease the value of the l index */
case 'l':
if (currentL > 0)
{
if (abs(currentM) == currentL)
{
if (currentM < 0) currentM++;
if (currentM > 0) currentM--;
}
currentL--;
}
printf(" l = %d - m = %d\t\r", currentL, currentM);
glutPostRedisplay();
break;
/* Increase the value of the l index */
case 'L':
currentL++;
printf(" l = %d - m = %d\t\r", currentL, currentM);
glutPostRedisplay();
break;
/* Decrease the value of the m index */
case 'm':
if ((currentM > 0) || ((currentM <= 0) && (-currentM < currentL))) currentM--;
printf(" l = %d - m = %d\t\r", currentL, currentM);
glutPostRedisplay();
break;
/* Increase the value of the m index */
case 'M':
if (currentM < currentL) currentM++;
printf(" l = %d - m = %d\t\r", currentL, currentM);
glutPostRedisplay();
break;
/* Decrease the resolution of the model. Minimum value is 16 */
case 'r':
if (iResolution > 16) iResolution--;
glutPostRedisplay();
break;
/* Increase the resolution of the model */
case 'R':
iResolution++;
glutPostRedisplay();
break;
/* Decrease the scale of the model. Minimum value is 0.1 */
case 'z':
if (scaleLevel > 0.1) scaleLevel -= 0.1;
glutPostRedisplay();
break;
/* Increase the scale of the model */
case 'Z':
scaleLevel += 0.1;
glutPostRedisplay();
break;
/* Toggle animation on/off */
case 'a':
case 'A':
if (iAnimatingStatus == 0)
{
glutIdleFunc(animateScene) ;
iAnimatingStatus = 1 ;
}
else
{
glutIdleFunc(NULL) ;
iAnimatingStatus = 0 ;
}
break ;
/* Switch between the different display modes: points, triangles and fill */
case 'y':
case 'Y':
if (g_ViewMode == GL_FILL) g_ViewMode = GL_LINE;
else if (g_ViewMode == GL_LINE) g_ViewMode = GL_POINT;
else if (g_ViewMode == GL_POINT) g_ViewMode = GL_FILL;
else g_ViewMode = GL_FILL;
glPolygonMode(GL_FRONT_AND_BACK, g_ViewMode);
glutPostRedisplay();
break ;
/* Quit */
case 'q':
case 'Q':
case '\27':
puts("");
exit(0);
break;
}
}
/* animates the scene */
void animateScene(void)
{
iYAngle += 2;
if (iYAngle >= 360) iYAngle = 0;
glutPostRedisplay();
}
/* displays the scene */
void displayScene(void)
{
/* The color of the SH */
GLfloat glfMaterialColor[] = { 0.0f, 0.0f, 1.0f, 1.0f };
/* Clear the screen */
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
/* Setup the view */
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(60, 1.3333, 0.01, 30);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0, -3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
/* Rotate the camera around the origin */
glRotatef(iXAngle, 1.f, 0.f, 0.f);
glRotatef(iYAngle, 0.f, 0.f, 1.f);
/* Assign material properties and render the SH for l and m indexes */
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, glfMaterialColor);
glEnable(GL_COLOR_MATERIAL);
renderSH(iResolution, 0.0, 0.0, 0.0, currentL, currentM, scaleLevel);
/* Swap the front and back buffers */
glutSwapBuffers();
}
/* Called when the window is reshaped */
void reshapeWindow(int iWidth, int iHeight)
{
/* Local variables */
GLdouble gldAspect;
gldAspect = (GLdouble) iWidth / (GLdouble) iHeight;
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(30.0, gldAspect, 0.01, 50.0);
glViewport(0, 0, (GLsizei) iWidth, (GLsizei) iHeight);
}
/* Handles mouse input */
void mouseInput(int iButton, int iState, int iX, int iY)
{
switch (iButton)
{
/* Left button: camera rotation around the origin */
case GLUT_LEFT_BUTTON:
if (iState == GLUT_DOWN)
{
iMovingStatus = 1;
iYBegin = iY;
iXBegin = iX;
}
break;
/* Right button: switch on/off animation */
case GLUT_RIGHT_BUTTON:
if (iState == GLUT_DOWN)
{
if (iAnimatingStatus == 0)
{
glutIdleFunc(animateScene) ;
iAnimatingStatus = 1 ;
}
else
{
glutIdleFunc(NULL) ;
iAnimatingStatus = 0 ;
}
}
break;
default:
break;
}
}
/* handles mouse input */
void mouseMotion(int iX, int iY)
{
/* Calculate the camera rotation angles based on the mouse movement */
if (iMovingStatus == 1)
{
iYAngle += (iX - iXBegin);
iXAngle += (iY - iYBegin);
iYBegin = iY;
iXBegin = iX;
glutPostRedisplay();
}
}
/* sets the OpenGL environment */
void initOpenGL(void)
{
/* Enable depth testing */
glEnable(GL_DEPTH_TEST);
/* Allow color */
glEnable(GL_COLOR_MATERIAL);
/* Enable Smooth Shading */
glShadeModel(GL_SMOOTH);
/* Enable rescale of normals after transform */
glEnable(GL_NORMALIZE);
/* Really Nice Perspective Calculations */
glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
/* Rendering mode */
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
/* Add a light to the scene */
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
}
/*
Render the SH for the given l and m indexes.
Scale is used to scale the SH.
iSteps is used asthe resolution of the SH.
*/
void renderSH(int iSteps, double x, double y, double z, int l, int m, double scale)
{
/* Loop indexes */
int i, j;
/* Angles in spherical coordinates */
double theta, phi;
/* Angles variations in spherical coordinates */
double step_theta, step_phi;
/* The radius is a function of phi and theta */
double R;
/* Initialize angles variations */
step_theta = (PI) / iSteps;
step_phi = (2 * PI) / iSteps;
glPushMatrix();
/* Translate the SH to (x,y,z) */
glTranslatef(x, y, z);
/* Uniform scalig of the SH */
glScalef(scale, scale, scale);
/* Loop around the unit sphere in spherical coordinates */
for (i = 0; i < iSteps; i++)
{
/* Current value of theta */
theta = step_theta * i;
/* Loop for phi */
for (j = 0; j < iSteps; j++)
{
/* Current value of theta */
phi = step_phi * j;
/* Draw quads */
glBegin(GL_QUADS);
/* First vertex */
R = evaluateSH(l, m, theta, phi);
glNormal3f(R * sin(phi) * cos(theta), R * sin(phi) * sin(theta), R * cos(phi));
glVertex3f(R * sin(phi) * cos(theta), R * sin(phi) * sin(theta), R * cos(phi));
/* Second vertex */
R = evaluateSH(l, m, theta + step_theta, phi);
glNormal3f(R * sin(phi) * cos(theta + step_theta), R * sin(phi) * sin(theta + step_theta), R * cos(phi));
glVertex3f(R * sin(phi) * cos(theta + step_theta), R * sin(phi) * sin(theta + step_theta), R * cos(phi));
/* Third vertex */
R = evaluateSH(l, m, theta + step_theta, phi + step_phi);
glNormal3f(R * sin(phi + step_phi) * cos(theta + step_theta), R * sin(phi + step_phi) * sin(theta + step_theta), R * cos(phi + step_phi));
glVertex3f(R * sin(phi + step_phi) * cos(theta + step_theta), R * sin(phi + step_phi) * sin(theta + step_theta), R * cos(phi + step_phi));
/* Fourth vertex */
R = evaluateSH(l, m, theta, phi + step_phi);
glNormal3f(R * sin(phi + step_phi) * cos(theta), R* sin(phi + step_phi) * sin(theta), R * cos(phi + step_phi));
glVertex3f(R * sin(phi + step_phi) * cos(theta), R* sin(phi + step_phi) * sin(theta), R * cos(phi + step_phi));
glEnd();
}
}
glPopMatrix();
}
/* Display the usage of the program */
void showUsage(void)
{
printf("\n Usage: SimpleSH [options]\n\n");
printf(" Example: SimpleSH -l 2 -m 2\n\n");
printf(" Options:\n");
printf(" -w/--width integer width of the display window\n");
printf(" -h/--height integer height of the display window\n");
printf(" -l/--l integer Band level\n");
printf(" -m/--m integer Secondary level\n");
printf(" -?/--help this help message\n\n");
}
/* Show the different credits */
void showCredits(void)
{
printf(VERSION_CREDITS);
puts("");
}
/* Show the help */
void help(void)
{
showUsage();
exit(-1);
}
/* Parse the command-line arguments and assign global variables accordingly */
void parseArgs(int argc, char **argv)
{
int i;
for (i=1; i<argc; ++i)
{
if (argv[i][0]=='-')
{
if ((strcmp(argv[i],"-?")==0) || (strcmp(argv[i],"--help")==0) || (strcmp(argv[i],"/?")==0))
{
help();
}
else if ((strcmp(argv[i],"-w")==0) || (strcmp(argv[i],"--width")==0))
{
++i;
if (i >= argc)
{
help();
}
if ((iDisplayWidth = atoi(argv[i])) <= 0)
{
printf(" Display width must be strictly positive.\n");
exit(-1);
}
}
else if ((strcmp(argv[i],"-l")==0) || (strcmp(argv[i],"--l")==0))
{
++i;
if (i >= argc)
{
help();
}
if ((currentL = atoi(argv[i])) < 0)
{
printf(" l must be positive.\n");
exit(-1);
}
}
else if ((strcmp(argv[i],"-m")==0) || (strcmp(argv[i],"--m")==0))
{
++i;
if (i >= argc)
{
help();
}
currentM = atoi(argv[i]);
if (abs(currentM) > currentL)
{
printf(" l must be greater or equal to m.\n");
exit(-1);
}
}
else if ((strcmp(argv[i],"-h")==0) || (strcmp(argv[i],"--height")==0))
{
++i;
if (i >= argc)
{
help();
}
if ((iDisplayHeight = atoi(argv[i])) <= 0)
{
printf(" Display height must be strictly positive.\n");
exit(-1);
}
}
}
}
}
/* Required include file */
#include "sh.h"
/* Calculate the double factorial of x */
int doubleFactorial(int x)
{
int result;
if (x == 0 || x == -1) return (1);
result = x;
while ((x -= 2) > 0) result *= x;
return (result);
}
/* Calculate the factorial of x */
int factorial(int x)
{
int result;
if (x == 0 || x == -1) return (1);
result = x;
while ((x -= 1) > 0) result *= x;
return (result);
}
/* Calculate the Associated Legendre Polynomial of index l and m at point x*/
double ALPStd(double x, int l, int m)
{
if (l == m) return (pow(-1, m) * doubleFactorial(2 * m - 1) * pow(sqrt(1 - x * x), m));
if (l == m + 1) return (x * (2 * m + 1) * ALPStd(x, m, m));
return ((x * (2 * l - 1) * ALPStd(x, l - 1, m) - (l + m - 1) * ALPStd(x, l - 2, m)) / (l - m));
}
/* Calculate the K constant used in the spherical harmonic calculation */
double evaluateK(int l, int m)
{
double result;
result = ((2.0 * l + 1.0) * factorial(l - m)) / (4 * PI * factorial(l + m));
return sqrt((result));
}
/* Calculate the spherical harmonic for l and m at angles theta and phi*/
double evaluateSH(int l, int m, double theta, double phi)
{
double SH = 0.0;
if (m == 0)
{
SH = evaluateK(l, 0) * ALPStd(cos(theta), l, 0);
}
else if (m > 0)
{
SH = SQRT2 * evaluateK(l, m) * cos(m * phi) * ALPStd(cos(theta), l, m);
}
else
{
SH = SQRT2 * evaluateK(l, -m) * sin(-m * phi) * ALPStd(cos(theta), l, -m);
}
return SH;
}
