Trying to insert an option to change the incidence angle for my simplest test case.

I don't think I've done it right. Essentially I need k = Abs[kx + i ky] and tk = Arg[kx + i ky] to be my input parameters for graphing cos they are the ones that stay consistent whereas the components kx and ky will change as k gets larger but the angle stays the same. I imagine that since Abs and Arg are functions that exist already in mathematica that this should be possible. But I can't find anything helpful to me in the documentation.

Currently this won't plot for tk = 0 which I need in order to check that the generalisation is working.

Any advice ye can give me on avoiding the singular matrices and 1/0 s?

Attachments:

This is interesting, but I am afraid it is too complicated to grasp in one sitting. One thing I can say is that sometimes I have been able to improve performance by many orders of magnitude by thinking about the problem differently. I mean, thinking of a different way to ask Mathematica the same question. But that is just generality. The only specific I can give right now is, as a first cut, to try using DiscretePlot:

DiscretePlot[f[x], {x, xmin, xmax, dx}]

This saves Mathematica the time to make decisions about where to evaluate, and it can give you a quick look at a complicated plot faster than you would get using Plot. Make sure you set Filling to False, otherwise the graph will look different than you expect.

I have a background in physics. If you send a reference to what you are trying to do written in human language, I may be able to digest it better. Best,

OL.

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...

<<Progress Update>>

Now I have successfully expanded this for two barriers (or a square potential well) and am attempting to expand fro arbitrarily many barriers. Getting errors I don't understand though. Also still can't delay evaluation without getting the wrong result.

Attachments:

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...

I'm not quite sure, but I guess generally it has to do with cases there the numbers become '0' and you have to do some division. So is it possible for your elements in your matrix to be zero? Just check what you expect, you can put some values in there and try. I put 3 random values for qx qy and l in there...

Also you can check the symbolic result, is that what you want?

Have a look at the following functions:

AbsArg
AngleVector
ReIm
ArcTan[x, y]

The first gives the magnitude and the 'angle' (argument) of a complex number at once. The second can construct a vector by a radius-angle pair The third gives you real and Imaginary part of a number as a pair. And perhaps use Arctan with two arguments; so it takes care of the 4 quadrants (the signs of x and y et cetera).

The problem is in your eigensystem function:

First symbolic finding the general eigenvalues and eigenvectors may involve certain assumptions, not valid for all cases! Take for example the simple case of this polynomial:

In the general case 3x^2/x is 3x but not for x=0 ! So 'underneath' the symbolic calculation there are always some assumptions (unavoidable in some sense). I'm not sure why the general solution of the eigensystem 'works' for your input but not the non-general case!

I don't see a difference!

I sped it up another factor 10 just now, but there is still lots to go. I'm sure.

Attachments:

I kinda figured you'd need more info. But I didn't want to spam everyone right away! Here is everything I have with annotations. In desperation I have also attached the file. Any tips at all would be appreciated.

P.S. I should say that the end goal of this will be to plot that value of T[n,a,d,kx,l,g] {k,-10,10} without it taking infinite time.

Remove["Global`*"];

$Assumptions = a > 0 && d > 0 && l > 0 && g > 0 && n \[Element] Integers;

Hamiltonian for region with SOI and Zeeman field:

Hsoz[qx_, l_, g_] = {{-g, 0, px + I py, -2 I l}, {0, g, 0, px + I py}, {px - I py, 0, -g, 
  0}, {2 I l, px - I py, 0, g}}

Evals[qx_, l_, g_] = Simplify[Eigenvalues[Hsoz[qx, l, g]]];

Ez1[qx_, l_, g_] = Evals[qx, l, g][[1]];
Ez2[qx_, l_, g_] = Evals[qx, l, g][[2]];
Ez3[qx_, l_, g_] = Evals[qx, l, g][[3]];
Ez4[qx_, l_, g_] = Evals[qx, l, g][[4]];

Evecs[qx_, l_, g_] = Simplify[Eigenvectors[Hsoz[qx, l, g]]];

Normalization below (I suspect there is a better way to do this)

n1[qx_, l_, g_] = g^2 + l^2 + Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];
n2[qx_, l_, g_] = g^2 + l^2 - Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];

Nz1[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez1[qx, l, g] + g)]^2)/
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez1[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz2[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez2[qx, l, g] + g)]^2)/
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez2[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz3[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez3[qx, l, g] + g)]^2)/
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez3[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];
Nz4[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez4[qx, l, g] + g)]^2)/
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez4[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];

wz1[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz1[qx, l, g]]) Evecs[qx, l, g][[1]]];
wz2[qx_, l_, g_] = 
  Simplify[(1./t[Nz2[qx, l, g]]) Evecs[qx, l, g][[2]]];
wz3[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz3[qx, l, g]]) Evecs[qx, l, g][[3]]];
wz4[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz4[qx, l, g]]) Evecs[qx, l, g][[4]]];

