Code covered by the BSD License

### Highlights from International Geomagnetic Reference Field (IGRF) Model

4.66667
4.7 | 12 ratings Rate this file 55 Downloads (last 30 days) File Size: 39.3 KB File ID: #34388 Version: 1.13

# International Geomagnetic Reference Field (IGRF) Model

### Drew Compston (view profile)

30 Dec 2011 (Updated )

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

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 12th IGRF generation (most recent as of 2015).
-loadigrfcoefs.m: Loads the proper IGRF coefficients at a given time (making the necessary interpolation).
-*grf*.dat: 12th 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 2015. 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)
MATLAB Search Path
```/
/datfiles```
14 Jul 2016 nikhil raykar

### nikhil raykar (view profile)

11 Mar 2016 Peter Gorham

### Peter Gorham (view profile)

05 Mar 2016 Yaguang Yang

### Yaguang Yang (view profile)

Thank you very much! I wrote some quick and dirty code yesterday for my application. Your info is very helpful.

Comment only
05 Mar 2016 Drew Compston

### Drew Compston (view profile)

I wrote some codes to do more or less what you're asking I think for a class, but because it was for a class it isn't really at a level worth posting on here. I think there are a lot of orbital mechanics toolboxes available you could use. A quick search on the exchange yields an Orbital Mechanics Library (File ID: #13439) that probably is similar to my codes (also done for a class) but likely more organized. Take a look at that and see if it will work for you.

Comment only
29 Feb 2016 Yaguang Yang

### Yaguang Yang (view profile)

Hi Drew,

Very nice work! I have an application that prevents me from directly using your code. I would like to see if you have a solution for my problem. Assuming a spacecraft orbit is circular, given the ascending node (assuming "now" the spacecraft pass the node), spacecraft altitude, and inclination, one should be able to calculate the coordinate of the spacecraft at any time after time "now" then use the coordinate and your code to find the magnetic field vector at the position of the spacecraft for any time after the time "now". Do you have a code or the formulas to compute the spacecraft coordinate given ascending node, altitude, and inclination so that one can use your code?

Comment only
24 Feb 2016 Drew Compston

### Drew Compston (view profile)

Note in general when having a question about a program, it always helps to say exactly what you input. For example:
>> igrf(now, 0, 0, 0)
ans =
27546 -2545.8 -15893
So yes, the units are definitely nT. As a guess, I wondered if you were inputting a radius in m. Apparently you are, as:
>> igrf(now, 0, 0, 6371e3)
ans =
2.934e-05 -4.7516e-06 2.9438e-06
But if you read the description, note that the fourth input is height from the Earth's surface in km when coord is 'geodetic' (the default). In geocentric, the fourth input is radius from the center of the Earth, but in km. So:
>> igrf(now, 0, 0, 6371, 'geoc')
ans =
27647 -2551.9 -15986

Comment only
24 Feb 2016 pedro nuñez

### pedro nuñez (view profile)

Hi,
Im not sure if the units given by the function igrf for the magnetic field are in Teslas or in nanoteslas, because the description says that is in nanoteslas but then the result given I think is too small to be in nanoteslas (something to the order of 10^-5). Could you help me with that?
Thanks
Thanks

Comment only
11 Feb 2016 Drew Compston

### Drew Compston (view profile)

I'm not sure why there was no igrf1925.dat. It looks like the program probably broke if you asked for a value between 1925 and 1930. I've uploaded a new version that should be more robust in that case (not assuming 5 year gaps between IGRF coefficient values) and included the igrf1925 coefficients.

Comment only
10 Feb 2016 Ryan Lewis

### Ryan Lewis (view profile)

Hi,
I am using this code for part of my final year dissertation project. I was wondering why there is no igrf1925.dat file? I have found the data and made my own .dat file however the code doesn't seem to recognize it. Any help would be much appreciated. Apart from that, working great and it's saved a lot of work for me so thanks!

13 Jan 2016 Huy Hoang

### Huy Hoang (view profile)

08 Oct 2015 Pat Welch

### Pat Welch (view profile)

Meltem Akan, please email me, pat at mousebrains.com, for an updated version of Drew's code which reads in the IGRF12 Excel file.

Comment only
10 Apr 2015 meltem akan

### meltem akan (view profile)

Pat Welch hello. I am working on the same project. I don't know how to contact with you , but I will be happy if we can work together.

Comment only
20 Mar 2015 Pat Welch

### Pat Welch (view profile)

Do you have an update to use the IGRF12 database for 2015-2020 yet? If not I'll take care of it, but please let me know.

Comment only
02 Dec 2014 John

### John (view profile)

Works pretty well so far.

09 Oct 2014 Drew Compston

### Drew Compston (view profile)

You should be able to get IGRF values at the website listed in the description above: http://ccmc.gsfc.nasa.gov/modelweb/models/igrf_vitmo.php

Comment only
06 Oct 2014 XAXRXTX

### XAXRXTX (view profile)

thanks for your response, in my case I used the formulas figured in wertz text book and they are the same as those figured in the paper that you indicated, and also I got the same results as yours.

now, if it's possible, could you show me how can I compare this model with the online model ?

Comment only
06 Oct 2014 Drew Compston

### Drew Compston (view profile)

I got the formulas for P(n,m) and dP(n,m) by looking at the direct Fortran translation of the IGRF done by Charles Rino (#28874, listed as inspiration for this file) and figuring out what that program was doing. The formulas in the paper here (http://hanspeterschaub.info/Papers/UnderGradStudents/MagneticField.pdf) are different, and I am not sure why.

Comment only
04 Oct 2014 XAXRXTX

### XAXRXTX (view profile)

Hi Drew

I have a question about Pn,m and dPn,m ,
where you get this formula to calculate them.

25 Apr 2014 Cesar Alberto Vazquez Garcia

### Cesar Alberto Vazquez Garcia (view profile)

Good work

Comment only
03 Apr 2014 Drew Compston

### Drew Compston (view profile)

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.

Comment only
03 Apr 2014 XAXRXTX

### XAXRXTX (view profile)

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

Comment only
12 Mar 2014 cao

09 Aug 2013 Paul

### Paul (view profile)

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

### Anna Kelbert (view profile)

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

### Charles Rino (view profile)

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

### Michael (view profile)

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.

Comment only
25 Jan 2012 Michael

### Michael (view profile)

31 Dec 2011 christie harper

### christie harper (view profile)

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

01 Jan 2012 1.3

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 1.4

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 1.5

Provided my own ECEF to geodetic conversion routine.

07 Aug 2012 1.6

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

02 Nov 2012 1.7

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

16 Oct 2014 1.8

Apparently, the earlier bug pointed out by Michael on 25 Jan 2012 resurfaced. Presumably the fix I had implemented before got undone, so this update includes that fix again.

29 Oct 2014 1.9

20 Mar 2015 1.10

Includes 2015 coefficients (12th generation IGRF).

20 Mar 2015 1.11

.

20 Apr 2015 1.12

Updated description text to reflect my most recent update: The newest 2015 coefficients are now included.

11 Feb 2016 1.13

Added igrf1925 coefficients and made loadigrfcoefs.m more robust to different than 5 year gaps between coefficient files.