def run_with_vertical_hseq(prefix, model, n_iter=10, mpi=False,
                           n_processes=multiprocessing.cpu_count(),
                           overwrite=False):

    from ..grid import CylindricalPolarGrid
    from .model_output import ModelOutput
    from ..util.integrate import integrate

    if not isinstance(model, AnalyticalYSOModel):
        raise TypeError("Can only run hydrostatic equilibrium for AnalyticalYSOModel instances")

    if model.grid['grid_type'] != 'cylindrical':
        raise TypeError("Can only run hydrostatic equilibrium for models with cylindrical polar grids")

    if model.star.mass is None:
        raise ValueError("Stellar mass needs to be defined for calculation of hydrostatic equilibrium")

    if len(model.disks) == 0:
        raise ValueError("Can only run hydrostatic equilibrium for models with disks")
    else:
        n_disks = len(model.disks)

    # Write out initial model
    model.write(prefix + '_00000.rtin', overwrite=overwrite, merge_if_possible=False)

    # Run the initial model
    mo = model.run(prefix + '_00000.rtout', overwrite=overwrite,
                   mpi=mpi, n_processes=n_processes)

    previous = prefix + '_00000.rtout'

    for iteration in range(1, n_iter + 1):

        # Read in output
        mo = ModelOutput(previous)

        # Extract the quantities
        g = mo.get_quantities()

        # Get wall positions
        rw, zw = g.w_wall, g.z_wall

        # Make a 2-d grid of wall positions
        R, Z = np.meshgrid(rw, zw)

        # Extract density and temperature
        density = g['density']
        temperature = g['temperature']

        # TODO: need to find a better way than just assuming the first n
        # density grids are disks

        for idisk in range(n_disks):

            # Vertically extrapolate temperatures
            for i in range(len(g.w)):
                for j in range(len(g.p)):
                    reset = temperature[idisk].array[j, :, i] < 1.
                    temperature[idisk].array[j, reset, i] = np.max(temperature[idisk].array[j, :, i])  # shouldn't be max, but will do for now

            # Compute new density
            for i in range(len(g.w)):
                for j in range(len(g.p)):
                    density[idisk].array[j, :, i] = hseq_profile(g.w[i], g.z, temperature[idisk].array[j, :, i], model.star.mass) * integrate(g.z, density[idisk].array[j, :, i])

        # Instantiate new model based on previous
        m = Model.read(previous)

        # Override the density
        m.grid['density'] = density

        # Write and run
        m.write('{0:s}_{1:05d}.rtin'.format(prefix, iteration), overwrite=overwrite)
        m.run('{0:s}_{1:05d}.rtout'.format(prefix, iteration),
              overwrite=overwrite, mpi=mpi, n_processes=n_processes)

        previous = '{0:s}_{1:05d}.rtout'.format(prefix, iteration)
Run a model with vertical hydrostatic equilibrium.
dataset/ETHPy150Open hyperion-rt/hyperion/hyperion/model/helpers.py/run_with_vertical_hseq