ENH: New method for calculating geomagnetic field vectors

[rstoneback]

A new local method for obtaining the magnetic unit vectors has been developed. The zonal unit vector is taken to be orthogonal to the magnetic field and point along a direction with no change in apex height. The accuracy of zonal vector is directly verifiable. New method includes direct accuracy/tolerance checks. Thus, most important pieces of code are now self-correcting.

To ensure smoothness in the unit vectors the numerical accuracy of the code had to be increased.

Some methods are no longer needed and where removed as this new method is better.

The old technique has been retained. It may be useful for characterizing how planar a given field line really is.

[rstoneback]

Review requested from @bharding512

[rstoneback]

Ready for review if anyone has the time, @jklenzing, @aburrell, @bharding512 .

[rstoneback]

unit_vector_grad_zonal.pdf

[rstoneback]

Here are some comparisons to the vectors provided by apexpy at 28.75 degrees latitude. Apex vectors are obtained via: apex = apexpy.Apex(date=date) apex_vecs = apex.basevectors_apex(data['glat'], data['glong'], data['alt'], coords='geo') apex_fa = apex_vecs[8] apex_mer = apex_vecs[7] apex_zon = apex_vecs[6]

The field-aligned and meridional vectors are very close here. The zonal vectors disagree, sometimes even on sign. I had to normalize the Richmond vectors, and there may be an issue there, but that wouldn't produce a sign difference. The Richmond zonal vector is based upon the gradient in longitude. My zonal vector lies along a constant apex height surface. The Richmond meridional based upon the gradient in magnetic latitude while my meridional is simply orthogonal to B and zonal. I think it is particularly interesting that both meridional vectors agree like they do, especially considering my meridional is driven by zonal and those don't agree (zonal Richmond and mine).

field_aligned_comparison.pdf meridional_comparison.pdf zonal_comparison.pdf

My calculation is different from the Richmond vectors and I have a different purpose, thus the vectors need not be the same. The Richmond vectors aren't orthonormal. Mine are. Included as a general comparison.

I should probably add the comparison code to the repo.

[rstoneback]

I've also performed a comparison neat the magnetic equator (1 MLAT). The pysatMagVect zonal vector has a nearly 0 upward component, consistent with the goal of my system; providing orthonormal basis vectors suitable for understanding equatorial ionospheric electrodynamics.

field_aligned_comparison.pdf meridional_comparison.pdf zonal_comparison.pdf

Dot product of basis vectors used in TIEGCM, apexpy and pysatMagVect, with local magnetic field vector dot_product.pdf

[bharding512]

I tested this with a couple orbits of simulated MIGHTI data I had laying around. (For context, the standard MIGHTI data product provides winds in geographic coordinates, but also uses pysatMagVect to provide the wind in magnetic coordinates).

Here is the difference between old and new:

• old = the SDC operational code as of now, presumably the "master" branch though I'm not sure.
• new = the "local_field_vector" branch

So the median change is zero, which is nice to see. There are some very large changes at low altitudes though... much larger than MIGHTI uncertainties, and large enough to make a significant difference in the dynamo calculations. I'm not sure what to make of this. Is the new version definitely more accurate, and why? Because it handles "twisted" flux tubes?

You mentioned the difference between the pysatMagVect vectors (which are orthonormal) and the Richmond vectors (which are not). Which should be used for dynamo calculations?

Given these large changes it will be important to make sure that when the SDC ingests the new magnetic coordinate file for creating ancillaries, they also update this new version of pysatMagVect to use for MIGHTI processing, so that drifts and winds are computed using the same coordinate systems. Unfortunately I had to change the MIGHTI code a little bit so it's no longer backwards compatible.

[bharding512]

I should have shown longitude on the above plots -- as you might expect, those large differences are near the SAA.

[rstoneback]

Thanks for the plots Brian. Yes, the new version is definitely more accurate, and needs to be incorporated into ICON science.

Suppose I give you the magnitude of a purely horizontal electric field directed along the magnetic equator. What is the electric field everywhere along those field lines? Using my new basis, the electric field is given directly by the specified input E magnitude dotted with my zonal basis vectors (the e_field only option from scalars_for_mapping_drifts method). There is no meridional E field component at all. Similarly, the same thing can be done for the meridional E field with an initial vertical E specification along the magnetic equator. Given an initial E at all altitudes at the magnetic equator I can specify the global E field distribution.

Let's consider the Richmond vectors for the same horizontal initial E field along equator. Describing this vector along the magnetic equator requires a combination of both the zonal and meridional vectors. The Richmond basis isn't as compact as my system. Second, since the Richmond vectors aren't orthogonal, taking the zonal and meridional amplitudes at the equator and applying them along the field lines won't produce the correct E fields along the field lines.

I may need to think about this last one a bit, but since my basis vectors are orthogonal they should conserve any curl or divergence properties found at the equator everywhere along the field line.

