text:
/* multifit_nlinear/trust.c
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>

#include "common.c"
#include "nielsen.c"

/*
 * This module contains a high level driver for a general trust
 * region nonlinear least squares solver. This container handles
 * the computation of all of the quantities relevant to all trust
 * region methods, including:
 *
 * residual vector: f_k = f(x_k)
 * Jacobian matrix: J_k = J(x_k)
 * gradient vector: g_k = J_k^T f_k
 * scaling matrix:  D_k
 */

typedef struct
{
  size_t n;                  /* number of observations */
  size_t p;                  /* number of parameters */
  double delta;              /* trust region radius */
  double mu;                 /* LM parameter */
  long nu;                   /* for updating LM parameter */
  gsl_vector *diag;          /* D = diag(J^T J) */
  gsl_vector *x_trial;       /* trial parameter vector */
  gsl_vector *f_trial;       /* trial function vector */
  gsl_vector *workp;         /* workspace, length p */
  gsl_vector *workn;         /* workspace, length n */

  void *trs_state;           /* workspace for trust region subproblem */
  void *solver_state;        /* workspace for linear least squares solver */

  double avratio;            /* current |a| / |v| */

  /* tunable parameters */
  gsl_multifit_nlinear_parameters params;
} trust_state_t;

static void * trust_alloc (const gsl_multifit_nlinear_parameters * params,
                           const size_t n, const size_t p);
static void trust_free(void *vstate);
static int trust_init(void *vstate, const gsl_vector * swts,
                      gsl_multifit_nlinear_fdf *fdf, const gsl_vector *x,
                      gsl_vector *f, gsl_matrix *J, gsl_vector *g);
static int trust_iterate(void *vstate, const gsl_vector *swts,
                         gsl_multifit_nlinear_fdf *fdf,
                         gsl_vector *x, gsl_vector *f, gsl_matrix *J,
                         gsl_vector *g, gsl_vector *dx);
static int trust_rcond(double *rcond, void *vstate);
static double trust_avratio(void *vstate);
static void trust_trial_step(const gsl_vector * x, const gsl_vector * dx,
                             gsl_vector * x_trial);
static double trust_calc_rho(const gsl_vector * f, const gsl_vector * f_trial,
                             const gsl_vector * g, const gsl_matrix * J,
                             const gsl_vector * dx, trust_state_t * state);
static int trust_eval_step(const gsl_vector * f, const gsl_vector * f_trial,
                           const gsl_vector * g, const gsl_matrix * J,
                           const gsl_vector * dx, double * rho, trust_state_t * state);
static double trust_scaled_norm(const gsl_vector *D, const gsl_vector *a);

static void *
trust_alloc (const gsl_multifit_nlinear_parameters * params,
             const size_t n, const size_t p)
{
  trust_state_t *state;

  state = calloc(1, sizeof(trust_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate lm state", GSL_ENOMEM);
    }

  state->diag = gsl_vector_alloc(p);
  if (state->diag == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for diag", GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workp", GSL_ENOMEM);
    }

  state->workn = gsl_vector_alloc(n);
  if (state->workn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workn", GSL_ENOMEM);
    }

  state->x_trial = gsl_vector_alloc(p);
  if (state->x_trial == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for x_trial", GSL_ENOMEM);
    }

  state->f_trial = gsl_vector_alloc(n);
  if (state->f_trial == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for f_trial", GSL_ENOMEM);
    }

  state->trs_state = (params->trs->alloc)(params, n, p);
  if (state->trs_state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for trs state", GSL_ENOMEM);
    }

  state->solver_state = (params->solver->alloc)(n, p);
  if (state->solver_state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for solver state", GSL_ENOMEM);
    }

  state->n = n;
  state->p = p;
  state->delta = 0.0;
  state->params = *params;

  return state;
}

static void
trust_free(void *vstate)
{
  trust_state_t *state = (trust_state_t *) vstate;
  const gsl_multifit_nlinear_parameters *params = &(state->params);

  if (state->diag)
    gsl_vector_free(state->diag);

  if (state->workp)
    gsl_vector_free(state->workp);

  if (state->workn)
    gsl_vector_free(state->workn);

  if (state->x_trial)
    gsl_vector_free(state->x_trial);

  if (state->f_trial)
    gsl_vector_free(state->f_trial);

  if (state->trs_state)
    (params->trs->free)(state->trs_state);

  if (state->solver_state)
    (params->solver->free)(state->solver_state);

  free(state);
}

/*
trust_init()
  Initialize trust region solver

Inputs: vstate - workspace
        swts   - sqrt(W) vector
        fdf    - user callback functions
        x      - initial parameter values
        f      - (output) f(x) vector
        J      - (output) J(x) matrix
        g      - (output) J(x)' f(x) vector

Return: success/error
*/

static int
trust_init(void *vstate, const gsl_vector *swts,
           gsl_multifit_nlinear_fdf *fdf, const gsl_vector *x,
           gsl_vector *f, gsl_matrix *J, gsl_vector *g)
{
  int status;
  trust_state_t *state = (trust_state_t *) vstate;
  const gsl_multifit_nlinear_parameters *params = &(state->params);
  double Dx;

  /* evaluate function and Jacobian at x and apply weight transform */
  status = gsl_multifit_nlinear_eval_f(fdf, x, swts, f);
  if (status)
   return status;

  status = gsl_multifit_nlinear_eval_df(x, f, swts, params->h_df,
                                        params->fdtype, fdf, J, state->workn);
  if (status)
    return status;

  /* compute g = J^T f */
  gsl_blas_dgemv(CblasTrans, 1.0, J, f, 0.0, g);

  /* initialize diagonal scaling matrix D */
  (params->scale->init)(J, state->diag);

  /* compute initial trust region radius */
  Dx = trust_scaled_norm(state->diag, x);
  state->delta = 0.3 * GSL_MAX(1.0, Dx);

  /* initialize LM parameter */
  status = nielsen_init(J, state->diag, &(state->mu), &(state->nu));
  if (status)
    return status;

  /* initialize trust region method solver */
  {
    gsl_multifit_nlinear_trust_state trust_state;

    trust_state.x = x;
    trust_state.f = f;
    trust_state.g = g;
    trust_state.J = J;
    trust_state.diag = state->diag;
    trust_state.sqrt_wts = swts;
    trust_state.mu = &(state->mu);
    trust_state.params = params;
    trust_state.solver_state = state->solver_state;
    trust_state.fdf = fdf;
    trust_state.avratio = &(state->avratio);

    status = (params->trs->init)(&trust_state, state->trs_state);

    if (status)
      return status;
  }

  /* set default parameters */

  state->avratio = 0.0;

  return GSL_SUCCESS;
}

