code: def example():
      
      # --- input data filename ---
      infile = os.path.join(os.path.dirname(__file__),
                            '../../ExampleDataFiles/MarsTopo719.shape')
      coeffs, lmax = shio.shread(infile)
  
      # --- plot grid ---
      grid = expand.MakeGridDH(coeffs, csphase=-1)
      fig_map = plt.figure()
      plt.imshow(grid)
  
      # ---- compute spectrum ----
      ls = np.arange(lmax + 1)
      pspectrum = spectralanalysis.spectrum(coeffs, unit='per_l')
      pdensity = spectralanalysis.spectrum(coeffs, unit='per_lm')
  
      # ---- plot spectrum ----
      fig_spectrum, ax = plt.subplots(1, 1)
      ax.set_xscale('log')
      ax.set_yscale('log')
      ax.set_xlabel('degree l')
      ax.grid(True, which='both')
  
      ax.plot(ls[1:], pspectrum[1:], label='power per degree l')
      ax.plot(ls[1:], pdensity[1:], label='power per degree l and order m')
  
      ax.legend()
  
      fig_map.savefig('SHRtopography_mars.png')
      fig_spectrum.savefig('SHRspectrum_mars.png')
      print('mars topography and spectrum saved')
docs: example that plots the power spectrum of Mars topography data