Merge pull request #108 from cpinte/Mueller_matrices

Updated infrastructure for Mueller matrices

src/Makefile

@@ -94,16 +94,16 @@ ifeq ($(SYSTEM), ifort)
         STATIC_FLAGS += -static-libstdc++
     endif
     FC= ifort
-    FFLAGS= -fpp -O3 -no-prec-div -fp-model fast=2 -traceback
+    FFLAGS= -fpp -O3 -fp-model fast=2 -traceback
     DBLFLAG= -r8
-    DEBUGFLAG= -check all -C -g -fpe0 -traceback -no-ftz \
+    DEBUGFLAG= -check all -g -fpe0 -traceback -no-ftz \
                -warn uninitialized -warn unused -warn truncated_source
     KNOWN_SYSTEM= yes
     FOMPFLAGS= -fopenmp
     IPOFLAGS= -ipo
     CC= icc
     CXX=icpc -std=c++11
-    CCFLAGS= -O3 -no-prec-div -fp-model fast=2 -traceback
+    CCFLAGS= -O3  -fp-model fast=2 -traceback
     COMPFLAGS= -fopenmp
     LIBCXX= -cxxlib
     LIBS= $(LIB)/ifort

src/dust_transfer.f90

@@ -217,13 +217,7 @@ subroutine transfert_poussiere()
      p_icell = icell_ref
      if (aniso_method==2) write(*,*) "g             ", tab_g_pos(p_icell,1)
      write(*,*) "albedo        ", tab_albedo_pos(p_icell,1)
-     if (lsepar_pola.and.(scattering_method == 2)) then
-        if (lmueller) then
-           write(*,*) "polarisability", maxval(-tab_mueller_pos(1,2,:,p_icell,1))
-        else
-           write(*,*) "polarisability", maxval(-tab_s12_o_s11_pos(:,p_icell,1))
-        endif
-     endif
+     if (lsepar_pola.and.(scattering_method == 2)) write(*,*) "polarisability", maxval(-tab_s12_o_s11_pos(:,p_icell,1))
      if (lopacite_only) call exit(0)
 
      if (l_em_disk_image) then ! le disque �met
@@ -482,11 +476,7 @@ subroutine transfert_poussiere()
         lambda = ind_etape - first_etape_obs + 1
 
         if (.not.lMueller_pos_multi .and. lscatt_ray_tracing) then
-           if (lmueller) then
-              call calc_local_scattering_matrices_mueller(lambda, p_lambda) ! Todo : this is not good, we compute this twice
-           else
-              call calc_local_scattering_matrices(lambda, p_lambda)
-           endif
+           call calc_local_scattering_matrices(lambda, p_lambda) ! Todo : this is not good, we compute this twice
         endif
 
         if (lspherical.or.l3D) then
@@ -986,6 +976,7 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
   logical, intent(inout) :: flag_star, flag_ISM, lpacket_alive
   logical, intent(out) :: flag_scatt
 
+  real(kind=dp), dimension(4,4) :: M
   real(kind=dp) :: u1,v1,w1, phi, cospsi, w02, srw02, argmt
   integer :: taille_grain, itheta
   integer :: n_iterations
@@ -1102,14 +1093,8 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
               ! direction de propagation apres diffusion
               call cdapres(cospsi, phi, u, v, w, u1, v1, w1)
               if (lsepar_pola) then
-                 ! Nouveaux param�tres de Stokes
-                 if (lmueller) then
-                    call new_stokes_mueller(lambda,itheta,rand2,taille_grain,u,v,w,u1,v1,w1,stokes)
-                 else if (laggregate) then
-                    call  new_stokes_gmm(lambda,itheta,rand2,taille_grain,u,v,w,u1,v1,w1,stokes)
-                 else
-                    call new_stokes(lambda,itheta,rand2,taille_grain,u,v,w,u1,v1,w1,stokes)
-                 endif
+                 call get_Mueller_matrix_per_grain(lambda,itheta,rand2,taille_grain, M)
+                 call update_Stokes(Stokes,u,v,w,u1,v1,w1,M)
               endif
            else ! fonction de phase HG
               call hg(tab_g(taille_grain,lambda),rand, itheta, COSPSI) !HG
