In [1]:
import yt
import os
import math
from yt import derived_field
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm
from yt.units import pc, kpc, second, Kelvin, gram, erg, cm
import matplotlib.patches as patches
from matplotlib.patches import Rectangle
import Reconstruct_UG as RUG

import copy
from matplotlib.colors import LogNorm

import timeit

start = timeit.default_timer()

get_ipython().magic('matplotlib inline')

# Define some constant parameters to be used.
mp      = 1.6726e-24  * gram # g
mu      = 1.2924
kb      = 1.3806e-16  *erg / Kelvin # erg K-1
GNewton = 6.6743e-8   * cm**3 / (gram * second**2 )# cm3 g-1 s-2
Msun    = 1.9884e33   * gram
mm      = mu*mp

ppc = 3.0856776e18

# Create a derived field.
@derived_field(name="numdens", units="1/cm**3", force_override=True)
def numdens(field, data):
    #mp      = 1.6726e-24 * gram
    #mu      = 1.2924
    #mm      = mu*mp
    dens_here = data["dens"].value
    dens_here = dens_here * gram / cm**3
    return dens_here/mm

yt.add_field('numdens', function=numdens, units="1/cm**3", force_override=True)
In [2]:
snapshot = 10
In [7]:
# Generalize for all the other clouds.

resolution = 0.06
pxM4 = 180
pyM4 = -25
pzM4 =  10
cloud_alphaM4 = 0.4

radM4 = 50

px_strM4 = str(pxM4)
py_strM4 = str(pyM4)
pz_strM4 = str(pzM4)

pzM4 = pzM4 - 3.

# Mean bulk velocity of the cloud.
bvx_meanM4 = -1.
bvy_meanM4 = -1.
bvz_meanM4 = -1.
    
txt_xoffsetM4 = 15
txt_yoffsetM4 = -25

cloud_dirM4 = "M4e3_" + "a%.2i" %(cloud_alphaM4*10) + "_x" + px_strM4 + "_y" + py_strM4 + "_z" + pz_strM4
data_dirM4  = "/data/manda/jcibanezm/StratBox/RealProd/1pc_and_AccClouds/AccPaper_Resims/Particles/" + cloud_dirM4 + "/006pc/"
#data_dirM4  = "/data/manda/jcibanezm/StratBox/RealProd/2pc/AccPaper_Resims/Particles/" + cloud_dirM4 + "/05pc/"

# Set the address and the name of the data set.
basenameM4 = "Strat_Box_M4e3_006pc_hdf5_"


######################################################################################

pxM3 = 450
pyM3 = -380
pzM3 = 25
cloud_alphaM3 = 0.4

radM3 = 50

px_strM3 = str(pxM3)
py_strM3 = str(pyM3)
pz_strM3 = str(pzM3)

pxM3 = pxM3 + 8
pyM3 = pyM3 
pzM3 = pzM3 - 8

# Mean bulk velocity of the cloud.
bvx_meanM3 = 0.
bvy_meanM3 = 3.0
bvz_meanM3 = -3.0
    
txt_xoffsetM3 = +30
txt_yoffsetM3 = -30


cloud_dirM3 = "M3e3_" + "a%.2i" %(cloud_alphaM3*10) + "_x" + px_strM3 + "_y" + py_strM3 + "_z" + pz_strM3
data_dirM3  = "/data/manda/jcibanezm/StratBox/RealProd/1pc_and_AccClouds/AccPaper_Resims/Particles/" + cloud_dirM3 + "/006pc/"

# Set the address and the name of the data set.
basenameM3 = "Strat_Box_M3e3_006pc_hdf5_"

######################################################################################

resolution = 0.1
pxM8 = 60
pyM8 = 370
pzM8 = 30
cloud_alphaM8 = 0.8

radM8 = 50

px_strM8 = str(pxM8)
py_strM8 = str(pyM8)
pz_strM8 = str(pzM8)

