docs: def closed_loop_edge_lengths_via_footpoint(glats, glons, alts, dates, direction,
                                             vector_direction, step_size=None, 
                                             max_steps=None, edge_length=25., 
                                             edge_steps=5):
      """
      Forms closed loop integration along mag field, satrting at input
      points and goes through footpoint. At footpoint, steps along vector direction
      in both positive and negative directions, then traces back to opposite
      footpoint. Back at input location, steps toward those new field lines 
      (edge_length) along vector direction until hitting distance of minimum
      approach. Loops don't always close. Returns total edge distance 
      that goes through input location, along with the distances of closest approach. 
      
      Note
      ----
          vector direction refers to the magnetic unit vector direction 
      
      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
      direction : string
          'north' or 'south' for tracing through northern or
          southern footpoint locations
      vector_direction : string
          'meridional' or 'zonal' unit vector directions
      step_size : float (km)
          Step size (km) used for field line integration
      max_steps : int
          Number of steps taken for field line integration
      edge_length : float (km)
          Half of total edge length (step) taken at footpoint location.
          edge_length step in both positive and negative directions.
      edge_steps : int
          Number of steps taken from footpoint towards new field line
          in a given direction (positive/negative) along unit vector
          
      Returns
      -------
      np.array, np.array, np.array
          A closed loop field line path through input location and footpoint in 
          northern/southern hemisphere and back is taken. The return edge length
          through input location is provided. 
          
          The distances of closest approach for the positive step along vector
          direction, and the negative step are returned.
  
      
      """
      
      if step_size is None:
          step_size = 100.
      if max_steps is None:
          max_steps = 1000
      steps = np.arange(max_steps)
  
      if direction == 'south':
          direct = -1
      elif direction == 'north':
          direct = 1
  
      # use spacecraft location to get ECEF
      ecef_xs, ecef_ys, ecef_zs = geodetic_to_ecef(glats, glons, alts)
  
      # prepare output
      full_local_step = []
      min_distance_plus = []
      min_distance_minus = []
  
      for ecef_x, ecef_y, ecef_z, glat, glon, alt, date in zip(ecef_xs, ecef_ys, ecef_zs, 
                                                               glats, glons, alts, 
                                                               dates):
          # going to try and form close loops via field line integration
          # start at location of interest, map down to northern or southern 
          # footpoints then take symmetric steps along meridional and zonal 
          # directions and trace back from location of interest, step along 
          # field line directions until we intersect or hit the distance of 
          # closest approach to the return field line with the known 
          # distances of footpoint steps, and the closet approach distance
          # we can determine the scalar mapping of one location to another
                      
          yr, doy = pysat.utils.getyrdoy(date)
          double_date = float(yr) + float(doy) / 366.
  
          # print (glat, glon, alt)
          # trace to footpoint, starting with input location
          sc_root = np.array([ecef_x, ecef_y, ecef_z])
          trace = field_line_trace(sc_root, double_date, direct, 120., 
                                   steps=steps,
                                   step_size=step_size, 
                                   max_steps=max_steps)
          # pull out footpoint location
          ftpnt = trace[-1, :]
          ft_glat, ft_glon, ft_alt = ecef_to_geodetic(*ftpnt)
          
          # take step from footpoint along + vector direction
          plus_step = step_along_mag_unit_vector(ftpnt[0], ftpnt[1], ftpnt[2], 
                                                 date, 
                                                 direction=vector_direction,
                                                 num_steps=edge_steps,
                                                 step_size=edge_length/edge_steps)
          # trace this back to other footpoint
          other_plus = field_line_trace(plus_step, double_date, -direct, 0., 
                                        steps=steps,
                                        step_size=step_size, 
                                        max_steps=max_steps)
          # take half step from first footpoint along - vector direction
          minus_step = step_along_mag_unit_vector(ftpnt[0], ftpnt[1], ftpnt[2], 
                                                 date, 
                                                 direction=vector_direction, 
                                                 scalar=-1,
                                                 num_steps=edge_steps,
                                                 step_size=edge_length/edge_steps)
          # trace this back to other footpoint
          other_minus = field_line_trace(minus_step, double_date, -direct, 0., 
                                         steps=steps,
                                         step_size=step_size, 
                                         max_steps=max_steps)
          # need to determine where the intersection of field line coming back from
          # footpoint through postive vector direction step and back
          # in relation to the vector direction from the s/c location. 
          pos_edge_length, _, mind_pos = step_until_intersect(sc_root,
                                          other_plus,
                                          1, date, 
                                          direction=vector_direction,
                                          field_step_size=step_size,
                                          step_size_goal=edge_length/edge_steps)        
          # take half step from S/C along - vector direction 
          minus_edge_length, _, mind_minus = step_until_intersect(sc_root,
                                          other_minus,
                                          -1, date, 
                                          direction=vector_direction,
                                          field_step_size=step_size,
                                          step_size_goal=edge_length/edge_steps)
          # collect outputs
          full_local_step.append(pos_edge_length + minus_edge_length)
          min_distance_plus.append(mind_pos)
          min_distance_minus.append(mind_minus)
          
      return np.array(full_local_step), np.array(min_distance_plus), np.array(min_distance_minus)