@@ -1122,7 +1107,6 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
               PHI = PI * ( 2.0 * rand - 1.0 )
               ! direction de propagation apres diffusion
               call cdapres(cospsi, phi, u, v, w, u1, v1, w1)
-              ! Param�tres de Stokes non modifi�s
            endif
 
         else ! methode 2 : diffusion sur la population de grains
@@ -1139,13 +1123,9 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
               PHI = PI * ( 2.0 * rand - 1.0 )
               ! direction de propagation apres diffusion
               call cdapres(cospsi, phi, u, v, w, u1, v1, w1)
-              ! Nouveaux param�tres de Stokes
               if (lsepar_pola) then
-                 if (lmueller) then
-                    call new_stokes_mueller_pos(p_lambda,itheta,rand2,p_icell,u,v,w,u1,v1,w1,Stokes)
-                 else
-                    call new_stokes_pos(p_lambda,itheta,rand2,p_icell,u,v,w,u1,v1,w1,Stokes)
-		 endif
+                 call get_Mueller_matrix_per_cell(lambda,itheta,rand2,p_icell, M)
+                 call update_Stokes(Stokes,u,v,w,u1,v1,w1,M)
 	      endif
            else ! fonction de phase HG
               call hg(tab_g_pos(p_icell,lambda),rand, itheta, cospsi) !HG
@@ -1158,7 +1138,6 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
               phi = pi * ( 2.0 * rand - 1.0 )
               ! direction de propagation apres diffusion
               call cdapres(cospsi, phi, u, v, w, u1, v1, w1)
-              ! Param�tres de Stokes non modifi�s
            endif
         endif
 

src/grains.f90

@@ -50,20 +50,18 @@ module grains
   real, dimension(:,:), allocatable :: amin, amax, aexp ! n_zones, n_especes
 
   ! Parametres de diffusion des grains
-  real, dimension(:,:,:), allocatable :: tab_s11, tab_s12, tab_s33, tab_s34, prob_s11 !n_lambda,n_grains,180
+  real, dimension(:,:,:), allocatable :: tab_s11, tab_s12, tab_s33, tab_s34, tab_s22, tab_s44, prob_s11 !n_lambda,n_grains,180
   real, dimension(:,:), allocatable :: tab_g, tab_albedo, C_ext, C_sca, C_abs, C_abs_norm !n_grains, n_lambda
   !real, dimension(:), allocatable :: q_geo ! n_grains section geometrique en m^2
 
   logical :: lforce_HG
   real :: forced_g
 
-  ! aggregats
-  real, dimension(:,:,:,:,:), allocatable :: tab_mueller !4,4, 180, n_grains,n_lambda
 
   ! Parametres de diffusion des cellules
   real, dimension(:,:), allocatable :: tab_albedo_pos, tab_g_pos ! n_cells,n_lambda
-  real, dimension(:,:,:), allocatable :: tab_s11_pos, tab_s12_o_s11_pos, tab_s33_o_s11_pos, tab_s34_o_s11_pos, prob_s11_pos ! 0:180, n_cells,n_lambda
-  real, dimension(:,:,:,:,:), allocatable :: tab_mueller_pos ! 4,4,0:180, n_cells,n_lambda
+  real, dimension(:,:,:), allocatable :: tab_s11_pos, tab_s12_o_s11_pos, tab_s33_o_s11_pos, tab_s34_o_s11_pos, &
+       tab_s22_o_s11_pos, tab_s44_o_s11_pos, prob_s11_pos ! 0:180, n_cells,n_lambda
 
   character(len=512) :: aggregate_file, mueller_aggregate_file, mueller_file
   real :: R_sph_same_M

src/init_mcfost.f90

@@ -46,8 +46,8 @@ subroutine set_default_variables()
   limg=.false.
   lorigine=.false.
   laggregate=.false.
-  lmueller=.false.
-  lper_size = .false.
+  lFresnel=.false.
+  lFresnel_per_size = .false.
   l3D=.false.
   lopacite_only=.false.
   lseed=.false.
@@ -522,27 +522,21 @@ subroutine initialisation_mcfost()
      case("-aggregate")
         laggregate=.true.
         i_arg = i_arg+1
