Wednesday, April 18, 2018

Studying discrete-time dynamical systems (VII): some updates and Python code

A kind reader of Hobbymaths asked me to publish the code to visualize the images that I showed at my former post (here). So here it is, it is a bit long, but I think the comments in the code are self-explanatory. You will need Python Anaconda in order to run properly the code and generate the images. Basically you will need to save the code into a text file, for instance "dynamical_system_test.py" and then run in command line (my testing environment is Windows 10 + Python) "python dynamical_system_test.py". Currently I am preparing a paper about this topic, but it will take some months to be finished, maybe there will be something ready for the end of $2018$.

The code as it is written will generate $100$ images, which are the visualizations of the evolution of the control parameter time $t$ of the following families of dynamical systems:

$$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{D}, Re(f(x_{n} + y_n i))\cdot\cos{D})\}$$

Where $f(z)=\frac{1}{(z+1)^t}$, and initially $t \in \Bbb [1,2]$ (it can be greater than $2$ but for the purpose of the code I have set $2$ as the maximum value to be calculated). It will generate the attractors for $t=1.01, t=1.02, \cdots , t=2.00$. For each $t$, the control parameters $D$ and $n$ will run through $D \in [0,2\pi]$ in $8 \cdot 10^3$ steps, and $n \in [0,8 \cdot 10^3]$, and the images are showing the region $z=[+/-5]+[+/-5]i$, for the initial condition $(x_0=0,y_0=0)$. 

Please be aware that generating the images takes time! each one of them takes around $10$ minutes in a good computer. And also you need to re-dimension the width of the images with a program like e.g. Gimp2 (the explanation is written in the code). So it will generate (between others) the following images one by one sequentially:

$f(z)=\frac{1}{(z+1)}$ (this will be the first one it generates):


$f(z)=\frac{1}{(z+1)^{1.1}}$:


$f(z)=\frac{1}{(z+1)^{1.2}}$:


$\cdots$

$f(z)=\frac{1}{(z+1)^{2}}$ (this image is the last one that will generate):


And if you use an application to join the images you will obtain the following video (it took me some days to render it):


So here is the code, enjoy it and, if you find new interesting variations of the model please let me know (please keep the reference to this post in the comments so other people can find the original information here, thanks!):

# Studying discrete-time dynamical systems (VI): mathematical models
# for symmetry and asymmetry and some biological coincidences
# keys: recurrence-relations, dynamical-systems, visualization,
# experimental-mathematics, chaotic-systems
# Done by Hobbymaths:
# http://hobbymaths.blogspot.jp/2018/01/studying-discrete-time-dynamical.html

