Code covered by the BSD License  

Highlights from
International Geomagnetic Reference Field (IGRF) Model

4.66667

4.7 | 6 ratings Rate this file 75 Downloads (last 30 days) File Size: 39 KB File ID: #34388
image thumbnail

International Geomagnetic Reference Field (IGRF) Model

by

 

30 Dec 2011 (Updated )

Compute the Earth's magnetic field at points in space according to the IGRF model.

| Watch this File

File Information
Description

The International Geomagnetic Reference Field (IGRF) is an internationally agreed upon mathematical model of the Earth's magnetic field. This program is a conversion of the FORTRAN subroutines that make the calculation into MATLAB. It does not use a compiled FORTRAN mex file, which probably makes it slower but at the advantage of being easier to use (as no compilation is necessary). In fact, my motivation in writing the program was to provide an IGRF implementation in MATLAB with minimal "fuss." Another motivation was a vectorized IGRF function, which this function is (with a separate routine adapted directly from the FORTRAN code that is faster for scalars implemented as well).

The following files are provided:

-igrf.m: Computes Earth's magnetic field at a point(s).
-igrfline.m: Gives the coordinates along a magnetic field line starting at a given point.
-getigrfcoefs.m: Extracts coefficients from the .dat files provided on the IGRF website and saves them to a .mat file.
-igrfcoefs.mat: IGRF coefficients of the 11th IGRF generation (most recent as of 2010).
-loadigrfcoefs.m: Loads the proper IGRF coefficients at a given time (making the necessary interpolation).
-*grf*.dat: 11th generation IGRF coefficient data files.
-plotbline: Plots a magnetic field line.
-plotbearth: Plots a number of magnetic field lines.
-geod2ecef: Converts geodetic coordinates to ECEF coordinates.
-ecef2geod: Converts ECEF coordinates to geodetic.

The only prerequisite to running either the function IGRF or the function IGRFLINE is to put the file igrfcoefs.mat in the MATLAB search path. The program is designed to be scalable with time: As new IGRF generations are released, simply replace the old .dat files with their newer versions in a subfolder called 'datfiles' within the same directory that the function getigrfcoefs.m is located and run getigrfcoefs, and then replace the file it generates (igrfcoefs.mat) with the old .mat file. Updates happen every five years, with the last update occurring in 2010. New .dat files will hopefully continue to be uploaded to the following ftp:
ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/

Finally, I have included two example scripts showing how the function IGRFLINE works: plotbline.m and plotbearth.m. These scripts both utilize the Mapping Toolbox to plot globes upon which magnetic field lines are plotted, but if the user does not have that package, a crude globe with just latitude and longitude lines is shown.

I've made some cursory comparisons with the online IGRF calculator at http://ccmc.gsfc.nasa.gov/modelweb/models/igrf_vitmo.php and found this function to be accurate to within 1 nT. I'm not sure why there is a discrepancy between the two, but my guess is round-off error.

Acknowledgements

Igrf Magnetic Field inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.5 (R2007b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (10)
25 Apr 2014 Cesar Alberto Vazquez Garcia

Good work

03 Apr 2014 Drew Compston

I do have an IRI and MSIS code that works by querying an online interface. So it's a lot slower than this but gets the job done. I don't foresee converting the FORTRAN source code for models like the IRI or Tsyganenko any time soon (if ever). The IGRF model is much simpler than those as it is really just one (vector) equation.

03 Apr 2014 XAXRXTX

hi
this is a very helpfull work for me,
do you have another mathematical model describing the atmospheric perturbations like MSIS or Jacchia models

12 Mar 2014 cao  
09 Aug 2013 Paul

Great work. Any plans do something similar with the Tsyganenko models?

http://ccmc.gsfc.nasa.gov/models/modelinfo.php?model=Tsyganenko%20Model

06 Jul 2012 Anna Kelbert

It's amazing how fast this thing is! Takes me much longer to compute the fields on a grid with a homemade spherical harmonic toolbox.

13 May 2012 Charles Rino

A significant improvement over my own direct translation of igrf11syn.
http://www.mathworks.com/matlabcentral/fileexchange/28874-igrf-magnetic-field
Easier to use, with cleaner interface with spherical harmonic coefficient files. Field tracing utilities are also very useful.

Users may find GPS transformation helpful:
http://www.mathworks.com/matlabcentral/fileexchange/28813-gps-coordinate-transformations

25 Jan 2012 Michael

I have noticed that results are inconsistant depending on the size if the inputs:

[Bx,By,Bz]=igrf(734305,[-10],0,0,'geodetic');
atand(Bz./sqrt(Bx.^2 + By.^2))

Gives an inclination of -48.9117, while:

[Bx,By,Bz]=igrf(734305,[-10,-5],0,0,'geodetic');
atand(Bz./sqrt(Bx.^2 + By.^2))

Gives an inclination of -48.6375, and other sized vectors produce slightly different results.

Otherwise it is well done.

25 Jan 2012 Michael  
31 Dec 2011 christie harper

nice job. few minor issues.

isleapyr function not included
here is a simple one
function [is_leap] = leap_year(year_in)
is_leap = (~mod(year_in,4) && mod(year_in,100)) || ~mod(year_in,400);
end

and the
referenced version of ecef2lla does not return a matrix but vectors
mods
change function signature to
function [lla] = ecef2lla(ecx)
change returns to
lla = zeros(length(glat),3);
lla(:,1) = glat;
lla(:,2) = glon;
lla(:,3) = eht;

also your spin logic is broken at the end
change while true to while spin

Updates
01 Jan 2012

Fixed the issues detailed in Christie Harper's comment (including some additional ones relating to the differences between MATLAB's built-in ecef2lla and the free file exchange ecef2lla) and made some minor changes to speed up the functions slightly.

25 Jan 2012

Fixed bug described in Michael's comment. For those interested, it resulted from mistakenly calculating dP using only the last latitude input rather than the vector.

25 Apr 2012

Provided my own ECEF to geodetic conversion routine.

07 Aug 2012

Added capability to plot globe without mapping toolbox and minor bug fix.

02 Nov 2012

Adjusted initial comment block in IGRF with correct definition of inclination using atan2 (rather than atan).

Contact us