version 1.14.0.0 (38.8 KB) by
Drew Compston

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

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 website (though it appears the files haven't been updated since 2007):

https://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.

Drew Compston (2020). International Geomagnetic Reference Field (IGRF) Model (https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field-igrf-model), MATLAB Central File Exchange. Retrieved .

Created with
R2007b

Compatible with any release

**Inspired by:**
IGRF Magnetic Field

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Huaixu GaoMitch HezelOkay, there is an issue with either how this code parses through the igrf2020.dat and igrf2020s.dat files. Firstly, the 2020 .dat files found at the link above are NOT in the form the matlab scripts expect. I found the correct forms here: https://ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code_error/. Although admittedly I do not like using a file coming from a directory with "error" in it's name. Then, you must include BOTH the igrf2020.dat and igrf2020s.dat, as those appear to set the date limits (2020-2025). However, the igrf2020s.dat file has an error. The number at the top, on the right should say "2025" instead of "0.0". Anyone wanting IGRF values after 2020 will have to make that edit. The igrf2020.dat file, however, does correctly have "2020.0".

Finally, to compare and make sure the file gave the correct values, I compared with the nasa.gov link above. I have no idea what in god's name start, stop, and step size do, although on first thought it might allow for the determination of the IGRF along some trajectory. Obviously to check the script solution I wanted just the IGRF at a few different points, not this interpolation or trajectory crap. Couldn't make the nasa.gov site agree. Using http://www.geomag.bgs.ac.uk/data_service/models_compass/igrf_calc.html led to an agreement, since it operates like I would expect an IGRF calculator to operate.

Very nice set of code, and works very well, but please either tell us about the changes we need to make for 2020-2025 or include the necessary parsing functionality in your code.

Mitch HezelI put both the igrf2020.dat and igrf2020s.dat files in the datfiles folder. Upon execution I receive an error on line 99: coefsvec(~h0) = thiscoefs;

"In an assignment A(:) = B, the number of elements in A and B must be the same."

What's the remedy?

zhaojin rongopen the code "igrf", and find line 154. enter a new line as line 155, and input i= abs(latitude)==90;sintheta(i)=0;, then this minor bug can be removed

zhaojin rongI found the calculated field at the south pole is not correct. the calculated field strength shows periodic variation against longitude nearby southpole.

fan zhangWayne Tsui好活

wenxin maChen LvLutz RastaetterModel download links including host name "hanna" (here: ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/) have become invalid.

Data and source code files are available instead at https://ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/.

Lutz Rastaetter - Community Coordinated Modeling Center

haha peteZhongqiu WangDavid KoronczayHi Drew,

apparently upgrading to the m36 (R2016b) version of matlab breaks this code:

>> loadigrfcoefs(now)

Error: File: /root/Documents/szakdoga/igrf/loadigrfcoefs.m Line: 55 Column: 36

The end operator must be used within an array index expression.

We found that in this matlab version, "years" is a function, and you're using a variable of the same name (in your code and also in igrfcoefs.mat).

(interestingly if I manually paste each line of the code into the matlab prompt, it seems to work, but invoking

it (see above) doesn't - maybe matlab thinks you're using a reserved word/function inappropriately).

Renaming the variable is one way to fix.

Changyong HeJonathan ParhamHey Drew,

The ftp server doesn't seem to respond anymore. Is there a new link to it? I found this link after some digging but cannot confirm, looks older though. https://ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/

Thanks

mary yangHi Drew,

Your work is excellent! I want to calculate the geodetic coordinate of the magnetic north pole .Could you tell me how can I do?

Drew CompstonSulav: If you're curious about the internals, I actually don't know why when inputting geodetic coordinates that the magnetic field is scaled like that in the last three lines. I got that bit of code from the original FORTRAN. I haven't check this myself, but you should get the same magnetic field at the same point in space regardless of whether the COORD input is 'geocentric' or 'geodetic'. Now that does NOT mean that these two statements will result in the same output:

>> igrf(now, 0, 0, 0, 'geoc')

>> igrf(now, 0, 0, 0, 'geod')

because (0, 0, 0) in geocentric coordinates is not the same point in space as (0, 0, 0) in geodetic coordinates. You'd actually have to convert between the two systems to determine what one point is in the other system. If you want to do this, you might reference, for example, https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Coordinate_system_conversion.

The program itself makes the computation of Earth's magnetic field in spherical coordinates. That is, when you provide a LATITUDE, LONGITUDE, and ALTITUDE where you want the magnetic field vector components, the program figures out where those are in a spherical coordinate system with an origin at Earth's center. Again, this is only relevant if you're interested in the internals.

Finally, I'll just say in case it is your cause for confusion: The output is a vector field, so it doesn't really make sense to convert it to ECEF, as that coordinate system defines a point in space, not a vector field. You could convert the coordinates of that point in space from geodetic latitude, longitude, and altitude that you input to ECEF using the provided function, but that wouldn't change the output magnetic vector field there.

SulavThanks a lot for this great work.

But, I am getting confused with all these reference frames.So, to obtain the output in geodetic coordinates we have to use the last three lines, is it?

And what to do, to convert the values to ECEF co-ordinates? Will the given geod2ecef work? To me, it seems like it works only for latitude, longitude values.

Minjie FanAlison MoraesIs there a way to get the geomagnetic longitude using your igrf function?

Regards

M J SchwartzSorry, I wasn't thinking... igrf.m is in Matlab code, so it is trival to look for myself. The answer to my question is that the output is in the same coordinate frame as the input, selected by COORD to be ['geodetic'] or 'geocentric'.

Drew CompstonThe output is the vector components of the magnetic field at each position given by the inputs. The output components specifically are those directed northward, eastward, and downward, respectively for Bx, By, and Bz. These are equivalent to the spherical coordinate components elevation (or -inclination), azimuth, and -radial. The requested input coordinate positions are interpreted as either geodetic or geocentric depending on the coord input, which in either case are converted to the ECEF frame.

CyberJogreat work. I am just wondering what is the reference frame of the output? Bx, By and Bz are in North, East and Down. Does that means the output is on NED frame?

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

Drew CompstonI 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.

Yaguang YangHi 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?

Drew CompstonNote 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

pedro nuñezHi,

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

Drew CompstonI'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.

Ryan LewisHi,

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!

Huy HoangHighly appreciate your nice work.

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

meltem akanPat 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.

Pat WelchDo 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.

JohnWorks pretty well so far.

Drew CompstonYou 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

XAXRXTXthanks 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 ?

thanks for your help and your time

Drew CompstonI 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.

XAXRXTXHi Drew

I have a question about Pn,m and dPn,m ,

where you get this formula to calculate them.

thanks in advance

Cesar Alberto Vazquez GarciaGood work

Drew CompstonI 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.

XAXRXTXhi

this is a very helpfull work for me,

do you have another mathematical model describing the atmospheric perturbations like MSIS or Jacchia models

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

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

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

Charles RinoA 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

MichaelI 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.

Michaelchristie harpernice 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