# Shearing sheet (Hill's approximation) (C)

This example simulates a small patch of Saturn's Rings in shearing sheet coordinates. If you have OpenGL enabled, you'll see one copy of the computational domain. Press g to see the ghost boxes which are used to calculate gravity and collisions. Particle properties resemble those found in Saturn's rings.

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "rebound.h"

double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v);
void heartbeat(struct reb_simulation* const r);

int main(int argc, char* argv[]) {
struct reb_simulation* r = reb_create_simulation();
// Setup constants
r->opening_angle2    = .5;                    // This determines the precission of the tree code gravity calculation.
r->integrator            = REB_INTEGRATOR_SEI;
r->boundary            = REB_BOUNDARY_SHEAR;
r->gravity            = REB_GRAVITY_TREE;
r->collision            = REB_COLLISION_TREE;
r->collision_resolve = reb_collision_resolve_hardsphere;
double OMEGA             = 0.00013143527;    // 1/s
r->ri_sei.OMEGA         = OMEGA;
r->G                 = 6.67428e-11;        // N / (1e-5 kg)^2 m^2
r->softening             = 0.1;            // m
r->dt                 = 1e-3*2.*M_PI/OMEGA;    // s
r->heartbeat            = heartbeat;    // function pointer for heartbeat
// This example uses two root boxes in the x and y direction.
// Although not necessary in this case, it allows for the parallelization using MPI.
// See Rein & Liu for a description of what a root box is in this context.
double surfacedensity         = 400;             // kg/m^2
double particle_density        = 400;            // kg/m^3
double particle_radius_min     = 1;            // m
double particle_radius_max     = 4;            // m
double boxsize             = 100;            // m
if (argc>1){                        // Try to read boxsize from command line
boxsize = atof(argv[1]);
}
reb_configure_box(r, boxsize, 2, 2, 1);
r->nghostx = 2;
r->nghosty = 2;
r->nghostz = 0;

// Initial conditions
printf("Toomre wavelength: %f\n",4.*M_PI*M_PI*surfacedensity/OMEGA/OMEGA*r->G);
// Use Bridges et al coefficient of restitution.
r->coefficient_of_restitution = coefficient_of_restitution_bridges;
// When two particles collide and the relative velocity is zero, the might sink into each other in the next time step.
// By adding a small repulsive velocity to each collision, we prevent this from happening.
r->minimum_collision_velocity = particle_radius_min*OMEGA*0.001;  // small fraction of the shear accross a particle

double total_mass = surfacedensity*r->boxsize.x*r->boxsize.y;
double mass = 0;
while(mass<total_mass){
struct reb_particle pt;
pt.x         = reb_random_uniform(r, -r->boxsize.x/2.,r->boxsize.x/2.);
pt.y         = reb_random_uniform(r, -r->boxsize.y/2.,r->boxsize.y/2.);
pt.z         = reb_random_normal(r, 1.);                    // m
pt.vx         = 0;
pt.vy         = -1.5*pt.x*OMEGA;
pt.vz         = 0;
pt.ax         = 0;
pt.ay         = 0;
pt.az         = 0;
pt.m         = particle_mass;     // kg
mass += particle_mass;
}
reb_integrate(r, INFINITY);
}

// This example is using a custom velocity dependend coefficient of restitution
double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v){
// assumes v in units of [m/s]
double eps = 0.32*pow(fabs(v)*100.,-0.234);
if (eps>1) eps=1;
if (eps<0) eps=0;
return eps;
}

void heartbeat(struct reb_simulation* const r){
if (reb_output_check(r, 1e-3*2.*M_PI/r->ri_sei.OMEGA)){
reb_output_timing(r, 0);
//reb_output_append_velocity_dispersion("veldisp.txt");
}
if (reb_output_check(r, 2.*M_PI/r->ri_sei.OMEGA)){
//reb_output_ascii("position.txt");
}
}


This example is located in the directory examples/shearing_sheet