3654. 陰関数表示プログラム。C言語ソース(2/2)

void draw1(void) {
int ix0, iy0, iz0;
int j;
double x0, y0, z0, x1, y1, z1;
double x2, y2, z2, x3, y3, z3;
double* p0,* q0;

// clear color buufer and z buffer
for (iy0 = 0; iy0 < HT; iy0++) {
j = iy0 * WT;
for (ix0 = 0; ix0 < WT; ix0++) {
zb[j + ix0] = -1;
b0[(j + ix0) * 3] = cl[9][0];
b0[(j + ix0) * 3 + 1] = cl[9][1];
b0[(j + ix0) * 3 + 2] = cl[9][2];
}
}

// v01[] initial value
iz0 = 0; z0 = ka * (double)(iz0 * 8 - DT / 2);
x2 = uz[0] * z0;
y2 = uz[1] * z0;
z2 = uz[2] * z0;
q0 = v01;
for (iy0 = 0; iy0 <= HT / 8; iy0++) {
y0 = ka * (double)(iy0 * 8 - HT / 2);
x3 = x2 + uy[0] * y0;
y3 = y2 + uy[1] * y0;
z3 = z2 + uy[2] * y0;
for (ix0 = 0; ix0 <= WT / 8; ix0++) {
x0 = ka * (double)(ix0 * 8 - WT / 2);
x1 = px + x3 + ux[0] * x0;
y1 = py + y3 + ux[1] * x0;
z1 = pz + z3 + ux[2] * x0;
*q0++ = func1(x1, y1, z1);
}
}
// outer loop to find v cross
for (iz0 = 0; iz0 < DT / 8; iz0++) {
// moving v00 <- v01
q0 = v00; p0 = v01;
for (iy0 = 0; iy0 <= HT / 8; iy0++) {
for (ix0 = 0; ix0 <= WT / 8; ix0++) {
*q0++ = *p0++;
}
}
// set the first position in the next plane
z0 = ka * (double)((iz0 + 1) * 8 - DT / 2);
x2 = uz[0] * z0;
y2 = uz[1] * z0;
z2 = uz[2] * z0;
p0 = v01;
// calculate
for (iy0 = 0; iy0 <= HT / 8; iy0++) {
y0 = ka * (double)(iy0 * 8 - HT / 2);
x3 = x2 + uy[0] * y0;
y3 = y2 + uy[1] * y0;
z3 = z2 + uy[2] * y0;
for (ix0 = 0; ix0 <= WT / 8; ix0++) {
x0 = ka * (double)(ix0 * 8 - WT / 2);
x1 = px + x3 + ux[0] * x0;
y1 = py + y3 + ux[1] * x0;
z1 = pz + z3 + ux[2] * x0;
*p0++ = func1(x1, y1, z1);
}
}
// compare
for (iy0 = 0; iy0 < HT / 8; iy0++) {
j = iy0 * W81;
for (ix0 = 0; ix0 < WT / 8; ix0++) {
if (v00[j + ix0] > 0) {
if (!((v00[j + ix0 + 1] > 0) &&
(v00[j + W81 + ix0] > 0) && (v00[j + W81 + ix0 + 1] > 0) &&
(v01[j + ix0] > 0) && (v01[j + ix0 + 1] > 0) &&
(v01[j + W81 + ix0] > 0) && (v01[j + W81 + ix0 + 1] > 0))
) draw2(ix0 * 8 - WT / 2, iy0 * 8 - HT / 2, iz0 * 8 - DT / 2);
} else { // (v00[j + ix0] <= 0)
if (!((v00[j + ix0 + 1] <= 0) &&
(v00[j + W81 + ix0] <= 0) && (v00[j + W81 + ix0 + 1] <= 0) &&
(v01[j + ix0] <= 0) && (v01[j + ix0 + 1] <= 0) &&
(v01[j + W81 + ix0] <= 0) && (v01[j + W81 + ix0 + 1] <= 0))
) draw2(ix0 * 8 - WT / 2, iy0 * 8 - HT / 2, iz0 * 8 - DT / 2);
}
}
}
}
}

