added raytracing scattering test

projects/test/raytracing_scattering/POLARIS.cmd

@@ -0,0 +1,57 @@
+<common>
+
+    <dust_component>    "input/dust_nk/silicate_d03.nk" "plaw" 1.0 3800.0 1e-06 1e-06 -3.5
+    <phase_function>    PH_MIE
+
+    <mass_fraction>    0.01
+
+    <nr_threads>    -1
+
+</common>
+
+<task> 1
+
+    <cmd>    CMD_DUST_EMISSION
+
+    <detector_dust nr_pixel = "255*255">    1e-6    1e-3    4    1    0.0    0.0    4.32e+18
+
+    <max_subpixel_lvl>    2
+
+    <path_grid>    "projects/test/raytracing_scattering/grid_3D_sphere_const_T_m1e-5.dat"
+    <path_out>    "projects/test/raytracing_scattering/dust/"
+
+</task>
+
+<task> 1
+
+    <cmd>    CMD_DUST_EMISSION
+
+    <detector_dust nr_pixel = "255*255">    1e-6    1e-3    4    1    0.0    0.0    4.32e+18
+
+    <source_star nr_photons = "1e6">    0    0    0    2    4500
+
+    <max_subpixel_lvl>    2
+
+    <rt_scattering> 1
+
+    <path_grid>    "projects/test/raytracing_scattering/grid_3D_sphere_const_T_m1e-5.dat"
+    <path_out>    "projects/test/raytracing_scattering/dust_rt/"
+
+</task>
+
+<task> 1
+
+    <cmd>    CMD_DUST_SCATTERING
+
+    <detector_dust_mc nr_pixel = "255*255">    1e-6    1e-3    4    0.00    0.00    4.32e+18
+
+    <source_star nr_photons = "1e6">    0    0    0    2    4500
+    <source_dust nr_photons = "2.1e6">
+
+    <path_grid>    "projects/test/raytracing_scattering/grid_3D_sphere_const_T_m1e-5.dat"
+    <path_out>    "projects/test/raytracing_scattering/dust_mc/"
+
+    <peel_off>    1
+    <enfsca>    1
+
+</task>

projects/test/raytracing_scattering/compare.py

@@ -0,0 +1,65 @@
+from astropy.io import fits
+from astropy import units as u
+import numpy as np
+
+
+"""
+----------------------------------------------
+Raytracing and scattering 
+----------------------------------------------
+
+Test whether scattering in a dust simulation
+(from radiation field) looks the same as
+scattering in a dust_mc simulation.
+
+In order to pass this test, all detectors
+(with different orientations) have to yield
+the same stellar + scattered + emitted flux 
+for each wavelength.
+(Should get better with more photons)
+"""
+
+
+def read_data(sed_fits_file):
+    fits_header = fits.getheader(sed_fits_file)
+    fits_data = fits.getdata(sed_fits_file)
+
+    nr_wave = np.shape(fits_data)[2]
+    sed_wavelengths = np.zeros(nr_wave) * u.m
+    for i_wave in range(nr_wave):
+        sed_wavelengths[i_wave] = float(fits_header[f'HIERARCH WAVELENGTH{i_wave+1}']) * u.m
+
+    _stokes = ['I', 'Q', 'U', 'V']
+    sed_data = {}
+    for i_s, i_stokes in enumerate(_stokes):
+        sed_data[i_stokes] = fits_data[i_s,0,:] * u.Jy
+
+    return sed_data
+
+
+def compare():
+    dust_sed_data = read_data('projects/test/raytracing_scattering/dust/data/polaris_detector_nr0001_sed.fits.gz')
+    dust_rt_sed_data = read_data('projects/test/raytracing_scattering/dust_rt/data/polaris_detector_nr0001_sed.fits.gz')
+    dust_mc_sed_data = read_data('projects/test/raytracing_scattering/dust_mc/data/polaris_detector_nr0001_sed.fits.gz')
+
+    max_rel_diff = np.max(np.abs( (dust_sed_data['I'] + dust_mc_sed_data['I']) / dust_rt_sed_data['I'] - 1.0 ))
+    if max_rel_diff > 1e-1:
+        raise Exception(f'Test failed: Stokes I does not match (max. relative difference = {max_rel_diff})')
+
+    mc_polarization = np.sqrt(dust_mc_sed_data['Q']**2 + dust_mc_sed_data['U']**2 + dust_mc_sed_data['V']**2) / (dust_mc_sed_data['I'] + dust_sed_data['I'])
+    rt_polarization = np.sqrt(dust_rt_sed_data['Q']**2 + dust_rt_sed_data['U']**2 + dust_rt_sed_data['V']**2) / dust_rt_sed_data['I']
+    max_abs_diff = np.max(np.abs( mc_polarization - rt_polarization ))
+    if max_abs_diff > 1e-2:
+        raise Exception(f'Test failed: Polarization does not match (max. absolute difference = {max_abs_diff})')
+
+    max_polarization = np.max(mc_polarization)
+    if max_polarization > 1e-2:
+        raise Exception(f'Test failed: Polarization is too large (max. value = {max_polarization})')
+
+    return True
+
+
+if __name__ == '__main__':
+    res = compare()
+    if res:
+        print('Test passed')

