Vesistön suunnittelu

Kun tälle ohjelmalle kertoo saarten muodon ja sijainnin, se näyttää lopputuloksen ja kirjoittaa tuloksen tiedostoon simulointiohjelman ymmärtämään muotoon

Ohjelmointiin kuuluu allaolevan kaltaisia mystisiä määrittelyjä. On melko helppo selvittää, mitä määrittelyjä ohjelma tarvitsee.

# kaaviokuvan piirtelyyn apuvälineitä
import matplotlib.patches as patch
from matplotlib.collections import PatchCollection
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
import re
from textwrap import fill, wrap
import params as par

Simulointia varten maailma jaetaan kuutioihin. Koska tietokoneeni teho ei riittänyt kolmiulotteisen maailman simulointiin, minulla on kuutioita vain yksi kerros kaksiulotteisessa ruudukossa.

Matematiikassa on tapana esittää pisteen paikka tasossa muodossa (x, y), missä x kertoo vaakasuoran etäisyyden origosta ja y pystysuoran etäisyyden.

Ruudukko, jota simuloinnissa käytän, on järkevä esittää taulukkona, joiden käsittelyyn on käteviä funktiota numpy-moduulissa. Harmi kyllä, taulukossa ruudun paikka ilmoitetaan muodossa taulukko[i, j], missä i ilmoittaa rivin ja j sarakkeen. Yleensä olen tottunut, että i viittaa x-akseliin ja j y-akseliin.

Geometrian kielellä ylläoleva tarkoittanee, että (x,y) tasosta [i,j] taulukkoon siirtyminen on peilaus suoran x=y suhteen. En tiedä, mutta virheitä tuli, kun en heti aluksi miettinyt loppuun saakka, mitä kaikkea tuosta seuraa.

DIM_i ja DIM_j kertovat, montako riviä ja saraketta ruudukossa on. Täydellisessä ohjelmassa saarten sijainti ja mittasuhteet ilmoitetaan niin, että maailma skaalautuu automaattisesti ruutujen määrää muutettaessa. Tämä ei ole ihan täydellinen, muttei huonokaan.

DIM_i = 240
DIM_j = 200

DI_2 = int(DIM_i/2.0)
DJ_2 = int(DIM_j/2.0)
DI_4 = int(DIM_i/4.0)
DJ_4 = int(DIM_j/4.0)
DI_5 = int(DIM_i/5.0)
DJ_5 = int(DIM_j/5.0)
DI_6 = int(DIM_i/6.0)
DJ_6 = int(DIM_j/6.0)

Funktioita saarten piirtämiseen

Piirretään saari

Alla aliohjelma, joka piirtää näytölle saaren sille määriteltyyn paikkaan.

'..add_patch(patch.Ellipse( ...))' näyttää sotkuiselta ja toisinaan onkin työlästä selvittää, millä komennolla saa näytölle haluamansa kuvion.

def piirra_saari(kaavio, saari, ds):
    color = (0.5, 1.0, 0.2)
    (name, (i0, j0), w, h, kulma) = saari
    alpha = np.degrees(kulma)
    kaavio.add_patch(patch.Ellipse(
        (j0*ds, i0*ds), 2.0*w*ds, 2.0*h*ds, 
        angle=-alpha, facecolor=color, edgecolor='black'))    

Piirretään kaikki saaret

Piirretään saaret käyttäen ylläolevaa funktiota. Mystisen näköistä. Onneksi webistä voi copy-pasteta esimerkkejä.

Olennainen tapahtuu riveillä:

    for saari in saaret:
        piirra_saari(kaavio, saari, ds)

def kaikki_saaret(saaret):
    fig_w = 8.0
    fig_h = 8.0
    ds = 1.0

#    ds = fig_w/150.0

#    fig, kaavio = plt.subplots(1, 1)
    fig, kaavio = plt.subplots(1, 1, figsize=(fig_w, fig_h))

    kaavio.add_patch(patch.Rectangle(
        (0, 0), DIM_j*ds, DIM_i*ds, angle=0.0, edgecolor='black', 
         facecolor='none', linewidth=2))

    kaavio.set_aspect('equal', adjustable='datalim')
    plt.axis('on')
    kaavio.grid(which='major')
    kaavio.grid(which='minor', linestyle='--')
    
# Major ticks every 10, minor ticks every 5
    major_ticks = np.arange(0, max(DIM_i, DIM_j)+1, 20)
    minor_ticks = np.arange(0, max(DIM_i, DIM_j)+1, 5)

    kaavio.set_xticks(major_ticks)
    kaavio.set_xticks(minor_ticks, minor=True)
    kaavio.set_yticks(major_ticks)
    kaavio.set_yticks(minor_ticks, minor=True)
    
    for saari in saaret:
        piirra_saari(kaavio, saari, ds)

    kaavio.plot()   #Causes an autoscale update.
    plt.show()
    plt.savefig('allas+saarekkeet.png', format='png')
 

Saarten sijainti ja muoto