def nds():
    import matplotlib.pyplot as plt
    from math import pi, e
    from random import randint
    from cmath import sqrt, sin , cos, tan
    from scipy.special import digamma
   
   
    # Each system is iterated up to [ntestlimit] times
    # value of n: It can be a higher value but for the
    # tests is a good enough number to arrive to the sink points
    ntestlimit = 8000
   
    # Complete list of points of each system, the first point
    # of the list is the (x,y,z) seed
    # each dynamical system will have up to [ntestlimit] points.
    # If we arrive to a sink point before the [ntestlimit]
    # iteration there is not need to continue, so the quantity
    # of elements will be smaller.
    st=[]
   
    # Splitted lists of x,y and z coordinates of the points obtained
    # on each iteration of the system stored at the [st] list
    lox = []
    loy = []
    loz = []
   
    # D parameter range
    MinD = 0
    MaxD = 8000
   
    for expo in range (100,201):
        t = expo/100
        print("Exp " + str(t))
       
        lox = []
        loy = []
        loz = []
        plt.clf()
   
        # We will generate a dynamical system for each D parameter
        # from [MinD] to [MaxD] included.
        for f in range (MinD,MaxD):
            D=(f/MaxD)*(2*pi)
           
            st=[]
           
            #Initial condition of the dynamical system
            icx=0
            icy=0
           
            st.append([icx,icy,0])
           
            # We will generate up to [ntestlimit] points of the current
            # dynamical systes for a fixed D.
            # if we arrive to a sink point we do not need to continue the
            # loop up to [ntestlimit] because the
            # point will be repited again (so we do not need to plot more
            # than once the sink point)
            for n in range (1,ntestlimit):
                # We read the former the previous iteration point
                # (oldx,oldy,oldz)
                oldx=st[len(st)-1][0]
                oldy=st[len(st)-1][1]
                oldz=st[len(st)-1][2]
               
                # Then we calculate the f(z) function for the current
                # (x_n-1,y_n-1) values.
                c=1/((oldx+(oldy*(1j)))+1)**(t)
               
                # calculate new (x_n,y_n)
                newx=(c.imag)*sin(D)
                newy=(c.real)*cos(D)
               
                # The z coordinate is just a trick to give a depth in
                # space to the points according to D
                st.append([newx.real,newy.real,D])
               
            # Once we have generated the points of the current D
            # dynamical system, we have a set of points S_D
            # (where D is the D number). E.g S_1, S_2, etc.
           
            # This loop below will detect where the cycle begins BEGIN
            tlox=[]
            tloy=[]
            tloz=[]
           
            there_is_cycle=False
           
            for idx in range(0,len(st)-1):
                if not(abs(st[idx][0])<5 and abs(st[idx][1])<5):
                    continue
                if len(tlox)==0:
                    tlox.append(st[idx][0])
                    tloy.append(st[idx][1])
                    tloz.append(st[idx][2])
                else:
                    tmpl=st[idx+1:]
                    if st[idx] in tmpl:
                        tmpval=tmpl.index(st[idx])+1
                        there_is_cycle=True
                        break
                       
                    tlox.append(st[idx][0])
                    tloy.append(st[idx][1])
                    tloz.append(st[idx][2])                       
                       
            if there_is_cycle:
                lox.extend(tlox)
                loy.extend(tloy)
                loz.extend(tloz)
                pass
            else:
                lox.extend(tlox)
                loy.extend(tloy)
                loz.extend(tloz)
                pass
            # This loop below will detect where the cycle begins END
           
            # we set black color to the background and white color to the points
            ax = plt.gca()
            #ax.set_axis_bgcolor((0, 0, 0)) # old Python code
            ax.set_facecolor((0, 0, 0))     # new Python code (change depending of Python version)
            figure = plt.gcf()
           
            print(str(f)+"\r", end='')
           
            if f==MaxD-1:
           
                #image resolution 32000x3200 pixels
                #you will need to redimension width to 3200 pixels to see the image
                #e.g. with the application Gimp2
                figure.set_size_inches(320, 32)
               
                # If we have generated all the family of dynamical systems,
                # then we will print them all together
                plt.plot(lox,loy,",w")
               
                # Uncomment if you just want to diplay them
                #plt.show()
               
                # This will save the image associated with the dynamical system in a file
                # You can create a video using VirtualDub.exe
                plt.savefig("dynamical_system_inv_1_exp_"+str(expo)+"_D"+str(MaxD)+"_n"+str(ntestlimit)+"_icx_icy_"+str(icx)+"_"+str(icy)+"_Acot5_5.png")
               
nds()



Update and caveat: I doubted about adding this to the post, because it is quite pareidolic and can be misunderstood, but just for fun! It is the following family of dynamical systems:

$$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot |\sin{D}|, Re(f(x_{n} + y_n i)) \cdot |\cos{D}|)\}$$

Where $f(z)=\frac{1}{(z+1)^2}$ (so $t=2$ is fixed). As usual, the control parameters $D$ and $n$ will run through $D \in [0,2\pi]$ in $8 \cdot 10^3$ steps, and $n \in [0,8 \cdot 10^3]$, and the images are showing the region $z=[+/-5]+[+/-5]i$, for the initial condition $(x_0=0,y_0=0)$.  The only difference from the inital model is the absolute value applied to the sine and cosine functions, just that.

So this is the generated image:

Now if we rotate left $90$ degrees, and zoom in the center of the image, we have a peculiar pareidolic version of the Kubrick's movie "2001: A Space Odyssey" starchild (well, sort of, at least in my case when I saw it I remembered the reference to the movie), so for your viewing pleasure:


Sunday, January 14, 2018

Studying discrete-time dynamical systems (VI): mathematical models for symmetry and asymmetry and some biological coincidences

Let us do today some experimental mathematics. In my former posts  I have shown some stages of stochastic systems that look very "organic". Indeed they seem to have organic-like symmetrical and asymmetrical properties: while there is a clear axis of symmetry, both sides are not exactly symmetrical. Verifying the content of both sides of the axis it is possible to see that the main structures are symmetrical but there are small sub-structures that are only in one of the sides. It is very similar to the asymmetry that can be found in Nature, e.g. the human body looks externally symmetrical in a global level, but it is not. The human face, limbs, etc. are not symmetrical indeed, they are not exact mirrored copies. 