projects/test/reemission_sphere/POLARIS.cmd

@@ -5,8 +5,6 @@
 
 	<mass_fraction>		0.01
 
-	<write_dust_files> 1
-
 	<nr_threads>		-1
 
 </common>

projects/test/reemission_sphere/compare.py

@@ -117,7 +117,7 @@ def calc_flux_ana():
     return flux_sphere_ana.to(u.Jy, equivalencies=u.spectral_density(wavelength))
 
 
-def read_data(sed_fits_file, map_fits_file, stokes='I'):
+def read_data(sed_fits_file, map_fits_file):
     sed_fits_header = fits.getheader(sed_fits_file)
     sed_fits_data = fits.getdata(sed_fits_file)
 
@@ -126,13 +126,10 @@ def read_data(sed_fits_file, map_fits_file, stokes='I'):
     for i_wave in range(nr_wave):
         sed_wavelengths[i_wave] = float(sed_fits_header[f'HIERARCH WAVELENGTH{i_wave+1}']) * u.m
 
-    _stokes = ['I', 'Q', 'U', 'V', 'TAU']
+    _stokes = ['I', 'Q', 'U', 'V']
     sed_data = {}
     for i_s, i_stokes in enumerate(_stokes):
-        if i_stokes == 'TAU':
-            sed_data[i_stokes] = sed_fits_data[i_s,0,:] * u.dimensionless_unscaled
-        else:
-            sed_data[i_stokes] = sed_fits_data[i_s,0,:] * u.Jy
+        sed_data[i_stokes] = sed_fits_data[i_s,0,:] * u.Jy
 
     map_fits_header = fits.getheader(map_fits_file)
     map_fits_data = fits.getdata(map_fits_file)
@@ -147,7 +144,7 @@ def read_data(sed_fits_file, map_fits_file, stokes='I'):
     for i_s, i_stokes in enumerate(_stokes):
         map_data[i_stokes] = map_fits_data[i_s,:,:,:] * u.Jy / u.pix
 
-    return sed_data[stokes], map_data[stokes]
+    return sed_data, map_data
 
 
 def compare():
@@ -159,19 +156,19 @@ def compare():
         'projects/test/reemission_sphere/dust/data/polaris_detector_nr0002.fits.gz')
     reference = calc_flux_ana()
 
-    max_rel_diff = np.max(np.abs( sed_data_pol / reference - 1.0 ))
+    max_rel_diff = np.max(np.abs( sed_data_pol['I'] / reference - 1.0 ))
     if max_rel_diff > 1e-3:
         raise Exception(f'Test failed: Polar detector and reference do not match (max. relative difference = {max_rel_diff})')
 
