＜空白ページ＞
ヤコビの楕円関数の複素平面表示のビットマップファイルを計算する、C言語のプログラムです。解説は次項でやります。(8月11日改訂)
// ConsoleApplication9.cpp
// ellptic4.cpp; Jacobi elliptic functions on a Gaussian plane
// based on bmp24.cpp
// 200828; making 24bit color bmp file
// 230711; 230722; 230729; 230811
#include <stdio.h>
#include <math.h>
const int WT = 1024; // width (mod 8 = 0) *** modify here
const int HT = 1024; // height (mod 8 = 0) *** modify here
const int BZ = WT * HT * 3; // buffer size
const int FZ = BZ + 0x36; // file size
unsigned char b1[0x36]; // header*2
unsigned char b2[BZ]; // bitmap sn
unsigned char b3[BZ]; // bitmap cn
unsigned char b4[BZ]; // bitmap dn
const double pi4 = atan(1.0);
// elliptic k (0.0 <= k <= 1.0)
double elk = 0.5; // k^2
// draw area definition
// UD,VD will be automatically modified with K'/K
double UC = 0.0001; // real base
double VC = 0.0001; // imaginary base
double UD = 0.0026; // real pitch // hard
double VD = 0.0026; // imaginary pitch // hard
// contour information
double DD = 15.0; // value height pitch
double DE = 0.125; // value edge width (0.01 - 0.99)
// pole area definition
double anp = 4.0; // pole area (infinity) in degree
double asp = 4.0; // zero point area in degree
// Riemann sphere magnitude (absolute value at equator)
// it dose not adjusted by 1/k or k'/k or k'. please rewrite manually
// const double alf = sqrt(1.0 - elk * elk) / elk; // example
double alf = 1.0; //
// color index
const unsigned char cl0[7][3] = { // *** modify here; {b, g, r} 0-255
{ 0x2e, 0x2e, 0x2e }, // background
{ 0x77, 0x77, 0xcc }, // pole (infinity)
{ 0xbb, 0xbb, 0x22 }, // zero point
{ 0x11, 0x11, 0x11 }, // real part = 0; pure imaginary
{ 0xff, 0xff, 0xff }, // imaginary part = 0; pure real
{ 0x66, 0x66, 0xaa }, // contour edge northern hemisphere
{ 0x99, 0x99, 0x66 } // contour edge southern hemisphere
};
const unsigned char cl1[24][3] = { // *** modify here; {b, g, r} 0-255
{ 0x99, 0xff, 0x33 },{ 0xff, 0xff, 0x33 },{ 0xff, 0xff, 0x33 },{ 0xff, 0x99, 0x33 },
{ 0xff, 0x99, 0x33 },{ 0xff, 0x33, 0x33 },{ 0xff, 0x33, 0x33 },{ 0xff, 0x33, 0x99 },
{ 0xff, 0x33, 0x99 },{ 0xff, 0x33, 0xff },{ 0xff, 0x33, 0xff },{ 0x99, 0x33, 0xff },
{ 0x99, 0x33, 0xff },{ 0x33, 0x33, 0xff },{ 0x33, 0x33, 0xff },{ 0x33, 0x99, 0xff },
{ 0x33, 0x99, 0xff },{ 0x33, 0xff, 0xff },{ 0x33, 0xff, 0xff },{ 0x33, 0xff, 0x99 },
{ 0x33, 0xff, 0x99 },{ 0x33, 0xff, 0x33 },{ 0x33, 0xff, 0x33 },{ 0x99, 0xff, 0x33 }
};
// elliptic function constants
double elq; // q
double ell, elqc; // l, q'
double elkc, elkck; // k', k'/k
double ellk, ellkc; // K, K'
double te2a, te3a, te0a, tf2a, tf3a, tf0a; // theta(0)s
// elliptic function lines
double te1[WT + 1], te2[WT + 1], te3[WT + 1], te0[WT + 1]; // using q
double tf1[HT + 1], tf2[HT + 1], tf3[HT + 1], tf0[HT + 1]; // using q'
// elliptic function buffers
double te11[WT + 1], te12[WT + 1], te13[WT + 1], te10[WT + 1];
double te23[WT + 1], te20[WT + 1], te30[WT + 1], te00[WT + 1];
// line value buffer
double sn1[WT + 1], sn2[WT + 1], cn1[WT + 1], cn2[WT + 1], dn1[WT + 1], dn2[WT + 1];
// theta functions series
void the10(double *t1, double *t2, double *t3, double *t0, int i, double q, double v) {
t1[i] = 2 * (pow(q, 0.25) * sin(pi4 * 4.0 * v) - pow(q, 2.25) * sin(pi4 * 4.0 * 3.0 * v)
+ pow(q, 6.25) * sin(pi4 * 4.0 * 5.0 * v) - pow(q, 12.25) * sin(pi4 * 4.0 * 7.0 * v)
+ pow(q, 20.25) * sin(pi4 * 4.0 * 9.0 * v) - pow(q, 30.25) * sin(pi4 * 4.0 * 11.0 * v)
+ pow(q, 42.25) * sin(pi4 * 4.0 * 13.0 * v) - pow(q, 56.25) * sin(pi4 * 4.0 * 15.0 * v)
+ pow(q, 72.25) * sin(pi4 * 4.0 * 17.0 * v) - pow(q, 90.25) * sin(pi4 * 4.0 * 19.0 * v));
t2[i] = 2 * (pow(q, 0.25) * cos(pi4 * 4.0 * v) + pow(q, 2.25) * cos(pi4 * 4.0 * 3.0 * v)
+ pow(q, 6.25) * cos(pi4 * 4.0 * 5.0 * v) + pow(q, 12.25) * cos(pi4 * 4.0 * 7.0 * v)
+ pow(q, 20.25) * cos(pi4 * 4.0 * 9.0 * v) + pow(q, 30.25) * cos(pi4 * 4.0 * 11.0 * v)
+ pow(q, 42.25) * cos(pi4 * 4.0 * 13.0 * v) + pow(q, 56.25) * cos(pi4 * 4.0 * 15.0 * v)
+ pow(q, 72.25) * cos(pi4 * 4.0 * 17.0 * v) + pow(q, 90.25) * cos(pi4 * 4.0 * 19.0 * v));
t3[i] = 1.0 + 2.0 * (q * cos(pi4 * 8.0 * v) + pow(q, 4.0) * cos(pi4 * 8.0 * 2.0 * v)
+ pow(q, 9.0) * cos(pi4 * 8.0 * 3.0 * v) + pow(q, 16.0) * cos(pi4 * 8.0 * 4.0 * v)
+ pow(q, 25.0) * cos(pi4 * 8.0 * 5.0 * v) + pow(q, 36.0) * cos(pi4 * 8.0 * 6.0 * v)
+ pow(q, 49.0) * cos(pi4 * 8.0 * 7.0 * v) + pow(q, 64.0) * cos(pi4 * 8.0 * 8.0 * v)
+ pow(q, 81.0) * cos(pi4 * 8.0 * 9.0 * v) + pow(q, 100.0) * cos(pi4 * 8.0 * 10.0 * v));
t0[i] = 1.0 + 2.0 * (-q * cos(pi4 * 8.0 * v) + pow(q, 4.0) * cos(pi4 * 8.0 * 2.0 * v)
- pow(q, 9.0) * cos(pi4 * 8.0 * 3.0 * v) + pow(q, 16.0) * cos(pi4 * 8.0 * 4.0 * v)
- pow(q, 25.0) * cos(pi4 * 8.0 * 5.0 * v) + pow(q, 36.0) * cos(pi4 * 8.0 * 6.0 * v)
- pow(q, 49.0) * cos(pi4 * 8.0 * 7.0 * v) + pow(q, 64.0) * cos(pi4 * 8.0 * 8.0 * v)
- pow(q, 81.0) * cos(pi4 * 8.0 * 9.0 * v) + pow(q, 100.0) * cos(pi4 * 8.0 * 10.0 * v));
}
void init1(void) {
int i;
double u;
// prepare constants
te2a = 2.0 * (pow(elq, 0.25) + pow(elq, 2.25) + pow(elq, 6.25) + pow(elq, 12.25) + pow(elq, 20.25)
+ pow(elq, 30.25) + pow(elq, 42.25) + pow(elq, 56.25) + pow(elq, 72.25) + pow(elq, 90.25));
te3a = 1.0 + 2.0 * (elq + pow(elq, 4.0) + pow(elq, 9.0) + pow(elq, 16.0) + pow(elq, 25.0)
+ pow(elq, 36.0) + pow(elq, 49.0) + pow(elq, 64.0) + pow(elq, 81.0) + pow(elq, 100.0));
te0a = 1.0 + 2.0 * (-elq + pow(elq, 4.0) - pow(elq, 9.0) + pow(elq, 16.0) - pow(elq, 25.0)
+ pow(elq, 36.0) - pow(elq, 49.0) + pow(elq, 64.0) - pow(elq, 81.0) + pow(elq, 100.0));
elqc = exp(pi4 * pi4 * 16.0 / log(elq));
tf2a = 2.0 * (pow(elqc, 0.25) + pow(elqc, 2.25) + pow(elqc, 6.25) + pow(elqc, 12.25) + pow(elqc, 20.25)
+ pow(elqc, 30.25) + pow(elqc, 42.25) + pow(elqc, 56.25) + pow(elqc, 72.25) + pow(elqc, 90.25));
tf3a = 1.0 + 2.0 * (elqc + pow(elqc, 4.0) + pow(elqc, 9.0) + pow(elqc, 16.0) + pow(elqc, 25.0)
+ pow(elqc, 36.0) + pow(elqc, 49.0) + pow(elqc, 64.0) + pow(elqc, 81.0) + pow(elqc, 100.0));
tf0a = 1.0 + 2.0 * (-elqc + pow(elqc, 4.0) - pow(elqc, 9.0) + pow(elqc, 16.0) - pow(elqc, 25.0)
+ pow(elqc, 36.0) - pow(elqc, 49.0) + pow(elqc, 64.0) - pow(elqc, 81.0) + pow(elqc, 100.0));
elk = pow(te2a / te3a, 2.0); elkc = pow(te0a / te3a, 2.0); elkck = pow(te0a / te2a, 2.0);
ellk = 2.0 * pi4 * pow(te3a, 2.0); // K
ellkc = ellk / (4.0 * pi4) * log(1 / elq); // K'
// automatic expansion
VD *= ellk / pi4 / 2.0; UD *= ellkc / pi4 / 2.0;
// prepare lines
for (i = 0; i <= WT; i++) {
u = UC + UD * (double)(i - WT / 2);
the10(te1, te2, te3, te0, i, elq, u);
}
for (i = 0; i <= HT; i++) {
u = VC + VD * (double)(i - HT / 2);
the10(tf1, tf2, tf3, tf0, i, elqc, u);
}
// prepare horizontal values
for (i = 0; i <= WT; i++) {
te11[i] = te1[i] * te1[i]; te12[i] = te1[i] * te2[i];
te13[i] = te1[i] * te3[i]; te10[i] = te1[i] * te0[i];
te23[i] = te2[i] * te3[i]; te20[i] = te2[i] * te0[i];
te30[i] = te3[i] * te0[i]; te00[i] = te0[i] * te0[i];
}
}
void drawxn1(unsigned char** pp, int i, double* xn1, double* xn2, double xn1a, double xn2a, double xn0) {
double z;
int k;
z = sqrt(xn1[i] * xn1[i] + xn2[i] * xn2[i]) / xn0;
if (z < asp) { // south pole area
*(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
return;
}
if (z > anp) { // north pole area
*(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
return;
}
if (xn1[i] * xn1[i + 1] <= 0) { // pure imaginary
*(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
return;
}
if (xn1[i] * xn1a <= 0) { // pure imaginary
*(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
return;
}
if (xn2[i] * xn2[i + 1] <= 0) { // pure real
*(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
return;
}
if (xn2[i] * xn2a <= 0) { // pure real
*(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
return;
}
z = atan(z / alf) - pi4; // -pi/4 <= z < pi/4
if (fmod(fabs(z), DD) <= DE * DD) { // in a contour
if (z > 0.0) { // north edge
*(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
return;
}
else { // south edge
*(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
return;
}
}
k = (int)((atan2(xn2[i], xn1[i]) / pi4) * 2.999 + 12.0); // hard
*(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];
}
void draw1(void) {
int i, j;
unsigned char* p2, *p3, *p4;
double xn0;
double tf11, tf22, tf12, tf13, tf10, tf23, tf20, tf30;
double sn1a, sn2a, cn1a, cn2a, dn1a, dn2a;
// the 1st line
tf12 = elkck * tf1[0] * tf2[0]; tf13 = -elkck * tf1[0] * tf3[0]; tf10 = -elkc * tf1[0] * tf0[0];
tf23 = elkc * tf2[0] * tf3[0]; tf20 = elkck * tf2[0] * tf0[0]; tf30 = elkck * tf3[0] * tf0[0];
tf11 = tf1[0] * tf1[0]; tf22 = tf2[0] * tf2[0];
for (i = 0; i < WT; i++) {
sn1[i] = te10[i] * tf30; sn2[i] = te23[i] * tf12;
cn1[i] = te20[i] * tf20; cn2[i] = te13[i] * tf13;
dn1[i] = te30[i] * tf23; dn2[i] = te12[i] * tf10;
}
p2 = b2; p3 = b3; p4 = b4;
for (j = 1; j <= HT; j++) {
// sentinel
sn1[WT] = te10[WT] * tf30; sn2[WT] = te23[WT] * tf12;
cn1[WT] = te20[WT] * tf20; cn2[WT] = te13[WT] * tf13;
dn1[WT] = te30[WT] * tf23; dn2[WT] = te12[WT] * tf10;
// calculate tne next upper line value
tf12 = elkck * tf1[j] * tf2[j]; tf13 = -elkck * tf1[j] * tf3[j]; tf10 = -elkc * tf1[j] * tf0[j];
tf23 = elkc * tf2[j] * tf3[j]; tf20 = elkck * tf2[j] * tf0[j]; tf30 = elkck * tf3[j] * tf0[j];
for (i = 0; i < WT; i++) {
xn0 = te11[i] * tf11 + te00[i] * tf22; // this point
sn1a = te10[i] * tf30; sn2a = te23[i] * tf12;
cn1a = te20[i] * tf20; cn2a = te13[i] * tf13;
dn1a = te30[i] * tf23; dn2a = te12[i] * tf10;
if (xn0 < 1e-10) { // pole
*p2++ = cl0[1][0]; *p2++ = cl0[1][1]; *p2++ = cl0[1][2];
*p3++ = cl0[1][0]; *p3++ = cl0[1][1]; *p3++ = cl0[1][2];
*p4++ = cl0[1][0]; *p4++ = cl0[1][1]; *p4++ = cl0[1][2];
}
else { // calculable
drawxn1(&p2, i, sn1, sn2, sn1a, sn2a, xn0);
drawxn1(&p3, i, cn1, cn2, cn1a, cn2a, xn0);
drawxn1(&p4, i, dn1, dn2, dn1a, dn2a, xn0);
}
sn1[i] = sn1a; sn2[i] = sn2a;
cn1[i] = cn1a; cn2[i] = cn2a;
dn1[i] = dn1a; dn2[i] = dn2a;
}
tf11 = tf1[j] * tf1[j]; tf22 = tf2[j] * tf2[j];
}
}
int main()
{
FILE* ftp;
errno_t err;
unsigned int ui;
double z;
// reading setting file
err = fopen_s(&ftp, "jac.txt", "r");
if (err != NULL) {
err = fopen_s(&ftp, "jac.txt", "w");
if (err != NULL) return err;
fprintf(ftp, "0.5; 0; 0.0001; 0.0001; 0.0026; 0.0026; 15; 0.125; 4; 4; 1.0\n");
fprintf(ftp, "\n");
fprintf(ftp, "// 0.5; 0; 0.0001; 0.0001; 0.0026; 0.0026; 15; 0.125; 4; 4; 1.0\n");
fprintf(ftp, "// k^2; q; UC; VC; UD; VD; DD; DE; anp; asp; alf\n");
fprintf(ftp, "\n");
fprintf(ftp, "// if q is void, k^2 is starting value (0 < k, q < 1)\n");
fprintf(ftp, "// UC; VC: drawing area center (0.00125)\n");
fprintf(ftp, "// UD; VD: next pixel value pitch (0.0025)\n");
fprintf(ftp, "// DD: hemisphere parallels pitch in degree (15.0)\n");
fprintf(ftp, "// DE: contour edge width (0.01 - 0.99)\n");
fprintf(ftp, "// anp; asp: pole area definition northern, southern in degree (2.0)\n");
fprintf(ftp, "// alf: Riemann sphere magnitude (absolute value at equator) (1.0)\n");
fprintf(ftp, "// sn: 1/k, cn: k'/k, dn: k' for example\n");
err = fclose(ftp);
if (err != NULL) return err;
err = fopen_s(&ftp, "jac.txt", "r");
if (err != NULL) return err;
}
fscanf_s(ftp, "%le;%le;%le;%le;%le;%le;%le;%le;%le;%le;%le\n", &elk, &elq, &UC, &VC, &UD, &VD, &DD, &DE, &anp, &asp, &alf);
err = fclose(ftp);
if (err != NULL) return err;
if (elk <= 1e-10) elk = 1e-10;
if (elk >= 1.0 - 1e-10) elk = 1.0 - 1e-10;
elk = sqrt(elk);
if (elq <= 0.0) {
ell = (1.0 - pow(1.0 - elk * elk, 0.25)) / 2.0 / (1.0 + pow(1.0 - elk * elk, 0.25));
elq = ell + 2.0 * pow(ell, 5.0) + 15.0 * pow(ell, 9.0) + 150.0 * pow(ell, 13.0) + 1707.0 * pow(ell, 17.0);
}
DD = pi4 * DD / 90.0; // value height pitch
anp = 1.0 / (tan(anp * pi4 / 90.0)) * alf; // pole area (infinity)
asp = tan(asp * pi4 / 90.0) * alf; // zero point area
// writing header
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
init1();
// making bitmap data
draw1();
// output file
err = fopen_s(&ftp, "sn.bmp", "wb");
if (err != NULL) return err;
for (ui = 0; ui < 0x36; ui++) {
fputc(b1[ui], ftp);
}
for (ui = 0; ui < BZ; ui++) {
fputc(b2[ui], ftp);
}
err = fclose(ftp);
if (err != NULL) return err;
err = fopen_s(&ftp, "cn.bmp", "wb");
if (err != NULL) return err;
for (ui = 0; ui < 0x36; ui++) {
fputc(b1[ui], ftp);
}
for (ui = 0; ui < BZ; ui++) {
fputc(b3[ui], ftp);
}
err = fclose(ftp);
if (err != NULL) return err;
err = fopen_s(&ftp, "dn.bmp", "wb");
if (err != NULL) return err;
for (ui = 0; ui < 0x36; ui++) {
fputc(b1[ui], ftp);
}
for (ui = 0; ui < BZ; ui++) {
fputc(b4[ui], ftp);
}
err = fclose(ftp);
if (err != NULL) return err;
err = fopen_s(&ftp, "el.txt", "w");
if (err != NULL) return err;
fprintf(ftp, "Jacobi elliptic function constants\n");
fprintf(ftp, "k = %10.5lg, k^2 = %10.5lg, k' = %10.5lg, k'^2 = %10.5lg\n", elk, elk * elk, elkc, elkc * elkc);
fprintf(ftp, "1/k = %10.5lg, k'/k = %10.5lg\n", 1 / elk, elkck);
fprintf(ftp, "K = %10.5lg, K' = %10.5lg, K/K' = %10.5lg\n", ellk, ellkc, ellk / ellkc);
fprintf(ftp, "K/(pi/2) = %10.5lg, K'/(pi/2) = %10.5lg\n", ellk / pi4 / 2.0, ellkc / pi4 / 2.0);
fprintf(ftp, "q = %10.5lg, te2(0) = %10.5lg, te3(0) = %10.5lg, te0(0) = %10.5lg\n", elq, te2a, te3a, te0a);
fprintf(ftp, "q' = %10.5lg, tf2(0) = %10.5lg, tf3(0) = %10.5lg, tf0(0) = %10.5lg\n", elqc, tf2a, tf3a, tf0a);
fprintf(ftp, "horizontal: WT = %6i, CT = %10.5lg, D = %10.5lg\n", WT, UC, UD);
fprintf(ftp, "vertical: HT = %6i, CT = %10.5lg, D = %10.5lg\n", WT, VC, VD);
fprintf(ftp, "northern hemisphere parallels: (degree) value\n");
for (z = 0; z < pi4 - 0.001; z += DD) {
fprintf(ftp, "(%5.2lf) %10.5lg\n", z / pi4 * 90.0, tan(pi4 + z) * alf);
}
fprintf(ftp, "north polar (infinity) area = %10.5lg(deg)\n", atan(1.0 / (anp / alf)) / pi4 * 90.0);
fprintf(ftp, "southern hemisphere parallels: (degree) value\n");
for (z = 0; z < pi4 - 0.001; z += DD) {
fprintf(ftp, "(%5.2lf) %10.5lg\n", z / pi4 * 90.0, tan(pi4 - z) * alf);
}
fprintf(ftp, "south polar (zero point) area = %10.5lg(deg)\n", atan(asp / alf) / pi4 * 90.0);
err = fclose(ftp);
if (err != NULL) return err;
printf("complete\n");
return err;
}
ヤコビの楕円関数の複素平面表示の、同じような絵が出てくる3回目のプログラムで恐縮です。設定用のテキストファイルでk値などが指定できるようになりました。動作は4146.のプログラムと全く同一で、係数一つ変えるにも再コンパイルの方式ではなく、外部のテキストファイルの一部を書き換えるだけで済みます。
私の世代だとこれが便利の範疇に入ります。
もちろん、最初のコンパイルは必要です。このブログでバイナリを配るのは適切では無いと考えます。Microsoft Visual Studio等をご用意下さい。コンパイルのやり方は前回同様です。
作成された実行形式プログラムには「ConsoleApplication9.exe」のような名前が付いているはずです。.exeの前の所は適当に書き換えても(たとえばCA9.exeなど)同じ動作をします。
他のパソコンで動作させるには、この26KB程度(リリース版)の実行ファイルを持って行くだけでOKです。しかしもしかしたら、C++の実行時ライブラリが無い、とのアナウンスが出るかもしれません。このライブラリは公式サイトからダウンロード出来ます。ネット検索するとやり方が出てくるはずで、導入は簡単です。
専用のフォルダを用意し、ConsoleApplication9.exeを入れ、クリックすると一瞬、コンソール画面が出てきてすぐに消えます。コンソール画面用プログラムなので、スタートメニューから「ターミナル」を起動して、文字画面で起動する方が正式な感じですが、クリックで構いません。
出てくるファイルは楕円関数の複素平面表示、sn.bmp, cn.bmp, dn.bmpとグラフを読むのに必要な情報のテキストファイル、el.txtとともに、設定ファイルのjac.txtです。
作成されるjac.txtの中身は、
0.5; 0; 0.0001; 0.0001; 0.0026; 0.0026; 15; 0.125; 4; 4; 1.0
// 0.5; 0; 0.0001; 0.0001; 0.0026; 0.0026; 15; 0.125; 4; 4; 1.0
// k^2; q; UC; VC; UD; VD; DD; DE; anp; asp; alf
// if q is void, k^2 is starting value (0 < k, q < 1)
// UC; VC: drawing area center (0.00125)
// UD; VD: next pixel value pitch (0.0025)
// DD: hemisphere parallels pitch in degree (15.0)
// DE: contour edge width (0.01 - 0.99)
// anp; asp: pole area definition northern, southern in degree (2.0)
// alf: Riemann sphere magnitude (absolute value at equator) (1.0)
// sn: 1/k, cn: k'/k, dn: k' for example
で、この1行目を書き換えると出てくる絵が変化します。元のソースファイルを見ると分かるように、1行目しか読んでいません。しかも、後ろは適当に省略できて、極端に言うと、
0.9755
の中身が6文字のファイル(jac.txt)でもOKなのです。最初の数値はk値の自乗の値の指定です。jac.txtをこれにして、実行形式を動作させるとグラフと共に出てくる情報ファイル、el.txtの中身は、
Jacobi elliptic function constants
k = 0.98767, k^2 = 0.9755, k' = 0.15652, k'^2 = 0.0245
1/k = 1.0125, k'/k = 0.15848
K = 3.2547, K' = 1.5806, K/K' = 2.0592
K/(pi/2) = 2.072, K'/(pi/2) = 1.0062
q = 0.21749, te2(0) = 1.4306, te3(0) = 1.4395, te0(0) = 0.56949
q' = 0.0015503, tf2(0) = 0.39686, tf3(0) = 1.0031, tf0(0) = 0.9969
horizontal: WT = 1024, CT = 0.0001, D = 0.0026161
vertical: HT = 1024, CT = 0.0001, D = 0.0053873
northern hemisphere parallels: (degree) value
( 0.00) 1
(15.00) 1.3032
(30.00) 1.7321
(45.00) 2.4142
(60.00) 3.7321
(75.00) 7.5958
north polar (infinity) area = 4(deg)
southern hemisphere parallels: (degree) value
( 0.00) 1
(15.00) 0.76733
(30.00) 0.57735
(45.00) 0.41421
(60.00) 0.26795
(75.00) 0.13165
south polar (zero point) area = 4(deg)
などとなります。
読み込みにはC言語に付属する標準の関数fscanf_s()を使用しているので、若干の注意が必要です。原則的には「;」で数値を区切れば良いのですが、必ず数値を書かないと、それ以降は評価されません。なので、デフォルト値をわざわざ最後まで書いてあるのです。この;で区切られた11個の数値を適当に変更するのが無難です。負数(-3)や指数表示(1e-4)などは受け付けます。
2番目の数値を0.0以外にすると、2番目の数値がq値と解釈されます。たとえば、
0.3; 0.04321
だと、k値0.3の指定は無効となり、q値0.04321から計算がスタートします。ただし、計算の精度上、q = 0.7程度よりも1に近い数値は推奨しません。
3番目と4番目のUCとVCはグラフ中央の座標です。単位は周期2Kと2K'です。原点の場合は(0K+0iK')でも良いのですが、デフォルト値のようにわずかにずらせると良好な絵になるようです。他の場合(0.2K+0.4iK'など)は普通に指定できると思います。
5番目と6番目のUDとVDは右ドットとの値の差、上ドットとの値の差、の初期値です。単位は周期2Kと2K'です。特に必要が無ければ、同じ数値を指定します。たとえば、縦横2倍の領域を見るには、0.0052と0.0052とします。
k値やq値などから周期KとK'が算出でき、その割合で縦横の比率が修正されます。その結果はel.txtに書かれます(horizontalとverticalのD値)。これはヤコビの楕円関数の入力値sn(u + iv)のuとvの変化の比率を画面上で1:1にするための措置です。初期値を別々の数値に変えると、1:1で無くなります。
7番目の値は等高線の間隔指定です。リーマン球面の緯度で指定します。赤道が0度で、15と指示すると15度ずつの間隔で等高線が描かれます。角度単位では無く具体的な数値はel.txtに出ています(等高線の外側(赤道に近い方)の縁の値)。
8番目の数値は等高線の太さの割合です。0.1だと間隔の1/10の幅、0.5だと間隔の1/2の幅になります。
9番目と10番目の数値は北極域(極の近傍)と南極域(零点の近傍)をリーマン球面の角度で指定します。値を上げると極または零点のドットが大きくなります。
11番目はリーマン球面の大きさ指定です。1.0だと単位球(ガウス平面とは単位円で交わる)なので、赤道、グラフでは赤系と青系の緯度の線が隣接しているところ、が関数値の絶対値が1の所です。この赤道の位置を1.0から移動させる数値です。例えば、関数値の絶対値の1.4142がどのあたりにあるのかを知りたい時などに、この数値を変更します。
sn, cn, dnとも、純実数と純虚数の線はKとK'の間隔で引かれているので、el.txtのK値とK'値の値から(ガウス平面の)座標が分かります。
再コンパイルが必要なのは、グラフの縦横のドット数の指定(画面の大きさ)、極点や等高線などの色、関数値の偏角を表す虹色の個々の色、です。これらをファイルから読み込ませるのは同様にすれば出来ますが、うかつに変えると訳が分からなくなりますから慎重に扱うことにしました。
昨日は普通に出勤日で普通にお仕事でした。スタッフが幾分夏バテ気味で、多少の何とかがあったくらいで、概ね順調です。
日本国内では特に大きなニュースは無かったと思います。
楕円関数の複素平面表示のプログラムはテキスト入力のついでに計算量の最適化というか結構余分な計算(ソース上)を減らしたバージョンを作りました。
しかし思いのほか、処理時間はあまり変わらないような感じです。おそらく入出力のハードやOS回りの処理時間が目立っている気がします。今の計算機には簡単な計算の部類に入るようです。
PS5とswitchに新機種の噂が立っていて、しかし発売はあるとしても、まだずっと先の話です。内容的には改良の範囲に入るとのことです。もしも驚きがあったら意外だ、の雰囲気です。
任天堂switchはソフトの関係で私は使っていません。発売後6年目でしたか、今も人気は抜群だそうです。
PS5の方は通常型の改良版に加えて、proと呼ばれる能力向上版が噂されているようです。こちらはやっと普通のゲーム機の感じになりました。初期は品不足で、私も手に入りませんでした。今は普通に購入できます。ソフトは順調に伸びている感じ。
今現在は経済の世界動向が流動的な感じなので、しばらくはこんな感じが続くのかな、と思います。VR等の新規の方向の研究はいろいろ進んでいるとは思います。
