function [a,b,c,d,e] = stripes(xsize,ysize,iterations,speed)
%
% [a,b,c,d,e] = stripes(xsize,ysize,iterations,speed)

p1 = 0.04;  p2 = 0.06;  p3 = 0.04;
diff1 = 0.009;  diff2 = 0.2;
arand = 0.02;
a = zeros(200,200); b = a; c = a; d = a; e = a;
da = a; db = a; dc = a; dd = a; de = a; ai = a;
ainit = 0;  binit = 0;  cinit = 0;  dinit = 0; einit = 0;

for row = 1:xsize
   ainit = p2/(2*p1);
   binit = ainit;
   cinit = 0.02*ainit*ainit*ainit/p2;
   dinit = ainit;
   einit = ainit;
   for col = 1:ysize
      a(row,col) = ainit;
      b(row,col) = binit;
      c(row,col) = cinit;
      d(row,col) = dinit;
      e(row,col) = einit;
      ai(row,col) = 1 + frand(-0.5*arand,0.5*arand);
   end
end

for k = 0:(iterations-1)
   for iindex = 0:(xsize-1)
      ka = -p1 - 4*diff1;
      kc = -p2;
      kd = -p3 - 4*diff2;
      iprev = rem((iindex+xsize-1),xsize);
      rowprev = iprev+1;
      inext = rem((iindex+1),xsize);
      rownext = inext+1;
      row = iindex+1;
      for jindex = 0:(ysize-1)
         jprev = rem((jindex+ysize-1),ysize);
         colprev = jprev+1;
         jnext = rem((jindex+1),ysize);
         colnext = jnext+1;
         col = jindex+1;
         aval = a(row,col);
         bval = b(row,col);
         cval = c(row,col);
         dval = d(row,col);
         eval = e(row,col);
         temp1 = 0.01*aval*aval*eval*ai(row,col);
         temp2 = 0.01*bval*bval*dval;
         
         dda = a(row,col) + a(row,colnext) + a(rowprev,col) + a(rownext,col);
         ddb = b(row,col) + b(row,colnext) + b(rowprev,col) + b(rownext,col);
         ddd = d(row,col) + d(row,colnext) + d(rowprev,col) + d(rownext,col);
         dde = e(row,col) + e(row,colnext) + e(rowprev,col) + e(rownext,col);
         
         da(row,col) = aval*ka + diff1*dda + temp1/cval;
         db(row,col) = bval*ka + diff1*ddb + temp2/cval;
         dc(row,col) = cval*kc + temp1 + temp2;
         de(row,col) = dval*kd + diff2*ddd + p3*aval;
         de(row,col) = eval*kd + diff2*dde + p3*bval;
      end
   end
   
   for row = 1:xsize
      for col = 1:ysize
         a(row,col) = a(row,col) + speed*da(row,col);
         b(row,col) = b(row,col) + speed*db(row,col);
         c(row,col) = c(row,col) + speed*dc(row,col);
         d(row,col) = d(row,col) + speed*dd(row,col);
         e(row,col) = e(row,col) + speed*de(row,col);
      end
   end
end