-    max_rel_diff = np.max(np.abs( sed_data_car / sed_data_pol - 1.0 ))
+    max_rel_diff = np.max(np.abs( sed_data_car['I'] / sed_data_pol['I'] - 1.0 ))
     if max_rel_diff > 1e-3:
         raise Exception(f'Test failed: Cartesian and polar detector do not match (max. relative difference = {max_rel_diff})')
 
-    max_rel_diff = np.max(np.abs( np.sum(map_data_pol, axis=(1,2)) * u.pix / reference - 1.0 ))
+    max_rel_diff = np.max(np.abs( np.sum(map_data_pol['I'], axis=(1,2)) * u.pix / reference - 1.0 ))
     if max_rel_diff > 1e-3:
         raise Exception(f'Test failed: Sum of polar map detector and reference do not match (max. relative difference = {max_rel_diff})')
 
-    max_rel_diff = np.max(np.abs( np.sum(map_data_car, axis=(1,2)) / np.sum(map_data_pol, axis=(1,2)) - 1.0 ))
+    max_rel_diff = np.max(np.abs( np.sum(map_data_car['I'], axis=(1,2)) / np.sum(map_data_pol['I'], axis=(1,2)) - 1.0 ))
     if max_rel_diff > 1e-3:
         raise Exception(f'Test failed: Sum of cartesian and polar map detector do not match (max. relative difference = {max_rel_diff})')
 

projects/test/stellar_scattering_sphere/POLARIS.cmd

@@ -5,8 +5,6 @@
 
 	<mass_fraction>		0.01
 
-	<write_dust_files> 1
-
 	<nr_threads>		-1
 
 </common>

projects/test/stellar_scattering_sphere/compare.py

@@ -36,7 +36,7 @@ def read_data(sed_fits_file):
     _stokes = ['I', 'Q', 'U', 'V']
     sed_data = {}
     for i_s, i_stokes in enumerate(_stokes):
-        sed_data[i_stokes] = fits_data[i_s,:,:] * u.Jy
+        sed_data[i_stokes] = fits_data[i_s,0,:] * u.Jy
 
     return sed_data
 

projects/test/stellar_sed/POLARIS.cmd

@@ -5,8 +5,6 @@
 
 	<mass_fraction>		0.01
 
-	<write_dust_files> 1
-
 	<nr_threads>		-1
 
 </common>

projects/test/stellar_sed/compare.py

@@ -23,7 +23,7 @@ def flux_planck_law(wavelength, temperature=4500*u.K, radius=2*u.R_sun, distance
     return flux.to(u.Jy, equivalencies=u.spectral_density(wavelength))
 
 
-def read_data(sed_fits_file, stokes='I'):
+def read_data(sed_fits_file):
     fits_header = fits.getheader(sed_fits_file)
     fits_data = fits.getdata(sed_fits_file)
 
@@ -35,16 +35,16 @@ def read_data(sed_fits_file, stokes='I'):
     _stokes = ['I', 'Q', 'U', 'V']
     sed_data = {}
     for i_s, i_stokes in enumerate(_stokes):
-        sed_data[i_stokes] = fits_data[i_s,:,:] * u.Jy
+        sed_data[i_stokes] = fits_data[i_s,0,:] * u.Jy
 
-    return sed_wavelengths, sed_data[stokes]
+    return sed_wavelengths, sed_data
 
 
 def compare():
     sed_wavelengths, sed_data = read_data('projects/test/stellar_sed/dust_mc/data/polaris_detector_nr0001_sed.fits.gz')
     reference = flux_planck_law(sed_wavelengths)
 
-    max_rel_diff = np.max( sed_data[0] / reference - 1.0 )
+    max_rel_diff = np.max( sed_data['I'] / reference - 1.0 )
     if max_rel_diff > 1e-4:
         raise Exception(f'Test failed: POLARIS and reference do not match (max. relative difference = {max_rel_diff})')