code: def sunzen_corr_cos(data, cos_zen, limit=88., max_sza=95.):
      
  
      # Convert the zenith angle limit to cosine of zenith angle
      limit_rad = np.deg2rad(limit)
      limit_cos = np.cos(limit_rad)
      max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza
  
      # Cosine correction
      corr = 1. / cos_zen
      if max_sza is not None:
          # gradually fall off for larger zenith angle
          grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad)
          # invert the factor so maximum correction is done at `limit` and falls off later
          grad_factor = 1. - np.log(grad_factor + 1) / np.log(2)
          # make sure we don't make anything negative
          grad_factor = grad_factor.clip(0.)
      else:
          # Use constant value (the limit) for larger zenith angles
          grad_factor = 1.
      corr = corr.where(cos_zen > limit_cos, grad_factor / limit_cos)
      # Force "night" pixels to 0 (where SZA is invalid)
      corr = corr.where(cos_zen.notnull(), 0)
  
      return data * corr
docs: Perform Sun zenith angle correction.
  
      The correction is based on the provided cosine of the zenith
      angle (``cos_zen``).  The correction is limited
      to ``limit`` degrees (default: 88.0 degrees).  For larger zenith
      angles, the correction is the same as at the ``limit`` if ``max_sza``
      is `None`. The default behavior is to gradually reduce the correction
      past ``limit`` degrees up to ``max_sza`` where the correction becomes
      0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same shape.