Contravariant and covariant basis vectors

[rstoneback]

Implemented the two mathematical forms of complete basis vectors (contravariant and covariant), in addition to the unit vectors already returned. This puts the system on firm mathematical ground and should make it easy to translate any system that currently uses vectors per Richmond over to pysatMagVect.

I also implemented a centered difference when calculating the change in apex height as part of the unit vector determination algorithm. A zonal vector with no apex height change is a fundamental definition of this basis, thus worth spending some cycles on for a more precise and accurate derivative algorithm.

Items to come

• Unit test plots for the E and D vectors (done)

• There are three numerical paths to the meridional vector. I'll implement the next hardest one (easiest in the code) to provide an independent check on the meridional unit vector (this is done. May be more appropriate to work on zonal magnitude instead since apex height calculations are very accurate)

• Test that a proper combination of E and D vectors scales electric fields and ion drifts through comparison to existing routines. These older routines may be removed, or shifted to a test support role.

• Think about and make plots on what these vectors mean for scaling field-aligned ion drifts, if that makes any sense (feedback will be requested from @aburrell)

• Incorporate isolated apexpy and pysatMagVect comparison code into repo (done)

• Future pull - implement a new mlat derivation based upon the meridional vector

Pinging @jklenzing and @bharding512 per the ICON science connection

[rstoneback]

The impact to field aligned drifts may be felt more for flux than anything. I've been a bit negligent (well, more than a bit) on the Mars work @aburrell but looks like the stuff here works in that area as well.

[bharding512]

Are the new E and D vectors non-orthogonal? If so, I think I understand and agree this is a nice way to connect with the Richmond apex vectors. I liked the interpretation you had in another thread, that these could be considered as the multi-pole generalization of those vectors.

Regarding ICON, I expect this will be a topic of conversation. I'm beginning to think we should use a non-orthogonal system, for the following reason:

1. Imagine there is a purely zonal (i.e., horizontal) electric field at the mag equator. This maps down the field line to an electric field at 150 km that let's call e1.

2. Now do the same for a purely meridional (i.e, vertical) electric field at the mag equator. Call this e2 at 150km.

3. Unless the magnetic field topology is exactly ideal (e.g., a dipole), then e1 is not perpendicular to e2.

4. Conversely, a purely zonal wind at 120km will create an electric field f1 that IVM would measure at the equator. A purely meridional wind would create f2. There is no guarantee f1 is perp to f2.

You have all thought about this much longer than I have, so let me know if this was already obvious or if I'm missing something important.

[rstoneback]

Are the new E and D vectors non-orthogonal? If so, I think I understand and agree this is a nice way to connect with the Richmond apex vectors. I liked the interpretation you had in another thread, that these could be considered as the multi-pole generalization of those vectors.

I will produce plots to demonstrate this, but my E and D should be orthogonal by construction. Orthonormal even. Thanks for the feedback on calling it a multi-pole generalization.

Regarding ICON, I expect this will be a topic of conversation. I'm beginning to think we should use a non-orthogonal system, for the following reason:

1. Imagine there is a purely zonal (i.e., horizontal) electric field at the mag equator. This maps down the field line to an electric field at 150 km that let's call e1.
2. Now do the same for a purely meridional (i.e, vertical) electric field at the mag equator. Call this e2 at 150km.
3. Unless the magnetic field topology is exactly ideal (e.g., a dipole), then e1 is not perpendicular to e2.
4. Conversely, a purely zonal wind at 120km will create an electric field f1 that IVM would measure at the equator. A purely meridional wind would create f2. There is no guarantee f1 is perp to f2.

You have all thought about this much longer than I have, so let me know if this was already obvious or if I'm missing something important.

Your description is correct for the Richmond vectors since they are non-orthogonal. His E and D vectors even pointed in different directions for the Earth. My E and D vectors should point the same directions as my unit vectors. E and D only differ in magnitude. I say should, but by code construction I know the d vectors are scaled unit vectors, and the e vectors follow from the d. Since unit vectors are orthogonal, d is orthogonal, thus so is e as e vectors point along d.

The following should be true, per my understanding. Soon to be demonstrated. My only real worry is a small calculation typo at the moment but the proof will be in the pudding.

1. Imagine there is a purely zonal (i.e., horizontal) electric field at the mag equator. This maps down the field line to an electric field at 150 km that let's call e1.
2. Now do the same for a purely meridional (i.e, vertical) electric field at the mag equator. Call this e2 at 150km.
3. e1 is perpendicular to e2 at 150 km and all other locations along the field line.
4. Conversely, a purely zonal wind at 120km will create an electric field f1 that IVM would measure at the equator. A purely meridional wind would create f2. f1 is perp to f2 at all locations along the field line.