We should probably bring this up to the larger ICON group, but yes, I believe pysatMagVect vectors should be used for dynamo calculations.

My new vectors are better than my old vectors since: 1) New vectors respond to local changes in field line 2) Zonal vector is purely horizontal - cleanest link to understanding how horizontal and vertical forces/fields at the magnetic equator map throughout the system 3) Algorithm is self-validating to ensure more consistent numerical performance 4) I can actually demonstrate accuracy, not just precision 5) Fundamental measure for basis is unambiguous, the maximum height of field line above geodetic Earth. 6) I can also calculate the new basis with higher precision

[bharding512]

Well this all sounds great. I am a bit confused about the curl and divergence preserving properties. There may be also a distinction here between the velocity mapping and the E-field mapping.

Laundal and Richmond (2017) Section 5.2 seems especially relevant here. I'm curious if pysatMagVect produces a similar mapping as their Fig 9. Based on your comparisons with apexpy, I'd guess it will be close. Maybe the pysatMagVect velocity mapping scalars would also need to be included?

In any case I think the important point for ICON is that we avoid the problem they bring up with local magnetic coordinate systems. I think this new way (and maybe, the old way too) addresses their concerns, but I'd need to think a lot more about it:

Laundal, K. M., & Richmond, A. D. (2017). Magnetic Coordinate Systems. Space Science Reviews, 206(1–4), 27–59. https://doi.org/10.1007/s11214-016-0275-y

[rstoneback]

Well this all sounds great.

Thanks!

I am a bit confused about the curl and divergence preserving properties.

  I’m still thinking about this myself. I looked at the equations that justify claims in the area but I need to spend more time.

There may be also a distinction here between the velocity mapping and the E-field mapping.

Could be. I need to start writing this up so all the details can be discussed and/or confirmed.

Laundal and Richmond (2017) Section 5.2 seems especially relevant here. I'm curious if pysatMagVect produces a similar mapping as their Fig 9.

Yes. I’ll post it once I get back to a computer.

Based on your comparisons with apexpy, I'd guess it will be close. Maybe the pysatMagVect velocity mapping scalars would also need to be included?

I’ll double check the coordinates routine but I believe I include both the e field and velocity mappings in the big file.

In any case I think the important point for ICON is that we avoid the problem they bring up with local magnetic coordinate systems. I think this new way (and maybe, the old way too) addresses their concerns, but I'd need to think a lot more about it:

[image]https://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F8883135%2F67113155-a7d5cc00-f18d-11e9-84f7-c15264d2646f.png&data=02%7C01%7Crstoneba%40utdallas.edu%7Cfc1d5d86c90344629a2b08d753ecaabf%7C8d281d1d9c4d4bf7b16e032d15de9f6c%7C0%7C0%7C637070148252763427&sdata=jDVvcfUZ6LQmVFjTe5XnrStVOacwW9pX8ff4rAOX5%2BY%3D&reserved=0

Laundal, K. M., & Richmond, A. D. (2017). Magnetic Coordinate Systems. Space Science Reviews, 206(1–4), 27–59. https://doi.org/10.1007/s11214-016-0275-yhttps://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1007%2Fs11214-016-0275-y&data=02%7C01%7Crstoneba%40utdallas.edu%7Cfc1d5d86c90344629a2b08d753ecaabf%7C8d281d1d9c4d4bf7b16e032d15de9f6c%7C0%7C0%7C637070148252768432&sdata=uByBbvQgp9JyIfx7QTH0LB6ZvEtBBhTUzpzrhcD%2BgmY%3D&reserved=0

Thanks for the text and link. I agree we should think about the details more. Since my system is orthogonal we should avoid the concerns listed above.

I am certain about my zonal vector given I directly test it. The meridional is the correct orthogonal vector when also given the local mag field. That being said I can imagine a few questions about its properties that I don’t have a plot for yet. I’ll put that on my list.

It may also make sense for me to implement a centered difference in one area of the code. I don’t see any more large changes in outputs but this may make a small refinement to the precision.

[rstoneback]

Looking at the 2017 reference,

In a dipole 𝐝1 would point horizontally eastward.

The Richmond vectors are built on equations that start with a dipole field. If the Earth's field were a dipole, centered, and probably aligned, then it looks like my vectors and the Richmond vectors would be the same. It could be appropriate to state that pysatMagVect is a generalization of Richmond to a multipolel magnetic field.

[rstoneback]

This should work for other planets then as well, including the E-field and velocity mapping.

[rstoneback]

Drift Scalars eq_mer_drift.pdf eq_zonal_drift.pdf

E-Field Scalars eq_mer_field.pdf eq_zon_field.pdf

I have similar maps for the footpoints.

[rstoneback]

I spent a bit of time over the weekend and cast my vectors in the general covariant and contravariant forms. Rather than add to this already large pull, I'm going to approve both of these, based upon the test outputs and feedback, and create a new pull for the updated outputs that produce d and e vectors (from the apex Richmond paper notation).