docs: def _compute_VesPoly(R=2.4, r=1., elong=0., Dshape=0.,
                      divlow=True, divup=True, nP=200):
      """ Utility to compute three 2D (R,Z) polygons
  
      One represents a vacuum vessel, one an outer bumper, one a baffle
  
      The vessel polygon is centered on (R,0.), with minor radius r
      It can have a vertical (>0) or horizontal(<0) elongation in [-1;1]
      It can be D-shaped (Dshape in [0.,1.], typically 0.2)
      It can be non-convex, with:
          * a lower divertor-like shape
          * a upper divertor-like shape
      The elongation also affects the outer bumper and baffle
  
      Parameters
      ----------
      R:          int / float
          Major radius used as a center of the vessel
      r :         int / float
          Minor radius of the vessel
      elong:      int / float
          Dimensionless elongation parameter in [-1;1]
      Dshape:     int / float
          Dimensionless parameter for the D-shape (in-out asymmetry) in [0;1]
      divlow:     bool
          Flag indicating whether to incude a lower divertor-like shape
      divup:      bool
          Flag indicating whether to incude an upper divertor-like shape
      nP :        int
          Parameter specifying approximately the number of points of the vessel
  
      Return
      ------
      poly:       np.ndarray
          Closed (2,nP) polygon of the vacuum vessel, optionnally with divertors
      pbump:      np.ndarray
          Closed (2,N) polygon defining the outer bumper
      pbaffle:    np.ndarray
          Closed (2,N) polygon defining the lower baffle
      """
  
      # Basics (center, theta, unit vectors)
      cent = np.r_[R,0.]
      theta = np.linspace(-np.pi,np.pi,nP)
      poly = np.array([np.cos(theta), np.sin(theta)])
  
      # Divertors
      pdivR = np.r_[-0.1,0.,0.1]
      pdivZ = np.r_[-0.1,0.,-0.1]
      if divlow:
          ind = (np.sin(theta)<-0.85).nonzero()[0]
          pinsert = np.array([pdivR, -1.+pdivZ])
          poly = np.concatenate((poly[:,:ind[0]], pinsert, poly[:,ind[-1]+1:]),
                                axis=1)
  
      if divup:
          theta = np.arctan2(poly[1,:], poly[0,:])
          ind = (np.sin(theta)>0.85).nonzero()[0]
          pinsert = np.array([pdivR[::-1], 1.-pdivZ])
          poly = np.concatenate((poly[:,:ind[0]], pinsert, poly[:,ind[-1]+1:]),
                                axis=1)
  
      # Modified radius (by elongation and Dshape)
      rbis = r*np.hypot(poly[0,:],poly[1,:])
      theta = np.arctan2(poly[1,:],poly[0,:])
      rbis = rbis*(1+elong*0.15*np.sin(2.*theta-np.pi/2.))
      if Dshape>0.:
          ind = np.cos(theta)<0.
          coef = 1 + Dshape*(np.sin(theta[ind])**2-1.)
          rbis[ind] = rbis[ind]*coef
  
      er = np.array([np.cos(theta), np.sin(theta)])
      poly = cent[:,np.newaxis] + rbis[np.newaxis,:]*er
  
      # Outer bumper
      Dbeta = 2.*np.pi/6.
      beta = np.linspace(-Dbeta/2.,Dbeta/2., 20)
      pbRin = 0.85*np.array([np.cos(beta), np.sin(beta)])
      pbRout = 0.95*np.array([np.cos(beta), np.sin(beta)])[:,::-1]
      pinsert = np.array([[0.95,1.05,1.05,0.95],
                          [0.05,0.05,-0.05,-0.05]])
  
      ind = (np.abs(pbRout[1,:])<0.05).nonzero()[0]
      pbump = (pbRin, pbRout[:,:ind[0]], pinsert,
               pbRout[:,ind[-1]+1:], pbRin[:,0:1])
      pbump = np.concatenate(pbump, axis=1)
      theta = np.arctan2(pbump[1,:],pbump[0,:])
      er = np.array([np.cos(theta), np.sin(theta)])
      rbis = r*(np.hypot(pbump[0,:],pbump[1,:])
                *(1.+elong*0.15*np.sin(2.*theta-np.pi/2.)))
      pbump = cent[:,np.newaxis] + rbis[np.newaxis,:]*er
  
      # Baffle
      offR, offZ = 0.1, -0.85
      wR, wZ = 0.2, 0.05
      pbaffle = np.array([offR + wR*np.r_[-1,1,1,-1,-1],
                          offZ + wZ*np.r_[1,1,-1,-1,1]])
      theta = np.arctan2(pbaffle[1,:],pbaffle[0,:])
      er = np.array([np.cos(theta), np.sin(theta)])
      rbis = r*(np.hypot(pbaffle[0,:],pbaffle[1,:])
                *(1.+elong*0.15*np.sin(2.*theta-np.pi/2.)))
      pbaffle = cent[:,np.newaxis] + rbis[np.newaxis,:]*er
  
      return poly, pbump, pbaffle