" Format x into nice Latex rounding to n"

power = np.floor(np.log10(round_to_n(x, 0)))

f_SF = round_to_n(x, n) * pow(10, -power)

if f_SF != 1.0:

mystr = "$"+format(f_SF,'.1f')+r"\times 10^{"+format(power,'.0f')+"}$"

return mystr

else:

mystr = r"$10^{"+format(power,'.0f')+"}$"

# first import global variables

import par

colored_cmap=par.mycolormap

# how many color bars

two_color_bars = False

if par.showdust == 'Yes':

two_color_bars = True

# several output numbers

if isinstance(par.on, int) == True:

on = [par.on]

else:

on = par.on

# if several directories specified in parameter file, only deal with the first one

if isinstance(par.directory, str) == False:

directory = par.directory[0]

else:

directory = par.directory

# if movie requested, then compute fields for all output numbers between range

# specified in parameter file via on (onmin,onmax)

if par.movie == 'Yes':

on = range(par.on[0],par.on[1]+1,par.take_one_point_every)

for k in range(len(on)):

if par.movie == 'Yes':

print('animation: output number '+str(k)+' / '+str(len(on)-1),end='\r')

if par.allfluids == 'Yes':

# prepare figure for each output number

if par.nbfluids > 6: # could be refined...

fig, axes = plt.subplots(2,int(np.ceil(par.nbfluids/2)), figsize=(24.,12.))

else:

fig, axes = plt.subplots(2,int(np.ceil(par.nbfluids/2)), figsize=(18.,12.))

if par.fieldofview == 'cart':

plt.subplots_adjust(left=0.25, right=0.75, top=0.90, bottom=0.10, wspace=0.6, hspace=0.6)

else:

plt.subplots_adjust(left=0.05, bottom=0.08, right=0.95, top=0.92, wspace=0.4, hspace=0.4)

# loop over fluids

for f in range(len(par.fluids)):

if par.allfluids == 'Yes':

if f < int(np.ceil(par.nbfluids/2)):

rowindex = 0

colindex = f

else:

rowindex = 1

colindex = f-int(np.ceil(par.nbfluids/2))

ax = axes[rowindex,colindex]

# case we have an odd number of 'windows' to be displayed

if len(par.fluids) % 2 != 0 and f == len(par.fluids)-1:

fig.delaxes(axes[1,int(np.ceil(par.nbfluids/2))-1])

myfield = Field(field=par.whatfield, fluid=par.fluids[f], on=on[k], directory=directory, physical_units=par.physical_units, nodiff=par.nodiff, fieldofview=par.fieldofview, slice=par.slice, z_average=par.z_average, onedprofile='No', override_units=par.override_units)

# case we plot the relative difference of a same field between two different directories

if ('subtract_directory' in open('paramsf2p.dat').read()) and (par.subtract_directory != '#'):

print('subtract_directory set in params.dat file: I understand you want to compare a same field from two different directories')

myfield2 = Field(field=par.whatfield, fluid=par.fluids[f], on=on[k], directory=par.subtract_directory, physical_units=par.physical_units, nodiff=par.nodiff, fieldofview=par.fieldofview, slice=par.slice, z_average=par.z_average, onedprofile='No', override_units=par.override_units)

array2 = myfield2.data

array = (myfield.data-myfield2.data)/myfield2.data

par.log_colorscale = 'No'

else:

array = myfield.data

# plot relative difference wrt initial field

if par.nodiff == 'No':

myfield0 = Field(field=par.whatfield, fluid=par.fluids[f], on=0, directory=directory, physical_units=par.physical_units, nodiff=par.nodiff, fieldofview=par.fieldofview, slice=par.slice, z_average=par.z_average, onedprofile='No', override_units=par.override_units)

array0 = myfield0.data