/*
trust_iterate()
  This function performs 1 iteration of the trust region algorithm.
It calls a user-specified method for computing the next step
(LM or dogleg), then tests if the computed step is acceptable.

Args: vstate - trust workspace
      swts   - data weights (NULL if unweighted)
      fdf    - function and Jacobian pointers
      x      - on input, current parameter vector
               on output, new parameter vector x + dx
      f      - on input, f(x)
               on output, f(x + dx)
      J      - on input, J(x)
               on output, J(x + dx)
      g      - on input, g(x) = J(x)' f(x)
               on output, g(x + dx) = J(x + dx)' f(x + dx)
      dx     - (output only) parameter step vector

Return:
1) GSL_SUCCESS if we found a step which reduces the cost
function

2) GSL_ENOPROG if 15 successive attempts were to made to
find a good step without success

3) If a scaling matrix D is used, inputs and outputs are
set to the unscaled quantities (ie: J and g)
*/

static int
trust_iterate(void *vstate, const gsl_vector *swts,
              gsl_multifit_nlinear_fdf *fdf, gsl_vector *x,
              gsl_vector *f, gsl_matrix *J, gsl_vector *g,
              gsl_vector *dx)
{
  int status;
  trust_state_t *state = (trust_state_t *) vstate;
  const gsl_multifit_nlinear_parameters *params = &(state->params);
  const gsl_multifit_nlinear_trs *trs = params->trs;
  gsl_multifit_nlinear_trust_state trust_state;
  gsl_vector *x_trial = state->x_trial;       /* trial x + dx */
  gsl_vector *f_trial = state->f_trial;       /* trial f(x + dx) */
  gsl_vector *diag = state->diag;             /* diag(D) */
  double rho;                                 /* ratio actual_reduction/predicted_reduction */
  int foundstep = 0;                          /* found step dx */
  int bad_steps = 0;                          /* consecutive rejected steps */

  /* store all state parameters needed by low level methods */
  trust_state.x = x;
  trust_state.f = f;
  trust_state.g = g;
  trust_state.J = J;
  trust_state.diag = state->diag;
  trust_state.sqrt_wts = swts;
  trust_state.mu = &(state->mu);
  trust_state.params = params;
  trust_state.solver_state = state->solver_state;
  trust_state.fdf = fdf;
  trust_state.avratio = &(state->avratio);

  /* initialize trust region subproblem with this Jacobian */
  status = (trs->preloop)(&trust_state, state->trs_state);
  if (status)
    return status;

  /* loop until we find an acceptable step dx */
  while (!foundstep)
    {
      /* calculate new step */
      status = (trs->step)(&trust_state, state->delta, dx, state->trs_state);

      /* occasionally the iterative methods (ie: CG Steihaug) can fail to find a step,
       * so in this case skip rho calculation and count it as a rejected step */

      if (status == GSL_SUCCESS)
        {
          /* compute x_trial = x + dx */
          trust_trial_step(x, dx, x_trial);

          /* compute f_trial = f(x + dx) */
          status = gsl_multifit_nlinear_eval_f(fdf, x_trial, swts, f_trial);
          if (status)
            return status;

          /* check if step should be accepted or rejected */
          status = trust_eval_step(f, f_trial, g, J, dx, &rho, state);
          if (status == GSL_SUCCESS)
            foundstep = 1;
        }
      else
        {
          /* an iterative TRS method failed to find a step vector */
          rho = -1.0;
        }

      /*
       * update trust region radius: if rho is large,
       * then the quadratic model is a good approximation
       * to the objective function, enlarge trust region.
       * If rho is small (or negative), the model function
       * is a poor approximation so decrease trust region. This
       * can happen even if the step is accepted.
       */
      if (rho > 0.75)
        state->delta *= params->factor_up;
      else if (rho < 0.25)
        state->delta /= params->factor_down;

      if (foundstep)
        {
          /* step was accepted */

          /* compute J <- J(x + dx) */
          status = gsl_multifit_nlinear_eval_df(x_trial, f_trial, swts,
                                                params->h_df, params->fdtype,
                                                fdf, J, state->workn);
          if (status)
            return status;

          /* update x <- x + dx */
          gsl_vector_memcpy(x, x_trial);

          /* update f <- f(x + dx) */
          gsl_vector_memcpy(f, f_trial);

          /* compute new g = J^T f */
          gsl_blas_dgemv(CblasTrans, 1.0, J, f, 0.0, g);

          /* update scaling matrix D */
          (params->scale->update)(J, diag);

          /* step accepted, decrease LM parameter */
          status = nielsen_accept(rho, &(state->mu), &(state->nu));
          if (status)
            return status;

          bad_steps = 0;
        }
      else
        {
          /* step rejected, increase LM parameter */
          status = nielsen_reject(&(state->mu), &(state->nu));
          if (status)
            return status;

          if (++bad_steps > 15)
            {
              /* if more than 15 consecutive rejected steps, report no progress */
              return GSL_ENOPROG;
            }
        }
    }

  return GSL_SUCCESS;
} /* trust_iterate() */

static int
trust_rcond(double *rcond, void *vstate)
{
  int status;
  trust_state_t *state = (trust_state_t *) vstate;
  const gsl_multifit_nlinear_parameters *params = &(state->params);

  status = (params->solver->rcond)(rcond, state->solver_state);

  return status;
}

static double
trust_avratio(void *vstate)
{
  trust_state_t *state = (trust_state_t *) vstate;
  return state->avratio;
}

/* compute x_trial = x + dx */
static void
trust_trial_step(const gsl_vector * x, const gsl_vector * dx,
                 gsl_vector * x_trial)
{
  size_t i, N = x->size;

  for (i = 0; i < N; i++)
    {
      double dxi = gsl_vector_get (dx, i);
      double xi = gsl_vector_get (x, i);
      gsl_vector_set (x_trial, i, xi + dxi);
    }
}

