#!/usr/bin/env sh
# -*- wisp -*-
guile -L $(dirname $(dirname $(realpath "$0"))) -c '(import (language wisp spec))'
exec guile -L $(dirname $(dirname $(realpath "$0"))) --language=wisp -e '(@@ (examples d20world) main)' -s "$0" "$@"
; !#
; A world projected on a d20 (20-sided die, ikosaeder)
; For this we need a vector with 20 elements, a vector which shows the
; neighboring elements and accessor functions which give us the
; relevant elements for any set of longitude and latitude as well as
; its inverse (element-id to lon+lat). For further subdivisions, just
; elevate the center of each edge and connect these centers.
; Advection: Give each field a wind direction: target fields with an
; advection fraction: The fraction of the value which will be
; transported into the other field. Basic system: Follow the numbers.
define-module : examples d20world
. #:export : world neighbors d20-as-text d20-diffuse
use-modules : ice-9 format
use-modules : srfi srfi-1
use-modules
: ice-9 popen
. #:select : open-output-pipe close-pipe
define world : make-vector 20 0
define neighbors : make-vector 20
; count from the top
; Contains the numbers instead of the indexes, to make it easier for
; me to think about them.
;
; 7 8
; 3 4
; 1
; 6 2 9
; 5 10
;
; 14 13
; 18 17
; 20
; 15 19 12
; 16 11
;
define neighbors-helper
' : 1 2 3 4
2 1 5 10
3 1 6 7
4 1 8 9
5 2 6 14
6 3 5 15
7 3 8 16
8 4 7 11
9 4 10 12
10 1 9 13
20 19 18 17
19 20 16 11
18 20 15 14
17 20 13 12
16 19 17 7
15 18 16 6
14 18 13 5
13 17 14 10
12 17 11 9
11 19 12 8
let loop : : relationships neighbors-helper
cond
: null? relationships
. neighbors
else
let*
: cur : car relationships
idx : 1- : car cur
vec : cdr cur
vector-set! world idx : 1+ idx
vector-set! neighbors idx : make-vector 3
let setidx : : idxtoset '(0 1 2)
cond
: null? idxtoset
; the outer loop continues here
loop : cdr relationships
else
vector-set!
vector-ref neighbors idx
car idxtoset
1- : list-ref vec : car idxtoset
setidx : cdr idxtoset
define advection-directions
make-vector 20
let loop : : index 20
cond
: = 0 index
. advection-directions
: = 20 index
vector-set! advection-directions (1- index) (1- index)
loop : 1- index
else
vector-set! advection-directions (1- index) index
loop : 1- index
define : d20-value-ascii-color-string letter value
. "Create an ascii color string for d20."
let
: csi "["
color : inexact->exact : max 17 : min 230 : floor : * 12 value
format #f "~A38;5;~dm~A~Am" csi color letter csi
define : d20-value-ascii-color-string-show-values letter value
. "Create an ascii color string for d20."
let
: csi "["
color : inexact->exact : max 17 : min 230 : floor : * 12 value
int : inexact->exact : floor : * 12 value
format #f "~A38;5;~dm~A~Am" csi color int csi
define : d20-as-text-base world-vector function
. "show the given d20 world as text"
let
: template "
~A ~A
~A ~A
~A
~A ~A ~A
~A ~A
~A ~A
~A ~A
~A
~A ~A ~A
~A ~A
"
indexes ' : 7 8 3 4 1 6 2 9 5 10 14 13 18 17 20 15 19 12 16 11
apply format : append (list #f template) : map function indexes : map (lambda (x) (vector-ref world (1- x))) indexes
define : d20-as-text world-vector
. "show the given d20 world as text"
d20-as-text-base world-vector d20-value-ascii-color-string-show-values
define : d20-cursor-up-text world-vector
. "Kill each line of the text of the world vector in a terminal."
let*
: text : d20-as-text-base world-vector d20-value-ascii-color-string-show-values
lines : string-split text #\newline
format #t "[~AA" : 1- : length lines
define : d20-diffuse world neighbors D
. "Diffuse the values on the d20 using the diffusion constant D. Step 1: Simply iterative."
let leapfrog : : targets '(0 1 2)
if : null? targets
. world
let loop : : neighbors-to-diffuse : iota : vector-length neighbors
cond
: null? neighbors-to-diffuse
leapfrog : cdr targets
else
let*
: originidx : car neighbors-to-diffuse ; index in world and in neighbors
targetleap : car targets
targetidx : vector-ref (vector-ref neighbors originidx) targetleap
originval : vector-ref world originidx
targetval : vector-ref world targetidx
diff : * (/ D 3) : - targetval originval
vector-set! world originidx : + originval diff
vector-set! world targetidx : - targetval diff
loop : cdr neighbors-to-diffuse
define : d20-advect world advection-directions A
. "Advect the values on the d20 using the advection constant A."
let loop : : neighbors-to-advect : iota : vector-length advection-directions
cond
: null? neighbors-to-advect
. world
else
let*
: source : car neighbors-to-advect
target : vector-ref advection-directions source
source-value : vector-ref world source
target-value : vector-ref world target
change : * A source-value
source-new : - source-value change
target-new : + target-value change
; format #t "target: ~A, source: ~A, change: ~A\n" target source change
when : not : = source target
vector-set! world source source-new
vector-set! world target target-new
loop : cdr neighbors-to-advect
define d20numbers '(1 14 10 6
19 18 4 8 9 16
2 3 17 13 12 5
11 15 7 20)
define : cellidx->dienumber idx
list-ref d20numbers idx
define : dienumber->cellidx number
list-index (λ(x)(= x number)) d20numbers
define : latlonsixthslabidx latfromtop lonfrac
. "calculate the index in a sixth longitude slab of the icosaeder"
; TODO: use shortest surface distance from center of triangle as faster metric.
let*
: triangleheight : / (sqrt 3) 2
length-top-to-bottom-at-lon0 : + 1 (* 2 triangleheight)
height-deg : * 180 : / triangleheight length-top-to-bottom-at-lon0
side-deg : * 180 : / 1 length-top-to-bottom-at-lon0
; in one sixth of the icosaeder, there are 6 reachable
; fields. I am indexing them from top to bottom.
; format #t "latfromtop: ~a, lonfrac: ~a, height-deg/3: ~a, side-deg: ~a\n" latfromtop lonfrac (/ height-deg 3) side-deg
cond
: < latfromtop : / height-deg 3
. 0
: < latfromtop : - (* 2 (/ height-deg 3)) (* lonfrac (/ height-deg 3))
. 0
: < latfromtop : * 2 : / height-deg 3
. 1
: < latfromtop : + (* 2 (/ height-deg 3)) (* lonfrac (* 2 (/ height-deg 3)))
. 1
: < latfromtop : * 4 : / height-deg 3
. 2
: < latfromtop : - (+ side-deg (* 2 (/ height-deg 3))) (* lonfrac (- (+ side-deg (* 2 (/ height-deg 3))) (* 4 (/ height-deg 3))))
. 2
: < latfromtop : + side-deg : * 2 : / height-deg 3
. 3
: < latfromtop : + (+ side-deg (* 2 (/ height-deg 3))) (* lonfrac (- (+ side-deg (* 4 (/ height-deg 3))) (+ side-deg (* 2 (/ height-deg 3)))))
. 3
: < latfromtop : - (+ side-deg (* 5 (/ height-deg 3))) (* lonfrac (- (+ side-deg (* 5 (/ height-deg 3))) (+ side-deg (* 4 (/ height-deg 3)))))
. 4
else
. 5
define : latlon2cellidx lat lon
. "Convert a position given as latitude (-90 .. 90) and
longitude (0 .. 360) into the correct cell index.
This uses heavy linear approximation."
; FIXME: there is still a bug, as shown in the map plot.
; cell 1 (index 0) is on top, cell 20 at the bottom. The left
; border of cell 2 is situated at longitude 0. We can cleanly
; divide the upper part of the icosaeder into 3 regions by
; longitude. Let's do that.
let* ; the sector number is defined by the uppermost triangle
: sector : if (< lon 120) 2 (if (< lon 240) 4 3)
; we start by calculating the fraction inside the sector
lonsectorfraction : modulo lon 120
; we can further subdivide the sector by longitude into two subsectors
subsector : if (< lonsectorfraction 60) 0 1
subseclon : if (= subsector 0) lonsectorfraction (- 120 lonsectorfraction)
lonfrac : / subseclon 60
latfromtop : - 90 lat
sixthslab : latlonsixthslabidx latfromtop lonfrac
; for each sector and subsector, set the dienumber
slabsec->index '((2 . ((1 14 19 13 15 20) (1 14 16 17 11 20)))
(4 . ((1 6 9 3 11 20) (1 6 8 2 7 20)))
(3 . ((1 10 4 5 7 20) (1 10 18 12 15 20))))
dienumber->cellidx
list-ref
list-ref
assoc-ref slabsec->index sector
. subsector
. sixthslab
define : main args
. "Test the code"
if : > 2 (length args)
set! args : append args '("88") ; lat
if : > 3 (length args)
set! args : append args '("45") ; lon
display : latlon2cellidx (string->number (first (take-right args 2))) (string->number (last args))
newline
display : d20-as-text world
newline
; format #t "Diffuse ~A\n" 0.01
; d20-diffuse world neighbors 0.01
; display : d20-as-text world
; newline
; format #t "Advect ~A\n" 0.1
; d20-advect world advection-directions 0.1
; display : d20-as-text world
; newline
; format #t "Diffuse ~A\n" 0.1
; d20-diffuse world neighbors 0.1
; display : d20-as-text world
; newline
format #t "Diffuse: ~A*(~A)\n" 100 0.1
let loop : : steps 100
cond
: = 0 steps
. world
else
d20-diffuse world neighbors 0.1
loop : 1- steps
display : d20-as-text world
newline
let
: number 20
val 1
format #t "disturb: ~A to ~A\n" number val
vector-set! world (1- number) val
display : d20-as-text world
newline
; format #t "Diffuse ~A\n" 0.1
; d20-diffuse world neighbors 0.1
; display : d20-as-text world
; newline
;
; format #t "Advect: ~A*(~A)\n" 1000 0.001
; let loop : : steps 1000
; cond
; : = 0 steps
; . world
; else
; d20-advect world advection-directions 0.001
; display : d20-as-text world
; d20-cursor-up-text world
; loop : 1- steps
; display : d20-as-text world
; newline
; format #t "Diffuse: ~A*(~A)\n" 1000 0.004
; let loop : : steps 1000
; cond
; : = 0 steps
; . world
; else
; d20-diffuse world neighbors 0.004
; display : d20-as-text world
; d20-cursor-up-text world
; loop : 1- steps
display : d20-as-text world
newline
format #t "Diffuse+Advect: ~A*(~A+~A)\n" 500 0.002 0.003
let loop : : steps 500
cond
: = 0 steps
. world
else
d20-diffuse world neighbors 0.002
d20-advect world advection-directions 0.003
display : d20-as-text world
d20-cursor-up-text world
loop : 1- steps
display : d20-as-text world
newline
let
:
v
let loop
: lon 359
lat 89
map '()
zone '()
cond
: and (= lat -90) (= lon 0)
cons : cons (vector-ref world (latlon2cellidx lat lon)) zone
. map
: = lon 0
loop
. 359
- lat 1
cons : cons (vector-ref world (latlon2cellidx lat lon)) zone
. map
. '()
else
loop
- lon 1
. lat
. map
cons : vector-ref world (latlon2cellidx lat lon)
. zone
port : open-output-pipe "python"
display "a = \"" port
write v port
display "\"" port
newline port
display "a = eval(a.replace('(', '[').replace(')', ']').replace(' ',', '))" port
newline port
display "import numpy as np
import pylab as pl
import mpl_toolkits.basemap as bm
arr = np.array(a)
m = bm.Basemap(projection='cea', resolution='l', lat_ts=37.5)
m.drawcoastlines(color='k',linewidth=0.3)
m.drawmeridians(np.arange(-120.0, 180.0, 60), labels=[0,0,0,1], linewidth=0.15) # , yoffset=6) # labels = [left,right,top,bottom]
m.drawparallels(np.arange(-60.0, 90.0, 30), labels=[1,0,0,0], linewidth=0.15)
ny, nx = arr.shape
lons, lats = pl.meshgrid(range(-nx/2, nx/2 + nx%2),
range(-ny/2, ny/2 + ny%2))
x, y = m(lons, lats)
m.pcolormesh(x, y, arr, cmap=pl.get_cmap('Paired'))
pl.colorbar()
pl.show()
" port
close-pipe port
newline
; ; now plot the result
; let : : port : open-output-pipe "python"
; format port "from mpl_toolkits.mplot3d import Axes3D, art3d
; import numpy as np
; import scipy as sp
; from matplotlib import cm
; import matplotlib.pyplot as plt
; from scipy.spatial import Delaunay
;
; def Icosahedron():
; h = 0.5*(1+np.sqrt(5))
; p1 = np.array([[0,1,h],[0,1,-h],[0,-1,h],[0,-1,-h]])
; p2 = p1[:,[1,2,0]]
; p3 = p1[:,[2,0,1]]
; return np.vstack((p1,p2,p3))
;
; Ico = Icosahedron()
; tri = Delaunay(Ico)
; CH = tri.convex_hull
; points = tri.points
;
; fig = plt.figure(figsize=(4.0,4.0))
; ax = fig.add_subplot(111, projection='3d')
;
; print points
; for i in range(points.shape[0]):
; neighbors = tri.neighbors[i,:]
; for n in range(points.shape[0]):
; pts = []
; for u in range(points.shape[0]):
; pt = np.zeros((3,3))
; pt[0,:] = points[(i),:]
; pt[1,:] = points[(n),:]
; pt[2,:] = points[(u),:]
; # print pt
; pt *= 0.5
; pt += 0.5
; pts.append(pt)
; tr = art3d.Poly3DCollection(pts)
; tr.set_color([(0.9*i)/points.shape[0]] + [(0.9*n)/points.shape[0]]*3)
; ax.add_collection3d(tr)
; # ax.plot_surface(x, y, z, color='g')
;
; plt.show()
;
; exit()\n"
; close-pipe port