/* Independently authored restricted-disk force for optional REBOUND reference runs.
 * Compile against the separately installed, pinned REBOUND header. No package
 * code is copied or distributed here. The compiled library stays temporary.
 */
#include <math.h>
#include "rebound.h"

static double surface_density = 1e-4;
static double aspect_at_one = 0.055;
static double temperature_slope = 0.5;
static double viscosity = 1e-6;
static int use_transition = 0;

void ac15_configure(double density, double aspect, double slope, double diffusion, int transition) {
    surface_density = density;
    aspect_at_one = aspect;
    temperature_slope = slope;
    viscosity = diffusion;
    use_transition = transition;
}

static double cutoff(double parameter, double turnover) {
    const double threshold = sqrt(4.0*turnover/(45.0*M_PI));
    if (parameter < threshold) {
        return 16.0/25.0*pow(45.0*M_PI/(4.0*turnover),0.75)*pow(parameter,1.5);
    }
    return 1.0-9.0/25.0*pow(4.0*turnover/(45.0*M_PI),4.0/3.0)*pow(parameter,-8.0/3.0);
}

static void disk_acceleration(double mass, double x, double y, double vx, double vy, double* ax, double* ay) {
    *ax = 0.0;
    *ay = 0.0;
    if (surface_density == 0.0) return;
    const double radius = hypot(x,y);
    const double gravity = 1.0+mass;
    const double velocity_squared = vx*vx+vy*vy;
    const double radial_product = x*vx+y*vy;
    const double energy = velocity_squared/2.0-gravity/radius;
    if (energy >= 0.0) return;
    const double semimajor = -gravity/(2.0*energy);
    const double ex = ((velocity_squared-gravity/radius)*x-radial_product*vx)/gravity;
    const double ey = ((velocity_squared-gravity/radius)*y-radial_product*vy)/gravity;
    const double eccentricity = hypot(ex,ey);
    double slope = temperature_slope;
    double log_temperature = -slope*log(semimajor);
    if (use_transition) {
        const double coordinate = log(semimajor/1.2)/0.5;
        const double reference = log(1.0/1.2)/0.5;
        slope = 0.9-0.5*tanh(coordinate);
        log_temperature = -0.9*log(semimajor)+0.25*(log(cosh(coordinate))-log(cosh(reference)));
    }
    const double aspect = aspect_at_one*sqrt(semimajor*exp(log_temperature));
    const double density = surface_density*pow(semimajor,-0.5);
    const double frequency = pow(semimajor,-1.5);
    const double gamma = 1.4;
    const double diffusivity = 2.0*viscosity;
    const double parameter = 2.0*diffusivity/(3.0*pow(aspect,3.0)*semimajor*semimajor*frequency);
    const double squared = parameter*parameter;
    const double gamma_squared = gamma*gamma;
    const double radical = sqrt(pow(gamma_squared*squared+1.0,2.0)-16.0*squared*(gamma-1.0));
    const double divided_inner = 2.0*((2.0*gamma_squared-16.0*(gamma-1.0)+gamma_squared*gamma_squared*squared)/(radical+1.0)+gamma_squared);
    const double gamma_effective = 2.0*gamma/(gamma+0.5*sqrt(fmax(0.0,divided_inner)));
    const double width = 1.1/pow(gamma_effective,0.25)*sqrt(mass/aspect);
    const double viscous_parameter = 2.0/3.0*sqrt(semimajor*semimajor*frequency*pow(width,3.0)/(2.0*M_PI*viscosity));
    const double thermal_parameter = sqrt(semimajor*semimajor*frequency*pow(width,3.0)/(2.0*M_PI*diffusivity));
    const double viscous_f = 1.0/(1.0+pow(viscous_parameter/1.3,2.0));
    const double thermal_f = 1.0/(1.0+pow(thermal_parameter/1.3,2.0));
    const double viscous_g = cutoff(viscous_parameter,2.0);
    const double thermal_g = cutoff(thermal_parameter,2.0);
    const double viscous_k = cutoff(viscous_parameter,7.0);
    const double thermal_k = cutoff(thermal_parameter,7.0);
    const double entropy_slope = slope-(gamma-1.0)*0.5;
    double wave = (-2.5-1.7*slope+0.05)/gamma_effective;
    const double barotropic = (1.1*viscous_f*viscous_g+0.7*(1.0-viscous_k))/gamma_effective;
    double entropy = 7.9*entropy_slope/pow(gamma_effective,2.0)*viscous_f*thermal_f*sqrt(viscous_g*thermal_g);
    entropy += (2.2-1.4/gamma_effective)*entropy_slope/gamma_effective*sqrt((1.0-viscous_k)*(1.0-thermal_k));
    const double scaled_eccentricity = eccentricity/aspect;
    wave *= (1.0-pow(scaled_eccentricity/2.02,4.0))/(1.0+pow(scaled_eccentricity/2.25,1.2)+pow(scaled_eccentricity/2.84,6.0));
    const double corotation = (barotropic+entropy)*exp(-eccentricity/(0.5*aspect+0.01));
    const double normalization = pow(mass/aspect,2.0)*density*pow(semimajor,4.0)*frequency*frequency;
    const double wave_time = pow(aspect,4.0)/(mass*density*semimajor*semimajor*frequency);
    const double eccentricity_time = wave_time/0.780*(1.0-0.14*scaled_eccentricity*scaled_eccentricity+0.06*pow(scaled_eccentricity,3.0));
    const double signed_momentum = mass*(x*vy-y*vx);
    const double rate = (wave+corotation)*normalization/signed_momentum;
    const double damping = -2.0*radial_product/(radius*radius*eccentricity_time);
    *ax = rate*vx+damping*x;
    *ay = rate*vy+damping*y;
}

void ac15_acceleration(double mass, double x, double y, double vx, double vy, double* ax, double* ay) {
    disk_acceleration(mass,x,y,vx,vy,ax,ay);
}

void ac15_forces(struct reb_simulation* simulation) {
    struct reb_particle* star = &simulation->particles[0];
    for (size_t index=1; index<simulation->N; index++) {
        struct reb_particle* planet = &simulation->particles[index];
        double ax,ay;
        disk_acceleration(planet->m,planet->x-star->x,planet->y-star->y,planet->vx-star->vx,planet->vy-star->vy,&ax,&ay);
        planet->ax += ax;
        planet->ay += ay;
    }
}