Saari = (saaren nimi, (keskipisteen x-koordinaatti, keskipisteen y-koordinaatti), ellipsin leveys, ellipsin korkeus, kallistuskulma)

E_I = DI_5
E_J = DJ_4
EW = 30.0
EH = 20.0

Y_I = DI_6 - 15
Y_J = DJ_5 - 15

e_kulma = np.arctan(float(DIM_i/DIM_j))

saaret = [
    # keskiympyrä
    ('kk', (DI_2, DJ_2), 20.0, 20.0, 0.0),
    # pikkuellipsit
    ('vk', (DI_2, Y_J), 17.0, 25.0, 0.0),
    ('ok', (DI_2, DIM_j - Y_J), 17.0, 25.0, 0.0),
    ('ak', (Y_I, DJ_2), 18.0, 15.0, np.pi/2.0),
    ('yk', (DIM_i - Y_I, DJ_2), 18.0, 15.0, np.pi/2.0),
    # isot ellipsit
    ('va', (E_I, E_J), EW, EH, -e_kulma),
    ('oa', (E_I, DIM_j-E_J), EW, EH, e_kulma),
    ('vy',(DIM_i - E_I, E_J), EW, EH, e_kulma),
    ('oy', (DIM_i - E_I, DIM_j - E_J), EW, EH, -e_kulma)
    ]

kaikki_saaret(saaret)

<Figure size 432x288 with 0 Axes>

Kerrotaan saarten paikat simulointiohjelmalle

# Tämmöisiä kummallisen näköisiä määrittelyjäkin tarvitaan
x, x_0, x_a, x_b, y, y_0, y_a, y_b, w, h, alpha, i, j, i_0, j_0, ds = \
    sp.symbols('x x_0 x_a x_b y y_0 y_a y_b w h alpha i j  i_0 j_0 ds')

Ellipsin yhtälö

Saaret ovat ellipsejä. Ruudukon ruutu[i, j] on saarella, mikäli i ja j ovat jonkin ellipsin sisällä.

Origossa vaakasuorassa oleva ellipsi

Piste (x, y) on ellipsin kehällä, mikäli allaoleva yhtälö toteutuu. Jos korvaamme =-merkin <-merkillä, saamme epäyhtälön, joka kertoo, onko piste ellipsin sisällä, siis saarella.

saarella = sp.expand((x/w)**2 + (y/h)**2 - 1)
display(sp.Eq(saarella, 0))

$\displaystyle -1 + \frac{x^{2}}{w^{2}} + \frac{y^{2}}{h^{2}} = 0$

Pyöräytetään ellipsiä

Yhtälö ellipsille, jota on pyöräytetty kulman alpha verran.

Luulin pyöräyttymisen kaavan helpoksi johtaa, mutta lopulta jouduin kopiomaan wikipediasta.

Tuntuu, että koordinaatistot (x,y) ja (x_a, y_a) tuli nimettyä harhaanjohtavasti nurinperin. Koska ei haittaa lopputulosta, en jää jahkailemaan.

#beta = alpha-sp.pi/2
rotate = [
    (x, x_a*sp.cos(alpha)-y_a*sp.sin(alpha)), 
    (y, x_a*sp.sin(alpha)+y_a*sp.cos(alpha))]
rotated_vector = (x*i+y*j).subs(rotate)
display(rotated_vector)

$\displaystyle i \left(x_{a} \cos{\left(\alpha \right)} - y_{a} \sin{\left(\alpha \right)}\right) + j \left(x_{a} \sin{\left(\alpha \right)} + y_{a} \cos{\left(\alpha \right)}\right)$

saarella_rotated = saarella.subs(rotate)
display(sp.Eq(saarella_rotated, 0))

$\displaystyle -1 + \frac{\left(x_{a} \cos{\left(\alpha \right)} - y_{a} \sin{\left(\alpha \right)}\right)^{2}}{w^{2}} + \frac{\left(x_{a} \sin{\left(\alpha \right)} + y_{a} \cos{\left(\alpha \right)}\right)^{2}}{h^{2}} = 0$

Siirrytään pois origosta

trans = [(x_a, x_b-x_0), (y_a, y_b-y_0)]
display(trans)

[(x_a, -x_0 + x_b), (y_a, -y_0 + y_b)]

saarella_rot_tr = saarella_rotated.subs(trans)
display(sp.Eq(saarella_rot_tr, 0))

$\displaystyle -1 + \frac{\left(\left(- x_{0} + x_{b}\right) \cos{\left(\alpha \right)} - \left(- y_{0} + y_{b}\right) \sin{\left(\alpha \right)}\right)^{2}}{w^{2}} + \frac{\left(\left(- x_{0} + x_{b}\right) \sin{\left(\alpha \right)} + \left(- y_{0} + y_{b}\right) \cos{\left(\alpha \right)}\right)^{2}}{h^{2}} = 0$

Kirjoitetaan yleinen ellipsin yhtälö python funktioksi