/*
trust_calc_rho()
  Calculate ratio of actual reduction to predicted
reduction.

rho = actual_reduction / predicted_reduction

actual_reduction = 1 - ( ||f+|| / ||f|| )^2
predicted_reduction = -2 g^T dx / ||f||^2 - ( ||J*dx|| / ||f|| )^2
                    = -2 fhat . beta - ||beta||^2

where: beta = J*dx / ||f||

Inputs: f        - f(x)
        f_trial  - f(x + dx)
        g        - gradient J^T f
        J        - Jacobian
        dx       - proposed step, size p
        state    - workspace

Return: rho = actual_reduction / predicted_reduction
If actual_reduction is < 0, return rho = -1
*/

static double
trust_calc_rho(const gsl_vector * f, const gsl_vector * f_trial,
               const gsl_vector * g, const gsl_matrix * J,
               const gsl_vector * dx, trust_state_t * state)
{
  int status;
  const gsl_multifit_nlinear_parameters *params = &(state->params);
  const gsl_multifit_nlinear_trs *trs = params->trs;
  const double normf = gsl_blas_dnrm2(f);
  const double normf_trial = gsl_blas_dnrm2(f_trial);
  gsl_multifit_nlinear_trust_state trust_state;
  double rho;
  double actual_reduction;
  double pred_reduction;
  double u;

  /* if ||f(x+dx)|| > ||f(x)|| reject step immediately */
  if (normf_trial >= normf)
    return -1.0;

  trust_state.x = NULL;
  trust_state.f = f;
  trust_state.g = g;
  trust_state.J = J;
  trust_state.diag = state->diag;
  trust_state.sqrt_wts = NULL;
  trust_state.mu = &(state->mu);
  trust_state.params = params;
  trust_state.solver_state = state->solver_state;
  trust_state.fdf = NULL;
  trust_state.avratio = &(state->avratio);

  /* compute numerator of rho (actual reduction) */
  u = normf_trial / normf;
  actual_reduction = 1.0 - u*u;

  /*
   * compute denominator of rho (predicted reduction); this is calculated
   * inside each trust region subproblem, since it depends on the local
   * model used, which can vary according to each TRS
   */
  status = (trs->preduction)(&trust_state, dx, &pred_reduction, state->trs_state);
  if (status)
    return -1.0;

  if (pred_reduction > 0.0)
    rho = actual_reduction / pred_reduction;
  else
    rho = -1.0;

  return rho;
}

/*
trust_eval_step()
  Evaluate proposed step to determine if it should be
accepted or rejected
*/

static int
trust_eval_step(const gsl_vector * f, const gsl_vector * f_trial,
                const gsl_vector * g, const gsl_matrix * J,
                const gsl_vector * dx, double * rho, trust_state_t * state)
{
  int status = GSL_SUCCESS;
  const gsl_multifit_nlinear_parameters *params = &(state->params);

  if (params->trs == gsl_multifit_nlinear_trs_lmaccel)
    {
      /* reject step if acceleration is too large compared to velocity */
      if (state->avratio > params->avmax)
        status = GSL_FAILURE;
    }

  /* compute rho */
  *rho = trust_calc_rho(f, f_trial, g, J, dx, state);
  if (*rho <= 0.0)
    status = GSL_FAILURE;

  return status;
}

/* compute || diag(D) a || */
static double
trust_scaled_norm(const gsl_vector *D, const gsl_vector *a)
{
  const size_t n = a->size;
  double e2 = 0.0;
  size_t i;

  for (i = 0; i < n; ++i)
    {
      double Di = gsl_vector_get(D, i);
      double ai = gsl_vector_get(a, i);
      double u = Di * ai;

      e2 += u * u;
    }

  return sqrt (e2);
}

static const gsl_multifit_nlinear_type trust_type =
{
  "trust-region",
  trust_alloc,
  trust_init,
  trust_iterate,
  trust_rcond,
  trust_avratio,
  trust_free
};

const gsl_multifit_nlinear_type *gsl_multifit_nlinear_trust = &trust_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
 * Monte Carlo.cpp
 *
 * Copyright 2017 Benjamin Church <ben@U430>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.
 *
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_integration.h>
#include <time.h>

#define ADIABATIC_CUTOFF 1

#define INTEGRATION_POINTS 81
#define INTEGRAL_FUDGE_FACTOR 1.0125
#define INTEGRATION_SECTIONS 8
#define CUTOFF_SCALE 1E3
#define min(x,y) (x < y? x : y)
#define sign(x) (x >= 0? 1.0 : -1.0)
#define dot_macro(A, B) ((A).x * (B).x + (A).y * (B).y + (A).z * (B).z)
#define mag_sq(A) dot_macro(A, A)
#define mag(A) sqrt(dot_macro(A, A))
#define assign_vec(A, x_val, y_val, z_val) {(A).x = x_val; (A).y = y_val; (A).z = z_val;}
#define assign_vec_diff(d, A, B) {(d).x = (A).x - (B).x; (d).y = (A).y - (B).y; (d).z = (A).z - (B).z;}
#define make_unit(A) {double len = mag(A); assign_vec(A, (A.x)/len, (A.y)/len, (A.z)/len);}
#define rho_profile(x, a, b) (pow(x, 2 - a)*pow((1 + x), a - b))

gsl_rng *RNG;
gsl_integration_workspace * w;

const int default_log_num_trials = 6;
int num_trials, num_points = 100;

const double pi = 3.14159265;
double crit_density = 1.3211775E-7,
	f = 0.1,
	g = 1E-4,
	p = 1.9,
	G = 0.0045,
	k = 2,
	M_prim = 1E12,
	T_age = 1E4;

double m_max, m_min, expected_num, R_core_prim, R_max_prim, c_prim;

double MaxRadius(double M)
{
    return pow(3.0*M/(4.0 * pi * 200.0 * crit_density), 1.0/3.0);
}

typedef struct
{
	double x, y, z;
} vector;

typedef struct
{
	double m, s, max_val;
}data_cell;

void update_cell(data_cell *ptr, double new_data, int k)
{
	double old_m = ptr->m;
	ptr->m += (new_data - old_m)/(double) k;
	ptr->s += (new_data - old_m)*(new_data - ptr->m);
	if(ptr->max_val < new_data) ptr->max_val = new_data;
}

