rot = function(theta) matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)),2,2,byrow=T)
mag = function(mag) matrix(c(mag,0,0,mag),2,2)

x = -1*c(9.00E-02, 0.165272, 0.2259414, 0.3012552, 0.3619247, 0.416318, 0.460251, 0.5313808, 0.5711297, 0.6213389, 0.6569037, 0.6569037, 0.6297071, 0.6004184, 0.6150628, 0.6569037, 0.6966527, 0.7322176, 0.7719665, 0.7866109, 0.7782427, 0.8117155, 0.8870293, 0.9497908, 0.9748954, 0.9476988, 0.8640167, 0.8514645, 0.8410042, 0.8117155, 0.7468619, 0.7259414, 0.7133891, 0.6464435, 0.5962343, 0.539749, 0.4979079, 0.5041841, 0.4916318, 0.4497908, 0.4225941, 0.3200837, 0.251046, 0.2217573, 0.1882845, 0.1192469, 7.95E-02, 8.79E-02, 7.74E-02, 8.79E-02, NA, 0.1171548, 0.1338912, 0.1569038, 0.1778243, 0.1882845, 0.1757322, 0.1527197, 0.1317992, 0.1150628, 0.1087866, 0.1150628)
y = c(0.9113082, 0.902439, 0.8691796, 0.8159645, 0.7738359, 0.7317073, 0.6984479, 0.6784922, 0.6651884, 0.6385809, 0.6031042, 0.578714, 0.5631929, 0.5498891, 0.5277162, 0.4678493, 0.4101996, 0.345898, 0.3436807, 0.3126386, 0.2993348, 0.2815965, 0.2594235, 0.2439024, 0.2195122, 0.1906874, 0.1685144, 0.1263858, 4.88E-02, 3.10E-02, 0.1064302, 0.1818182, 0.2062084, 0.2084257, 0.2727273, 0.2949002, 0.3281596, 0.286031, 0.286031, 0.323725, 0.3880266, 0.4634146, 0.5410199, 0.5321507, 0.6274945, 0.7250555, 0.8226164, 0.8403547, 0.8957871, 0.9068736, NA, 0.8558758, 0.8736142, 0.8780488, 0.864745, 0.842572, 0.8181818, 0.8048781, 0.8070953, 0.8181818, 0.8403547, 0.8536585)
x = x-mean(x[!is.na(x)])
y = y-mean(y[!is.na(x)])
px = as.matrix(cbind(x,y))
px = px%*%rot(.8)%*%mag(3)

nfish=30
nlines=40
xmin=0.1
xmax=5.9
xvec = seq(xmin,xmax,length=nfish)
liney = seq(7.9,0.1,length=nlines)
linerot = seq(0,1,length=nlines)
px = px/20

pdfON=0
if(pdfON==1){pdf('c:/GurarieFish.pdf',12,16)}else{windows(6,8)}

par(mar=rep(0,4))
plot(0,type='n',axes=F,xlim=c(0,6),ylim=c(0,8),xlab='',ylab='')
for(iline in 1:nlines)
{
  for(ifish in 1:nfish)
  {
    xval = xvec[ifish]+rnorm(1,0,linerot[iline/5])
    yval = sin(xval*3)/10 + liney[iline]+rnorm(1,0,linerot[iline]/5)
    rot1 = rot(-cos(xval*3)/8 + rnorm(1,0,linerot[iline]/2))
    mag1 = mag(exp(rnorm(1,0,linerot[iline]/2)))
    newvals = px%*%mag1%*%rot1
    colval = floor(runif(1,60-linerot[iline]*15, 60+linerot[iline]*15))
    polygon(xval+newvals[,1], yval+newvals[,2], border=paste('grey',colval,sep=''),col=c(paste('grey',colval+20,sep=''),'white'))
  }
}
if(pdfON==1){dev.off()}