'''

array = (myfield.data-array0)/array0

else:

if par.nodiff != 'normnao':

# conversion in physical units

if par.physical_units == 'Yes':

array *= myfield.unit

# rotate field by user-defined angle specified in degree in .dat file

if ('rotate_angle' in open('paramsf2p.dat').read()) and (par.rotate_angle != '#'):

shift = int(myfield.nsec*par.rotate_angle/360.0)

array = np.roll(array,shift=shift,axis=1)

# roll array by nsec/2 along azimuthal axis if grid's

# azimuthal extent is smaller than 2pi:

if abs(myfield.pedge.max()-2.0*np.pi) > 0.1:

array = np.roll(array,shift=myfield.nsec//2,axis=1)

# get field's full name

strfield = myfield.strname

# -------------------

# Stuff we do only at first output number

# -------------------

if k == 0 and f==0:

# get computational grid coordinates

# start with cylindrical radius (Fargo2D) or spherical radius (Fargo3D)

R = myfield.redge # myfield.rmed

# radius conversion in physical units

if par.physical_units == 'Yes':

R *= (myfield.culength / 1.5e11) # in au

# define visualisation params

if (('rbox' in open('paramsf2p.dat').read()) and (par.rbox != '#')) or (('box_around_planet' in open('paramsf2p.dat').read()) and (par.box_around_planet != '#')):

planet_to_track = 0

if ('planet_to_track' in open('paramsf2p.dat').read()) and (par.planet_to_track != '#'):

planet_to_track = par.planet_to_track

# start by finding planet's orbita radius

if myfield.fargo3d == 'Yes':

f1, xpla, ypla, f4, f5, f6, f7, f8, date, omega = np.loadtxt(directory+"/planet"+str(planet_to_track)+".dat",unpack=True)

else:

f1, xpla, ypla, f4, f5, f6, f7, date, omega, f10, f11 = np.loadtxt(directory+"/planet"+str(planet_to_track)+".dat",unpack=True)

if par.physical_units == 'Yes':

xpla *= (myfield.culength / 1.5e11) # in au

ypla *= (myfield.culength / 1.5e11) # in au

else:

if (par.myrmin != '#'):

myrmin = par.myrmin

else:

myrmin = R.min()

if (par.myrmax != '#'):

myrmax = par.myrmax

else:

myrmax = R.max()

imin = np.argmin(np.abs(R-myrmin))

imax = np.argmin(np.abs(R-myrmax))

# VISUALISATION IN MIDPLANE: get azimuth

if par.fargo3d == 'No' or (par.fieldofview != 'latitudinal' and par.fieldofview != 'vertical'):

T = myfield.pedge # myfield.pmed

myphimin = par.myphimin

if (par.myphimin == '#'):

myphimin = T.min()

jmin = np.argmin(np.abs(T-myphimin))

myphimax = par.myphimax

if (myphimax == '#'):

myphimax = T.max()

jmax = np.argmin(np.abs(T-myphimax))

if ('flip_xaxis' in open('paramsf2p.dat').read()) and (par.flip_xaxis == 'Yes'):

bufmin = myphimin

bufmax = myphimax

myphimax = bufmin

myphimin = bufmax

# VISUALISATION IN VERTICAL (LATITUDINAL) PLANE

if par.fargo3d == 'Yes' and par.fieldofview == 'latitudinal':

#T = myfield.tmed # latitude

T = myfield.tedge # latitude

# number of grid cells in the radial and azimuthal directions

jmin = np.argmin(np.abs(T-T.min()))

jmax = np.argmin(np.abs(T-T.max()))

if par.fargo3d == 'Yes' and par.fieldofview == 'vertical':

T = myfield.zedge # altitude above midplane

# number of grid cells in the radial and azimuthal directions

jmin = np.argmin(np.abs(T-T.min()))

jmax = np.argmin(np.abs(T-T.max()))

# color bar for background field

mycolormap = colored_cmap

'''

# end of stuff done only at first output number

# -------------------

# case where fields are displayed with a fixed radial range about planet's orbital radius

# -------------------

if ('rbox' in open('paramsf2p.dat').read()) and (par.rbox != '#'):

rpla = np.sqrt( xpla[int(on[k])]*xpla[int(on[k])] + ypla[int(on[k])]*ypla[int(on[k])] )

myrmin = rpla-par.rbox

myrmax = rpla+par.rbox

imin = np.argmin(np.abs(R-myrmin))

imax = np.argmin(np.abs(R-myrmax))

# -------------------

# case where fields are displayed over a squared box around one (tracked) planet

# -------------------

if ('box_around_planet' in open('paramsf2p.dat').read()) and (par.box_around_planet != '#'):

rpla = np.sqrt( xpla[int(on[k])]*xpla[int(on[k])] + ypla[int(on[k])]*ypla[int(on[k])] )

myrmin = rpla-par.box_around_planet

myrmax = rpla+par.box_around_planet

imin = np.argmin(np.abs(R-myrmin))

imax = np.argmin(np.abs(R-myrmax))

tpla = math.atan2(ypla[int(on[k])],xpla[int(on[k])]) # between 0 and 2pi

if par.fieldofview == 'polar':

rotate_angle_box = -tpla

shift_box = int(myfield.nsec*tpla/2.0/np.pi)

array = np.roll(array,shift=-shift_box,axis=1)

myphimin = np.pi-2.0*par.box_around_planet

myphimax = np.pi+2.0*par.box_around_planet

