code: def _chand(phi, muv, mus, taur):
      # FROM FUNCTION CHAND
      # phi: azimuthal difference between sun and observation in degree
      #      (phi=0 in backscattering direction)
      # mus: cosine of the sun zenith angle
      # muv: cosine of the observation zenith angle
      # taur: molecular optical depth
      # rhoray: molecular path reflectance
      # constant xdep: depolarization factor (0.0279)
      #          xfd = (1-xdep/(2-xdep)) / (1 + 2*xdep/(2-xdep)) = 2 * (1 - xdep) / (2 + xdep) = 0.958725775
      # */
      xfd = 0.958725775
      xbeta2 = 0.5
      #         float pl[5];
      #         double fs01, fs02, fs0, fs1, fs2;
      as0 = [
          0.33243832,
          0.16285370,
          -0.30924818,
          -0.10324388,
          0.11493334,
          -6.777104e-02,
          1.577425e-03,
          -1.240906e-02,
          3.241678e-02,
          -3.503695e-02,
      ]
      as1 = [0.19666292, -5.439061e-02]
      as2 = [0.14545937, -2.910845e-02]
      #         float phios, xcos1, xcos2, xcos3;
      #         float xph1, xph2, xph3, xitm1, xitm2;
      #         float xlntaur, xitot1, xitot2, xitot3;
      #         int i, ib;
  
      xph1 = (
          1.0
          + (3.0 * mus * mus - 1.0)
          * (3.0 * muv * muv - 1.0)
          * xfd
          / 8.0
      )
      xph2 = (
          -xfd
          * xbeta2
          * 1.5
          * mus
          * muv
          * np.sqrt(1.0 - mus * mus)
          * np.sqrt(1.0 - muv * muv)
      )
      xph3 = (
          xfd * xbeta2 * 0.375 * (1.0 - mus * mus) * (1.0 - muv * muv)
      )
  
      # pl[0] = 1.0
      # pl[1] = mus + muv
      # pl[2] = mus * muv
      # pl[3] = mus * mus + muv * muv
      # pl[4] = mus * mus * muv * muv
  
      fs01 = (
          as0[0]
          + (mus + muv) * as0[1]
          + (mus * muv) * as0[2]
          + (mus * mus + muv * muv) * as0[3]
          + (mus * mus * muv * muv) * as0[4]
      )
      fs02 = (
          as0[5]
          + (mus + muv) * as0[6]
          + (mus * muv) * as0[7]
          + (mus * mus + muv * muv) * as0[8]
          + (mus * mus * muv * muv) * as0[9]
      )
      #         for (i = 0; i < 5; i++) {
      #                 fs01 += (double) (pl[i] * as0[i]);
      #                 fs02 += (double) (pl[i] * as0[5 + i]);
      #         }
  
      # for refl, (ah2o, bh2o, ao3, tau) in zip(reflectance_bands, coefficients):
  
      # ib = _find_coefficient_index(center_wl)
      # if ib is None:
      #     raise ValueError("Can't handle band with wavelength '{}'".format(center_wl))
  
      xlntaur = np.log(taur)
  
      fs0 = fs01 + fs02 * xlntaur
      fs1 = as1[0] + xlntaur * as1[1]
      fs2 = as2[0] + xlntaur * as2[1]
      del xlntaur, fs01, fs02
  
      trdown = np.exp(-taur / mus)
      trup = np.exp(-taur / muv)
  
      xitm1 = (1.0 - trdown * trup) / 4.0 / (mus + muv)
      xitm2 = (1.0 - trdown) * (1.0 - trup)
      xitot1 = xph1 * (xitm1 + xitm2 * fs0)
      xitot2 = xph2 * (xitm1 + xitm2 * fs1)
      xitot3 = xph3 * (xitm1 + xitm2 * fs2)
      del xph1, xph2, xph3, xitm1, xitm2, fs0, fs1, fs2
  
      phios = np.deg2rad(phi + 180.0)
      xcos1 = 1.0
      xcos2 = np.cos(phios)
      xcos3 = np.cos(2.0 * phios)
      del phios
  
      rhoray = (
          xitot1 * xcos1 + xitot2 * xcos2 * 2.0 + xitot3 * xcos3 * 2.0
      )
      return rhoray, trdown, trup
  
imports: from __future__ import annotations
  import logging
  from typing import Optional, Type, Union
  import dask.array
  import numpy
  import xarray
  from satpy.dataset.dataid import WavelengthRange
label_desc: Use image node_id 9 for calling a global function with example usage: _chand(phi, muv, mus, taur) and returns: rhoray,  trdown,  trup