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.