# Mean bulk velocity of the cloud.
bvx_meanM8 =  1.
bvy_meanM8 = -1.
bvz_meanM8 = -1.

pxM8 = pxM8 + 5
pyM8 = pyM8 - 10
pzM8 = pzM8 - 5
    
txt_xoffsetM8 = 30
txt_yoffsetM8 = -50


cloud_dirM8 = "M8e3_" + "a%.2i" %(cloud_alphaM8*10) + "_x" + px_strM8 + "_y" + py_strM8 + "_z" + pz_strM8
data_dirM8  = "/data/manda/jcibanezm/StratBox/RealProd/1pc_and_AccClouds/AccPaper_Resims/Particles/" + cloud_dirM8 + "/01pc/"

# Set the address and the name of the data set.
basenameM8 = "Strat_Box_M8e3_01pc_hdf5_"
In [8]:
plt_fileM3  = data_dirM3 + basenameM3 + "plt_cnt_%.4i" %snapshot
part_fileM3 = data_dirM3 + basenameM3 + "part_%.4i"    %snapshot

# Load data set
pfM3 = yt.load(plt_fileM3 , particle_filename=part_fileM3)

# Move the center of the frame.
px_nowM3 = pxM3 + bvx_meanM3 * snapshot/10.
py_nowM3 = pyM3 + bvy_meanM3 * snapshot/10.
pz_nowM3 = pzM3 + bvz_meanM3 * snapshot/10.

cM3 = [px_nowM3*ppc, py_nowM3*ppc, pz_nowM3*ppc]
#rad = 60

leM3 = [cM3[0] - radM3*ppc, cM3[1]-radM3*ppc, cM3[2]-radM3*ppc]
reM3 = [cM3[0] + radM3*ppc, cM3[1]+radM3*ppc, cM3[2]+radM3*ppc]