else:

myphimin = tpla-2.0*par.box_around_planet

myphimax = tpla+2.0*par.box_around_planet

# -------------------

# read information on the dust particles

# -------------------

if par.showdust == 'Yes':

(rd, td, vrd, vtd, Stokes, sizedust) = np.loadtxt(directory+'/dustsystat'+str(on[k])+'.dat', unpack=True)

if par.physical_units == 'Yes':

rd *= (myfield.culength / 1.5e11) # in au

sizemin = par.sizemin

if (par.sizemin == '#'):

sizemin = sizedust.min()

sizemax = par.sizemax

if (par.sizemax == '#'):

sizemax = sizedust.max()

# rotate particles azimuth by user-defined angle specified in degree in .dat file

if ('rotate_angle' in open('paramsf2p.dat').read()) and (par.rotate_angle != '#'):

pc_shift_angle = np.pi*par.rotate_angle/180 # in radian

for i in range(len(td)):

td[i] += pc_shift_angle

if td[i] > 2.0*np.pi:

td[i] -= 2.0*np.pi

if td[i] < 0.0:

td[i] += 2.0*np.pi

# Only select particles with size between sizemin and sizemax

mydustsize = sizedust.compress((sizedust<sizemax)&(sizedust>sizemin).flat)

mydustsize2 = mydustsize[mydustsize.ravel().argsort()]

myrd = rd.compress((sizedust<sizemax)&(sizedust>sizemin).flat)

myrd2 = myrd[mydustsize.ravel().argsort()]

mytd = td.compress((sizedust<sizemax)&(sizedust>sizemin).flat)

mytd2 = mytd[mydustsize.ravel().argsort()]

rd = myrd2

td = mytd2

sizedust = mydustsize2

# -------------------

# POLAR FIELD OF VIEW

# -------------------

if par.fieldofview == 'polar':

# shift by nsec/2 along azimuthal direction, same for dust

if par.fargo3d == 'No':