void reset_cell(data_cell *ptr)
{
	ptr->m = 0;
	ptr->s = 0;
	ptr->max_val = 0;
}

typedef struct
{
	double R, theta, phi;
	double cos_theta;
	double M, r_core, r_max, alpha, beta, c, normalization;
	double integs[INTEGRATION_SECTIONS];
	vector v, position;
} halo;

double Power(double x, int y)
{
	if(y > 0)
	{
		double result = 1;
		while(y > 0)
		{
			if(y % 2 == 0)
			{
				x = x*x;
				y = y/2;
			}
			else
			{
				result = x*result;
				y--;
			}
		}
		return result;
	}
	else if(y < 0)
		return 1/Power(x, -y);
	else
		return 1;
}

void print_vector(vector *A)
{
	printf("(%.3f, %.3f, %.3f)\n", A->x, A->y, A->z);
}

double perp_mag_sq(vector *A, vector *u)
{
	double dot_prod = dot_macro(*A, *u);
	return mag_sq(*A) - dot_prod*dot_prod/mag_sq(*u);
}

double MFreeNFW(double r)
{
    if(r < R_max_prim)
        return M_prim*(log(1 + r/R_core_prim)-r/(r + R_core_prim))/(log(1 + c_prim) - c_prim/( 1 + c_prim));
    else
        return M_prim;
}

double DFreeNFW(double r)
{
    if(r < R_max_prim)
        return 200.0/3.0 * crit_density / (log(1 + c_prim) - c_prim/(1 + c_prim))*pow(c_prim, 3) * 1.0/(r/R_core_prim*pow(1+r/R_core_prim, 2.0));
    else
        return 0.0;
}

double PhiFreeNFW(double r)
{
    if(r < R_max_prim)
        return -M_prim*G*((R_max_prim/r * log(1 + r/R_core_prim) - log(1 + c_prim))/(log(1 + c_prim) - c_prim/(1 + c_prim)) + 1)/R_max_prim;
    else
        return -M_prim*G/r;
}

double TidalRadius(double M, double R)
{
	return R*pow(M/(2*MFreeNFW(R)), 1.0/3.0);
}

double NFW_func(double x, double r)
{
    return (log(1.0 + x) - x/(1.0 + x))/(log(1.0 + c_prim) - c_prim/(1.0 + c_prim)) - r;
}

double df (double x)
{
    return x/((1.0 + x)*(1.0 + x)) * 1.0/(log(1.0 + c_prim) - c_prim/(1.0 + c_prim));
}

double newton(double r)
{
    int itr, maxmitr = 100;
    double h, x0 = 0.1, x1, allerr = pow(10, -6);

    for (itr = 0; itr < maxmitr; itr++)
    {
        h = NFW_func(x0, r)/df(x0);
        x1 = x0 - h;

        if (fabs(h) < allerr)
        {
            return x1;
        }
        x0 = x1;
    }
    return x1;
}

double rho_profile_func(double x, void *params)
{
	double *param_ptr = (double *)params;
	return pow(x, 2 - param_ptr[0])*pow((1 + x), param_ptr[0] - param_ptr[1]);
}

double integ_profile(double alpha, double beta, double r)
{
	/*if(alpha > 3)
		return 1;
	double dist = 0, step = r/INTEGRATION_POINTS, sum = 0;
	int i;

	for(i = 0; i < INTEGRATION_POINTS; i++)
	{
		if(i == 0 && alpha > 2 && alpha < 3)
			sum += 8.0/(3.0*step) * 1/(3.0 - alpha) * pow(step, 3 - alpha);
		else if(i == 0 || i == INTEGRATION_POINTS - 1)
			sum += rho_profile(dist, alpha, beta);
		else if(i%3 == 0)
			sum += 2*rho_profile(dist, alpha, beta);
		else
			sum += 3*rho_profile(dist, alpha, beta);
		dist += step;
	}

	return 3.0*step/8.0*sum * INTEGRAL_FUDGE_FACTOR;*/
	//printf("int %f\n", 3*step/8*sum);

	double result, error;
	double arr[2];
	arr[0] = alpha;
	arr[1] = beta;

	gsl_function F;
	F.function = &rho_profile_func;
  	F.params = arr;

  	gsl_integration_qags (&F, 0, r, 0, 1e-7, 1000, w, &result, &error);
	if(error/result > 1e-7)printf("LARGE INTEG ERROR WARNING %f for result %f", error, result);

	//printf("alpha: %.3f beta: %.3f r: %.3f err %f\n", alpha, beta, r, (result - quick_result)/result);

	return result;
}