The mathematical models I was able to find seem to have some of these characteristics, so I am showing in the current post some biological coincidences of some stages of them and real life organic structures at different scale levels. For instance, jellyfish umbrella-shaped bells, sand wasp bodies and tardigrade limbs can be reflected at some extent with this type of complex systems. The pictures I am using at the right side of the images are just for the sake of completeness (they belong to their respective owners, I do not own them, if there is any problem I will remove them, so just please let me know):

1. The mathematical model (left) explained in this previous post shows similarities with (right) jellyfish structures (e.g. umbrella-shaped bell, including bioluminiscence patterns), the sand wasps and other similar insect bodies and the limbs of tardigrades. Click to enlarge:


2. The mathematical model (left) explained in this previous post shows similarities with some families of moths with rounded wings (right). Click to enlarge:


It would be great if these complex systems were able to be humbly used at some extent for some field of biology. For now they are just coincidences, I hope the readers will enjoy the samples.

3. A slight variation of the mathematical model explained at point 1 shows similarities with the bodies of some families of moths (e.g. Acherontia atropos). Click to enlarge:


This is the original image in full version. Click to enlarge:


I will try to expand this chimerical bestiary as much as possible.

Wednesday, December 20, 2017

Studying discrete-time dynamical systems (V): generalization of dynamical systems able to produce spirals and clusters of points

Following my previous post, here is a new beauty I was able to find recently. The following is an animation of the families: $$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(f(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ Where $f(z)=\frac{1}{(z^{1.5}+1)^t}, t$ from $1$ to $2$ by $0.01$ steps. Each frame is a complete family plot, calculated using $D \in [0,8 \cdot 10^3], n \in [0,8 \cdot 10^3]$ zoom into region $z=[+/-5]+[+/-5]i$. It took me some days to render it as well:


And this is the family of dynamical systems of $f(z)=\frac{1}{(z^{1.5}+1)^0.9}$, Rorschach test anyone?





I am trying to find some academic group working on this kind of complex maps to share these "stochastic" models. If I am able to make some advances on this topic I will post the news here. Crossing fingers!

Thursday, November 30, 2017

Studying discrete-time dynamical systems (IV): generalization of dynamical systems able to produce spirals and clusters of points

Here is an animation of the families mentioned in my previous post: $$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(f(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ Where $f(z)=\frac{1}{(z+1)^t}, t$ from $1$ to $2$ by $0.01$ steps. Each frame is a complete family plot, calculated using $D \in [0,8 \cdot 10^3], n \in [0,8 \cdot 10^3]$ zoom into region $z=[+/-5]+[+/-5]i$. It took me some days to render it:


It is visually astonishing! Still I did not receive any feedback regarding my question at Mathematics Stack Exchange about the dynamical systems that make these beauties. As far as I know it is the first time they are shown.

Monday, November 27, 2017

Studying discrete-time dynamical systems (III): generalization of dynamical systems able to produce spirals and clusters of points

In a previous post I wondered if there were any other dynamical systems able to generate spirals and clusters of points as the found expression:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$
 To answer this, first I had a look to the image of the complete plot of the family of spirals: 


We can see that there is a quasi-horizontal symmetry there. Both the "upper body" ($z=a+bi , b \ge 0$) and "lower body" ($z=a+bi , b \lt 0$) are quite symmetrical, but not exactly equal. The structure also has three vertical bodies / structures, a right side arc-like structure, a central body and a left-side body. These are the regions where the density of points is bigger.

The density of points make this structures to look like they are "glowing" (e.g. they remind a bioluminiscent jellyfish). The areas where there is more density of points seems to "glow" compared with other areas of lower density.

First question: who is responsible for making the spirals and who is responsible for deciding the location and distortion of the spirals?

To answer this, first I need to amend the expression:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ And make it more general, it is: $$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(f(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ So indeed: the family of dynamical systems is even bigger. Depending on the complex function that is used we will have a different structure. All of them have similar transformations but the results depend obviously on the selected function. So basically my original question was about the above more general family of dynamical systems when $f(z)=\psi(z)$ is specifically the digamma function.

Now, going back to the digamma expression:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ Let us split the expression into two versions:

Version 1: No sine or cosine are applied:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)), Re(\psi(x_{n} + y_n i)))\}$$and

Version 2: No function is applied, only the crossing of $Re$ and $Im$ (the meaning of "crossing" is as explained in the question):
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(x_{n} + y_n i)\cdot\cos{\frac{D}{f}})\}$$ Version 1 will provide the following result:



Version 2 will provide the following result (for the same $n$ iterations and $f$ and $D$ parameters):



So the conclusion is: the sine and cosine expressions provide the observed horizontal quasi-symmetry of the dynamical systems, and the digamma function was providing the spiral-like structures. The combination of both expressions provide a final structure in which the spirals of the digamma function are transformed, displaced and distorted by the sine and cosine transformations.
Second question: what do we obtain if we replace the digamma function by other functions?

My observations are that we obtain different "glowing" systems, but they share similar characteristics: horizontal quasi-symmetry, and three vertical structures (left, center, right), in other words, regions in which the density of points is bigger.

These are some of those beauties:

Scaled complementary error function, $f(z)=exp(z^2) \cdot erfc(z)$. (see Python scipy.special functions list, Steven G. Johnson, Faddeeva W function implementation). $D \in [0,10^4], n \in [0,10^4]$ zoom into region $z=[+/-5]+[+/-5]i$:



$f(z) = $ Exponential integral $E_1$ of complex argument $z$. (see Python scipy.special functions list). $D \in [0,10^3], n \in [0,10^4]$ zoom into region $z=[+/-1]+[+/-4]i$:



$f(z)=$ Exponential integral $E_i$. (see Python scipy.special functions list). $D \in [0,10^3], n \in [0,10^4]$ zoom into region $z=[+/-5]+[+/-5]i$:



Finally, my favorite families so far, $f(z)=\frac{1}{(z+1)^t}, t \in \Bbb R$. $D \in [0,8 \cdot 10^3], n \in [0,8 \cdot 10^3]$ zoom into region $z=[+/-5]+[+/-5]i$::

$f(z)=\frac{1}{(z+1)}$:



$f(z)=\frac{1}{(z+1)^{1.1}}$:



$f(z)=\frac{1}{(z+1)^{1.2}}$:



$\cdots$

$f(z)=\frac{1}{(z+1)^{2}}$:



$\cdots$

$f(z)=\frac{1}{(z+1)^{3}}$:



I just can say that these structures are visually beautiful. I will keep my question regarding these families of dynamical systems open at MSE in hope that somebody will be able to provide more insights about why this transformation / mapping works in this way.

Wednesday, November 15, 2017

Trying to understand the maths behind the movements of a two-legged air dancer (aka skydancer, tube man)

I am trying to understand the maths behind the movements of a two-legged air dancer, aka skydancer aka tube man. Well, I mean these cheery friends here:

 

(This is a mirror of my question on this topic at MSE) The following discrete time algorithm is my devised solution, but I am sure that other much better approaches are possible. Indeed it is not very realistic and I am not very happy. I would like to learn other solutions more closer to reality and related with dynamical systems, complex dynamics or fluid dynamics:

1. Take for instance a unit square, same width and height. Then generate two partially (explanation of the "partially" wording is done below) random walks with these rules:

2. The first partially random walk $R1$ starts at $(0,0)$ and finishes at $(1,1)$

3. The second random walk $R2$ starts at $(0,1)$ and finishes at $(1,0)$

4. We will generate $n$ random points sequentially for $R1$ as follows. The first point is $(x_0,y_0)=(0,0)$, then we will take a random point from the intervals $([0,1],[0,1])$, it will be $(x_1,y_1)$. Then $(x_2,y_2)$ will be a random point from the intervals $([x_2,1],[y_2,1])$, then next point $([x_3,1],[y_3,1])$, etc. So we are approaching $(1,1)$ on each iteration. Usually $n=20$ iterations is good enough to be visually very close to $(1,1)$ so the simulation is good enough. The reason of the "partially" wording is that on every iteration we reduce the intervals, so it is not always a random point of the square unit, so it is not totally random.

5. In the same fashion we will generate points for $R2$ but the intervals will approach on each iteration to $(1,0)$.

6. Now we will create segments joining the sequential points of $R1$ and same for $R2$.

7. As both paths are crossing completely each diagonal of the square unit, there is an intersection point between them. Mark with a disk the intersection of both random paths and display the unit square including the segments and the intersection.

8. Execution: the algorithm can be looped for an infinite time of $t$ units in a loop, and each iteration $t$ will generate the aspect of the air dancer for that unit of time. So in my case the algorithm makes a kind of discrete system.

And this is how looks my air dancer simulation:



And another one with two dancers:



Some facts and trivia:

1. Due to the way that the segments are generated, the lengths of them are not totally random, see the square line picking problem.