array = np.roll(array, shift=int(myfield.nsec//2), axis=1)

#array = np.roll(array, shift=int(myfield.nsec//2), axis=1)

array_orig = array

if par.showdust == 'Yes':

td += np.pi

for i in range(len(td)):

if td[i] > 2.0*np.pi:

td[i] -= 2.0*np.pi

# plot radius in y-axis, azimuth in x-axis

if par.rvsphi == 'Yes':

X = T

Y = R

if par.showdust == 'Yes': # particles

xdust = td

ydust = rd

else:

# plot azimuth in y-axis, radius in x-axis

array = np.transpose(array)

X = R

Y = T

if par.showdust == 'Yes': # particles

xdust = rd

ydust = td

# figure

if par.allfluids == 'No':

if two_color_bars == True:

fig = plt.figure(figsize=(8.5,8.))

plt.subplots_adjust(left=0.12, right=0.86, top=0.88, bottom=0.11)

else:

fig = plt.figure(figsize=(8.,8.))

plt.subplots_adjust(left=0.16, right=0.94, top=0.88, bottom=0.11)

ax = plt.gca()

if par.physical_units == 'No':

ax.xaxis.set_minor_locator(MultipleLocator(0.1))

ax.yaxis.set_minor_locator(MultipleLocator(0.1))

if par.rvsphi == 'Yes':

ax.set_xlabel('Azimuth [rad]')

if par.physical_units == 'Yes':

ax.set_ylabel('Radius [au]')

else:

ax.set_ylabel(r'Radius [$R_0$]')

xlim_min = myphimin

xlim_max = myphimax

ylim_min = myrmin

ylim_max = myrmax

#ax.set_xlim(myphimin,myphimax)

#ax.set_ylim(myrmin,myrmax)

else:

ax.set_ylabel('Azimuth [rad]')

if par.physical_units == 'Yes':

ax.set_xlabel('Radius [au]')

else:

ax.set_xlabel(r'Radius [$R_0$]')

xlim_min = myrmin

xlim_max = myrmax

ylim_min = myphimin

ylim_max = myphimax

#ax.set_ylim(myphimin,myphimax)

#ax.set_xlim(myrmin,myrmax)

#

# -------------------------

# LATITUDINAL FIELD OF VIEW

# -------------------------

if par.fieldofview == 'latitudinal':

X = R

Y = T

array_orig = array

array = np.transpose(array)

if par.allfluids == 'No':

if two_color_bars == True:

fig = plt.figure(figsize=(8.5,8.))

plt.subplots_adjust(left=0.19, right=0.86, top=0.88, bottom=0.11)

else:

fig = plt.figure(figsize=(8.,8.))

plt.subplots_adjust(left=0.19, right=0.94, top=0.88, bottom=0.11)

ax = plt.gca()

if par.physical_units == 'No':

ax.xaxis.set_minor_locator(MultipleLocator(0.1))

ax.yaxis.set_minor_locator(MultipleLocator(0.1))

ax.set_ylabel('Latitude [rad]')

ax.set_xlabel(r'Radius [$R_0$]')

else:

ax.set_xlabel('Radius [au]')

ax.set_ylabel('Latitude [rad]')

xlim_min = myrmin

xlim_max = myrmax

ylim_min = T.min()

ylim_max = T.max()

#ax.set_ylim(T.min(),T.max())

#ax.set_xlim(myrmin,myrmax)

#

# -----------------------

# CARTESIAN FIELD OF VIEW

# -----------------------

if par.fieldofview == 'cart':

if par.fargo3d == 'Yes' and myfield.cartesian_grid == 'No':

array = np.roll(array, shift=int(myfield.nsec//2), axis=1)

array_orig = array

if myfield.cartesian_grid == 'No':

radius_matrix, theta_matrix = np.meshgrid(R,myfield.pedge)

X = radius_matrix * np.cos(theta_matrix)

if ('flip_xaxis' in open('paramsf2p.dat').read()) and (par.flip_xaxis == 'Yes'):

X = -X

Y = radius_matrix * np.sin(theta_matrix)

# roll array by nsec/2 along azimuthal axis if grid's

# azimuthal extent is smaller than 2pi:

if abs(myfield.pedge.max()-2.0*np.pi) > 0.1:

array = np.roll(array,shift=myfield.nsec//2,axis=1)

#

array = np.transpose(array)

#

if par.showdust == 'Yes': # particles

xdust = rd*np.cos(td)

if ('flip_xaxis' in open('paramsf2p.dat').read()) and (par.flip_xaxis == 'Yes'):

xdust = -xdust

ydust = rd*np.sin(td)

if myfield.cartesian_grid == 'Yes':

X = myfield.xmed

Y = myfield.ymed

# figure

if par.allfluids == 'No':

if two_color_bars == True:

fig = plt.figure(figsize=(8.8,8.))

# do not edit subplot position below!

plt.subplots_adjust(left=0.14, right=0.85, top=0.88, bottom=0.1)

else:

fig = plt.figure(figsize=(8.,8.))

# do not edit subplot position below!

plt.subplots_adjust(left=0.18, right=0.93, top=0.89, bottom=0.11)

ax = plt.gca()

xlim_min = -myrmax

xlim_max = myrmax

ylim_min = -myrmax

ylim_max = myrmax

#ax.set_xlim(-myrmax,myrmax)

#ax.set_ylim(-myrmax,myrmax)

# special case where a disc wedge is simulated (=not full 2pi in azimuth)

if np.abs(myfield.pmed[-1]-myfield.pmed[0]) < 3.2:

ylim_min = 0.0

fig = plt.figure(figsize=(16.,8.))

# do not edit subplot position below!

plt.subplots_adjust(left=0.17, right=0.92, top=0.88, bottom=0.1)

ax = plt.gca()

if par.physical_units == 'Yes':

ax.set_xlabel('x [au]')

ax.set_ylabel('y [au]')

else:

ax.set_xlabel(r'x [$R_0$]')

ax.set_ylabel(r'y [$R_0$]')

#

# -----------------------

# VERTICAL FIELD OF VIEW

# -----------------------

if par.fieldofview == 'vertical':

array_orig = array

array = np.transpose(array)

if myfield.cylindrical_grid == 'No':

radius_matrix, theta_matrix = np.meshgrid(R,myfield.tedge)

#radius_matrix, theta_matrix = np.meshgrid(myfield.rmed,myfield.tmed)

X = radius_matrix * np.cos(theta_matrix)

Y = radius_matrix * np.sin(theta_matrix)

else:

X = R

Y = T

# figure

if par.allfluids == 'No':

if two_color_bars == True:

fig = plt.figure(figsize=(8.5,8.))

plt.subplots_adjust(left=0.15, right=0.86, top=0.88, bottom=0.11)

else:

fig = plt.figure(figsize=(8.,8.))

plt.subplots_adjust(left=0.18, right=0.96, top=0.88, bottom=0.11)

ax = plt.gca()

xlim_min = myrmin #X.min()

xlim_max = myrmax #X.max()

if (par.myzmin != '#'):

ylim_min = par.myzmin

else:

ylim_min = Y.min()

if (par.myzmax != '#'):

ylim_max = par.myzmax

else:

ylim_max = Y.max()

if par.physical_units == 'No':

ax.xaxis.set_minor_locator(MultipleLocator(0.1))

ax.yaxis.set_minor_locator(MultipleLocator(0.1))

ax.set_ylabel(r'Altitude [$R_0$]')

ax.set_xlabel(r'Radius [$R_0$]')

else:

ax.set_xlabel('Radius [au]')

ax.set_ylabel('Altitude [au]')

# common display

ax.set_xlim(xlim_min,xlim_max)

ax.set_ylim(ylim_min,ylim_max)

ax.tick_params(top='on', right='on', length = 5, width=1.0, direction='out')

ax.tick_params(axis='x', which='minor', top=True)

ax.tick_params(axis='y', which='minor', right=True)

#ax.xaxis.set_major_locator(plt.MaxNLocator(7))

#ax.yaxis.set_major_locator(plt.MaxNLocator(7))

ax.xaxis.set_major_locator(plt.MaxNLocator(5))

ax.yaxis.set_major_locator(plt.MaxNLocator(5))

# -------------------

# adjust min and max in color bar according to field of view

# -------------------

#print(array_orig.shape,imin,imax,jmin,jmax)

myfieldmin = par.fieldmin

if (par.fieldmin == '#'):

myfieldmin = array_orig[imin:imax+1,jmin:jmax+1].min()

if par.fieldofview == 'latitudinal':

myfieldmin = array_orig[imin:imax+1,jmin:jmax+1].min()

myfieldmax = par.fieldmax

if (par.fieldmax == '#'):

myfieldmax = array_orig[imin:imax+1,jmin:jmax+1].max()

if par.fieldofview == 'latitudinal':

myfieldmax = array_orig[imin:imax+1,jmin:jmax+1].max()

if ( ('symmetrical_colorscale' in open('paramsf2p.dat').read()) and (par.symmetrical_colorscale == 'Yes') ):

if myfieldmax > np.abs(myfieldmin):

myfieldmin = -myfieldmax

else:

myfieldmax = -myfieldmin

if par.log_colorscale == 'Yes':

if (par.fieldmin == 'auto' or par.fieldmax == 'auto'):

minarray = array.min() #1e-3*array.max()?

maxarray = array.max()

array = np.log(array/minarray)/np.log(maxarray/minarray)

#print(par.fluids[f],minarray,maxarray,array.min(),array.max())

#print(array)

myfieldmin = 1e-5

myfieldmax = 1.0

mynorm = matplotlib.colors.Normalize(vmin=myfieldmin,vmax=myfieldmax)

strfield = r' $\log(\rho/\rho_{min})~/~\log(\rho_{max}/\rho_{min})$'+' at '+myfield.strtime

else:

if (myfieldmax/myfieldmin > 1e3 and (par.fieldmin == '#')):

myfieldmin = 1e-3*myfieldmax

# Normalization for colorbar: linear or logarithmic scale

if par.log_colorscale == 'Yes':

if (par.fieldmin != 'auto' and par.fieldmax != 'auto'):

mynorm = matplotlib.colors.LogNorm(vmin=myfieldmin,vmax=myfieldmax)

else:

mynorm = matplotlib.colors.Normalize(vmin=myfieldmin,vmax=myfieldmax)

# -----------------------

# display contour field

# -----------------------

CF = ax.pcolormesh(X,Y,array,cmap=mycolormap,norm=mynorm,rasterized=True)

#CF = ax.imshow(array, origin='lower', cmap=mycolormap, interpolation='bilinear', vmin=myfieldmin, vmax=myfieldmax, aspect='auto', extent=[X.min(),X.max(),Y.min(),Y.max()])

# ------------------

# overlay streamlines

# ------------------

if par.streamlines == 'Yes':

for s in range(par.nstreamlines):

# initial radius and azimuth of streamlines

#myR0 = myrmin + np.random.rand()*(myrmax-myrmin)

myR0 = myrmin + s*(myrmax-myrmin)/par.nstreamlines

myT0 = myphimin + np.random.rand()*(myphimax-myphimin)

# forward integration of streamlines

xstr,ystr = myfield.compute_streamline(niterations=10000,R0=myR0,T0=myT0,rmin=myrmin,rmax=myrmax,pmin=myphimin,pmax=myphimax,forward=True,fieldofview=par.fieldofview,slice=par.slice, z_average=par.z_average)

ax.scatter(xstr,ystr,s=1,marker='.',color='white')

# backward integration of streamlines

xstr,ystr = myfield.compute_streamline(niterations=10000,R0=myR0,T0=myT0,rmin=myrmin,rmax=myrmax,pmin=myphimin,pmax=myphimax,forward=False,fieldofview=par.fieldofview,slice=par.slice, z_average=par.z_average)

ax.scatter(xstr,ystr,s=1,marker='.',color='white')

# ------------------

# overlay particles

# ------------------

if par.showdust == 'Yes':

mynorm = matplotlib.colors.LogNorm(vmin=sizemin,vmax=sizemax)

CD = ax.scatter(xdust,ydust,s=1,c=sizedust,cmap='nipy_spectral',alpha=0.3,norm=mynorm)

#CD = ax.scatter(xdust,ydust,s=1,c=sizedust,cmap='nipy_spectral',alpha=0.3,vmin=sizemin,vmax=sizemax,norm=matplotlib.colors.LogNorm())

# ------------------

# overlay planets

# ------------------

if par.showplanet == 'Yes':

# Find how many 'planets' there are

nbplanets = len(fnmatch.filter(os.listdir(directory), 'planet*.dat'))

xp = np.zeros(nbplanets)

yp = np.zeros(nbplanets)

for l in range(nbplanets):

# read information on the planets (inner one so far)

if par.fargo3d == 'Yes':

f1, xplabis, yplabis, f4, f5, f6, f7, f8, date, omega = np.loadtxt(directory+"/planet"+str(l)+".dat",unpack=True)

else:

if par.fargo_orig == 'No':

f1, xplabis, yplabis, f4, f5, f6, f7, date, omega, f10, f11 = np.loadtxt(directory+"/planet"+str(l)+".dat",unpack=True)

else:

f1, xplabis, yplabis, f4, f5, f6, f7, date, omega = np.loadtxt(directory+"/planet"+str(l)+".dat",unpack=True)

xp[l] = xplabis[on[k]]

yp[l] = yplabis[on[k]]

if par.physical_units == 'Yes':

xp[l] *= (myfield.culength / 1.5e11) # in au

yp[l] *= (myfield.culength / 1.5e11) # in au

rp = np.sqrt(xp[l]*xp[l] + yp[l]*yp[l]) # planet's orbital radius

tp = math.atan2(yp[l],xp[l]) # planet's azimuthal angle

# rotate particles azimuth by user-defined angle specified in degree in .dat file

if ('rotate_angle' in open('paramsf2p.dat').read()) and (par.rotate_angle != '#'):

planet_shift_angle = np.pi*par.rotate_angle/180 # in radian

tp += planet_shift_angle

if tp > 2.0*np.pi:

tp -= 2.0*np.pi

if tp < 0.0:

tp += 2.0*np.pi

xp[l] = rp*np.cos(tp)

yp[l] = rp*np.sin(tp)

if ('box_around_planet' in open('paramsf2p.dat').read()) and (par.box_around_planet != '#'):

planet_shift_angle = rotate_angle_box

tp += planet_shift_angle

if tp > 2.0*np.pi:

tp -= 2.0*np.pi

if tp < 0.0:

tp += 2.0*np.pi

xp[l] = rp*np.cos(tp)

yp[l] = rp*np.sin(tp)

if par.fieldofview == 'polar':

tp += np.pi

if tp > 2.0*np.pi:

tp -= 2.0*np.pi

if par.rvsphi == 'Yes':

xp[l] = tp

yp[l] = rp

else:

xp[l] = rp

yp[l] = tp

if par.fieldofview == 'cart' and myfield.cartesian_grid == 'No':

if ('flip_xaxis' in open('paramsf2p.dat').read()) and (par.flip_xaxis == 'Yes'):

xp[l] = -xp[l]

# add planets via scatter plot

if par.verbose == 'Yes':

print('xp = ', xp, ' , yp = ', yp)

CP = ax.scatter(xp,yp,s=10,c='lightpink',cmap=colored_cmap,alpha=1)

# ------------------

# overlay CPU

# ------------------

if ('showcpus' in open('paramsf2p.dat').read()) and par.showcpus == 'Yes':

if par.fargo3d == 'No':

cpunb, cpurmin, cpurmax = np.loadtxt(directory+"/minmaxradii.dat",unpack=True)

if par.physical_units == 'Yes':

cpurmin *= (myfield.culength / 1.5e11) # in au

cpurmax *= (myfield.culength / 1.5e11) # in au

for i in range(len(cpurmin)):

if par.fieldofview == 'polar':

if par.rvsphi == 'No':

ax.plot([cpurmin[i],cpurmin[i]],[Y.min(),Y.max()],'-',linewidth=1,color='grey')

else:

ax.plot([X.min(),X.max()],[cpurmin[i],cpurmin[i]],'-',linewidth=1,color='grey')

# ------------------

# overlay wave-killing zones

# ------------------

if ('showwkzones' in open('paramsf2p.dat').read()) and par.showwkzones == 'Yes':

if par.fargo3d == 'No':

buf = subprocess.getoutput(par.awk_command+' " /^Rmin / " '+directory+'/*.par')

grid_rmin = float(buf.split()[1])

buf = subprocess.getoutput(par.awk_command+' " /^WKZRmin/ " '+directory+'/*.par')

grid_wkzrmin = float(buf.split()[1])

buf = subprocess.getoutput(par.awk_command+' " /^Rmax / " '+directory+'/*.par')

grid_rmax = float(buf.split()[1])

buf = subprocess.getoutput(par.awk_command+' " /^WKZRmax/ " '+directory+'/*.par')

grid_wkzrmax = float(buf.split()[1])

else:

command = par.awk_command+' " /^YMIN/ " '+directory+'/variables.par'

if sys.version_info[0] < 3: # python 2.X

buf = subprocess.check_output(command, shell=True)

else: # python 3.X

buf = subprocess.getoutput(command)

grid_rmin = float(buf.split()[1])

command = par.awk_command+' " /^YMAX/ " '+directory+'/variables.par'

if sys.version_info[0] < 3: # python 2.X

buf = subprocess.check_output(command, shell=True)

else: # python 3.X

buf = subprocess.getoutput(command)

grid_rmax = float(buf.split()[1])

command = par.awk_command+' " /^DAMPINGZONE/ " '+directory+'/variables.par'

if sys.version_info[0] < 3: # python 2.X

buf = subprocess.check_output(command, shell=True)

else: # python 3.X

buf = subprocess.getoutput(command)

dampingzone = float(buf.split()[1])

grid_wkzrmin = grid_rmin*(dampingzone**(2./3))

grid_wkzrmax = grid_rmax*(dampingzone**(-2./3))

if par.physical_units == 'Yes':

grid_rmin *= (myfield.culength / 1.5e11) # in au

grid_wkzrmin *= (myfield.culength / 1.5e11) # in au

grid_rmax *= (myfield.culength / 1.5e11) # in au

grid_wkzrmax *= (myfield.culength / 1.5e11) # in au

from matplotlib.patches import Ellipse, Polygon, Rectangle

if par.fieldofview == 'polar':

if par.rvsphi == 'No':

ax.plot([grid_wkzrmin,grid_wkzrmin],[Y.min(),Y.max()],'-',linewidth=1,color='grey')

ax.fill_between([grid_rmin,grid_wkzrmin],[Y.min(),Y.min()], [Y.max(),Y.max()], hatch="/", linewidth=0.0, alpha=0.0, color='grey')

ax.plot([grid_wkzrmax,grid_wkzrmax],[Y.min(),Y.max()],'-',linewidth=1,color='grey')

ax.fill_between([grid_wkzrmax,grid_rmax],[Y.min(),Y.min()], [Y.max(),Y.max()], hatch="/", linewidth=0.0, alpha=0.0, color='grey')

else:

ax.plot([X.min(),X.max()],[grid_wkzrmin,grid_wkzrmin],'-',linewidth=1,color='grey')

ax.fill_between([X.min(),X.max()],[grid_rmin,grid_rmin], [grid_wkzrmin,grid_wkzrmin], hatch="/", linewidth=0.0, alpha=0.0, color='grey')

ax.plot([X.min(),X.max()],[grid_wkzrmax,grid_wkzrmax],'-',linewidth=1,color='grey')

ax.fill_between([X.min(),X.max()],[grid_wkzrmax,grid_wkzrmax], [grid_rmax,grid_rmax], hatch="/", linewidth=0.0, alpha=0.0, color='grey')

else:

circle_rmin = plt.Circle(xy=(0,0), radius=grid_rmin, color='gray', fill=False, linestyle='-')

ax.add_patch(circle_rmin)

circle_wkrmin = plt.Circle(xy=(0,0), radius=grid_wkzrmin, color='gray', fill=False, linestyle='--')

ax.add_patch(circle_wkrmin)

circle_wkrrmaz = plt.Circle(xy=(0,0), radius=grid_wkzrmax, color='gray', fill=False, linestyle='--')

ax.add_patch(circle_wkrrmaz)

circle_rmax = plt.Circle(xy=(0,0), radius=grid_rmax, color='gray', fill=False, linestyle='-')

ax.add_patch(circle_rmax)

# ------------------

# Add user-defined string in top right corner

# ------------------

if ( ('display_label' in open('paramsf2p.dat').read()) and (par.display_label != '#') ):

mylabel = par.display_label.replace('_',' ')

xlabel = xlim_max - 0.05*(xlim_max-xlim_min)

ylabel = ylim_max - 0.05*(ylim_max-ylim_min)

#print('xlabel=',xlabel, ' ylabel=', ylabel)

if par.mycolormap == 'RdBu_r' or par.mycolormap == 'seismic':

color_label = 'black'

else:

color_label = 'white'

ax.text(xlabel, ylabel, mylabel, fontsize=20, color = color_label,weight='bold', horizontalalignment='right')

# ----------------

# special case all fluids

# ----------------

if par.allfluids == 'Yes':

colorstr = 'lightpink' # 'white'

xmin,xmax = ax.get_xlim()

ymin,ymax = ax.get_ylim()

if f > 0:

if par.dust_size[f-1] > 0.0:

strsize = str_fmt(par.dust_size[f-1])

strsize += ' m'

ax.text(xmax,ymin,strsize,fontsize=20,color=colorstr,horizontalalignment='right',verticalalignment='bottom')

elif par.dust_stokes[f-1] > 0.0:

strstnb = 'St = '+str_fmt(par.dust_stokes[f-1])

ax.text(xmax,ymin,strstnb,fontsize=20,color=colorstr,horizontalalignment='right',verticalalignment='bottom')

else:

ax.text(xmax,ymin,'gas',fontsize=20,color=colorstr,horizontalalignment='right',verticalalignment='bottom')

# ----------------

# plot color-bars

# ----------------

divider = make_axes_locatable(ax)

cax = divider.append_axes("top", size="2.5%", pad=0.12)

cb = plt.colorbar(CF, cax=cax, orientation='horizontal')

cax.xaxis.tick_top()

cax.xaxis.set_tick_params(direction='out')

#cax.xaxis.set_major_locator(plt.MaxNLocator(3))

cax.xaxis.set_major_locator(plt.MaxNLocator(4))

# title on top

cax.xaxis.set_label_position('top')

cax.set_xlabel(strfield)

cax.xaxis.labelpad = 8

if par.log_colorscale == 'Yes' and par.fieldmin != 'auto' and par.fieldmax != 'auto':

cax.xaxis.set_major_locator(ticker.LogLocator(base=10.0,numticks=8))

if par.log_colorscale == 'Yes' and (par.fieldmin == 'auto' or par.fieldmax == 'auto'):

cax.xaxis.set_tick_params(direction='out')

if two_color_bars == True:

cax2 = divider.append_axes("right", size="2.5%", pad=0.12)

cb2 = plt.colorbar(CD, cax=cax2, orientation='vertical')

cax2.yaxis.tick_right()

cb2.set_alpha(0.7)

from packaging import version

# if matplotlib version >= 3.6, draw_all no longer

# exists and becomes _draw_all

if version.parse(matplotlib.__version__) < version.parse("3.6"):

cb2.draw_all()

else:

cb2._draw_all()

cax2.yaxis.set_tick_params(direction='out')

# title on right-hand side

cax2.set_ylabel('dust size [meter]')

cax2.yaxis.labelpad = 8

# ------------------

# save in pdf or png files

# ------------------

if ('filename' in open('paramsf2p.dat').read()) and (par.filename != '#'):

outfile = par.filename

else:

outfile = par.fluid+'_'+par.whatfield+'_'+directory+'_'+par.fieldofview+'_'+par.slice+'_'+par.z_average+'_'+str(on[k]).zfill(4)

if par.movie == 'Yes' and par.take_one_point_every != 1:

outfile = par.fluid+'_'+par.whatfield+'_'+directory+'_'+par.fieldofview+'_'+par.slice+'_'+par.z_average+'_'+str(k).zfill(4)

if par.showdust == 'Yes':

outfile += '_dust'

fileout = outfile+'.pdf'

if par.saveaspdf == 'Yes':

plt.savefig('./'+fileout, dpi=80)

if par.saveaspng == 'Yes':

plt.savefig('./'+re.sub('.pdf', '.png', fileout), dpi=80) #80 is optimum

if par.movie == 'Yes':

plt.close(fig) # close figure as we reopen figure at every output number

# ------------------

# finally concatenate png if movie requested

# ------------------

if par.movie == 'Yes':

# png files that have been created above

allpngfiles = [par.fluid+'_'+par.whatfield+'_'+directory+'_'+par.fieldofview+'_'+par.slice+'_'+par.z_average+'_'+str(on[x]).zfill(4)+'.png' for x in range(len(on))]