[bharding512]

Ok thanks that helps clear things up. So you and I have written two different versions of steps 1-4. If I understand it, one is correct for the Richmond vectors, one is correct for the pysatMagVect vectors. Which is correct for nature? My inclination is to say that (insofar as IGRF is a good representation of the real B field) e1 and e2 are not perpendicular in reality.

[rstoneback]

Ok thanks that helps clear things up. So you and I have written two different versions of steps 1-4. If I understand it, one is correct for the Richmond vectors, one is correct for the pysatMagVect vectors. Which is correct for nature? My inclination is to say that (insofar as IGRF is a good representation of the real B field) e1 and e2 are not perpendicular in reality.

The 2017 paper directly states the Richmond vectors are only correct for a dipole field and a spherical Earth. Once I get my plots, and assuming they confirm the expected math properties, then my vectors are consistent with reality (observations), presuming IGRF and geodetic model.

My claim can be independently verified with a ionosphere model that incorporates IGRF, assumes equipotential field lines, and fully solves the EM equations in 3D space without resorting to a Richmond like basis. Is that SAMI? If that model were driven by purely horizontal or vertical electric field at the equator then we could test if the electric field distribution was the same as my basis or not.

[bharding512]

Ok thanks that helps clear things up. So you and I have written two different versions of steps 1-4. If I understand it, one is correct for the Richmond vectors, one is correct for the pysatMagVect vectors. Which is correct for nature? My inclination is to say that (insofar as IGRF is a good representation of the real B field) e1 and e2 are not perpendicular in reality.

The 2017 paper directly states the Richmond vectors are only correct for a dipole field and a spherical Earth. Once I get my plots, and assuming they confirm the expected math properties, then my vectors are consistent with reality (observations), presuming IGRF and geodetic model.

What I don't understand is why IGRF would produce an orthogonal electric field mapping, while an even simpler quasi-dipole model produces a non-orthogonal mapping (according to my understanding of Richmond's coordinates)

My claim can be independently verified with a ionosphere model that incorporates IGRF, assumes equipotential field lines, and fully solves the EM equations in 3D space without resorting to a Richmond like basis. Is that SAMI? If that model were driven by purely horizontal or vertical electric field at the equator then we could test if the electric field distribution was the same as my basis or not.

That would certainly put this question to rest. Impose a horizontal drift at the equator, see what the drift is at 150 km. Impose a vertical drift is at the equator, see what the drift is at 150 km. Are those drifts orthogonal? SAMI assumes equipotential field lines, so that won't work, but I know Dave Hysell has the type of model you describe. What I don't know is whether it has a realistic enough magnetic field geometry to see the effect we're after.

But maybe that's all too complicated. I keep coming back to the simple circuit analogy. At the risk of stating the obvious, an electric field means there is a potential drop between two adjacent field lines. This potential drop is constant no matter how those field lines twist or converge (as long as the specific conductivity is larger than Pedersen, of course). The field line that is adjacent to me in the zonal direction and the field line adjacent in the meridional direction are in orthogonal directions from me at first, but that's not necessarily true if I move some distance along my field line.

What I don't have a good feel for is how much this matters. Perhaps it's negligible, but some of the plots in the Laundal and Richmond paper show alarmingly large non-orthogonalities if I understand it correctly.

[rstoneback]

Ok thanks that helps clear things up. So you and I have written two different versions of steps 1-4. If I understand it, one is correct for the Richmond vectors, one is correct for the pysatMagVect vectors. Which is correct for nature? My inclination is to say that (insofar as IGRF is a good representation of the real B field) e1 and e2 are not perpendicular in reality.

The 2017 paper directly states the Richmond vectors are only correct for a dipole field and a spherical Earth. Once I get my plots, and assuming they confirm the expected math properties, then my vectors are consistent with reality (observations), presuming IGRF and geodetic model.

What I don't understand is why IGRF would produce an orthogonal electric field mapping, while an even simpler quasi-dipole model produces a non-orthogonal mapping (according to my understanding of Richmond's coordinates)

Richmond does produce orthogonal coordinates for a dipole and spherical Earth. When this same model gets translated to the Earth the dipole model has to be stretched, compressed, etc. to fit the geomagnetic field. Thats the part that distorts the vector orthogonality.

The non-orthogonality of the vectors isn't a property of the magnetic field. EM fields are linear. I believe that implies that an orthogonal basis is fundamentally there.

My claim can be independently verified with a ionosphere model that incorporates IGRF, assumes equipotential field lines, and fully solves the EM equations in 3D space without resorting to a Richmond like basis. Is that SAMI? If that model were driven by purely horizontal or vertical electric field at the equator then we could test if the electric field distribution was the same as my basis or not.

