code: def rotfreq_to_omega(rotfreq, scale=c.R_sun.si.value, solar_units=False):
      
      if solar_units:
          omega = rotfreq / (2*np.pi) / np.sqrt(c.GM_sun.to(u.solRad**3/u.d**2).value/scale**3)
      else:
          # then SI units
          omega = rotfreq / (2*np.pi) / np.sqrt(c.GM_sun.value/scale**3)
  
      # print "*** rotstar.rotfreq_to_omega", rotfreq, scale, solar_units, omega
  
      return omega
docs: TODO: add documentation
  
      NOTE: everything MUST be in consistent units according to solar_units bool