sort berremann qi's according to (12)

[dschick]

Hi @pyMatJ ,

I am currently checking the GTM core code I stumbled a bit over this code here:

pyGTM/GTM/GTMcore.py

Lines 461 to 477 in 1f62af2

	## sort berremann qi's according to (12)
	if any(np.abs(np.imag(qsunsorted))):
	for km in range(0,4):
	if np.imag(qsunsorted[km])>=0 :
	transmode[kt] = km
	kt = kt + 1
	else:
	reflmode[kr] = km
	kr = kr +1
	else:
	for km in range(0,4):
	if np.real(qsunsorted[km])>0 :
	transmode[kt] = km
	kt = kt + 1
	else:
	reflmode[kr] = km
	kr = kr +1

Shouldn't the if-else-case be within the for loop .
As I understood it, you need to sort according to the four cases:

1. complex number with imaginary part >= 0: transmode[0]
2. real number with >= 0: transmode[1]
3. complex number with imaginary part < 0: reflmode[0]
4. real number with < 0: reflmode[1]

Maybe you can check and tell me the idea of the current code.

Thanks

Daniel

[dschick]

okay, I just noticed, that in case of a real epsilon of a layer, the qsunsorted is also real
I guess this is why you are checking for if any(np.abs(np.imag(qsunsorted))):
However, I think the sorting in this for loop is not correct and I would be interested what is the sorting parameters for a fully real qsunsorted.

[dschick]

oh, I think I missunderstood eq. (12) which means that q_ij are either real or imag and depending on that fact you simply sort according to trans and refl and within the two refl and trans you do the sorting for s- and p-polarization in the next step.

[pyMatJ]

Hi Daniel,
You beat me to it. This is indeed why we check for if any(np.abs(np.imag(qsunsorted))):
So with respect to your first question, we could swap the for loop with the main if statement but I don't think this is very useful, and the current form is easier (at lest for me) to read.
What makes you say that the sorting is wrong?

Usually, eigenvalues are returned in pairs (of complex conjugates in case of complex eigenvalues). So you need to sort them (reflected and transmitted) but also pair them properly (e.g. transmode[0] and reflmode[0] should refer to a pair of complex conjugate). So far, with the sorting output from scipy.linalg.eig this works (after sorting with the Cp criterion if I recall correctly)

[dschick]

What makes you say that the sorting is wrong?

because I did not understood the sorting correctly.

I think I am now fine with it

[dschick]

is there a reason why the indices for in the if structures are swapped?
Don't you always want to check that _2 > _1 and if so, flip it?

pyGTM/GTM/GTMcore.py

Lines 506 to 521 in 1f62af2

	if Cp_t2>Cp_t1:
	transmode = np.flip(transmode,0) ## flip the two values
	## then calculate for reflected waves if necessary
	Cp_r1 = np.abs(self.Py[0,reflmode[1]])**2/(np.abs(self.Py[0,reflmode[1]])**2+np.abs(self.Py[1,reflmode[1]])**2)
	Cp_r2 = np.abs(self.Py[0,reflmode[0]])**2/(np.abs(self.Py[0,reflmode[0]])**2+np.abs(self.Py[1,reflmode[0]])**2)
	if Cp_r1>Cp_r2:
	reflmode = np.flip(reflmode,0) ## flip the two values
	else: ### No birefringence, use the Electric field s-pol/p-pol
	Cp_te1 = np.abs(psiunsorted[0,transmode[1]])**2/(np.abs(psiunsorted[0,transmode[1]])**2+np.abs(psiunsorted[2,transmode[1]])**2)
	Cp_te2 = np.abs(psiunsorted[0,transmode[0]])**2/(np.abs(psiunsorted[0,transmode[0]])**2+np.abs(psiunsorted[2,transmode[0]])**2)
	if Cp_te1>Cp_te2:
	transmode = np.flip(transmode,0) ## flip the two values
	Cp_re1 = np.abs(psiunsorted[0,reflmode[1]])**2/(np.abs(psiunsorted[0,reflmode[1]])**2+np.abs(psiunsorted[2,reflmode[1]])**2)
	Cp_re2 = np.abs(psiunsorted[0,reflmode[0]])**2/(np.abs(psiunsorted[0,reflmode[0]])**2+np.abs(psiunsorted[2,reflmode[0]])**2)
	if Cp_re1>Cp_re2:

[pyMatJ]

Ah, I see. This is a direct translation of the matlab code, but there are so many indices that you end up lost...
The idea is directly from https://doi.org/10.1364/AO.27.001334 just below eq. (8):

The eigenvectors of the transmitting modes with the large and small Cp occupy the first and third columns, respectively. Similarly, the eigenvectors of the reflecting modes with the large and small Cp will occupy the second and fourth columns, respectively.

So you are right, but in order to stick as much as possible to the matlab code, I kept the same ordering also in indexing transmode and reflmode. As you can see in lines 501-502:

pyGTM/GTM/GTMcore.py

Lines 501 to 502 in 1f62af2

	Cp_t1 = np.abs(self.Py[0,transmode[0]])**2/(np.abs(self.Py[0,transmode[0]])**2+np.abs(self.Py[1,transmode[0]])**2)
	Cp_t2 = np.abs(self.Py[0,transmode[1]])**2/(np.abs(self.Py[0,transmode[1]])**2+np.abs(self.Py[1,transmode[1]])**2)

Cp_t1 deals with transmode[0] and Cp_t2 with transmode[1]
but then

pyGTM/GTM/GTMcore.py

Lines 509 to 510 in 1f62af2

	Cp_r1 = np.abs(self.Py[0,reflmode[1]])**2/(np.abs(self.Py[0,reflmode[1]])**2+np.abs(self.Py[1,reflmode[1]])**2)
	Cp_r2 = np.abs(self.Py[0,reflmode[0]])**2/(np.abs(self.Py[0,reflmode[0]])**2+np.abs(self.Py[1,reflmode[0]])**2)

Cp_r1 with reflmode[1] and Cp_r2 with reflmode[0].

I will correct this. (EDIT: I should probably wait until I can test a bit more and accept your PR so this goes along)