code: def igrf_variation(lat, lon, alt=0., year=2005):
      
      X, Y, Z, F = calculate.igrf12syn(0, year, 1, alt, lat, lon)
      H = np.sqrt(X * X + Y * Y)
      DX, DY, DZ, DF = calculate.igrf12syn(1, year, 1, alt, lat, lon)
      DD = (60.0 * FACT * (X * DY - Y * DX)) / (H * H)
      DH = (X * DX + Y * DY) / H
      DS = (60.0 * FACT * (H * DZ - Z * DH)) / (F * F)
      DF = (H * DH + Z * DZ) / F
      return DD, DS, DH, DX, DY, DZ, DF
docs: Annual variation
           D is declination (+ve east)
           I is inclination (+ve down)
           H is horizontal intensity
           X is north component
           Y is east component
           Z is vertical component (+ve down)
           F is total intensity