code: def b_vdm(B, lat):
      
      # changed radius of the earth from 3.367e6 3/12/2010
      fact = ((6.371e6)**3) * 1e7
      colat = np.radians(90. - lat)
      return fact * B / (np.sqrt(1 + 3 * (np.cos(colat)**2)))
docs: Converts a magnetic field value (input in units of tesla) to a virtual
      dipole moment (VDM) or a virtual axial dipole moment (VADM); output
      in units of Am^2)
  
      Parameters
      ----------
      B: local magnetic field strength in tesla
      lat: latitude of site in degrees
  
      Returns
      ----------
      V(A)DM in units of Am^2
  
      Examples
      --------
      >>> pmag.b_vdm(33e-6,22)*1e-21
  
      71.58815974511788