I define a piecewise function `Lens(x_,y_)`

, which contains `Subscript(Y, i)`

. `Y()`

is a different function for different values of i.

But some Y() calculations will produce a small imaginary part of 0.000000I, which will tell me when I run my segmented function `Less::nord: tried to use 37344.8-3240.47i for an invalid comparison.`

and the stack trace `30102.9 + 0. I <= 31825 < 37344.8 - 3240.47 I`

.

And then I take my piecewise function `Lens(x_,y_)`

and convert it to an array `LensArray`

and plot it using the `MatrixPlot`

function. The `MatrixPlot`

result is correct.

The imaginary part is supposed to be extremely small, so if I put `Abs()`

in front of `Y()`

or `Chop()`

in front of `Y()`

it shouldn’t matter, but it does make a big difference when I MatrixPlot.

So I guess at some points the imaginary part greatly affects the compare of the piecewise function, but if I Plot Y() alone, the function can be plotted, which means that the imaginary part of Y() is small enough to be ignored. I’m not quite sure why the piecewise function is making this kind of comparison: `Yrange((i)) <= y < Yrange((i + 1))`

indicates it’s an invalid comparison.

And then I check the matrix of `LensArray`

, and I see that the matrix is fine, but when I take the discrete Fourier transform of the matrix of `LensArray`

it tells me that `Fourier::fftl : LensArray is not a non-empty list, it's not a rectangular array of numeric quantities.`

I checked the stack trace, but I see the matrix `LensArray`

is really an array of numbers.

So how can I take the `Fourier`

to the `LensArray`

correctly?

The defining function of `LensArray`

and `MatrixPlot`

