From 429f0890d6b02a7aed8d42c763048e5de47441f0 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 15 Jun 2022 12:20:45 +0200 Subject: [PATCH] generate beautiful pictures of the complex Anosov set --- complex_anosov.c | 242 +++++++++++++++++++++++++++++---------------- compute_picture.sh | 9 ++ enumerate.c | 134 ++++++++++++++++++------- enumerate.h | 4 +- generators.c | 62 ++++++------ make_picture.py | 67 +++++++++++++ 6 files changed, 362 insertions(+), 156 deletions(-) create mode 100755 compute_picture.sh create mode 100755 make_picture.py diff --git a/complex_anosov.c b/complex_anosov.c index cc1a8a4..1c61928 100644 --- a/complex_anosov.c +++ b/complex_anosov.c @@ -8,35 +8,39 @@ #include #include #include +#include #define LOOP(i,n) for(int i = 0; i < (n); i++) #define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0); +//#define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__) +#define INFO(msg, ...) /* - Elements up to length 0: 1 - Elements up to length 1: 4 - Elements up to length 2: 10 - Elements up to length 3: 22 - Elements up to length 4: 46 - Elements up to length 5: 91 - Elements up to length 6: 175 - Elements up to length 7: 334 - Elements up to length 8: 634 - Elements up to length 9: 1198 - Elements up to length 10: 2260 - Elements up to length 11: 4261 - Elements up to length 12: 8029 - Elements up to length 13: 15124 - Elements up to length 14: 28486 - Elements up to length 15: 53650 - Elements up to length 16: 101038 - Elements up to length 17: 190279 - Elements up to length 18: 358339 - Elements up to length 19: 674830 - Elements up to length 20: 1270846 - Elements up to length 21: 2393266 - Elements up to length 22: 4507012 - Elements up to length 23: 8487625 + Number of elements + up to length 0: 1 + up to length 1: 4 + up to length 2: 10 + up to length 3: 22 + up to length 4: 46 + up to length 5: 91 + up to length 6: 175 + up to length 7: 334 + up to length 8: 634 + up to length 9: 1198 + up to length 10: 2260 + up to length 11: 4261 + up to length 12: 8029 + up to length 13: 15124 + up to length 14: 28486 + up to length 15: 53650 + up to length 16: 101038 + up to length 17: 190279 + up to length 18: 358339 + up to length 19: 674830 + up to length 20: 1270846 + up to length 21: 2393266 + up to length 22: 4507012 + up to length 23: 8487625 */ static double gaussian_sqrt5_real(NUMBER x) @@ -82,9 +86,56 @@ static double gaussian_sqrt5_imag(NUMBER x) return result; } +static int read_idlist(const char *filename, int *list) +{ + FILE *f = fopen(filename, "r"); + if(f == NULL) { + fprintf(stderr, "Could not open %s\n", filename); + exit(1); + } + + char *line = NULL; + size_t len = 0; + ssize_t read; + int n = 0; + + while ((read = getline(&line, &len, f)) != -1) + list[n++] = atoi(line); + + if (line) + free(line); + + fclose(f); + + return n; +} + +static int compare_int(const void *a, const void *b) +{ + const int *aa = a; + const int *bb = b; + return (*aa > *bb) - (*aa < *bb); +} + +static double runtime() +{ + static struct timespec starttime; + static int started = 0; + + if(!started) { + clock_gettime(CLOCK_MONOTONIC, &starttime); + started = 1; + } + + struct timespec curtime; + double diff; + clock_gettime(CLOCK_MONOTONIC, &curtime); + return (curtime.tv_sec - starttime.tv_sec) + (curtime.tv_nsec - starttime.tv_nsec) / 1e9; +} + enum mode { MODE_HELP, - MODE_TRACES, + MODE_EIGENVALUES, MODE_SUMMARY, MODE_TRACE_IDS }; @@ -94,10 +145,13 @@ int main(int argc, char *argv[]) char buf[100]; int mode; + runtime(); // start timer + + // parse arguments if(argc < 2 || strcmp(argv[1], "help") == 0) mode = MODE_HELP; - else if(strcmp(argv[1], "traces") == 0) - mode = MODE_TRACES; + else if(strcmp(argv[1], "evs") == 0) + mode = MODE_EIGENVALUES; else if(strcmp(argv[1], "summary") == 0) mode = MODE_SUMMARY; else if(strcmp(argv[1], "trace_ids") == 0) @@ -106,72 +160,78 @@ int main(int argc, char *argv[]) mode = MODE_HELP; if(mode == MODE_HELP) { - fprintf(stderr, "Usage: %s [arguments]\n", argv[0]); + fprintf(stderr, "Usage: %s [arguments]\n", argv[0]); fprintf(stderr, "%s help display this page\n", argv[0]); - fprintf(stderr, "%s traces enumerate group and output unique trace/trace inverse pairs\n", argv[0]); - fprintf(stderr, "%s summary only output max slope etc.\n", argv[0]); - fprintf(stderr, "%s trace_ids list of ids of unique traces\n", argv[0]); + fprintf(stderr, "%s evs enumerate group and output unique (log) eigenvalue triples\n", argv[0]); + fprintf(stderr, "%s summary only output max slope etc.\n", argv[0]); + fprintf(stderr, "%s trace_ids list of ids of unique traces\n", argv[0]); return 0; } - if(argc < 5) { + if(argc < 8) { fprintf(stderr, "Not enough arguments!\n"); return 0; } int n = atoi(argv[2]); + int nlist, nuniq; - mpq_t qreal, qimag; - mpq_inits(qreal, qimag, NULL); - mpq_set_str(qreal, argv[3], 10); - mpq_set_str(qimag, argv[4], 10); -// mpq_set_si(qreal, 50, 10); -// mpq_set_si(qimag, 1, 10); - mpq_canonicalize(qreal); - mpq_canonicalize(qimag); - - /* - int length = 0; - LOOP(i, n) { - if(group->elements[i].length > length) { - printf("Elements up to length %d: %d\n", length, i); - length = group->elements[i].length; - } - } - return 0; - */ - - /* - LOOP(i, n) { - groupelement_t *cur = &group->elements[i]; - groupelement_t *other; - - cur->conjugacy_class = cur; // start with itself and reduce if possible - - LOOP(j, 3) { - if(cur->left[j] && cur->left[j]->right[j]) { - other = cur->left[j]->right[j]; - if(other->id < cur->id) - cur->conjugacy_class = other->conjugacy_class; - } - if(cur->right[j] && cur->right[j]->left[j]) { - other = cur->right[j]->left[j]; - if(other->id < cur->id) - cur->conjugacy_class = other->conjugacy_class; - } - } - } - */ + int q[3]; + mpq_t treal, timag; + q[0] = atoi(argv[3]); + q[1] = atoi(argv[4]); + q[2] = atoi(argv[5]); + mpq_inits(treal, timag, NULL); + mpq_set_str(treal, argv[6], 10); + mpq_set_str(timag, argv[7], 10); + mpq_canonicalize(treal); + mpq_canonicalize(timag); + // enumerate group + INFO("generate group\n"); group_t *group = coxeter_init_triangle(5, 5, 5, n); + // read list of elements we need to compute, or just fill it with all elements + INFO("prepare list\n"); + int *idlist = malloc(n*sizeof(int)); + char *id_file = getenv("IDLIST"); + if(id_file != NULL) { + nlist = read_idlist(id_file, idlist); + // sort and symmetrize the list + qsort(idlist, nlist, sizeof(int), compare_int); + int ninverses = 0; + LOOP(i, nlist) { + int id = idlist[i]; + int invid = group->elements[id].inverse->id; + if(!bsearch(&invid, idlist, nlist, sizeof(int), compare_int)) { + idlist[nlist+ninverses] = invid; + ninverses++; + } + } + nlist += ninverses; + qsort(idlist, nlist, sizeof(int), compare_int); + } else { + // just list all elements which have inverses + nlist = 0; + LOOP(i, n) { + if(group->elements[i].inverse) + idlist[nlist++] = i; + } + } + + // get generator matrices + INFO("make generators\n"); mat gen[3]; LOOP(i, 3) mat_init(gen[i], 3, QT_GAUSS_SQRT5); - generators_triangle_reflection_group_555_complex(gen, 2, 2, 2, qreal, qimag); + generators_triangle_reflection_group_555_complex(gen, q[0], q[1], q[2], treal, timag); + // compute the traces of all elements in idlist + INFO("enumerate traces\n"); struct tracedata *traces; - int nuniq = enumerate_coxeter_group_traces(group, gen, &traces); + nuniq = enumerate_coxeter_group_traces(group, gen, &traces, idlist, nlist, 1); + // compute eigenvalues out of traces + INFO("compute eigenvalues\n"); mps_context *solver = mps_context_new(); mps_monomial_poly *poly = mps_monomial_poly_new(solver, 3); mps_context_set_output_prec(solver, 20); // relative precision @@ -198,7 +258,14 @@ int main(int argc, char *argv[]) if(fabs(ev_abs2[0]) < fabs(ev_abs2[1])) SWAP(double, ev_abs2[0], ev_abs2[1]); - if(log(ev_abs2[0]) < 1e-3) // we regard this as a finite order element + if(mode == MODE_TRACE_IDS) { + // we only want to record unordered pairs here + if(CMP(traces[i].tr, traces[i].trinv) >= 0) + printf("%d\n", traces[i].id); + continue; + } + + if(log(ev_abs2[0]) < 1e-6) // we regard this as a finite order element continue; double slope = - log(ev_abs2[0]) / log(ev_abs2[2]); @@ -207,28 +274,31 @@ int main(int argc, char *argv[]) max_slope_id = traces[i].id; } - if(mode == MODE_TRACES) { + if(mode == MODE_EIGENVALUES) { printf("%d %f %f %f\n", - traces[i].id, log(ev_abs2[0]), log(ev_abs2[1]), log(ev_abs2[2])); - } else if(mode == MODE_TRACE_IDS) { - printf("%d\n", traces[i].id); + traces[i].id, log(ev_abs2[0])/2, log(ev_abs2[1])/2, log(ev_abs2[2])/2); } } + // output summary coxeter_snprint(buf, sizeof(buf), &group->elements[max_slope_id]); - if(mode == MODE_SUMMARY) { - gmp_fprintf(stdout, "%f %f %d %d %f %s\n", mpq_get_d(qreal), mpq_get_d(qimag), n, nuniq, max_slope, buf); - } else { - gmp_fprintf(stderr, "q = %Qd + i*%Qd\tElements: %d\tTraces: %d\tMaximal slope: %f at %s\n", qreal, qimag, n, nuniq, max_slope, buf); + gmp_fprintf(stdout, "%f %f %d %d %d %f %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf); + } else if(mode == MODE_EIGENVALUES) { + gmp_fprintf(stderr, "q = %.2f + i*%.2f\tn = %d\tnlist = %d\tnuniq = %d\tMaximal slope: %f at %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf); + } else if(mode == MODE_TRACE_IDS) { + gmp_fprintf(stderr, "q = %.2f + i*%.2f\tElements: %d\tTraces: %d\n", mpq_get_d(treal), mpq_get_d(timag), n, nuniq); } - mps_monomial_poly_free(solver, MPS_POLYNOMIAL(poly)); + // clean up + INFO("clean up\n"); +// mps_monomial_poly_free(solver, MPS_POLYNOMIAL(poly)); mps_context_free(solver); enumerate_tracedata_clear(traces, nuniq); + free(idlist); LOOP(i, 3) mat_clear(gen[i]); coxeter_clear(group); - mpq_clears(qreal, qimag, NULL); + mpq_clears(treal, timag, NULL); return 0; } diff --git a/compute_picture.sh b/compute_picture.sh new file mode 100755 index 0000000..226905a --- /dev/null +++ b/compute_picture.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +for i in $(seq -50 50); do + for j in $(seq 0 50); do + if [ $i -ne 0 ] || [ $j -ne 0 ]; then + IDLIST=output/idlist_13 ./complex_anosov summary 15124 1 1 1 $i/50 $j/50 + fi + done +done diff --git a/enumerate.c b/enumerate.c index bd847e9..dec5be7 100644 --- a/enumerate.c +++ b/enumerate.c @@ -3,25 +3,84 @@ #include #define LOOP(i,n) for(int i = 0; i < (n); i++) +// #define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__) +#define INFO(msg, ...) -void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices) +static int compare_int(const void *a, const void *b) { - TYPE t = GETTYPE(gen[0]->x[0]); + const int *aa = a; + const int *bb = b; + return (*aa > *bb) - (*aa < *bb); +} + +static int index_in_list(const int *list, int n, int key) +{ + int *ptr = bsearch(&key, list, n, sizeof(int), compare_int); + if(ptr) + return ptr - list; + else + return -1; +} + +// idlist is assumed to be sorted and symmetric +void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist) +{ + TYPE type = GETTYPE(gen[0]->x[0]); mat_workspace *ws; - ws = mat_workspace_init(3, t); + // mark elements which we need to compute + int ncompute = 0; + LOOP(i, group->size) group->elements[i].need_to_compute = 0; + LOOP(i, nlist) { + groupelement_t *cur = &group->elements[idlist[i]]; + while(cur && cur->need_to_compute == 0) { + cur->need_to_compute = 1; + ncompute++; + cur = cur->parent; + } + } + INFO("need to compute %d elements for %d elements in list\n", ncompute, nlist); - mat_identity(matrices[0]); - for(int i = 1; i < group->size; i++) { - if(!group->elements[i].inverse) - continue; + // compute them depth first + int maxlen = group->elements[group->size-1].length; - int parent = group->elements[i].parent->id; - int letter = group->elements[i].letter; + groupelement_t **stack = malloc((maxlen+1)*sizeof(groupelement_t)); + int *letter_stack = malloc((maxlen+1)*sizeof(int)); + mat *matrix_stack = malloc((maxlen+1)*sizeof(mat)); + int level = 0; - mat_multiply(ws, matrices[i], matrices[parent], gen[letter]); + LOOP(i, maxlen+1) mat_init(matrix_stack[i], 3, type); + ws = mat_workspace_init(3, type); + + stack[0] = &group->elements[0]; + mat_identity(matrix_stack[0]); + letter_stack[0] = 0; + + while(level >= 0) { + int letter = letter_stack[level]; + groupelement_t *next = stack[level]->left[letter]; + if(next && next->length > level && next->need_to_compute) { + mat_multiply(ws, matrix_stack[level+1], gen[letter], matrix_stack[level]); + int id = next->id; + int listid = index_in_list(idlist, nlist, id); + if(listid != -1) + mat_copy(matrices[listid], matrix_stack[level+1]); + next->need_to_compute = 0; + stack[level+1] = next; + letter_stack[level+1] = 0; + level++; + } else { + letter = ++letter_stack[level]; + while(letter >= group->rank) { // done with this level, go back down + level--; + if(level < 0) + break; + letter = ++letter_stack[level]; + } + } } + LOOP(i, maxlen+1) mat_clear(matrix_stack[i]); mat_workspace_clear(ws); } @@ -48,44 +107,45 @@ static int compare_tracedata_id(const void *a, const void *b) return ida > idb ? 1 : ida < idb ? -1 : 0; } -int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out) +int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique) { TYPE t = GETTYPE(gen[0]->x[0]); - int n = group->size; + int nuniq; - mat *matrices = malloc(n*sizeof(mat)); - struct tracedata *traces = malloc(n*sizeof(struct tracedata)); - struct tracedata **distinct_traces = malloc(n*sizeof(struct tracedata*)); + mat *matrices = malloc(nlist*sizeof(mat)); + struct tracedata *traces = malloc(nlist*sizeof(struct tracedata)); + struct tracedata **distinct_traces = malloc(nlist*sizeof(struct tracedata*)); - LOOP(i, n) mat_init(matrices[i], 3, t); - enumerate_coxeter_group(group, gen, matrices); + LOOP(i, nlist) mat_init(matrices[i], 3, t); + enumerate_coxeter_group(group, gen, matrices, idlist, nlist); - LOOP(i, n) { + LOOP(i, nlist) { INIT(traces[i].tr, t); INIT(traces[i].trinv, t); - } - LOOP(i, n) { - if(!group->elements[i].inverse) - continue; + int inv_id_in_list = index_in_list(idlist, nlist, group->elements[idlist[i]].inverse->id); mat_trace(traces[i].tr, matrices[i]); - mat_trace(traces[i].trinv, matrices[group->elements[i].inverse->id]); - traces[i].id = i; + mat_trace(traces[i].trinv, matrices[inv_id_in_list]); + traces[i].id = idlist[i]; distinct_traces[i] = &traces[i]; } - qsort(distinct_traces, n, sizeof(struct tracedata*), compare_tracedata); + if(unique) { + qsort(distinct_traces, nlist, sizeof(struct tracedata*), compare_tracedata); - int nuniq = 0; - LOOP(i, n) { - if(i == 0 || compare_tracedata(&distinct_traces[i], &distinct_traces[nuniq-1]) != 0) { - distinct_traces[nuniq] = distinct_traces[i]; - nuniq++; - } else { - int oldlength = group->elements[distinct_traces[nuniq-1]->id].length; - int newlength = group->elements[distinct_traces[i]->id].length; - if(newlength < oldlength) - distinct_traces[nuniq-1]->id = distinct_traces[i]->id; + nuniq = 0; + LOOP(i, nlist) { + if(i == 0 || compare_tracedata(&distinct_traces[i], &distinct_traces[nuniq-1]) != 0) { + distinct_traces[nuniq] = distinct_traces[i]; + nuniq++; + } else { + int oldlength = group->elements[distinct_traces[nuniq-1]->id].length; + int newlength = group->elements[distinct_traces[i]->id].length; + if(newlength < oldlength) + distinct_traces[nuniq-1]->id = distinct_traces[i]->id; + } } + } else { + nuniq = nlist; } qsort(distinct_traces, nuniq, sizeof(struct tracedata*), compare_tracedata_id); @@ -99,8 +159,8 @@ int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata ** unique_traces[i].id = distinct_traces[i]->id; } - LOOP(i, n) mat_clear(matrices[i]); - LOOP(i, n) { + LOOP(i, nlist) mat_clear(matrices[i]); + LOOP(i, nlist) { CLEAR(traces[i].tr); CLEAR(traces[i].trinv); } diff --git a/enumerate.h b/enumerate.h index 0e3b171..ea208ac 100644 --- a/enumerate.h +++ b/enumerate.h @@ -12,8 +12,8 @@ struct tracedata { NUMBER trinv; }; -void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices); -int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out); +void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist); +int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique); void enumerate_tracedata_clear(struct tracedata *traces, int n); diff --git a/generators.c b/generators.c index da958af..fe773b7 100644 --- a/generators.c +++ b/generators.c @@ -3,39 +3,39 @@ #define LOOP(i,n) for(int i = 0; i < (n); i++) -void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER q) +void generators_triangle_reflection_generic(mat *gen, NUMBER sqrho1, NUMBER sqrho2, NUMBER sqrho3, NUMBER t) { NUMBER tmp; - NUMBER qinv; - TYPE t = GETTYPE(rho1); + NUMBER tinv; + TYPE type = GETTYPE(sqrho1); - INIT(tmp, t); - INIT(qinv, t); - INV(qinv, q); + INIT(tmp, type); + INIT(tinv, type); + INV(tinv, t); LOOP(i, 3) { - mat_init(gen[i], 3, t); + mat_init(gen[i], 3, type); mat_zero(gen[i]); LOOP(j, 3) { SET_INT(*mat_ref(gen[i], j, j), i == j ? 1 : -1); } } - MUL(tmp, q, rho1); + MUL(tmp, t, sqrho1); NEG(*mat_ref(gen[2], 1, 2), tmp); - MUL(tmp, q, rho2); - NEG(*mat_ref(gen[0], 2, 0), tmp); - MUL(tmp, q, rho3); - NEG(*mat_ref(gen[1], 0, 1), tmp); + NEG(*mat_ref(gen[0], 2, 0), sqrho2); + NEG(*mat_ref(gen[1], 0, 1), sqrho3); - NEG(*mat_ref(gen[1], 2, 1), qinv); - NEG(*mat_ref(gen[2], 0, 2), qinv); - NEG(*mat_ref(gen[0], 1, 0), qinv); + MUL(tmp, tinv, sqrho1); + NEG(*mat_ref(gen[1], 2, 1), tmp); + NEG(*mat_ref(gen[2], 0, 2), sqrho2); + NEG(*mat_ref(gen[0], 1, 0), sqrho3); CLEAR(tmp); - CLEAR(qinv); + CLEAR(tinv); } +// warning: this is not compatible anymore! int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t t) { int p[3] = {p1, p2, p3}; @@ -107,23 +107,23 @@ int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, i int q[3] = {q1, q2, q3}; - NUMBER rho[3]; - LOOP(i, 3) INIT(rho[i], QT_GAUSS_SQRT5); + NUMBER sqrho[3]; + LOOP(i, 3) INIT(sqrho[i], QT_GAUSS_SQRT5); - // compute rho1, rho2, rho3 + // compute sqrho1, sqrho2, sqrho3 LOOP(i, 3) { if(q[i] == 1) { - // 4cos(pi/5)^2 = 3/2 + 1/2*sqrt(5) - mpq_set_si(rho[i]->a[0], 3, 2); - mpq_set_si(rho[i]->a[1], 7, 12); - mpq_set_si(rho[i]->a[2], 0, 1); - mpq_set_si(rho[i]->a[3], -1, 24); + // 2cos(pi/5) = 1/2 + 1/2*sqrt(5) = 1/2 + 7/12 a - 1/24 a^3 + mpq_set_si(sqrho[i]->a[0], 1, 2); + mpq_set_si(sqrho[i]->a[1], 7, 12); + mpq_set_si(sqrho[i]->a[2], 0, 1); + mpq_set_si(sqrho[i]->a[3], -1, 24); } else if(q[i] == 2) { - // 4cos(pi/5)^2 = 3/2 - 1/2*sqrt(5) - mpq_set_si(rho[i]->a[0], 3, 2); - mpq_set_si(rho[i]->a[1], -7, 12); - mpq_set_si(rho[i]->a[2], 0, 1); - mpq_set_si(rho[i]->a[3], 1, 24); + // 2cos(pi/5) = -1/2 + 1/2*sqrt(5) + mpq_set_si(sqrho[i]->a[0], -1, 2); + mpq_set_si(sqrho[i]->a[1], 7, 12); + mpq_set_si(sqrho[i]->a[2], 0, 1); + mpq_set_si(sqrho[i]->a[3], -1, 24); } else { return 0; } @@ -138,9 +138,9 @@ int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, i mpq_set_si(param->a[3], 1, 12); mpq_mul(param->a[3], param->a[3], timag); - generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], param); + generators_triangle_reflection_generic(gen, sqrho[0], sqrho[1], sqrho[2], param); - LOOP(i, 3) CLEAR(rho[i]); + LOOP(i, 3) CLEAR(sqrho[i]); CLEAR(param); return 1; diff --git a/make_picture.py b/make_picture.py new file mode 100755 index 0000000..5fc877f --- /dev/null +++ b/make_picture.py @@ -0,0 +1,67 @@ +#!/usr/bin/python + +import sys +import math + +def sharpen(x): + alpha = 1 + y = abs(2*x-1)**alpha # between 0 and 1 + if x > 0.499 and x < 0.501: + return 0.5 + elif x > 0.5: + return (y+1)/2 + else: + return (1-y)/2 + +f = open(sys.argv[1]) +lines = [x.split() for x in f] +f.close() + +res1 = int(sys.argv[2]) # 400 pixel per unit +res2 = int(sys.argv[3]) # 50 pixel in picture (acutally one quadrant) + +data = {(round(float(l[0])*res1), round(float(l[1])*res1)) : float(l[5]) for l in lines} +data[(0,0)] = 2.0 + +print("P3") +print("1000 1000") +print("255") + +for i in range (-500,500): + for j in range(-500,500): + x = j/500*res2 + y = i/500*res2 + + if y < 0: + y = -y + + x0 = math.floor(x) + y0 = math.floor(y) + x1 = math.ceil(x) + y1 = math.ceil(y) + tx = sharpen(x-x0) + ty = sharpen(y-y0) + + if not (x0,y0) in data or not (x0,y1) in data or not (x1,y0) in data or not (x1,y1) in data: + value = 0 + else: + value = 0.0 + value += data[(x0,y0)]*(1-tx)*(1-ty) + value += data[(x0,y1)]*(1-tx)*ty + value += data[(x1,y0)]*tx*(1-ty) + value += data[(x1,y1)]*tx*ty + value -= 1 + + value = 0 if value < 0 else 1 if value > 1 else value + + r = round(255*math.sqrt(value)) + g = round(255*(value)**3) + b = round(255*math.sin(2*math.pi*(value))) + + r = 0 if r < 0 else 255 if r > 255 else r + g = 0 if g < 0 else 255 if g > 255 else g + b = 0 if b < 0 else 255 if b > 255 else b + + print("%d %d %d" % (r,g,b)) + +#print(data)