Conservation of energy:

qz3x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 - 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];
qz4x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 + 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];

Generic wave equation constructed:

M1[x_, kx_, l_, g_] = 
  Simplify[wz3[qz3x[kx, l, g], l, g]*E^(I qz3x[kx, l, g] x)];
M2[x_, kx_, l_, g_] = 
  Simplify[wz4[qz4x[kx, l, g], l, g]*E^(I qz4x[kx, l, g] x)];
M3[x_, kx_, l_, g_] = 
  Simplify[wz3[-qz3x[kx, l, g], l, g]*E^(-I qz3x[kx, l, g] x)];
M4[x_, kx_, l_, g_] = 
  Simplify[wz4[-qz4x[kx, l, g], l, g]*E^(-I qz4x[kx, l, g] x)];

And in matrix form:

Mz[x_, kx_, l_, g_] = ({
    {M1[x, kx, l, g][[1]], M2[x, kx, l, g][[1]], M3[x, kx, l, g][[1]],
      M4[x, kx, l, g][[1]]},
    {M1[x, kx, l, g][[2]], M2[x, kx, l, g][[2]], M3[x, kx, l, g][[2]],
      M4[x, kx, l, g][[2]]},
    {M1[x, kx, l, g][[3]], M2[x, kx, l, g][[3]], M3[x, kx, l, g][[3]],
      M4[x, kx, l, g][[3]]},
    {M1[x, kx, l, g][[4]], M2[x, kx, l, g][[4]], M3[x, kx, l, g][[4]],
      M4[x, kx, l, g][[4]]}
   });

Hamiltonian for region without SOI and Zeeman field:

L[x_, kx_] = 1./Sqrt[2.] ({{E^(I kx x), 0, -E^(-I kx x), 0},
{0, E^(I kx x), 0, -E^(-I kx x)},
{E^(I kx x), 0, E^(-I kx x), 0},
{0, E^(I kx x), 0, E^(-I kx x)}});

Generalising for potential barriers of width "d" every "a" units of distance

Sl[a_, kx_] = Simplify[L[2 a + d, kx].Inverse[L[a + d, kx]]];

Where it starts to really slow down

Smz[a_, d_, kx_, l_, g_] = Mz[a + d, kx, l, g].Inverse[Mz[a, kx, l, g]];

Sz[a_, d_, kx_, l_, g_] = Sl[a, kx].Smz[a, d, kx, l, g];

Impossible to compute beyond here. Even takes forever if I define every variable:

Anz[n_, a_, d_, kx_, l_, g_] := Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

Tferz[n_, a_, d_, kx_, l_, g_] := ({
   {1, 0, -Anz[n, a, d, kx, l, g][[1, 3]], -Anz[n, a, d, kx, l, g][[1,4]]},
   {0, 1, -Anz[n, a, d, kx, l, g][[2, 3]], -Anz[n, a, d, kx, l, g][[2,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[3, 3]], -Anz[n, a, d, kx, l, g][[3,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[4, 3]], -Anz[n, a, d, kx, l, g][[4,4]]}})

Auz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1,1]]}, {Anz[n, a, d, kx, l, g][[2,1]]}, {Anz[n, a, d, kx, l, g][[3,1]]}, {Anz[n, a, d, kx, l, g][[4,1]]}}

Tuz[n_, a_, d_, kx_, l_, g_] := Abs[Auz[n, a, d, kx, l, g][[1]]]^2 + Abs[Auz[n, a, d, kx, l, g][[2]]]^2

Ruz[n_, a_, d_, kx_, l_, g_] := 1 - Tuz[n, a, d, kx, l, g]

Adz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1, 2]]}, {Anz[n, a, d, kx, l, g][[2, 2]]}, {Anz[n, a, d, kx, l, g][[3, 2]]}, {Anz[n, a, d, kx, l, g][[4, 2]]}}

Tdz[n_, a_, d_, kx_, l_, g_] := Abs[Adz[n, a, d, kx, l, g][[1]]]^2 + Abs[Adz[n, a, d, kx, l, g][[2]]]^2
Rdz[n_, a_, d_, kx_, l_, g_] := 1 - Tdz[n, a, d, kx, l, g]

Tz[n_, a_, d_, kx_, l_, g_] := (Tuz[n, a, d, kx, l, g] + Tdz[n, a, d, kx, l, g])/2.
Rz[n_, a_, d_, kx_, l_, g_] := 1 - Tz[n, a, d, kx, l, g]

Attachments:

Yes, I have tried to evaluate T for suitable numbers. It takes, I would say, over 30 minutes the way the code is written now. And that's for a single value of kx.

I know about putting periods in numbers. That sped things up a lot when I was trying this for a much simpler case.

Are you saying I should hold off the evaluation of everything in the book until the end?

Oh sure. Sorry. In this workbook it's Tz[n,a,d,kx,l,g]

I was thinking about an older workbook. I want to plot Tz over kx