-        if (lmueller) call error ("You can't use both -aggregate and -mueller options")
-        if (lper_size) call error ("You can't use both -aggregate and -mueller_size options")
         if (i_arg > nbr_arg) call error("GMM input file needed")
         call get_command_argument(i_arg,aggregate_file)
         i_arg = i_arg+1
         if (i_arg > nbr_arg) call error("GMM input file needed")
         call get_command_argument(i_arg,mueller_aggregate_file)
-     case("-mueller")
-        lmueller=.true.
+     case("-Fresnel")
+        lFresnel=.true.
         i_arg = i_arg+1
-        if (laggregate) call error ("You can't use both  -mueller and -aggregate options")
-        if (lper_size) call error ("You can't use both -mueller and -mueller_size options")
         if (i_arg > nbr_arg) call error("Mueller input file needed")
         call get_command_argument(i_arg,mueller_file)
         i_arg = i_arg+1
-     case("-mueller_size")
-        lper_size = .true.
+     case("-Fresnel_size")
+        lFresnel_per_size = .true.
         i_arg = i_arg+1
-        if (laggregate) call error ("You can't use both -mueller_size and -aggregate options")
-        if (lmueller) call error ("You can't use both -mueller_size and -mueller options")
-        lmueller=.true.
+        lFresnel=.true.
         if (i_arg > nbr_arg) call error("Mueller input pathfile needed")
         call get_command_argument(i_arg,mueller_file)
         i_arg = i_arg+1
@@ -1851,8 +1845,6 @@ end subroutine initialisation_mcfost
 !********************************************************************
 
 subroutine display_help()
-! Ajout du cas ou les matrices de Mueller sont donnees en entrees
-! 20/04/2023
 
   implicit none
 
@@ -1979,29 +1971,11 @@ subroutine display_help()
   write(*,*) "        : -op <wavelength> (microns) : computes dust properties at"
   write(*,*) "                                    specified wavelength and stops"
   write(*,*) "        : -aggregate <GMM_input_file> <GMM_output_file>"
-  write(*,*) "        : -mueller <Mueller_input_file> "
-  write(*,*) "     "
-  write(*,*) "        Mueller_input_file contain the mean mueller matrix averaged"
-  write(*,*) "        over the size distribution. Every element is divided by s11."
-  write(*,*) "        The format of the input file must be the following."
-  write(*,*) "     "
-  write(*,*) "             Qext        Qsca        <cos(theta)> "
-  write(*,*) "          Qext_value  Qsca_value  <cos(theta)>_value "
-  write(*,*) "     "
-  write(*,*) "     "
-  write(*,*) "                               Mueller Scattering Matrix "
-  write(*,*) "       an_value   s11_value   s12_value   s13_value    s14_value "
-  write(*,*) "                  s21_value   s22_value   s23_value    s24_value "
-  write(*,*) "                  s31_value   s32_value   s33_value    s34_value "
-  write(*,*) "                  s41_value   s42_value   s43_value    s44_value "
-  write(*,*) "       an_value   s11_value   s12_value   s13_value    s14_value "
-  write(*,*) "                  s21_value   s22_value   s23_value    s24_value "
-  write(*,*) "     ....... "
-  write(*,*) "     "
-  write(*,*) "        : -mueller_size <Mueller_input_pathfile> "
-  write(*,*) "                   Argument pathfile contain the size of each grain, and the path"
-  write(*,*) "                   for each associated matrix for every grain size, sorted"
-  write(*,*) "                   from the first to the last grain size considered."
+  write(*,*) "        : -Fresnel <Mueller_matrix_input_file>  Uses a Mueller matrix from the Fresnel database"
+  write(*,*) "        : -Fresnel_per_size <Mueller_matricesinput_pathfile> "
+!  write(*,*) "                   Argument pathfile contain the size of each grain, and the path"
+!  write(*,*) "                   for each associated matrix for every grain size, sorted"
+!  write(*,*) "                   from the first to the last grain size considered."
   write(*,*) "     "
   write(*,*) "        : -optical_depth_map ot -tau_map   : create an map of the optical depth"
   write(*,*) "        : -tau=1_surface : creates a map of the tau=1 surface"