2. The above solution might be defined as a kind of discrete map where each unit of time $t$ is independent of the status of the previous unit of time (but this can be modified to make movements depend on the previous unit of time at some extent so it would be a standard discrete-time map).

3. The algorithm works for any square, not only the unit square, it is just an example. They can be any desired width and height.

4. The legs of the air dancers are usually fixed to the ground, but not the arms, and that is a part of my simulation that is not very realistic (apart from others) although I think that it would be more or less easy to enhance this point specifically.

I would like to simulate this kind of objects with a more realistic approach, but initially I did not find a reference on Internet. My guessing is that it must be very related to dynamical systems and more specifically to fluid dynamics, so it might be possible to express the movement by a discrete (map) or continuous (based on derivatives) dynamical-system.

Below is the Python code for my simulation (also available at my account at repl.it). It will work in, for instance, a personal computer using Python Anaconda if the pygame library is downloaded as usual. Feel free to use and modify it and enjoy!

#!/usr/bin/env python
import pygame
import random
from collections import defaultdict
from time import sleep

# This script aims to simulate a tube man aka air dancer aka skydancer
# https://en.wikipedia.org/wiki/Tube_man

# The code will generate a TOTWIDTHxTOTHEIGHT pixels screen
# where HOBJECTSxVOBJECTS skydancers are generated
# each one with a similar WIDTHxHEIGHT
# the body of each skydancer is a disk of radius = RBIG pixels
# the simulation runs until ESC is pressed and each
# simulation time unit is shown during DELAYSEC seconds
# then the screen is wiped and the next time unit image is calculated

# Initialization constraints
HOBJECTS=2 #number of air dancers per line
VOBJECTS=1 #number of rows
WIDTH, HEIGHT = 200, 200
TOTWIDTH, TOTHEIGHT = WIDTH*HOBJECTS, HEIGHT*VOBJECTS
RBIG = 17#10#17
RSMALL = 8#8
DELAYSEC=0.3

