Dear yt users,

I am running a binary star simulation with enzo, and in my box there is a star made up of a gas envelope and a particle as core, and a secondary star represented by a single particle (no gas).

I would like to plot the equipotential surfaces of the system on a slice plot, hence, since the default Grav_Potential field does not include the particles potential, I created a custom made field as follows:

def _Total_Grav_Potential(field,data):

#G in cgs

grav_const_cgs = 6.674e-8

#data for gas

gas_grav_potential_cgs = data['Grav_Potential'] * (length_unit1/time_unit1)**2.0 #in units of velocity^2

gas_position_x_cgs = data['x'] * length_unit1

gas_position_y_cgs = data['y'] * length_unit1

gas_position_z_cgs = data['z'] * length_unit1

#initialize the array for the total grav potential to the gas grav potential

total_grav_potential_cgs = gas_grav_potential_cgs

#loops on the particles

for i in range(len(data['ParticleMassMsun'])):

#data of the current particle

current_particle_mass_cgs= data['ParticleMassMsun'][i] * mass_unit

current_particle_position_x_cgs = data['particle_position_x'][i] * length_unit1

current_particle_position_y_cgs = data['particle_position_y'][i] * length_unit1

current_particle_position_z_cgs = data['particle_position_z'][i] * length_unit1

#computes the array of grav potential for the current particle

current_particle_grav_potential_cgs = -(grav_const_cgs * current_particle_mass_cgs)/(((gas_position_x_cgs - current_particle_position_x_cgs)**2.0 + (gas_position_y_cgs - current_particle_position_y_cgs)**2.0 + (gas_position_z_cgs - current_particle_position_z_cgs)**2.0)**0.5)

#adds the current particle potential to the total potential

total_grav_potential_cgs = total_grav_potential_cgs + current_particle_grav_potential_cgs

return total_grav_potential_cgs

#adds the field to the available yt fields

add_field("Total_Grav_Potential", function=_Total_Grav_Potential, take_log=True ,units=r"\rm{erg}/\rm{g}")

where the length unit and time unit variables have the function to convert my quantities in cgs.

If I then call my "Total_Grav_Potential" field, I get reasonable values, but when I try to inspect the shape of the equipotential surfaces on the xy-plane by plotting the following:

sp = SlicePlot(pf, 'z', "Density", width = 1.0)

sp.annotate_velocity(factor=16, normalize=True)

sp.annotate_contour('Total_Grav_Potential',clim=(max_total_pot,min_total_pot))

sp.annotate_particles(1.0, p_size=50.0, marker='o', col='black')

sp.annotate_text((0.7,1.05), "time = "+str(pf.current_time)+"yr")

sp.save(plot_dir+"density_slice_z_"+str(index)+".png")

nothing happens, and the plot does not show any additional contour; while if I do the same with the default "Grav_Potential" field, the gas equipotential contours get plotted without problems.

I also tried to play with the "annotate_contour" options "ncont", "factor" and "clim" to refine my plot but the result is always the same.

Additionally, if I do not set the "clim" option the following error shows up:

ValueError: zero-size array to minimum.reduce without identity

I would like to ask if my custom definition of the Total_Grav_Potential field is correct from a yt coding point of view and if I should modify it (or if there is any improvement I can do).

Also, I would like to ask if you have any idea of what could be the problem with the plotting.

Thanks for the help,

Roberto