double worst_thing_ever(int lr, double a, double b, double r)
{
	switch(lr)
	{
		case -2: return pow(2,-2 + b)*pow(3,-4 + a - b)*(-20.25 - ((8*a*a*a + 12*a*a*(-2 + b) + b*(38 + (-15 + b)*b) +
					2*a*(8 + 3*(-7 + b)*b))*(-0.0078125 + Power(1 - 2*r,4)/8.))/2. +
					81*r - 27*(-6 + 2*a + b)*(0.1875 + (-1 + r)*r) + (3*(4*a*a + 4*a*(-4 + b) + (-9 + b)*(-2 + b))*(-1 + 4*r)*(7 + 4*r*(-5 + 4*r)))/32.);
		break;

		case -1: return (pow(2,-4 + a - b)*(-24 + 48*r - ((-2 + a + b)*(a*a + (-7 + b)*b + a*(-1 + 2*b))*(-3 + 2*r)*(-1 + 2*r)*(5 + 4*(-2 + r)*r))/64. -
					3*(-4 + a + b)*(3 - 8*r + 4*r*r) + 6*(8 + a*a - 7*b + b*b + a*(-5 + 2*b))*(-0.2916666666666667 + (r*(3 + (-3 + r)*r))/3.)))/3.0;
		break;

		case 0: return pow(2,-4 - a)*pow(3,-4 + a - b)*(-(a*a*a*(-1 + Power(-2 + r,4))) -
					16*(27 + (3*b)/4. + (b*b*b*(-1 + Power(-2 + r,4)))/2. + 40*b*r - 48*b*r*r - 27*r*r*r + 4*b*r*r*r +
					(13*b*r*r*r*r)/4. - 3*b*b*(1 + Power(-2 + r,3)*r)) - 3*a*a*(-1 + r)*(-41 + 2*b*(-3 + r)*(5 + (-4 + r)*r) - r*(-23 + r + r*r)) -
					2*a*(-1 + r)*(-75 + 6*b*b*(-3 + r)*(5 + (-4 + r)*r) + r*(-187 + r*(77 + r)) - 3*b*(37 + r*(5 + r*(-19 + 5*r)))));
		break;

		case 1: return  -(pow(2,-5 - 2*a)*pow(5,-3 + a - b)*(-2 + r)*(a*a*a*(-120 + 68*r - 14*r*r + r*r*r) +
					a*a*(-1880 + 596*r - 38*r*r - 3*r*r*r + 12*b*(-120 + 68*r - 14*r*r + r*r*r)) +
					2*a*(-2200 - 1932*r + 426*r*r + r*r*r + b*(-3920 + 344*r + 268*r*r - 42*r*r*r) + 24*b*b*(-120 + 68*r - 14*r*r + r*r*r)) +
					8*(-500*(4 + 2*r + r*r) + 8*b*b*b*(-120 + 68*r - 14*r*r + r*r*r) - 4*b*b*(40 + 212*r - 86*r*r + 9*r*r*r) +
					b*(-200 - 1892*r + 206*r*r + 31*r*r*r))))/3.0;
		break;

		case 2: return  pow(2,-4 - 3*a)*pow(3,-7 + 2*a - 2*b)*(-8957952 - ((-2 + a)*(-1 + a)*a + 8*(182 + 3*(-11 + a)*a)*b + 192*(-10 + a)*b*b + 512*b*b*b)*
					(-64 + Power(-8 + r,4)/4.) + 2239488*r - 15552*(-18 + a + 8*b)*(48 - 16*r + r*r) +
					72*(162 + a*a - 224*b + 64*b*b + a*(-19 + 16*b))*(-448 + r*(192 + (-24 + r)*r)));

		break;

		case 3: return (pow(2,-5 - 4*a)*pow(17,-3 + a - b)*(-965935104 - (a*a*a + a*a*(-3 + 48*b) + 32*b*(307 - 432*b + 128*b*b) + a*(2 - 912*b + 768*b*b))*
					(-1024 + Power(-16 + r,4)/4.) + 120741888*r - 221952*(-34 + a + 16*b)*(192 - 32*r + r*r) +
					272*(578 + a*a - 832*b + 256*b*b + a*(-35 + 32*b))*(-3584 + r*(768 + (-48 + r)*r))))/3.0;

		break;

		case 4: return pow(2,-6 - 5*a)*pow(33,-4 + a - b)*11*(-113048027136 -
				(a*a*a + a*a*(-3 + 96*b) + 64*b*(1123 - 1632*b + 512*b*b) + a*(2 - 3360*b + 3072*b*b))*(-16384 + Power(-32 + r,4)/4.) + 7065501696*r -
				3345408*(-66 + a + 32*b)*(768 - 64*r + r*r) + 1056*(2178 + a*a - 3200*b + 1024*b*b + a*(-67 + 64*b))*(-28672 + r*(3072 + (-96 + r)*r)));

		break;

		case 5: return (pow(2,-7 - 6*a)*pow(65,-3 + a - b)*
     			(-13822328832000 + 431947776000*r - (((-2 + a)*(-1 + a)*a + 64*(8582 + 3*(-67 + a)*a)*b +
            	12288*(-66 + a)*b*b + 262144*b*b*b)*(-96 + r)*(-32 + r)*
          		(5120 + (-128 + r)*r))/4. - 51916800*(-130 + a + 64*b)*(3072 - 128*r + r*r) + 4160*(8450 + a*a - 12544*b + 4096*b*b + a*(-131 + 128*b))*
        		(-229376 + r*(12288 + (-192 + r)*r))))/3.0;
		break;
	};

	return 1;
}

double sec_integ_profile(halo *ptr, double r)
{
	double a = ptr->alpha, b = ptr->beta;
	int lr = floor(log2(r)); //MAKE BETTER!
	if(lr < -2) return pow(r, 3-a)*(1/(3-a) + r*(a-b)/(4-a) + (a*a - a + b - 2*a*b + b*b)/(5-a)*r*r/2.0);
	if(lr >= INTEGRATION_SECTIONS - 2) return integ_profile(a, b, r); //MAKE BETTER
	return ptr->integs[lr + 2] + worst_thing_ever(lr, a, b, r);
}

double get_M()
{
	double r = gsl_rng_uniform(RNG);
	return M_prim*pow(pow(m_max/M_prim, 1.0-p)*r + pow(m_min/M_prim, 1.0-p)*(1.0-r), 1.0/(1.0-p));
}

double c_bar(double M)
{
	return pow(10.0, 1.0 - 0.1 * (log10(M) - 12));
}

void set_shape(halo *ptr, double M, double R)
{
	ptr->r_max = pow(3*M/(4 * pi * 200.0 * crit_density), 1.0/3.0);
	ptr->alpha = gsl_ran_gaussian_ziggurat(RNG, 0.2) + 1.25; //suggested by J. Ostriker
	ptr->beta = 3; //probalby need to change
	ptr->c = gsl_ran_lognormal(RNG, log(c_bar(M)), 0.25); //Ludlow el. al. (2013) and  Frank van den Bosch
	ptr->r_core = ptr->r_max/ptr->c;
	int i;
	double r = 1.0/4.0, a = ptr->alpha, b = ptr->beta;
	ptr->integs[0] = pow(r, 3-a)*(1/(3-a) + r*(a-b)/(4-a) + (a*a - a + b - 2*a*b + b*b)/(5-a)*r*r/2.0);
	r *= 2;
	for(i = 0; i < INTEGRATION_SECTIONS - 1; i++)
	{
		ptr->integs[i + 1] = ptr->integs[i]  + worst_thing_ever(i - 2, ptr->alpha, ptr->beta, r);
		r *= 2;
	}
	ptr->normalization = sec_integ_profile(ptr, ptr->c);

}

void set_velocity(halo *ptr, double R)
{
	double v_sigma = sqrt(1.0/3.0)*sqrt(-PhiFreeNFW(R));

	ptr->v.x = gsl_ran_gaussian_ziggurat(RNG, v_sigma);
	ptr->v.y = gsl_ran_gaussian_ziggurat(RNG, v_sigma);
	ptr->v.z = gsl_ran_gaussian_ziggurat(RNG, v_sigma);
}