def air_dancer_simulator():
    pygame.init()
    win=pygame.display.set_mode((TOTWIDTH, TOTHEIGHT))

    keysPressed = defaultdict(bool)

    def ScanKeyboard():
        while True:
            # Update the keysPressed state:
            evt = pygame.event.poll()
            if evt.type == pygame.NOEVENT:
                break
            elif evt.type in [pygame.KEYDOWN, pygame.KEYUP]:
                keysPressed[evt.key] = evt.type == pygame.KEYDOWN

    # detects the intersection point (x,y) of two lines, if exists.
    def line_intersection(line1, line2):
        xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
        ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])

        def det(a, b):
            return a[0] * b[1] - a[1] * b[0]

        div = det(xdiff, ydiff)
        if div == 0:
            raise Exception('lines do not intersect')

        d = (det(*line1), det(*line2))
        x = det(d, xdiff) / div
        y = det(d, ydiff) / div
        return x, y
   
    bClearScreen = True
    pygame.display.set_caption('Skydancer aka tube man aka air dancer simulator')
    loop=0
   
    #Seed points of the air dancer sequential segments to be generated are (0,0) and (0,1)
    x_currpos_incdiag = 0
    y_currpos_incdiag = 0
    x_currpos_decdiag = 0
    y_currpos_decdiag = HEIGHT-1
   
    test_points_inc=[[x_currpos_incdiag,y_currpos_incdiag]]
    test_points_dec=[[x_currpos_decdiag,y_currpos_decdiag]]
   
    frames=0
    # the simulation will run infinite discrete units of time (ESC to finish the simulation)
    while True:
        # draw biash
        for biash in range(0,HOBJECTS):
            for biasv in range(0,VOBJECTS):
                loop=1
                while True:
                    # each dancer is finished in 20 iteracions
                    if loop%20==0:
                       
                        found=False
                        for i in range(0,len(test_points_inc)-2):
                            for j in range(0,len(test_points_dec)-2):
                                try:                       
                                    x,y=line_intersection([test_points_inc[i],test_points_inc[i+1]],[test_points_dec[j],test_points_dec[j+1]])
                                    if x>=0 and x<=WIDTH and y>=0 and y<=HEIGHT:
                                        if (test_points_inc[i][0]<=x) and (x<=test_points_inc[i+1][0]):
                                            if (test_points_inc[i][1]<=y) and (y<=test_points_inc[i+1][1]):
                                                if (test_points_dec[j][0]<=x) and (x<=test_points_dec[j+1][0]):
                                                    if (test_points_dec[j][1]>=y) and (y>=test_points_dec[j+1][1]):
                                                        # Detected interection of two segments, will be the body of the air dancer
                                                        pygame.draw.circle(win, (255,255,255),(int(x)+(biash*WIDTH),int(y)+(biasv*HEIGHT)),RBIG)       
                                                        found=True
                                                        break
                                except:
                                    pass
                                   
                            if found==True:
                                break
                       
                        #pygame.display.flip()
                        #sleep(5)
       
                        #initial positions
                        x_currpos_incdiag = 0
                        y_currpos_incdiag = 0
                        x_currpos_decdiag = 0
                        y_currpos_decdiag = HEIGHT-1
                        test_points_inc=[[x_currpos_incdiag,y_currpos_incdiag]]
                        test_points_dec=[[x_currpos_decdiag,y_currpos_decdiag]]
                        break
                       
                    x_newpos_incdiag = random.randint(x_currpos_incdiag,WIDTH-1)
                    y_newpos_incdiag = random.randint(y_currpos_incdiag,HEIGHT-1)
                    x_newpos_decdiag = random.randint(x_currpos_decdiag,WIDTH-1)
                    y_newpos_decdiag = random.randint(0,y_currpos_decdiag)
                   
                    test_points_inc.append([x_newpos_incdiag,y_newpos_incdiag])
                    test_points_dec.append([x_newpos_decdiag,y_newpos_decdiag])
               
                    pygame.draw.line(win, (150,200,255),(x_currpos_incdiag+(biash*WIDTH),y_currpos_incdiag+(biasv*HEIGHT)),(x_newpos_incdiag+(biash*WIDTH),y_newpos_incdiag+(biasv*HEIGHT)))
                    pygame.draw.line(win, (150,200,255),(x_currpos_decdiag+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)),(x_newpos_decdiag+(biash*WIDTH),y_newpos_decdiag+(biasv*HEIGHT)))
                   
                    # Add horizontal and vertical lines too
                    #pygame.draw.line(win, (255,255,255),(x_currpos_incdiag+(biash*WIDTH),y_currpos_incdiag+(biasv*HEIGHT)),(x_currpos_incdiag+(biash*WIDTH),0+(biasv*HEIGHT)))
                    #pygame.draw.line(win, (255,255,255),(x_currpos_incdiag+(biash*WIDTH),y_currpos_incdiag+(biasv*HEIGHT)),(x_currpos_incdiag+(biash*WIDTH),HEIGHT+(biasv*HEIGHT)))
                    #pygame.draw.line(win, (255,255,255),(x_currpos_decdiag+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)),(0+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)))
                    #pygame.draw.line(win, (255,255,255),(x_currpos_decdiag+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)),(WIDTH+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)))
                   
                    # Add circles on each beginning and end of segment
                    #pygame.draw.circle(win, (150,200,255),(x_currpos_incdiag+(biash*WIDTH),y_currpos_incdiag+(biasv*HEIGHT)),RSMALL)
                    #pygame.draw.circle(win, (150,200,255),(x_currpos_decdiag+(biash*WIDTH),y_currpos_decdiag+(biasv*HEIGHT)),RSMALL)       
                       
                    x_currpos_incdiag = x_newpos_incdiag
                    y_currpos_incdiag = y_newpos_incdiag
                    x_currpos_decdiag = x_newpos_decdiag
                    y_currpos_decdiag = y_newpos_decdiag
               
                    win.lock()
                    win.unlock()
                    ScanKeyboard()

                    if keysPressed[pygame.K_ESCAPE]:
                        pygame.image.save(win, "intersections.png")
                        break
                   
                    loop = loop+1
                   
                    if keysPressed[pygame.K_ESCAPE]:
                        pygame.image.save(win, "intersections.png")
                        break
           
                if keysPressed[pygame.K_ESCAPE]:
                    pygame.image.save(win, "intersections.png")
                    break
           
            if keysPressed[pygame.K_ESCAPE]:
                pygame.image.save(win, "intersections.png")
                break
       
        if keysPressed[pygame.K_ESCAPE]:
            pygame.image.save(win, "intersections.png")
            break       
           
        pygame.display.flip()
        sleep(DELAYSEC)
       
        # Uncomment to save the frames, then you can create an animated gif or video with
        # your favorite program (e.g. VirtualDub)
        #pygame.image.save(win, "airdancer_" + str(frames) + ".png")
       
        frames=frames+1
       
        if bClearScreen:
            win.fill((0, 0, 0))
                   
air_dancer_simulator()