台形則による数値積分とOpenGLとアフィン変換

OpenGLで台形則を使った積分を描いてみた。

ところでOpenGLの視点設定が覚えられない。難しすぎる。いつも視点移動させずに、座標を移動させたり縮小させたりしてデフォルトの範囲内におさまるようにしている。

#include<stdio.h>
#include<GL/glut.h>

// 自前のアフィン変換に使う
double m[3][3] = {{0.7, 0, -.5}, {0, 0.7, -.5}, {0,0,1}};
double v[3] = {0, 0, 1};

double MulMatVec()
{
	int i, j, k;
	for(i=0; i<3; i++){
		double x=0;
		for(j=0; j<3; j++){
			x += m[i][j] * v[j];
		}
		v[i] = x;
	}
}
double vertex(double x, double y)
{
	v[0] = x;
	v[1] = y;
	// アフィン変換
	MulMatVec();
	glVertex2d(v[0], v[1]);
}

double f(double x)
{
	return x*x*x - 3./2*x*x +1;
}


double integration(double a, double b)
{
	// 下限a,上限bの積分
	// assert a < b
	int i;
	int d=4; // 分割数
	double h = (b-a)/d;

	glColor3f(1,0,0);
	glBegin(GL_LINES);
	for(i=0; i<d; i++){
		double an = a+ h*i;
		vertex(an, 0); vertex(an, f(an));
		vertex(an, f(an)); vertex(an+h, f(an+h));
	}
	vertex(b, 0); vertex(b, f(b));
	glEnd();
}

void display()
{
	double x;
	glClear(GL_COLOR_BUFFER_BIT);
		
	glColor3f(1,1,1);
	glBegin(GL_LINES);
	// 軸の描画
	vertex(-1000, 0); vertex(1000, 0);
	vertex(0, -1000); vertex(0, 1000);
	glEnd();

	// グラフの描画
	glBegin(GL_POINTS);
	for(x=-10; x<10; x+=1e-3)
		vertex(x, f(x));
	glEnd();

	integration(0, 1.5);
}

int main(int argc, char **argv)
{
	glutInit(&argc, argv);
	glutCreateWindow("Integration");
	glutDisplayFunc(display);
	glutMainLoop();
	return 0;
}

OpenGL(R) Programming Guide: The Official Guide to Learning OpenGL(R), Version 2 (5th Edition)