function  [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
% function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf);
% 
% compute magnetic deviation
% input:
%      flat    latitude (degree)
%      flon    longitude (degree)
%      elevkm  elevation above mean geoid (km)
%      nmax    number of harmonics
%      igrf    (need file, default 'IGRF95')
%
% output:
%     dev     mag. deviation (degree)
% 
% 
% based of FORTRAN ROUTINE GEOMAG.FOR
% more info under  http://fdd.gsfc.nasa.gov/IGRF.html
%
% M. Visbeck, LDEO FEB 2000
%

if nargin<5
 igrf='IGRF95';
 igrf='IGRF00';
end
eval(['gh=',igrf,';'])

if nargin<4
 nmax=10;
end

if nargin<3
 elevkm=0;
end



[x,y,z] = shval3(flat,flon,elevkm,nmax,gh);
[d,i,h,f]=dihf (x,y,z);

dev=d*180/pi;

if nargout<1
 disp([' model ',igrf,'  harmonics: ',int2str(nmax)])
 disp([' lat: ',num2str(flat)])
 disp([' lon: ',num2str(flon)])
 disp([' mag dev [deg]: ',num2str(d*180/pi)])
end

%=====================================================================

function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
% function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh)
% ================================================================
%
%       version 1.01
%
%       calculates field components from spherical harmonic (sh)
%       models.
%
%       input:
%           flat  - north latitude, in degrees
%           flon  - east longitude, in degrees
%           elevkm  - elevation above mean sea level 
%           nmax  - maximum degree and order of coefficients
%           gh    - schmidt quasi-normal internal spherical
%                   harmonic coefficients
%           iext  - external coefficients flag (= 0 if none)
%           ext   - the three 1st-degree external coefficients
%                   (not used if iext = 0)
%
%       output:
%           x     -  northward component
%           y     -  eastward component
%           z     -  vertically-downward component
%
%       based on subroutine 'igrf' by d. r. barraclough and
%       s. r. c. malin, report no. 71/1, institute of geological
%       sciences, u.k.
%
%       norman w. peddie, u.s. geological survey, mail stop 964,
%       federal center, box 25046, denver, colorado 80225
%
% ================================================================
%       the required sizes of the arrays used in this subroutine
%       depend on the value of nmax.  the minimum dimensions
%       needed are indicated in the table below.  (note that this
%       version is dimensioned for nmax of 14 or less).
%
% ================================================================
dtr = pi/180;
r=elevkm;
erad=6371.2;
a2=40680925.;
b2=40408588.;



slat = sin(flat*dtr);
aa = min(89.999,max(-89.999,flat));
clat = cos(aa*dtr);
sl(1) = sin(flon*dtr);
cl(1) = cos(flon*dtr);

x=0.0;
y=0.0;
z=0.0;
sd = 0.0;
cd = 1.0;
n=0;
l=1;
m=1;

npq = (nmax*(nmax+3))/2;
aa = a2*clat*clat;
bb = b2*slat*slat;
cc = aa+bb;
dd = sqrt(cc);
r=sqrt(elevkm*(elevkm+2.0*dd)+(a2*aa+b2*bb)/cc);
cd = (elevkm+dd)/r;
sd = (a2-b2)/dd*slat*clat/r;
aa = slat;
slat = slat*cd-clat*sd;
clat = clat*cd+aa*sd;
ratio = erad/r;
aa = sqrt(3.0);
p(1) = 2.0*slat;
p(2) = 2.0*clat;
p(3) = 4.5*slat*slat-1.5;
p(4) = 3.0*aa*clat*slat;
q(1) = -clat;
q(2) = slat;
q(3) = -3.0*clat*slat;
q(4) = aa*(slat*slat-clat*clat);


for k = 1: npq

if (n<m)
    m=0;
    n=n+1;
    rr = ratio^(n+2);
    fn=n;
end;

fm=m;

if (k>=5)
    if (m==n)
         aa = sqrt(1.0-.5/fm);
         j=k-n-1;
         p(k) = (1.0+1.0/fm)*aa*clat*p(j);
         q(k) = aa*(clat*q(j)+slat/fm*p(j));
         sl(m) = sl(m-1)*cl(1)+cl(m-1)*sl(1);
         cl(m) = cl(m-1)*cl(1)-sl(m-1)*sl(1);
    else 
         aa = sqrt(fn*fn-fm*fm);
         bb = sqrt((fn-1.0)^2-fm*fm)/aa;
         cc = (2.0*fn-1.0)/aa;
         i=k-n;
         j=k-2*n+1;
         p(k) = (fn+1.0)*(cc*slat/fn*p(i)-bb/(fn-1.0)*p(j));
         q(k) = cc*(slat*q(i)-clat/fn*p(i))-bb*q(j);
    end;
end;

aa = rr*gh(l);

if (m==0)
   x=x+aa*q(k);
   z=z-aa*p(k);
   l=l+1;
else 
   bb = rr*gh(l+1);
   cc = aa*cl(m)+bb*sl(m);
   x=x+cc*q(k);
   z=z-cc*p(k);
   if (clat>0.0)
     y=y+(aa*sl(m)-bb*cl(m))*fm*p(k)/((fn+1.0)*clat);
   else 
     y=y+(aa*sl(m)-bb*cl(m))*q(k)*slat;
   end;
   l=l+2;
end;

m=m+1;

end;

aa=x;
x=x*cd+z*sd;
z=z*cd-aa*sd;
          
return;

%==================================================================
function [d,i,h,f] = dihf(x,y,z)
% function [d,i,h,f] = dihf(x,y,z)
% ===============================================================
%
%       version 1.01
%
%       computes the geomagnetic elements d, i, h, and f from
%       x, y, and z.
%
%       input:
%           x   - northward component
%           y   - eastward component
%           z   - vertically-downward component
%
%       output:
%           d   - declination
%           i   - inclination
%           h   - horizontal intensity
%           f   - total intensity
%
%       a. zunde
%       usgs, ms 964, box 25046 federal center, denver, co  80225
%
% ===============================================================
sn=0.0001;

% ---------------------------------------------------------------
%       if d and i cannot be determined, set equal to NaN
% ---------------------------------------------------------------
h=sqrt(x*x+y*y);
f=sqrt(x*x+y*y+z*z);
if (f<sn)
  d=NaN;
  i=NaN;
else 
  i=atan2(z,h);
  if (h<sn)
    d=NaN;
  else 
    hpx = h+x;
    if (hpx<200)
      d=pi;
    else 
      d=2.0*atan2(y,hpx);
    end;
  end;
end;
return;