Tää on vähän sotkuista hommaa. Turha lukea, ellei aio tehdä samaa.

Python koodin tuottamiseen on valmiitakin komentoja, mutta ne eivät tee ihan sitä, mitä tarvitsen.

Valmis koodi alempana.

ij = [(x_b, j*ds), (y_b, i*ds), (x_0, j_0*ds), (y_0, i_0*ds)]
saarella3 = str(saarella_rot_tr.subs(ij).evalf())

iv, jv = sp.symbols('iv jv')
vec_a = x*iv + y*jv
vec_b = vec_a.subs(rotate).subs(trans).subs(ij)
display(vec_b)
print(vec_b)

$\displaystyle iv \left(- \left(ds i - ds i_{0}\right) \sin{\left(\alpha \right)} + \left(ds j - ds j_{0}\right) \cos{\left(\alpha \right)}\right) + jv \left(\left(ds i - ds i_{0}\right) \cos{\left(\alpha \right)} + \left(ds j - ds j_{0}\right) \sin{\left(\alpha \right)}\right)$

iv*(-(ds*i - ds*i_0)*sin(alpha) + (ds*j - ds*j_0)*cos(alpha)) + jv*((ds*i - ds*i_0)*cos(alpha) + (ds*j - ds*j_0)*sin(alpha))

code_subs = [('ds', 'par.ds'), ('sin', 'math.sin'), ('cos', 'math.cos')]
for sub in code_subs:
    saarella3 = re.sub(sub[0], sub[1], saarella3)

display(saarella3)

'-1.0 + (-(par.ds*i - par.ds*i_0)*math.sin(alpha) + (par.ds*j - par.ds*j_0)*math.cos(alpha))**2/w**2 + ((par.ds*i - par.ds*i_0)*math.cos(alpha) + (par.ds*j - par.ds*j_0)*math.sin(alpha))**2/h**2'

with open('saaret.py', 'w') as fil:
    fil.write("# -*- coding: utf-8 -*-\n\n")
    fil.write("import math\n")
    fil.write("import params as par\n")

    fil.write('DIM_i = '+str(DIM_i)+'\n')
    fil.write('DIM_j = '+str(DIM_j)+'\n')
    fil.write('DI_2 = '+str(DI_2)+'\n')
    fil.write('DJ_2 = '+str(DJ_2)+'\n')
    fil.write('DI_4 = '+str(DI_4)+'\n')
    fil.write('DJ_4 = '+str(DJ_4)+'\n')
    fil.write('DI_5 = '+str(DI_5)+'\n')
    fil.write('DJ_5 = '+str(DJ_5)+'\n')
#    fil.write('off_i = '+str(off_i)+'\n')
#    fil.write('off_j = '+str(off_j)+'\n')
    fil.write('E_I = '+str(E_I)+'\n')
    fil.write('E_J = '+str(E_J)+'\n')
    fil.write('Y_I = '+str(Y_I)+'\n')
    fil.write('Y_J = '+str(Y_J)+'\n')
 
    fil.write('saaret = \\\n')
    fil.write(fill(str(saaret), initial_indent='      ', subsequent_indent='       ', break_long_words=False))

    fil.write('\n\n\ndef saarella(i, j, i_0, j_0, w, h, alpha):\n')
    # fil.write('    alpha = math.radians(alpha)\n')
    fil.write('    return \\\n')
    fil.write(fill(str(saarella3), initial_indent='      (', subsequent_indent='       ', break_long_words=False))
    fil.write(')\n')

Valmis koodi

with open('saaret.py', 'r') as code:
    print(code.read())

# -*- coding: utf-8 -*-

import math
import params as par
DIM_i = 240
DIM_j = 200
DI_2 = 120
DJ_2 = 100
DI_4 = 60
DJ_4 = 50
DI_5 = 48
DJ_5 = 40
E_I = 48
E_J = 50
Y_I = 25
Y_J = 25
saaret = \
      [('kk', (120, 100), 20.0, 20.0, 0.0), ('vk', (120, 25), 17.0,
       25.0, 0.0), ('ok', (120, 175), 17.0, 25.0, 0.0), ('ak', (25,
       100), 18.0, 15.0, 1.5707963267948966), ('yk', (215, 100), 18.0,
       15.0, 1.5707963267948966), ('va', (48, 50), 30.0, 20.0,
       -0.8760580505981934), ('oa', (48, 150), 30.0, 20.0,
       0.8760580505981934), ('vy', (192, 50), 30.0, 20.0,
       0.8760580505981934), ('oy', (192, 150), 30.0, 20.0,
       -0.8760580505981934)]


def saarella(i, j, i_0, j_0, w, h, alpha):
    return \
      (-1.0 + (-(par.ds*i - par.ds*i_0)*math.sin(alpha) + (par.ds*j -
       par.ds*j_0)*math.cos(alpha))**2/w**2 + ((par.ds*i -
       par.ds*i_0)*math.cos(alpha) + (par.ds*j -
       par.ds*j_0)*math.sin(alpha))**2/h**2)