// A.4 marching-Square Program /* Generates Contours Using marching squares * // * region size */# define x_max 1.0 # define y_max 1.0 # define x_min-1.0 # define y_min-1.0/* Number of cells */ # define N_x 50 # define n_y 50/* contour value */# define threshold 0.0 # include <Gl/glut. h> void display () {double F (double, double); int cell (Double, double); void lines (INT, Int, Int, Double, double, double, double, double); double data [N_x] [n_y]; int I, j; int C; glclear (gl_color_buffer_bit ); /* form data array from function */for (I = 0; I <N_x; I ++) for (j = 0; j <n_y; j ++) data [I] [J] = f (x_min + I * (X_MAX-X_MIN)/(N_X-1.0), y_min + J * (Y_MAX-Y_MIN)/(N_Y-1.0 )); /* process each cell */for (I = 0; I <N_x; I ++) for (j = 0; j <n_y; j ++) {c = cell (data [I] [J], data [I + 1] [J], data [I + 1] [J + 1], data [I] [J + 1]); lines (C, I, j, data [I] [J], data [I + 1] [J], data [I + 1] [J + 1], data [I] [J + 1]);} glflush () ;}/ * define function f (x, y) */double F (Double X, Double Y) {double A = 0.49, B = 0.5; /* ovals of Cassini */Return (x * x + y * Y + A * A) * (x * x + y * Y + A *) -4 * a * x * X-B * B * B;}/* define cell vertices */INT cell (double A, double B, double C, double D) {int n = 0; if (a> threshold) N + = 1; if (B> threshold) N + = 8; If (C> threshold) N + = 4; If (D> threshold) N + = 2; return N;}/* draw line segments for each case */void lines (INT num, int I, int J, double A, double B, double C, double D) {void draw_one (INT, Int, Int, Double, double); void draw_adjacent (INT, int, Int, Double, double, double); void draw_opposite (INT, Int, Int, Double, double); Switch (Num) {Case 1: case 2: Case 4: Case 7: Case 8: Case 11: Case 13: Case 14: draw_one (Num, I, j, A, B, C, D); break; case 3: Case 6: Case 9: Case 12: draw_adjacent (Num, I, j, A, B, C, D); break; Case 5: Case 10: draw_opposite (Num, I, j, A, B, C, D); break; Case 0: Case 15: Break ;}} void draw_one (INT num, int I, int J, double A, double B, double C, double D) {double X1, Y1, X2, Y2; double ox, oy; double dx, Dy; dx = (X_MAX-X_MIN) /(N_X-1.0); DY = (Y_MAX-Y_MIN)/(N_Y-1.0); ox = x_min + I * (X_MAX-X_MIN)/(N_X-1.0); Oy = y_min + J * (Y_MAX-Y_MIN) /(N_Y-1.0); Switch (Num) {Case 1: Case 14: X1 = ox; Y1 = Oy + dy * (threshold-A)/(d-); x2 = ox + dx * (threshold-A)/(B-A); y2 = Oy; break; Case 2: Case 13: X1 = ox; y1 = Oy + dy * (threshold-A)/(d-A); x2 = ox + dx * (threshold-d)/(c-d ); y2 = Oy + dy; break; Case 4: Case 11: X1 = ox + dx * (threshold-d)/(c-d); Y1 = Oy + dy; x2 = ox + dx; y2 = Oy + dy * (threshold-B)/(c-B); break; Case 7: Case 8: x1 = ox + dx * (threshold-A)/(B-A); Y1 = Oy; x2 = ox + dx; y2 = Oy + dy * (threshold-B) /(c-B); break;} glbegin (gl_lines); glvertex2d (x1, Y1); glvertex2d (X2, Y2); glend ();} void draw_adjacent (INT num, int I, Int J, double A, double B, double C, double D) {double X1, Y1, X2, Y2; double ox, oy; double dx, Dy; DX = (X_MAX-X_MIN)/(N_X-1.0); DY = (Y_MAX-Y_MIN)/(N_Y-1.0); ox = x_min + I * (X_MAX-X_MIN)/(N_X-1.0 ); oy = y_min + J * (Y_MAX-Y_MIN)/(N_Y-1.0); Switch (Num) {Case 3: Case 12: X1 = ox + dx * (threshold-) /(B-A); Y1 = Oy; x2 = ox + dx * (threshold-d)/(c-d); y2 = Oy + dy; break; Case 6: case 9: X1 = ox; Y1 = Oy + dy * (threshold-A)/(d-A); x2 = ox + dx; y2 = Oy + dy * (threshold-B)/(c-B); break;} glbegin (gl_lines); glvertex2d (x1, Y1); glvertex2d (X2, Y2 ); glend ();} void draw_opposite (INT num, int I, Int J, double A, double B, double C, double D) {double X1, Y1, X2, Y2, x3, Y3, X4, Y4; double ox, oy; double dx, Dy; dx = (X_MAX-X_MIN)/(N_X-1.0); DY = (Y_MAX-Y_MIN)/(N_Y-1.0 ); ox = x_min + I * (X_MAX-X_MIN)/(N_X-1.0); Oy = y_min + J * (Y_MAX-Y_MIN)/(N_Y-1.0); Switch (Num) {Case 5: X1 = ox; y1 = Oy + dy * (threshold-A)/(d-A); x2 = ox + dx * (threshold-A)/(B-A); y2 = Oy; x3 = ox + dx * (threshold-d)/(c-d); Y3 = Oy + dy; X4 = ox + dx; y4 = Oy + dy * (threshold-B)/(c-B); break; case 10: X1 = ox; Y1 = Oy + dy * (threshold-) /(d-A); x2 = ox + dx * (threshold-d)/(c-d); y2 = Oy + dy; x3 = ox + dx * (threshold-d)/(c-d); Y3 = Oy; X4 = ox + dx; Y4 = Oy + dy * (threshold-B) /(c-B); break;} glbegin (gl_lines); glvertex2d (x1, Y1); glvertex2d (X2, Y2); glvertex2d (X3, Y3); glvertex2d (X4, y4); glend ();} void myreshape (int w, int h) {glviewport (0, 0, W, H); glmatrixmode (gl_projection); glloadidentity (); if (W <= h) gluortho2d (x_min, x_max, y_min * (glfloat) h/(glfloat) W, y_max * (glfloat) h/(glfloat) W ); else gluortho2d (x_min * (glfloat) W/(glfloat) h, x_max * (glfloat) W/(glfloat) h, y_min, y_max); glmatrixmode (gl_modelview );} void main (INT argc, char ** argv) {gluinit (& argc, argv); gluinitwindowsize (500,500); glucreatewindow ("Contour Plot"); glureshapefunc (myreshape ); gludisplayfunc (Display); glclearcolor (0.0, 0.0, 0.0, 1.0); glcolor3f (1.0, 1.0, 1.0); glumainloop ();}