test_spectroscopy2/ was included.

phoebe/backend/pyterpolmini.py

@@ -444,7 +444,7 @@ def truncate_spectrum(self, wmin=None, wmax=None):
         if (self.wave.min() > wmin) | (self.wave.max() < wmax):
             raise ValueError('The spectrum %s does not cover the whole spectral region <%s,%s>.' % (str(self).rstrip('\n'), str(wmin), str(wmax)))
 
-        ind = np.where((self.wave >= wmin) & (self.wave <= wmax))[0]
+        ind = np.where((self.wave >= wmin-ZERO_TOLERANCE) & (self.wave <= wmax+ZERO_TOLERANCE))[0]
         self.wave = self.wave[ind]
         self.intens = self.intens[ind]
 
@@ -506,8 +506,9 @@ def set_mode(self, mode):
         for i, d in enumerate(directories):
             gridlist = os.path.join(self.grid_directory, d, 'gridlist')
             directory = os.path.join(self.grid_directory, d)
+            family = families[i]
 
-            self.read_gridlist(gridlist, columns=columns, directory=directory, family=families)
+            self.read_gridlist(gridlist, columns=columns, directory=directory, family=family)
 
     def __str__(self):
         """String representation."""
@@ -779,7 +780,7 @@ def resolve_degeneracy(self, spectra):
         """
         ind = []
         for i in range(0, len(spectra)):
-            ind.append(self.grid_order.index(spectra[i]['family']))
+            ind.append(self.grid_order.index(spectra[i].props['family']))
 
         ind = np.array(ind)
         if np.any(ind < -1):

phoebe/backend/spectroscopy.py

@@ -144,11 +144,11 @@ def spe_integrate(b, system, wavelengths=None, info={}, k=None):
     loggs = meshes.get_column_flat('loggs', components)
     zs = 10.0**meshes.get_column_flat('abuns', components)
 
-    Lum = abs_intensities*areas*mus*visibilities		# J s^-1 m^-1
+    Lum = abs_intensities*areas*mus*visibilities	# J s^-1 m^-1
 
-    step = 0.01					# Ang
-    angstroms = wavelengths*1.0e10		# Ang
-    fluxes = np.zeros(len(wavelengths))		# 1
+    step = 0.01						# Ang
+    angstroms = wavelengths*1.0e10			# Ang
+    fluxes = np.zeros(len(wavelengths))			# 1
 
     for i in range(len(Lum)):
         if Lum[i] == 0.0:
@@ -203,7 +203,14 @@ def sed_integrate(b, system, wavelengths=None, bandwidths=None, info={}, k=None)
     loggs = meshes.get_column_flat('loggs', components)
     zs = 10.0**meshes.get_column_flat('abuns', components)
 
+    d = system.distance				# m
     Lum = areas*mus*visibilities		# m^2
+    Lum /= d**2					# 1
+
+    # Note: a factor 1/pi is needed to obtain the solar values:
+    # F_lambda ~ 2.e9 W m^-2 m^-1 (at Earth, 550 nm; Verbunt 2008)
+    # F = \int F_lambda dlambda = 1363 W m^-2 (Kopp & Lean 2011)
+    Lum /= np.pi				# 1
 
     step = 0.1					# Ang
     angstroms = wavelengths*1.0e10		# Ang

tests/nosetests/test_spectroscopy2/Sed.plt

@@ -0,0 +1,16 @@
+#!/usr/bin/gnuplot
+
+set xl "wave [A]"
+set yl "flux [1]"
+
+set zeroaxis
+
+p \
+  "Sed.dat" u ($2/1.e-10):3 w p,\
+
+pa -1
+
+set term png small
+set out "Sed.png"
+rep
+

tests/nosetests/test_spectroscopy2/Sed.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python3
+
+f = open("Sed.dat", "w")
+f.write("# time wavelength flux sigma\n")
+
+t1 = 0.0            # d
+t2 = 0.0            # d
+dt = 0.25           # d
+wave1 = 6500.0e-10  # m
+wave2 = 6600.0e-10  # m
+dwave = 0.1e-10     # m
+flux = 1.0          # 1
+sigma = 0.01        # 1
+
+t = t1
+while t <= t2+0.5*dt:
+
+    wave = wave1
+    while wave < wave2+0.5*dwave:
+       f.write("%.8f  %.8e  %.8f  %.8f\n" % (t, wave, flux, sigma))
+       wave += dwave
+
+    t += dt
+

tests/nosetests/test_spectroscopy2/grids_ABS/TEST/gridlist

@@ -0,0 +1,4 @@
+lte05700-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.vis.dat 05700  4.00 1.0
+lte05700-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.vis.dat 05700  4.50 1.0
+lte05800-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.vis.dat 05800  4.00 1.0
+lte05800-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.vis.dat 05800  4.50 1.0