That would certainly put this question to rest. Impose a horizontal drift at the equator, see what the drift is at 150 km. Impose a vertical drift is at the equator, see what the drift is at 150 km. Are those drifts orthogonal? SAMI assumes equipotential field lines, so that won't work, but I know Dave Hysell has the type of model you describe. What I don't know is whether it has a realistic enough magnetic field geometry to see the effect we're after.

But maybe that's all too complicated. I keep coming back to the simple circuit analogy. At the risk of stating the obvious, an electric field means there is a potential drop between two adjacent field lines. This potential drop is constant no matter how those field lines twist or converge (as long as the specific conductivity is larger than Pedersen, of course). The field line that is adjacent to me in the zonal direction and the field line adjacent in the meridional direction are in orthogonal directions from me at first, but that's not necessarily true if I move some distance along my field line.

It is true if the zonal and meridional vectors vary appropriately. My zonal vector is directly determined by asking the question, is this the field line that is horizontal at the magnetic equator? My zonal vector is along no apex height change, and the mag vector is along zero apex height change, thus the third vector that goes along with this system is always along the maximum change in apex height. Pointing to the same field line that is purely vertical wrt geodetic Earth at magnetic equator.

What I don't have a good feel for is how much this matters. Perhaps it's negligible, but some of the plots in the Laundal and Richmond paper show alarmingly large non-orthogonalities if I understand it correctly.

I don't have a calculation for that yet. I've confirmed that my E and D vectors check out (low res on laptop). My zonal vector has even gotten better with the centered difference.

That covers two of the items to come in summary above. On to the third, checking scaling of electric fields and ion drifts. I can do two checks now, one against my existing electric and ion drift mapping, the other against apexpy.

[rstoneback]

Basic comparison with apexpy at the magnetic equator.

comparison_field_aligned.pdf comparison_meridional.pdf comparison_zonal.pdf

[rstoneback]

Some updates:

I have two paths to the d and e vectors. They are consistent within about 99%, more in areas. Plot below is the 'worst'. One path uses the expansion of magnetic field along zonal directions at apex, the other is the same but for meridional. d_diff_zon_east_norm.pdf I should make one just of magnitude.

I've tried a few things to reduce this error. I think it is driven by gradients in the meridional vector, shown here in uncertainty in the obtained apex position when - taking a linear step - and taking multiple smaller steps along the vector. In both case the gradient I've chosen defines one vector, the mag field the other, then the third is the cross product of the two. meridional_step_diff_apex_height_z.pdf

Here are uncertainties for this process along meridional and zonal directions: meridional_step_diff_v_longitude.pdf zonal_step_diff_v_longitude.pdf

Zonal looks better. I've taken the d and e vectors generated from the zonal apex locations as the primary. Still thinking this over some as the quantity that goes into the d,e vectors is the apex height gradient, which looks pretty good. The zonal step error doesn't have any spatial variations that coincide with the d difference.

Additional check on d, e vectors. I haven't written up the general routine for mapping on these, but I can cheat for the mag equator. Results from my heritage technique, which is effectively the same but different code. eq_mer_drift.pdf From my calculated apex-height gradient. Essentially the d-vector. unit_vector_grad_meridional.pdf

While writing this up I thought of a couple things I can try to reduce the d-vector difference: 1) I currently use a zonal vector from centered step to get the linear step direction for the meridional. It may make sense for me to switch to a euler step here as at keeps me on the same apex height surface 2) There may be an additional term to include based upon the local magnitude gradient in the magnetic field. I tried reducing step sizes which didn't seem to help. If 2 was true, seems like a reduced step would reduce errors. A plot of magnetic field gradients would probably be helpful.

The places with a d difference are equivalent to saying the changes in distance I calculate aren't consistent with the change in the magnetic field magnitude.

My results aren't bad but I can't help but feel the loss of precision. 6-7+ digits on the unit vector directions, unit_vector_tol_zonal.pdf

[rstoneback]

I tried variants of option 1 above. Made things worse.

[rstoneback]

1) I've let this pull get a bit stale 2) The main thrust of the pull, an updated unit vector system plus D, E vectors is in there and complete. I have testing plots that characterize performance plus I've already incorporated some lessons. 3) Accuracy and sensitivity of the vectors are better than ever. Thus, I'm going to merge and close this pull.

4) Remaining issue, incorporating D, E vectors into the scaling methods will be handled by the next pull request. The current code already effectively uses these vectors, however it uses unique code. Needs to be replaced with the tested methods.