```
Clear("Global`*")
δ = 1.2131226805013586 10^-6;
β = 4.75349804 10^-8;
χ = -2 δ + 2 I β;
δ1 = -1.2131226805013586 10^-6;
δ2 = (1.2131226805013586 10^-6)/(1 - 1.2131226805013586 10^-6);
f = 5151736.70155;
xStart = -50;
xEnd = 50;
Δx = 1;
Nx = (xEnd - xStart)/Δx + 1;
Nz = 2000 + 1;
Maxz = 50000;
Table(Subscript(Axr,
i) = ((-1)/9) (q1((i))^2 - q1((i)) q2((i)) + q2((i))^2), {i, 1,
62});
Table(Subscript(Bxr, i) =
Module({δs}, δs =
If(OddQ(i), δ1, δ2); ((-(q1((i)) -
q2((i))))/(4δs)) x^2 + (1/54) (q1((i)) +
q2((i))) (9 q1 ((i)) q2((i)) - 2 (q1((i)) + q2((i)))^2)), {i,
1, 62});
Table(Subscript(β, i) =
ArcCos(Subscript(Bxr, i)/(-Subscript(Axr, i))^((3/2))), {i, 1, 62});
Table(Subscript(Y, i) =
Abs(-(2 Sqrt(-Subscript(Axr, i)) Cos(
Subscript(β, i)/3) - (1/3) (q1((i)) + q2((i))))), {i,
1, 62});
Len1 = {χ,
0 <= y < (Sqrt(x^2 + f^2) - f)/δ + 10 &&
0 < y < Yrange((1)) && -50 <= x <= 50};
Lens(x_, y_) =
Piecewise(
Join({Len1},
Table({χ,
Subscript(Y, i) + Apcie((i)) <= y <
Subscript(Y, i + 1) + Apcie((i + 1)) &&
Yrange((i)) <= y < Yrange((i + 1)) && -50 <= x <= 50}, {i, 1,
61, 2})));
Lensarray = Array(Lens, {Nx, Nz}, {{xStart, xEnd}, {0, Maxz}});
MatrixPlot(Abs(Lensarray), ColorFunction -> "Monochrome",
PlotLegends -> Automatic)
Fourier(Lensarray((1)))
```

The definition of the matrix:

```
q1={5151536.701550,2499990.000000,1614774.131776,1169990.000000,904742.171199,729990.000000,605711.526739,513990.000000,443182.277830,387990.000000,343258.780795,306590.000000,275934.940760,249990.000000,227506.115311,208190.000000,191077.508363,176290.000000,162848.416323,151190.000000,140369.167291,130990.000000,121982.805265,114190.000000,106694.176102,100290.000000,93923.837408,88490.000000,82991.123566,78390.000000,73613.061532,69670.000000,65476.568200,62140.000000,58472.599719,55590.000000,52312.925767,49820.000000,46892.253210,44720.000000,42084.946419,40180.000000,37775.450513,36102.000000,33901.020142,32420.000000,30386.185769,29060.000000,27166.524869,25970.000000,24198.228993,23120.000000,21453.974356,20470.000000,18886.440141,17985.000000,16481.195581,15650.000000,14213.422727,13445.000000,12071.827247,11355.000000};
q2={2500000.00,1615000.00,1170000.00,905000.00,730000.00,606000.00,514000.00,443500.00,388000.00,343600.00,306600.00,276300.00,250000.00,227900.00,208200.00,191500.00,176300.00,163300.00,151200.00,140850.00,131000.00,122500.00,114200.00,107240.00,100300.00,94500.00,88500.00,83600.00,78400.00,74250.00,69680.00,66140.00,62150.00,59150.00,55600.00,53010.00,49830.00,47600.00,44730.00,42800.00,40190.00,38500.00,36112.00,34630.00,32430.00,31120.00,29070.00,27905.00,25980.00,24940.00,23130.00,22198.00,20480.00,19637.00,17995.00,17234.00,15660.00,14970.00,13455.00,12830.00,11365.00,10798.00};
Apcie={422.115008,432.115008,900.361323,910.361323,1440.287177,1450.287177,2041.061643,2051.061643,2697.077993,2707.077993,3403.093596,3413.093596,4158.904889,4168.904889,4972.510513,4982.510513,5841.823636,5851.823636,6768.969999,6778.969999,7754.012562,7764.012562,8813.501181,8823.501181,9926.509448,9936.509448,11105.652528,11115.652528,12346.394885,12356.394885,13645.086416,13655.086416,14984.837936,14994.837936,16358.648861,16368.648861,17766.072000,17776.072000,19196.022557,19206.022557,20642.637485,20652.637485,22102.984088,22112.984088,23571.045417,23581.045417,25050.097410,25060.097410,26540.372202,26550.372202,28034.705468,28044.705468,29536.538260,29546.538260,31048.730592,31058.730592,32567.164839,32577.164839,34090.612687,34100.612687,35618.226849,35628.226849};
Yrange={210.000000,657.983232,657.983232,1168.190124,1168.190124,1738.760438,1738.760438,2368.783813,2368.783813,3048.297198,3048.297198,3778.152836,3778.152836,4562.789578,4562.789578,5405.002150,5405.002150,6303.407313,6303.407313,7259.802708,7259.802708,8281.207297,8281.207297,9369.325079,9369.325079,10512.672040,10512.672040,11724.528962,11724.528962,12993.333353,12993.333353,14318.518216,14318.518216,15672.238217,15672.238217,17065.723094,17065.723094,18483.818790,18483.818790,19921.076138,19921.076138,21377.186972,21377.186972,22841.963946,22841.963946,24314.859648,24314.859648,25798.572541,25798.572541,27292.143209,27292.143209,28788.731112,28788.731112,30297.098119,30297.098119,31811.535011,31811.535011,33333.742112,33333.742112,34858.785440,34858.785440,36389.282916};
```

If you want to look at the function `Subscript(Y, i)`

separately:

```
Clear("Global`*")
Print("Q1 = ", DecimalForm(q1 = 122500. - 517.194735, {20, 6}));
Print("Q2 = ", DecimalForm(q2 = 114200., {20, 6}));
Print("A24 = ",
DecimalForm(
A24 = 97.001875 (122500. - 517.194735)/122500., {20, 6}));
δ = -1.2131226805013586 10^-6;
Axr = -1/9 (q1^2 - q1 q2 + q2^2);
Bxr(x_) := -(q1 - q2)/(4 δ) x^2 +
1/54 (q1 + q2) (9 q1 q2 - 2 (q1 + q2)^2);
β(x_, bxr_) := ArcCos(bxr(x)/(-Axr)^(3/2));
Y(x_, β2_) :=
2 Sqrt(-Axr) Cos(β2(x, Bxr)/3) - 1/3 (q1 + q2);
Xd0 = 2 (-δ/(q1 - q2))^(1/2) (Bxr(0) + (-Axr)^(3/2))^(1/2);
DecimalForm(Y(A24/2, β), {12, 6})
Plot(-Y(x, β), {x, -Xd0, Xd0})
```

The imaginary part was answered,but I still can’t think of a solution to the above `Fourier`

problem.