double enclosed_mass(halo *ptr, double r)
{
	double x = r/(ptr->r_core);
	if(x > ptr->c)
		return ptr->M;
	else
	{
	/*	double sec = sec_integ_profile(ptr, x), integ = integ_profile(ptr->alpha, ptr->beta, x);
		double lerr = log(fabs(sec - integ));
		if(lerr > -4)
		{
			printf("alpha: %f beta: %f r: %f \n", ptr->alpha, ptr->beta, x);
			printf("lerr: %f sec: %f real: %f r: \n", lerr, sec, integ);
		}*/
		return ptr->M * sec_integ_profile(ptr, x)/ptr->normalization;
	}
}

void truncate(halo *ptr, double R)
{
	double l2 = R*R*perp_mag_sq(&(ptr->v), &(ptr->position));
	double r0 = l2/(M_prim * G);
	double E = 0.5*mag_sq(ptr->v) + PhiFreeNFW(R);
	double ecc =  sqrt(1.0 + 2.0*E*l2/(M_prim * G * M_prim * G));
	double R_min = min(fabs(r0/(1.0 + ecc)), fabs(r0/(1.0 - ecc)));

	/*printf("velocity: ");
	print_vector(ptr->v);
	printf("position ");
	print_vector(halo_pos(ptr));
	printf("R = %.3f theta = %.3f phi = %.3f\n l2 = %.3f e = %.3f \n", ptr->R, ptr->theta, ptr->phi, l2, ecc);*/

	if(E > 0 || R_min < 0.1)
	{
		ptr->M = 0;
	}
	else
	{
		double Rt = TidalRadius(ptr->M, R_min);
		double R_max = min(ptr->r_max, Rt);
		double new_M = enclosed_mass(ptr, R_max);

		ptr->r_max = R_max;
		ptr->c = R_max/ptr->r_core;
		ptr->normalization *= (new_M/ptr->M);
		ptr->M = new_M;
	}
}

halo *make_halo(halo *ptr)
{
	ptr->R = newton(gsl_rng_uniform(RNG)) * R_core_prim;
	ptr->cos_theta = 2*gsl_rng_uniform(RNG) - 1;
	ptr->theta = acos(ptr->cos_theta);
	ptr->phi = 2*pi*gsl_rng_uniform(RNG);
	assign_vec(ptr->position, ptr->R * sin(ptr->theta) * cos(ptr->phi), ptr->R * sin(ptr->theta) * sin(ptr->phi), ptr->R * ptr->cos_theta);

	ptr->M = get_M();
	set_shape(ptr, ptr->M, ptr->R);
	set_velocity(ptr, ptr->R);
	truncate(ptr, ptr->R);

	return ptr;
}

void print_halo_basic(halo *ptr)
{
	printf("\n R = %3f \n theta = %3f \n M = %3f \n", ptr->R, ptr->theta, ptr->M);
}

double Fluc(halo *halos, int num_halos, double D)
{
	double sum = 0;
	int i;
	vector my_pos = {0, 0, D}, diff;
	double natural_fluc = 2 * pi* pow(D, 3.0/2.0)/sqrt(MFreeNFW(D) * G) * 1/(PhiFreeNFW(D));
	for(i = 0; i < num_halos; i++)
	{
		double R = halos[i].R;
		double r = sqrt(R*R + D*D - 2.0*R*D*halos[i].cos_theta);

		assign_vec_diff(diff, halos[i].position, my_pos);
		make_unit(diff)
		double v_r = dot_macro(halos[i].v, diff);

		double produced_fluc = ((r > CUTOFF_SCALE) ? (enclosed_mass(halos + i, r) * G /(r*r) * v_r) : 0);
		//sum += produced_fluc; //KILLEM
		// kill heating which is adiabtic

		if(produced_fluc/natural_fluc > 1 || ! ADIABATIC_CUTOFF)
			sum += produced_fluc;
		/*printf("%.3f\n", D);
		printf("velocity: ");
		print_vector(halos[i]->v);
		printf("position ");
		print_vector(halo_pos(halos[i]));
		printf("R = %.3f theta = %.3f phi = %.3f\n v_r = %.3f \n", halos[i]->R, halos[i]->theta, halos[i]->phi, v_r);*/
	}
	return sum*sum;
}

double Fluc_mass_range(halo *halos, int num_halos, double D, double min_mass, double max_mass)
{
	double sum = 0;
	int i;
	vector my_pos = {0, 0, D}, diff;
	double natural_fluc = 2 * pi* pow(D, 3.0/2.0)/sqrt(MFreeNFW(D) * G) * 1/(PhiFreeNFW(D));
	for(i = 0; i < num_halos; i++)
	{
		double R = halos[i].R;
		double r = sqrt(R*R + D*D - 2.0*R*D*halos[i].cos_theta);
		if(halos[i].M < max_mass && halos[i].M > min_mass)
		{
			assign_vec_diff(diff, halos[i].position, my_pos);
			make_unit(diff)
			double v_r = dot_macro(halos[i].v, diff);

			double produced_fluc = ((r > CUTOFF_SCALE) ? (enclosed_mass(halos + i, r) * G /(r*r) * v_r) : 0);
			//sum += produced_fluc; //KILLEM
			// kill heating which is adiabtic

			if(produced_fluc/natural_fluc > 1 || ! ADIABATIC_CUTOFF)
				sum += produced_fluc;
			/*printf("%.3f\n", D);
			printf("velocity: ");
			print_vector(halos[i]->v);
			printf("position ");
			print_vector(halo_pos(halos[i]));
			printf("R = %.3f theta = %.3f phi = %.3f\n v_r = %.3f \n", halos[i]->R, halos[i]->theta, halos[i]->phi, v_r);*/
		}
	}
	return sum*sum;
}

double H_Density(halo *halos, int num_halos, double D, double dD)
{
	double sum = 0;
	int i;
	for(i = 0; i < num_halos; i++)
	{
		if(halos[i].R < D && halos[i].R > D - dD)
			sum += halos[i].M;
	}
	return sum/(4.0*pi*D*D*dD);
}

int print_out = 0;

