Valentin version

dust_temperature.py

@@ -41,7 +41,7 @@ def plot_dust_temperature(mystring):
     from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator, LogLocator, LogFormatter)
 
     matplotlib.rcParams.update({'font.size': 20})
-    matplotlib.rc('font', family='Arial')
+    #matplotlib.rc('font', family='Arial')
     fontcolor='white'
 
 

fargo2radmc3d.py

@@ -4,18 +4,16 @@
 # Program post-processes FARGO2D/3D outputs as inputs for RADMC-3D.
 # It can run with either Python 2.X or Python 3.X.
 # 
-# code written by Clement Baruteau (CB), Sebastian Perez (SP) and
-# Marcelo Barraza (MB) with substantial contributions from Simon
-# Casassus (SC) and Gaylor Wafflard-Fernandez (GWF)
+# code written by Clement Baruteau (CB), Sebastian Perez (SP), 
+# Marcelo Barraza (MB) and Valentin Oncle (VC), with substantial contributions 
+# from Simon Casassus (SC) and Gaylor Wafflard-Fernandez (GWF)
 #
 # =================================================================== 
 
 # =================================
 #            TO DO LIST
 # =================================
 # - add text banner with 'banner' command
-# - results from actual 3D simulations from Fargo3D (check
-# fargo2python): check latitude expression in mesh.py
 # - check again that all goes well without x-axisflip!
 # =================================
 
@@ -97,18 +95,19 @@
 if par.RTdust_or_gas == 'gas':
 
     # we start by computing gas temperature, density and velocity
-    if par.recalc_gas_quantities == 'Yes':
+    if (par.recalc_gas_quantities == 'Yes' or par.plot_gas_quantities == 'Yes'):
 
         from gas_density import *
         from gas_temperature import *
         from gas_velocity import *
-
-        print('--------- Setting gas temperature to temperature of hydro simulation ----------')
-        compute_hydro_temperature()
-        if par.plot_gas_quantities == 'Yes':
-            print('--------- Plotting gas temperature ----------')
-            plot_gas_temperature()
-
+
+        if (par.Tdust_eq_Thydro == 'Yes' or (par.Tdust_eq_Thydro == 'No' and par.Tdust_eq_Tgas == 'No') ):
+            print('--------- Computing temperature of hydro simulation ----------')
+            compute_hydro_temperature()
+            if par.plot_gas_quantities == 'Yes':
+                print('--------- Plotting hydro temperature ----------')
+                plot_gas_temperature()
+
         print('--------- Computing gas density ----------')
         compute_gas_mass_volume_density()
 
@@ -134,7 +133,7 @@
         tgas_eq_tdust = 0
 
     # start by computing gas temperature, density and velocity
-    if par.recalc_gas_quantities == 'Yes':
+    if (par.recalc_gas_quantities == 'Yes' or par.plot_gas_quantities == 'Yes'):
 
         from gas_density import *
         from gas_temperature import *

field.py

@@ -33,7 +33,16 @@ def __init__(self, field, staggered='c', directory='', dtype='float64'):
         nsec = int(nsec)
         self.nrad = nrad
         self.nsec = nsec
-        self.ncol = par.ncol              # colatitude (ncol is a global variable)
+        if par.hydro2D == 'Yes':
+            self.ncol = par.ncol              # colatitude (ncol is a global variable)
+        else:
+            command = par.awk_command+' " /^NZ/ " ' + directory + 'variables.par'
+            # check which version of python we're using
+            if sys.version_info[0] < 3:   # python 2.X
+                buf = subprocess.check_output(command, shell=True)
+            else:                         # python 3.X
+                buf = subprocess.getoutput(command)
+            self.ncol = int(buf.split()[1])
 
         Mesh.__init__(self, directory)    # all Mesh attributes inside Field
 
@@ -64,6 +73,7 @@ def __init__(self, field, staggered='c', directory='', dtype='float64'):
                 self.cumass = float(buf.split()[1])*2e30  #from Msol to kg        
                 # unit of temperature = mean molecular weight * 8.0841643e-15 * M / L;
                 self.cutemp = 2.35 * 8.0841643e-15 * self.cumass / self.culength
+                self.cutime = np.sqrt( self.culength**3 / 6.673e-11 / self.cumass)
 
         # override units used in Fargo simulations:
         else:
@@ -89,10 +99,14 @@ def __open_field(self, f, dtype):
         Reading the data
         """
         field = np.fromfile(f, dtype=dtype)
-        array = field.reshape(self.nrad,self.nsec) # 2D only!
-
-        # need to roll field by nsec/2 along azimuthal direction for FARGO3D runs!
-        if par.fargo3d == 'Yes':
-            array = np.roll(array, shift=int(self.nsec//2), axis=1)
+        if par.hydro2D == 'No': #3D simulation with fargo3d
+            array = field.reshape(self.ncol,self.nrad,self.nsec) # 3D
+            # need to roll field by nsec/2 along azimuthal direction for FARGO3D runs!
+            array = np.roll(array, shift=int(self.nsec//2), axis=2)
+        else:   # 2D simulation with dusty fargo adsg or fargo3d
+            array = field.reshape(self.nrad,self.nsec) # 2D
+            if par.fargo3d == 'Yes':
+                # need to roll field by nsec/2 along azimuthal direction for FARGO3D runs!
+                array = np.roll(array, shift=int(self.nsec//2), axis=1)
 
         return array