R script to make a flow field based on a rational function in complex space

flow.r

# Set random number seed for reproducibility
set.seed(2210)

# Number of starting points
N=5000

# Describe the flow function
flow <- function(z,bs) {    
    fz <- ((z+bs[1])*(z+bs[2])*(z+bs[3])*(z+bs[4])*(z+bs[8]))/
          ((z+bs[5])*(z+bs[6])*(z+bs[7]))
    fz / Mod(fz);
    }

# Choose random parameters for function (zeros for each polynomial)
bs <- runif(1,0.0,1)*(runif(15,-2,2) + 1i*runif(15,-2,2))

# Choose random starting points
z_start <- runif(N,-1,1) + 1i * runif(N,-1,1)

# Start the png file and set up the plot
png("flow.png",width=1600,height=2000,type="cairo")
par(mar=rep(10,4),bg="#dddddd")
plot(NA,asp=1,axes=F,ann=F,xlim=c(-0.8,0.8),ylim=c(-1,1))

# Draw forwards paths
z <- z_start
for(i in 1:100){
  points(z,pch=19,col=hsv(1,0,0,.1),cex=0.5)
  z <- z + .01*runif(N)*flow(z,bs)
  }

# Draw backwards paths
z <- z_start
for(i in 1:100){
  z <- z - .01*runif(N)*flow(z,bs)
  points(z,pch=19,col=hsv(1,0,0,.1),cex=0.5)
  }

# Close the graphics device
dev.off()


# For those not using R, the parameters for the flow function are:
# bs[1] = -0.3373520+0.6742209i  
# bs[2] = 1.1404262-0.9528159i 
# bs[3] = -0.0087374-1.3331986i
# bs[4] = -0.6954777-1.0817921i  
# bs[5] = 0.3198156+0.2273527i  
# bs[6] = 0.4125582+0.1261198i
# bs[7] = 1.3838665+1.2481734i 
# bs[8] = -0.8905171-0.5191771i
#
# The flow field is then calculated at each position 'z' in the complex plane by:
# 
# fz <- ((z+bs[1])*(z+bs[2])*(z+bs[3])*(z+bs[4])*(z+bs[8]))/((z+bs[5])*(z+bs[6])*(z+bs[7]))
# 
# Finally the modulus of the flow is set to 1 by diving fz by Mod(fz) to get the final flow.
# 
# N=5000 starting points are chosen at random.  Then each is iterated 100 times forwards (using z=z+0.01*fz) 
# then reset and run 100 times backwards (using z=z-0.01*fz) to create the final image.
#
#
#