void init(int argc, char **argv)
{
	const gsl_rng_type *T;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	RNG = gsl_rng_alloc (T);
	w = gsl_integration_workspace_alloc (1000);

	if(argc >= 2)
		num_trials = (int)pow(10, atoi(argv[1]));
	else
		num_trials = (int)pow(10, default_log_num_trials);

	if(argc > 2)
		print_out = 1;

	R_max_prim = MaxRadius(M_prim);
	c_prim = c_bar(M_prim);
	R_core_prim = R_max_prim/c_prim;
	m_max = f*M_prim;
	m_min = g*M_prim;
	expected_num = (2.0-p)/(1.0-p)*f*(pow(m_max/M_prim, 1.0-p) - pow(m_min/M_prim, 1.0-p))
		/ (pow(m_max/M_prim, 2.0-p) - pow(m_min/M_prim, 2.0-p));
}

double std_err_mean(double sum_of_squares)
{
	return sqrt(sum_of_squares/((double) num_trials*(num_trials - 1)));
}

void sq_root_data(data_cell *data, int num)
{
	int i;
	for(i = 0; i < num; i++)
	{
		data[i].m = sqrt(data[i].m);
		data[i].s *= pow(1.0/(2.0*data[i].m), 2);
		data[i].max_val = sqrt(data[i].max_val);
	}
}

void print_to_file(char *name, double *Ds, data_cell *data)
{
	int j;
	char path[100], filename[50], snum[10], addendum[10];
	strcpy(path, "data_files/");
	strcpy(filename, name);
	strcpy(addendum, ".txt");
	snprintf(snum, 10, "%d_%d", (int)gsl_rng_default_seed, (int)log10(num_trials));
	strcat(path, filename);
	strcat(path, snum);
	strcat(path, addendum);

	FILE *f = NULL;

	if(!print_out)
		f = fopen(path, "w");

	if(!print_out && f != NULL)
	{
		for(j = 0; j < num_points; j++)
		{
			double sig = std_err_mean(data[j].s)/data[j].m/log(10);
			printf("%f, %f, %f, %f\n", data[j].m, data[j].s, sig, log10(data[j].m));
			if(data[j].m > 0)
				fprintf(f, "%f %f %f %f\n", log10(Ds[j]), log10(data[j].m), sig, log10(data[j].max_val));
		}
		fclose(f);
	}
	else
	{
		for(j = 0; j < num_points; j++)
		{
			double sig = std_err_mean(data[j].s)/data[j].m/log(10);
			printf("%f, %f, %f, %f\n", data[j].m, data[j].s, sig, log10(data[j].m));
			if(data[j].m > 0)
				printf("%f : %f +/- %f\n", log10(Ds[j]), log10(data[j].m), sig);
		}
	}
}

int main(int argc, char **argv)
{
	init(argc, argv);
	unsigned int max_allowed_halos = (int)(10*expected_num);

	double Ds[num_points];
	data_cell Flucs[num_points], Dens[num_points];
	halo *halolist = malloc(max_allowed_halos*sizeof(halo));

	double mass = 0, avg_mass = 0;
	int i,j;

	for(i = 0; i < num_points; i++)
	{
		Ds[i] = pow(10, 5.0*i/(double) num_points);
		reset_cell(&Flucs[i]);
		reset_cell(&Dens[i]);
	}
	int hist_num = 10;

	data_cell Hist_Flucs[hist_num];
	for(int i = 0; i < hist_num; i++)
		reset_cell(Hist_Flucs + i);

	for(i = 0; i < num_trials; i++)
	{
		if(i % 100 == 0)printf("%d out of %d \n", i/100, num_trials/100);
		unsigned int num_halos = gsl_ran_poisson(RNG, expected_num);
		if(num_halos > max_allowed_halos)
		{
			max_allowed_halos *= 2;
			free(halolist);
			halolist = malloc(max_allowed_halos*sizeof(halo));
			printf("ALLOCATING\n");
		}
		mass = 0;
		//printf("Expect: %f   Have: %lu \n", expected_num, num_halos);

		for(j = 0; j < num_halos; j++)
		{
			make_halo(halolist + j);
			//print_halo_basic(halolist + j);
			mass += halolist[j].M;

			if(mass != mass)
			{
				print_halo_basic(halolist + j);
				printf("%f %f %d %f\n", mass, halolist[j].M, j, halolist[j].alpha);
				return -1;
			}
		}

		for(j = 0; j < num_points; j++)
		{
			update_cell(Flucs + j, Fluc(halolist, num_halos, Ds[j]), i + 1);
			update_cell(Dens + j, (j == 0 ? H_Density(halolist, num_halos, Ds[j], Ds[j]) : H_Density(halolist, num_halos, Ds[j], Ds[j] - Ds[j-1])), i + 1);
		}
		//printf("Mass frac: %f \n", mass/M_prim);
		avg_mass += mass/num_trials;

		for(int j = 0; j < hist_num; j++)
			update_cell(Hist_Flucs + j, Fluc_mass_range(halolist, num_halos, 1e4, M_prim*g*pow(f/g, (double) j/ (double) hist_num), M_prim*g*pow(f/g, (double) (j + 1)/ (double) hist_num)), i + 1);

		/*for(j = 0; j < num_points; j++)
		{
			//printf("%f : %f\n", pow(10, 5.0*j/num_points), log10(Flucs[j]));
			//printf("%f : %f\n", pow(10, 5.0*j/num_points), log10(Dens[j]));
			//printf("%f : %f\n", pow(10, 5.0*j/num_points), Hist[j]);
		}*/

	}
	print_to_file("Density", Ds, Dens);
	sq_root_data(Flucs, num_points);
	print_to_file("Flucs", Ds, Flucs);

	printf("Mass frac: %f \n", avg_mass/M_prim);
	sq_root_data(Hist_Flucs, hist_num);
	for(int i = 0; i < hist_num; i++)
		printf("Mass: %f Fluc: %f SDEV: %f\n", M_prim*g*pow(f/g, (double) i/ (double) hist_num), Hist_Flucs[i].m, std_err_mean(Hist_Flucs[i].s));
	gsl_rng_free (RNG);
	gsl_integration_workspace_free (w);
	return 0;
}

ext: c

max_stars_repo_name: benvchurch/project-eva

max_issues_repo_name: benvchurch/project-eva

