I have a program to draw a grid of ellipses with a uniform phase distribution. However, it is very slow.

I'd like my code to be faster, so that I can use, for example, `N = 150`

, `M = 150`

.
How can I speed up this code?

```
N = 10;
M = 10;
y = 1;
x = 1;
a = 1;
b = 2;
for k = 1:N
for m = 1:N
w = rand(1,1);
for l = 1:N
for s = 1:N
if(((l-x)*cos(w*pi)+(s-y)*sin(w*pi)).^2/a^2 + (-(l-x)*sin(w*pi) + (s-y)*cos(w*pi)).^2/b.^2 <= 1)
f(l,s) = 1*(cos(0.001)+i*sin(0.001));
end
end
end
y = y+4;
end
y = 1;
x = x+5;
end
image(arg(f),'CDataMapping','scaled');
```

This is what the code produces:

Updated:

```
N = 10;
M = 10;
y = 1;
x = 1;
a = 1;
b = 2;
for x = 1:5:N
for y = 1:4:N
w = rand(1);
for l = 1:N
for s = 1:N
if(((l-x).*cos(w.*pi)+(s-y).*sin(w.*pi)).^2/a.^2 + (-(l-x).*sin(w.*pi) + (s-y).*cos(w.*pi)).^2/b.^2 <= 1)
f(l,s) = cos(0.001)+i.*sin(0.001);
end
end
end
end
y = 1;
end
image(arg(f),'CDataMapping','scaled');
```

There are many things you can do to speed up the computation. One important one is to remove loops and replace them with vectorized code. Octave works much faster when doing many computations as once, as opposed to one at a time.

For example, instead of

one can write

However, here we still do too much work because we compute

`index`

at locations that we know will be outside the extend of the ellipse. Ideally we'd find a smaller region around each point`(x,y)`

that we know the ellipse fits into.Another important thing to do is preallocate the array.

`f`

grows within the loop iterations. Instead, one should create`f`

with its final size before the loop starts.There are also many redundant computations being made. For example

`w.*pi`

is computed multiple times, and the`cos`

and`sin`

of it too. You also assign`cos(0.001)+i.*sin(0.001)`

to output pixels over and over again, this could be a constant computed once.The following code runs in MATLAB in a tiny fraction of a second (though Octave will be a lot slower). I've also separated

`N`

and`M`

properly (so the output is not always square) and made the step sizes a variable for improved understanding of the code. I'm assigning`1`

to the ellipse elements, you can replace them by your constant by multiplying`f = f * (cos(0.001)+i*sin(0.001))`

.