int main()
{
unsigned int ui;
double a, x, y, z;

printf("start\n");
// writing hedder
b1[0x00] = 'B'; b1[0x01] = 'M';
b1[0x02] = FZ & 0xff; b1[0x03] = FZ >> 8 & 0xff; b1[0x04] = FZ >> 16 & 0xff; b1[0x05] = FZ >> 24;
b1[0x06] = 0; b1[0x07] = 0; b1[0x08] = 0; b1[0x09] = 0;
b1[0x0a] = 0x36; b1[0x0b] = 0; b1[0x0c] = 0; b1[0x0d] = 0;
b1[0x0e] = 0x28; b1[0x0f] = 0; b1[0x10] = 0; b1[0x11] = 0;
b1[0x12] = WT & 0xff; b1[0x13] = WT >> 8 & 0xff; b1[0x14] = WT >> 16 & 0xff; b1[0x15] = WT >> 24;
b1[0x16] = HT & 0xff; b1[0x17] = HT >> 8 & 0xff; b1[0x18] = HT >> 16 & 0xff; b1[0x19] = HT >> 24;
b1[0x1a] = 0x1; b1[0x1b] = 0; b1[0x1c] = 0x18; b1[0x1d] = 0;
b1[0x1e] = 0; b1[0x1f] = 0; b1[0x20] = 0; b1[0x21] = 0;
b1[0x22] = BZ & 0xff; b1[0x23] = BZ >> 8 & 0xff; b1[0x24] = BZ >> 16 & 0xff; b1[0x25] = BZ >> 24;
b1[0x26] = 0; b1[0x27] = 0; b1[0x28] = 0; b1[0x29] = 0;
b1[0x2a] = 0; b1[0x2b] = 0; b1[0x2c] = 0; b1[0x2d] = 0;
b1[0x2e] = 0; b1[0x2f] = 0; b1[0x30] = 0; b1[0x31] = 0;
b1[0x32] = 0; b1[0x33] = 0; b1[0x34] = 0; b1[0x35] = 0;

// initialize
ux[0] = 1.0; ux[1] = 0.0; ux[2] = 0.0;
uy[0] = 0.0; uy[1] = 1.0; uy[2] = 0.0;
uz[0] = 0.0; uz[1] = 0.0; uz[2] = 1.0;

// initial transposition *** modify here
// px += 0.2;
// py += 0.1;

// rotation *** modify here
a = -135.0 * pi2 / 90.0;
x = ux[0]; y = uy[0]; ux[0] = x * cos(a) - y * sin(a); uy[0] = x * sin(a) + y * cos(a);
x = ux[1]; y = uy[1]; ux[1] = x * cos(a) - y * sin(a); uy[1] = x * sin(a) + y * cos(a);
x = ux[2]; y = uy[2]; ux[2] = x * cos(a) - y * sin(a); uy[2] = x * sin(a) + y * cos(a);
a = -60.0 * pi2 / 90.0;
y = uy[0]; z = uz[0]; uy[0] = y * cos(a) - z * sin(a); uz[0] = y * sin(a) + z * cos(a);
y = uy[1]; z = uz[1]; uy[1] = y * cos(a) - z * sin(a); uz[1] = y * sin(a) + z * cos(a);
y = uy[2]; z = uz[2]; uy[2] = y * cos(a) - z * sin(a); uz[2] = y * sin(a) + z * cos(a);
a = 15.0 * pi2 / 90.0;
x = ux[0]; z = uz[0]; ux[0] = x * cos(a) + z * sin(a); uz[0] = x * -sin(a) + z * cos(a);
x = ux[1]; z = uz[1]; ux[1] = x * cos(a) + z * sin(a); uz[1] = x * -sin(a) + z * cos(a);
x = ux[2]; z = uz[2]; ux[2] = x * cos(a) + z * sin(a); uz[2] = x * -sin(a) + z * cos(a);

// making bitmap data b0
draw1();

// output file
err = fopen_s(&ftp, "t24.bmp", "wb");
if (err != NULL) return err;
for (ui = 0; ui < 0x36; ui++) {
fputc(b1[ui], ftp);
}
for (ui = 0; ui < BZ; ui++) {
fputc(b0[ui], ftp);
}
err = fclose(ftp);
printf("complete\n");
return err;
}