docs: def scalars_for_mapping_ion_drifts(glats, glons, alts, dates, step_size=None, 
                                     max_steps=None, e_field_scaling_only=False):
      """
      Calculates scalars for translating ion motions at position
      glat, glon, and alt, for date, to the footpoints of the field line
      as well as at the magnetic equator.
      
      All inputs are assumed to be 1D arrays.
      
      Note
      ----
          Directions refer to the ion motion direction e.g. the zonal
          scalar applies to zonal ion motions (meridional E field assuming ExB ion motion)
      
      Parameters
      ----------
      glats : list-like of floats (degrees)
          Geodetic (WGS84) latitude
      glons : list-like of floats (degrees)
          Geodetic (WGS84) longitude 
      alts : list-like of floats (km)
          Geodetic (WGS84) altitude, height above surface
      dates : list-like of datetimes
          Date and time for determination of scalars
      e_field_scaling_only : boolean (False)
          If True, method only calculates the electric field scalar, ignoring 
          changes in magnitude of B. Note ion velocity related to E/B.
          
      Returns
      -------
      dict
          array-like of scalars for translating ion drifts. Keys are,
          'north_zonal_drifts_scalar', 'north_mer_drifts_scalar', and similarly
          for southern locations. 'equator_mer_drifts_scalar' and 
          'equator_zonal_drifts_scalar' cover the mappings to the equator.
      
      """
      
      if step_size is None:
          step_size = 100.
      if max_steps is None:
          max_steps = 1000
      steps = np.arange(max_steps)
  
      # use spacecraft location to get ECEF
      ecef_xs, ecef_ys, ecef_zs = geodetic_to_ecef(glats, glons, alts)
  
      # prepare output
      eq_zon_drifts_scalar = []
      eq_mer_drifts_scalar = []
      # magnetic field info
      north_mag_scalar = []
      south_mag_scalar = []
      eq_mag_scalar = []
      out = {}
      # meridional e-field scalar map, can also be
      # zonal ion drift scalar map
      # print ('Starting Northern')
      north_zon_drifts_scalar, mind_plus, mind_minus = closed_loop_edge_lengths_via_footpoint(glats,
                                                          glons, alts, dates, 'north',
                                                          'meridional',
                                                          step_size=step_size,
                                                          max_steps=max_steps, 
                                                          edge_length=25.,
                                                          edge_steps=5)
  
      north_mer_drifts_scalar, mind_plus, mind_minus = closed_loop_edge_lengths_via_footpoint(glats,
                                                          glons, alts, dates, 'north',
                                                          'zonal',
                                                          step_size=step_size,
                                                          max_steps=max_steps, 
                                                          edge_length=25.,
                                                          edge_steps=5)
  
      # print ('Starting Southern')
      south_zon_drifts_scalar, mind_plus, mind_minus = closed_loop_edge_lengths_via_footpoint(glats,
                                                          glons, alts, dates, 'south',
                                                          'meridional',
                                                          step_size=step_size,
                                                          max_steps=max_steps, 
                                                          edge_length=25.,
                                                          edge_steps=5)
  
      south_mer_drifts_scalar, mind_plus, mind_minus = closed_loop_edge_lengths_via_footpoint(glats,
                                                          glons, alts, dates, 'south',
                                                          'zonal',
                                                          step_size=step_size,
                                                          max_steps=max_steps, 
                                                          edge_length=25.,
                                                          edge_steps=5)
      # print ('Starting Equatorial')                
      # , step_zon_apex2, mind_plus, mind_minus                                 
      eq_zon_drifts_scalar = closed_loop_edge_lengths_via_equator(glats, glons, alts, dates,
                                                          'meridional',
                                                          edge_length=25., 
                                                          edge_steps=5)
      # , step_mer_apex2, mind_plus, mind_minus                                                      
      eq_mer_drifts_scalar = closed_loop_edge_lengths_via_equator(glats, glons, alts, dates,
                                                          'zonal',
                                                          edge_length=25., 
                                                          edge_steps=5)
      # print ('Done with core')
      north_zon_drifts_scalar = north_zon_drifts_scalar/50. 
      south_zon_drifts_scalar = south_zon_drifts_scalar/50. 
      north_mer_drifts_scalar = north_mer_drifts_scalar/50. 
      south_mer_drifts_scalar = south_mer_drifts_scalar/50. 
      # equatorial 
      eq_zon_drifts_scalar = 50./eq_zon_drifts_scalar
      eq_mer_drifts_scalar = 50./eq_mer_drifts_scalar
  
      if e_field_scaling_only:
          # prepare output
          out['north_mer_fields_scalar'] = north_zon_drifts_scalar
          out['south_mer_fields_scalar'] = south_zon_drifts_scalar
          out['north_zon_fields_scalar'] = north_mer_drifts_scalar
          out['south_zon_fields_scalar'] = south_mer_drifts_scalar
          out['equator_mer_fields_scalar'] = eq_zon_drifts_scalar
          out['equator_zon_fields_scalar'] = eq_mer_drifts_scalar
      
      else:
          # figure out scaling for drifts based upon change in magnetic field
          # strength
          for ecef_x, ecef_y, ecef_z, glat, glon, alt, date in zip(ecef_xs, ecef_ys, ecef_zs, 
                                                                  glats, glons, alts, 
                                                                  dates):            
              yr, doy = pysat.utils.getyrdoy(date)
              double_date = float(yr) + float(doy) / 366.
              # get location of apex for s/c field line
              apex_x, apex_y, apex_z, apex_lat, apex_lon, apex_alt = apex_location_info(
                                                                          [glat], [glon], 
                                                                          [alt], [date])    
              # trace to northern footpoint
              sc_root = np.array([ecef_x, ecef_y, ecef_z])
              trace_north = field_line_trace(sc_root, double_date, 1., 120., 
                                          steps=steps,
                                          step_size=step_size, 
                                          max_steps=max_steps)
              # southern tracing
              trace_south = field_line_trace(sc_root, double_date, -1., 120., 
                                          steps=steps,
                                          step_size=step_size, 
                                          max_steps=max_steps)
              # footpoint location
              north_ftpnt = trace_north[-1, :]
              nft_glat, nft_glon, nft_alt = ecef_to_geodetic(*north_ftpnt)
              south_ftpnt = trace_south[-1, :]
              sft_glat, sft_glon, sft_alt = ecef_to_geodetic(*south_ftpnt)
      
              # scalar for the northern footpoint electric field based on distances
              # for drift also need to include the magnetic field, drift = E/B
              tbn, tbe, tbd, b_sc = igrf.igrf12syn(0, double_date, 1, alt, 
                                                  np.deg2rad(90.-glat), 
                                                  np.deg2rad(glon))
              # get mag field and scalar for northern footpoint
              tbn, tbe, tbd, b_nft = igrf.igrf12syn(0, double_date, 1, nft_alt, 
                                                  np.deg2rad(90.-nft_glat), 
                                                  np.deg2rad(nft_glon))
              north_mag_scalar.append(b_sc/b_nft)            
              # equatorial values
              tbn, tbe, tbd, b_eq = igrf.igrf12syn(0, double_date, 1, apex_alt, 
                                                   np.deg2rad(90.-apex_lat), 
                                                   np.deg2rad(apex_lon))
              eq_mag_scalar.append(b_sc/b_eq)        
              # scalar for the southern footpoint
              tbn, tbe, tbd, b_sft = igrf.igrf12syn(0, double_date, 1, sft_alt, 
                                                    np.deg2rad(90.-sft_glat), 
                                                    np.deg2rad(sft_glon))
              south_mag_scalar.append(b_sc/b_sft)
          
          # make E-Field scalars to drifts
          # lists to arrays
          north_mag_scalar = np.array(north_mag_scalar)
          south_mag_scalar = np.array(south_mag_scalar)
          eq_mag_scalar = np.array(eq_mag_scalar)
          # apply to electric field scaling to get ion drift values
          north_zon_drifts_scalar = north_zon_drifts_scalar*north_mag_scalar
          south_zon_drifts_scalar = south_zon_drifts_scalar*south_mag_scalar 
          north_mer_drifts_scalar = north_mer_drifts_scalar*north_mag_scalar
          south_mer_drifts_scalar = south_mer_drifts_scalar*south_mag_scalar
          # equatorial 
          eq_zon_drifts_scalar = eq_zon_drifts_scalar*eq_mag_scalar
          eq_mer_drifts_scalar = eq_mer_drifts_scalar*eq_mag_scalar
          # output
          out['north_zonal_drifts_scalar'] = north_zon_drifts_scalar
          out['south_zonal_drifts_scalar'] = south_zon_drifts_scalar
          out['north_mer_drifts_scalar'] = north_mer_drifts_scalar
          out['south_mer_drifts_scalar'] = south_mer_drifts_scalar
          out['equator_zonal_drifts_scalar'] = eq_zon_drifts_scalar
          out['equator_mer_drifts_scalar'] = eq_mer_drifts_scalar
  
      return out