# Generate regions to access the particle and fluid data.
boxM3  = pfM3.region(cM3, leM3, reM3)
sphM3  = pfM3.sphere(cM3, radM3*pc)
yt : [INFO     ] 2017-01-17 12:08:08,432 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:08,432 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:08,433 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:08,465 Parameters: current_time              = 7.58294610627e+15
yt : [INFO     ] 2017-01-17 12:08:08,466 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-01-17 12:08:08,467 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:08,468 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:08,468 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-01-17 12:08:09,166 Loading field plugins.
yt : [INFO     ] 2017-01-17 12:08:09,168 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,169 Loaded astro (15 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,170 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,171 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,173 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,175 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,176 Loaded local (111 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,177 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,178 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:09,179 Loaded species (121 new fields)
In [9]:
plt_fileM4  = data_dirM4 + basenameM4 + "plt_cnt_%.4i" %snapshot
part_fileM4 = data_dirM4 + basenameM4 + "part_%.4i"    %snapshot

# Load data set
pfM4 = yt.load(plt_fileM4 , particle_filename=part_fileM4)

# Move the center of the frame.
px_nowM4 = pxM4 + bvx_meanM4 * snapshot/10.
py_nowM4 = pyM4 + bvy_meanM4 * snapshot/10.
pz_nowM4 = pzM4 + bvz_meanM4 * snapshot/10.

cM4 = [px_nowM4*ppc, py_nowM4*ppc, pz_nowM4*ppc]
#rad = 60

leM4 = [cM4[0] - radM4*ppc, cM4[1]-radM4*ppc, cM4[2]-radM4*ppc]
reM4 = [cM4[0] + radM4*ppc, cM4[1]+radM4*ppc, cM4[2]+radM4*ppc]

# Generate regions to access the particle and fluid data.
boxM4  = pfM4.region(cM4, leM4, reM4)
sphM4  = pfM4.sphere(cM4, radM4*pc)
yt : [INFO     ] 2017-01-17 12:08:17,648 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:17,649 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:17,650 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:17,667 Parameters: current_time              = 7.59874825137e+15
yt : [INFO     ] 2017-01-17 12:08:17,668 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-01-17 12:08:17,670 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:17,671 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:17,671 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-01-17 12:08:19,013 Loading field plugins.
yt : [INFO     ] 2017-01-17 12:08:19,017 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,018 Loaded astro (15 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,019 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,021 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,023 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,024 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,025 Loaded local (111 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,026 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,027 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:19,028 Loaded species (121 new fields)
In [10]:
plt_fileM8  = data_dirM8 + basenameM8 + "plt_cnt_%.4i" %snapshot
part_fileM8 = data_dirM8 + basenameM8 + "part_%.4i"    %snapshot

# Load data set
pfM8 = yt.load(plt_fileM8 , particle_filename=part_fileM8)

# Move the center of the frame.
px_nowM8 = pxM8 + bvx_meanM8 * snapshot/10.
py_nowM8 = pyM8 + bvy_meanM8 * snapshot/10.
pz_nowM8 = pzM8 + bvz_meanM8 * snapshot/10.

cM8 = [px_nowM8*ppc, py_nowM8*ppc, pz_nowM8*ppc]
#rad = 60

leM8 = [cM8[0] - radM8*ppc, cM8[1]-radM8*ppc, cM8[2]-radM8*ppc]
reM8 = [cM8[0] + radM8*ppc, cM8[1]+radM8*ppc, cM8[2]+radM8*ppc]

# Generate regions to access the particle and fluid data.
boxM8  = pfM8.region(cM8, leM8, reM8)
sphM8  = pfM8.sphere(cM8, radM8*pc)
yt : [INFO     ] 2017-01-17 12:08:21,235 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:21,236 integer runtime parameter forcedplotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:21,236 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:21,237 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:21,237 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:21,251 Parameters: current_time              = 7.58294474853e+15
yt : [INFO     ] 2017-01-17 12:08:21,252 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-01-17 12:08:21,253 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:21,253 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:21,254 Parameters: cosmological_simulation   = 0.0
yt : [INFO     ] 2017-01-17 12:08:22,436 Loading field plugins.
yt : [INFO     ] 2017-01-17 12:08:22,438 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,439 Loaded astro (15 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,440 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,442 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,444 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,445 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,446 Loaded local (111 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,447 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,448 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-01-17 12:08:22,449 Loaded species (121 new fields)
In [12]:
# Midplane
# I think Im loading a higher resolution dataset for this plot.

data_dir = "/data/manda/jcibanezm/StratBox/RealProd/1pc_and_AccClouds/AccPaper_Resims/Particles/Midplane"

#filename = data_dir + "/Strat_Box_05pc_hdf5_plt_cnt_0001"

filename = data_dir + "/Strat_Box_05pc_forced_hdf5_plt_cnt_0000"

pfMP = yt.load(filename)

c  = [40*ppc, 0, 0]
le = [-500*ppc, -500*ppc, -50*ppc]
re = [ 500*ppc,  500*ppc,  50*ppc]

boxMP = pfMP.region(c, le, re)
yt : [INFO     ] 2017-01-17 12:08:39,503 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:39,504 integer runtime parameter nbegin overwrites a simulation scalar of the same name
yt : [INFO     ] 2017-01-17 12:08:39,521 Parameters: current_time              = 7.41867344844e+15
yt : [INFO     ] 2017-01-17 12:08:39,521 Parameters: domain_dimensions         = [ 16  16 640]
yt : [INFO     ] 2017-01-17 12:08:39,522 Parameters: domain_left_edge          = [ -1.50000000e+21  -1.50000000e+21  -6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:39,523 Parameters: domain_right_edge         = [  1.50000000e+21   1.50000000e+21   6.00000000e+22]
yt : [INFO     ] 2017-01-17 12:08:39,523 Parameters: cosmological_simulation   = 0.0

In [13]:
prjMP = pfMP.proj("numdens", "z", center=c, data_source=boxMP)
frbMP = prjMP.to_frb(972*pc, 2048, center=c)["numdens"].value
yt : [INFO     ] 2017-01-17 12:09:09,629 Loading field plugins.
yt : [INFO     ] 2017-01-17 12:09:09,631 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,632 Loaded astro (15 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,632 Loaded cosmology (22 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,634 Loaded fluid (62 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,636 Loaded fluid_vector (94 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,637 Loaded geometric (110 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,638 Loaded local (111 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,639 Loaded magnetic_field (119 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,639 Loaded my_plugins (119 new fields)
yt : [INFO     ] 2017-01-17 12:09:09,640 Loaded species (121 new fields)
yt : [INFO     ] 2017-01-17 12:12:20,191 Projection completed
yt : [INFO     ] 2017-01-17 12:12:20,197 Making a fixed resolution buffer of (numdens) 2048 by 2048
In [14]:
prjM3 = pfM3.proj("numdens", "z", center=cM3, data_source=boxM3)
frbM3 = prjM3.to_frb(2*radM3*pc, 2048, center=cM3)["numdens"].value

prjM4 = pfM4.proj("numdens", "z", center=cM4, data_source=boxM4)
frbM4 = prjM4.to_frb(2*radM4*pc, 2048, center=cM4)["numdens"].value

prjM8 = pfM8.proj("numdens", "z", center=cM8, data_source=boxM8)
frbM8 = prjM8.to_frb(2*radM8*pc, 2048, center=cM8)["numdens"].value
yt : [INFO     ] 2017-01-17 12:12:24,229 Projection completed
yt : [INFO     ] 2017-01-17 12:12:24,233 Making a fixed resolution buffer of (numdens) 2048 by 2048
yt : [INFO     ] 2017-01-17 12:12:28,422 Projection completed
yt : [INFO     ] 2017-01-17 12:12:28,428 Making a fixed resolution buffer of (numdens) 2048 by 2048
yt : [INFO     ] 2017-01-17 12:12:32,067 Projection completed
yt : [INFO     ] 2017-01-17 12:12:32,072 Making a fixed resolution buffer of (numdens) 2048 by 2048
In [15]:
# Is it working ?

#part_range = [:]

pxM3 = np.array(boxM3["particle_posx"])/ppc
pyM3 = np.array(boxM3["particle_posy"])/ppc

pxM4 = np.array(boxM4["particle_posx"])/ppc
pyM4 = np.array(boxM4["particle_posy"])/ppc

pxM8 = np.array(boxM8["particle_posx"])/ppc
pyM8 = np.array(boxM8["particle_posy"])/ppc
In [16]:
# Get the cloud data, reconstruct it in a Uniform Grid and get the contour of the cloud.

#------------------------------------------------------------------------------------------------------------------
# Cloud M3
#------------------------------------------------------------------------------------------------------------------

UGM3 = RUG.Reconstruct_UG("M3e3", 0.06, snapshot)

NxUGM3 = np.shape(UGM3["x"])[0]
NyUGM3 = np.shape(UGM3["x"])[1]

xminUGM3, xmaxUGM3 = np.min(UGM3["x"]), np.max(UGM3["x"])
yminUGM3, ymaxUGM3 = np.min(UGM3["y"]), np.max(UGM3["y"])

surface_zM3 = np.zeros((NxUGM3, NyUGM3))

for i in range(NxUGM3):
    for j in range(NyUGM3):
        surface_zM3[i,j] = np.max(UGM3["numdens"][i, j, :])

del UGM3
del NxUGM3
del NyUGM3
  
#------------------------------------------------------------------------------------------------------------------
# Cloud M4
#------------------------------------------------------------------------------------------------------------------

UGM4 = RUG.Reconstruct_UG("M4e3", 0.06, snapshot)

NxUGM4 = np.shape(UGM4["x"])[0]
NyUGM4 = np.shape(UGM4["x"])[1]

xminUGM4, xmaxUGM4 = np.min(UGM4["x"]), np.max(UGM4["x"])
yminUGM4, ymaxUGM4 = np.min(UGM4["y"]), np.max(UGM4["y"])

surface_zM4 = np.zeros((NxUGM4, NyUGM4))

for i in range(NxUGM4):
    for j in range(NyUGM4):
        surface_zM4[i,j] = np.max(UGM4["numdens"][i, j, :])

del UGM4
del NxUGM4
del NyUGM4
        
#------------------------------------------------------------------------------------------------------------------
# Cloud M8
#------------------------------------------------------------------------------------------------------------------

UGM8 = RUG.Reconstruct_UG("M8e3", 0.1, snapshot)

NxUGM8 = np.shape(UGM8["x"])[0]
NyUGM8 = np.shape(UGM8["x"])[1]

xminUGM8, xmaxUGM8 = np.min(UGM8["x"]), np.max(UGM8["x"])
yminUGM8, ymaxUGM8 = np.min(UGM8["y"]), np.max(UGM8["y"])

surface_zM8 = np.zeros((NxUGM8, NyUGM8))

for i in range(NxUGM8):
    for j in range(NyUGM8):
        surface_zM8[i,j] = np.max(UGM8["numdens"][i, j, :])
        
del UGM8
del NxUGM8
del NyUGM8
Reading the File M3e3_CloudObject_snp010.dat
===============================================
Reconstructing a uniform grid for this cloud
Uniform Grid properties:
Lz = 37.98   Ly = 37.62   Lz = 30.26
dx = 0.06 pc
Nx = 641,    Ny = 636,    Nz = 512
xmin = 433.42      xmax = 471.40
ymin = -403.04      ymax = -365.42
zmin = 2.49      zmax = 32.76
===============================================
Reading the File M4e3_CloudObject_snp010.dat
===============================================
Reconstructing a uniform grid for this cloud
Uniform Grid properties:
Lz = 28.60   Ly = 24.21   Lz = 35.60
dx = 0.24 pc
Nx = 122,    Ny = 104,    Nz = 152
xmin = 167.81      xmax = 196.42
ymin = -39.64      ymax = -15.43
zmin = -12.58      zmax = 23.02
===============================================
Reading the File M8e3_CloudObject_snp010.dat
===============================================
Reconstructing a uniform grid for this cloud
Uniform Grid properties:
Lz = 31.69   Ly = 41.18   Lz = 28.96
dx = 0.12 pc
Nx = 269,    Ny = 349,    Nz = 246
xmin = 47.83      xmax = 79.52
ymin = 336.34      ymax = 377.52
zmin = 16.38      zmax = 45.34
===============================================
In [22]:
xsize = 10  # 18
ysize = 11 # 20

cmap = "RdYlBu_r"

part_size = 0.005
part_alpha = 0.1

fig = plt.figure(figsize=(xsize, ysize))

wi = 0.48
he = wi*xsize/ysize
ax = fig.add_axes([0.08, 0.15, wi, he])

extent = [-486, 486, -486, 486]
cax = ax.imshow(np.log10(frbMP), cmap=cmap, vmin=18, vmax=23, extent=extent, origin="lower")

ax.add_patch(patches.Rectangle((cM3[0]/ppc-radM3-40, cM3[1]/ppc-radM3), 2*radM3, 2*radM3, fill=False, linewidth=3))
ax.add_patch(patches.Rectangle((cM4[0]/ppc-radM4-40, cM4[1]/ppc-radM4), 2*radM4, 2*radM4, fill=False, linewidth=3))
ax.add_patch(patches.Rectangle((cM8[0]/ppc-radM8-40, cM8[1]/ppc-radM8), 2*radM8, 2*radM8, fill=False, linewidth=3))

# We change the fontsize of minor ticks label 
ax.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)


#ax.scatter(pxM3, pyM3, s=part_size, alpha=part_alpha, c="k")
#ax.scatter(pxM4, pyM4, s=part_size, alpha=part_alpha, c="k")
#ax.scatter(pxM8, pyM8, s=part_size, alpha=part_alpha, c="k")


for tick in ax.xaxis.get_major_ticks():
    tick.label.set_fontsize(13) 
for tick in ax.yaxis.get_major_ticks():
    tick.label.set_fontsize(13) 

ax.set_xlabel("position x [pc]", fontsize=15)
ax.set_ylabel("position y [pc]", fontsize=15, labelpad=0.0)

ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])

####################################################################

ax3 = fig.add_axes([0.0, 0.0, 1, 1])
ax3.set_axis_bgcolor('none')

ax3.axes.get_xaxis().set_visible(False)
ax3.axes.get_yaxis().set_visible(False)
ax3.axis('off')

ax3.set_xlim((0, 720))
ax3.set_ylim((0, 792))

#ax3.invert_xaxis()
ax3.invert_yaxis()

# M8e3
ax3.plot([352.8, 257.42],[21.6, 354.98], "-k", linewidth=2)
ax3.plot([57.6, 221.87],[21.6, 354.98], "-k", linewidth=2)
ax3.plot([57.6, 221.87],[316.8, 390.53], "-k", linewidth=2)
ax3.plot([352.8, 257.42],[316.8, 390.53], "-k", linewidth=2)

# M4e3
ax3.plot([417.6, 262.04],[21.6, 491.87], "-k", linewidth=2)
ax3.plot([297.6, 712.8], [491.87, 21.6], "-k", linewidth=2)
ax3.plot([417.6, 262.04],[316.8, 527.42], "-k", linewidth=2)
ax3.plot([297.6, 712.8], [527.42, 316.8], "-k", linewidth=2)

# M3e3
ax3.plot([361.24, 417.6],[616.67, 378], "-k", linewidth=2)
ax3.plot([396.8, 712.8], [616.67, 378], "-k", linewidth=2)
ax3.plot([361.24, 417.6],[652.22, 673.2], "-k", linewidth=2)
ax3.plot([396.8, 712.8], [652.22, 673.2], "-k", linewidth=2)


#############################################################################
#                        M3e3 Cloud
#############################################################################

wi = 0.41
he = wi*xsize/ysize
ax0 = fig.add_axes([0.58, 0.15, wi, he])

extent = [cM3[0]/ppc-radM3, cM3[0]/ppc+radM3, cM3[1]/ppc-radM3, cM3[1]/ppc+radM3]
ax0.imshow(np.log10(frbM3), cmap=cmap, vmin=18, vmax=23, extent=extent, origin="lower")

ax0.contour(np.transpose(surface_zM3),1, levels=[100], extent=[xminUGM3/ppc, xmaxUGM3/ppc, yminUGM3/ppc, ymaxUGM3/ppc], colors='k', linewidth=2)


#ax0.scatter(pxM3, pyM3, s=part_size, alpha=part_alpha, c="k")

ax0.xaxis.set_ticks_position("both")
ax0.yaxis.set_ticks_position("both")

ax0.axes.xaxis.set_ticklabels([])
ax0.axes.yaxis.set_ticklabels([])

ax0.set_xlim(extent[0], extent[1])
ax0.set_ylim(extent[2], extent[3])

for axis in ['top','bottom','left','right']:
  ax0.spines[axis].set_linewidth(2)

ax0.set_title("M3", fontsize=18)

# We change the fontsize of minor ticks label 
ax0.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax0.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)


#############################################################################
#                        M4e3 Cloud
#############################################################################

#wi = 0.3
#he = wi*xsize/ysize
ax1 = fig.add_axes([0.58, 0.6, wi, he])

extent = [cM4[0]/ppc-radM4, cM4[0]/ppc+radM4, cM4[1]/ppc-radM4, cM4[1]/ppc+radM4]
ax1.imshow(np.log10(frbM4), cmap=cmap, vmin=18, vmax=23, extent=extent, origin="lower")

ax1.contour(np.transpose(surface_zM4),1, levels=[100], extent=[xminUGM4/ppc, xmaxUGM4/ppc, yminUGM4/ppc, ymaxUGM4/ppc], colors='k', linewidth=2)

#ax1.scatter(pxM4, pyM4, s=part_size, alpha=part_alpha, c="k")

ax1.xaxis.set_ticks_position("both")
ax1.yaxis.set_ticks_position("both")

ax1.axes.xaxis.set_ticklabels([])
ax1.axes.yaxis.set_ticklabels([])

for axis in ['top','bottom','left','right']:
  ax1.spines[axis].set_linewidth(2)

ax1.set_xlim(extent[0], extent[1])
ax1.set_ylim(extent[2], extent[3])

ax1.set_title("M4", fontsize=18)

# We change the fontsize of minor ticks label 
ax1.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax1.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)


#############################################################################
#                        M8e3 Cloud
#############################################################################

#wi = 0.3
#he = wi*xsize/ysize
ax2 = fig.add_axes([0.08, 0.6, wi, he])

extent = [cM8[0]/ppc-radM8, cM8[0]/ppc+radM8, cM8[1]/ppc-radM8, cM8[1]/ppc+radM8]
ax2.imshow(np.log10(frbM8), cmap=cmap, vmin=18, vmax=23, extent=extent, origin="lower")

ax2.contour(np.transpose(surface_zM8),1, levels=[100], extent=[xminUGM8/ppc, xmaxUGM8/ppc, yminUGM8/ppc, ymaxUGM8/ppc], colors='k', linewidth=2)

#ax2.scatter(pxM8, pyM8, s=part_size, alpha=part_alpha, c="k")

for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(13) 
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(13) 

ax2.xaxis.set_ticks_position("both")
ax2.yaxis.set_ticks_position("both")

ax2.axes.xaxis.set_ticklabels([])
ax2.axes.yaxis.set_ticklabels([])

for axis in ['top','bottom','left','right']:
  ax2.spines[axis].set_linewidth(2)

ax2.set_xlim(extent[0], extent[1])
ax2.set_ylim(extent[2], extent[3])

ax2.set_title("M8", fontsize=18)

# We change the fontsize of minor ticks label 
ax2.tick_params(axis='both', which='major', length=10, width=2, labelsize=20)
ax2.tick_params(axis='both', which='minor', length=5, width=1.5, labelsize=12)



# ------------------------------------------------------------------------------------    
# Colorbar
cbar_ax = fig.add_axes([0.03, 0.055, 0.95, 0.04])
cbar = fig.colorbar(cax, cax=cbar_ax, orientation='horizontal')

cbar.set_label("Column Density [cm$^{-2}$]", fontsize=17)
cbar.ax.tick_params(labelsize=15) 
cbar.ax.tick_params('both', length=10, width=1.5, which='major')
#cbar.ax.tick_params('both', length=10, width=1.5, which='minor')
cbar.ax.xaxis.set_ticks_position('bottom')
cbar.set_ticks([16, 17, 18, 19, 20, 21, 22, 23])
cbar.set_ticklabels(["10$^{16}$", "10$^{17}$", "10$^{18}$", "10$^{19}$", "10$^{20}$", "10$^{21}$", "10$^{22}$", "10$^{23}$"])
cbar.ax.xaxis.set_label_position("bottom")



fig.show()

#fig.savefig("/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/Figure1.pdf", format="pdf", dpi=200)
fig.savefig("/home/jcibanezm/codes/StratBox/AccretingClouds_Paper/Plots/Figure1_new.pdf", format="pdf", dpi=100)