max_forks_repo_name: benvchurch/project-eva


text:
#include <ceed.h>
#include <petsc.h>
#include "../problems/neo-hookean.h"

// Build libCEED context object
PetscErrorCode PhysicsContext_NH(MPI_Comm comm, Ceed ceed, Units *units,
                                 CeedQFunctionContext *ctx) {
  PetscErrorCode ierr;
  Physics_NH phys;

  PetscFunctionBegin;

  ierr = PetscMalloc1(1, units); CHKERRQ(ierr);
  ierr = PetscMalloc1(1, &phys); CHKERRQ(ierr);
  ierr = ProcessPhysics_NH(comm, phys, *units); CHKERRQ(ierr);
  CeedQFunctionContextCreate(ceed, ctx);
  CeedQFunctionContextSetData(*ctx, CEED_MEM_HOST, CEED_COPY_VALUES,
                              sizeof(*phys), phys);
  ierr = PetscFree(phys); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Build libCEED smoother context object
PetscErrorCode PhysicsSmootherContext_NH(MPI_Comm comm, Ceed ceed,
    CeedQFunctionContext ctx, CeedQFunctionContext *ctx_smoother) {
  PetscErrorCode ierr;
  PetscScalar nu_smoother = 0;
  PetscBool nu_flag = PETSC_FALSE;
  Physics_NH phys, phys_smoother;

  PetscFunctionBegin;

  ierr = PetscOptionsBegin(comm, NULL,
                           "Neo-Hookean physical parameters for smoother", NULL);
  CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-nu_smoother", "Poisson's ratio for smoother",
                            NULL, nu_smoother, &nu_smoother, &nu_flag);
  CHKERRQ(ierr);

  ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics

  if (nu_flag) {
    // Copy context
    CeedQFunctionContextGetData(ctx, CEED_MEM_HOST, &phys);
    ierr = PetscMalloc1(1, &phys_smoother); CHKERRQ(ierr);
    ierr = PetscMemcpy(phys_smoother, phys, sizeof(*phys)); CHKERRQ(ierr);
    CeedQFunctionContextRestoreData(ctx, &phys);
    // Create smoother context
    CeedQFunctionContextCreate(ceed, ctx_smoother);
    phys_smoother->nu = nu_smoother;
    CeedQFunctionContextSetData(*ctx_smoother, CEED_MEM_HOST, CEED_COPY_VALUES,
                                sizeof(*phys_smoother), phys_smoother);
    ierr = PetscFree(phys_smoother); CHKERRQ(ierr);
  } else {
    *ctx_smoother = NULL;
  }

  PetscFunctionReturn(0);
}

// Process physics options - Neo-Hookean
PetscErrorCode ProcessPhysics_NH(MPI_Comm comm, Physics_NH phys, Units units) {
  PetscErrorCode ierr;
  PetscBool nu_flag = PETSC_FALSE;
  PetscBool Young_flag = PETSC_FALSE;
  phys->nu = 0;
  phys->E = 0;
  units->meter     = 1;        // 1 meter in scaled length units
  units->second    = 1;        // 1 second in scaled time units
  units->kilogram  = 1;        // 1 kilogram in scaled mass units

  PetscFunctionBeginUser;

  ierr = PetscOptionsBegin(comm, NULL, "Neo-Hookean physical parameters", NULL);
  CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-nu", "Poisson's ratio", NULL, phys->nu, &phys->nu,
                            &nu_flag); CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-E", "Young's Modulus", NULL, phys->E, &phys->E,
                            &Young_flag); CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units",
                            NULL, units->meter, &units->meter, NULL);
  CHKERRQ(ierr);
  units->meter = fabs(units->meter);

  ierr = PetscOptionsScalar("-units_second", "1 second in scaled time units",
                            NULL, units->second, &units->second, NULL);
  CHKERRQ(ierr);
  units->second = fabs(units->second);

  ierr = PetscOptionsScalar("-units_kilogram", "1 kilogram in scaled mass units",
                            NULL, units->kilogram, &units->kilogram, NULL);
  CHKERRQ(ierr);
  units->kilogram = fabs(units->kilogram);

  ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics

  // Check for all required options to be set
  if (!nu_flag) {
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "-nu option needed");
  }
  if (!Young_flag) {
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "-E option needed");
  }

  // Define derived units
  units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second));

  // Scale E to Pa
  phys->E *= units->Pascal;

  PetscFunctionReturn(0);
};

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#ifdef USING_FFTW

#include <fftw.h>

#elif defined USING_FFTW3

#include <fftw3.h>

#else

typedef double fftw_real;
typedef struct {
     fftw_real re, im;
} fftw_complex_orig;
typedef fftw_real HPCC_Complex[2];
typedef HPCC_Complex fftw_complex;
#endif
typedef enum {
     FFTW2_FORWARD = -1, FFTW2_BACKWARD = 1
} fftw_direction;

struct hpcc_fftw_plan_struct {
  fftw_complex *w1, *w2, *ww1, *ww2, *ww3, *ww4, *c, *d;
  int n, c_size, d_size;
  int flags;
  fftw_direction dir;
};
typedef struct hpcc_fftw_plan_struct *hpcc_fftw_plan;

extern hpcc_fftw_plan HPCC_fftw_create_plan(int n, fftw_direction dir, int flags);
extern void HPCC_fftw_destroy_plan(hpcc_fftw_plan plan);
extern void HPCC_fftw_one(hpcc_fftw_plan plan, fftw_complex *in, fftw_complex *out);

#ifdef USING_FFTW

#elif defined USING_FFTW3

#define c_re(c)  ((c)[0])
#define c_im(c)  ((c)[1])

#else

typedef struct hpcc_fftw_plan_struct *fftw_plan;

#define c_re(c)  ((c)[0])
#define c_im(c)  ((c)[1])

#define fftw_malloc malloc
#define fftw_free free
/* flags for the planner */
#define  FFTW_ESTIMATE (0)
#define  FFTW_MEASURE  (1)

#define FFTW_OUT_OF_PLACE (0)
#define FFTW_IN_PLACE (8)
#define FFTW_USE_WISDOM (16)

#define fftw_create_plan HPCC_fftw_create_plan
#define fftw_destroy_plan HPCC_fftw_destroy_plan
#define fftw_one HPCC_fftw_one

#endif

ext: h