code: def get_total_spatial_integral(self, z=None):  
          
  
          dL= self.l_max.value-self.l_min.value if self.l_max.value > self.l_min.value else 360 + self.l_max.value - self.l_max.value
  
          #integral -inf to inf exp(-b**2 / 2*sigma_b**2 ) db = sqrt(2pi)*sigma_b 
          #Note that K refers to the peak diffuse flux (at b = 0) per square degree.
          integral = np.sqrt( 2*np.pi ) * self.sigma_b.value * self.K.value * dL 
  
          if isinstance( z, u.Quantity):
              z = z.value
          return integral * np.power( 180. / np.pi, -2 ) * np.ones_like( z )
docs: Returns the total integral (for 2D functions) or the integral over the spatial components (for 3D functions).
          needs to be implemented in subclasses.
  
          :return: an array of values of the integral (same dimension as z).