text:
/* multifit/gsl_multifit.h
 *
 * Copyright (C) 2000, 2007, 2010 Brian Gough
 * Copyright (C) 2013, 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.
 */

#ifndef __GSL_MULTIFIT_H__
#define __GSL_MULTIFIT_H__

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t nmax;         /* maximum number of observations */
  size_t pmax;         /* maximum number of parameters */
  size_t n;            /* number of observations in current SVD decomposition */
  size_t p;            /* number of parameters in current SVD decomposition */
  gsl_matrix * A;      /* least squares matrix for SVD, n-by-p */
  gsl_matrix * Q;
  gsl_matrix * QSI;
  gsl_vector * S;
  gsl_vector * t;
  gsl_vector * xt;
  gsl_vector * D;
  double rcond;        /* reciprocal condition number */
}
gsl_multifit_linear_workspace;

gsl_multifit_linear_workspace *
gsl_multifit_linear_alloc (const size_t n, const size_t p);

void
gsl_multifit_linear_free (gsl_multifit_linear_workspace * w);

int
gsl_multifit_linear (const gsl_matrix * X,
                     const gsl_vector * y,
                     gsl_vector * c,
                     gsl_matrix * cov,
                     double * chisq,
                     gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_tsvd (const gsl_matrix * X,
                          const gsl_vector * y,
                          const double tol,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double * chisq,
                          size_t * rank,
                          gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_svd (const gsl_matrix * X,
                         gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_bsvd (const gsl_matrix * X,
                          gsl_multifit_linear_workspace * work);

size_t
gsl_multifit_linear_rank(const double tol, const gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_solve (const double lambda,
                           const gsl_matrix * X,
                           const gsl_vector * y,
                           gsl_vector * c,
                           double *rnorm,
                           double *snorm,
                           gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_applyW(const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           gsl_matrix * WX,
                           gsl_vector * Wy);

int
gsl_multifit_linear_stdform1 (const gsl_vector * L,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_wstdform1 (const gsl_vector * L,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau);

int
gsl_multifit_linear_stdform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_matrix * M,
                              gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_matrix * M,
                               gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_genform1 (const gsl_vector * L,
                              const gsl_vector * cs,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_genform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              const gsl_vector * cs,
                              const gsl_matrix * M,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               const gsl_vector * cs,
                               const gsl_matrix * M,
                               gsl_vector * c,
                               gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_lreg (const double smin, const double smax,
                          gsl_vector * reg_param);

int
gsl_multifit_linear_lcurve (const gsl_vector * y,
                            gsl_vector * reg_param,
                            gsl_vector * rho, gsl_vector * eta,
                            gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_lcorner(const gsl_vector *rho,
                            const gsl_vector *eta,
                            size_t *idx);

int
gsl_multifit_linear_lcorner2(const gsl_vector *reg_param,
                             const gsl_vector *eta,
                             size_t *idx);

int
gsl_multifit_linear_Lk(const size_t p, const size_t k, gsl_matrix *L);

int
gsl_multifit_linear_Lsobolev(const size_t p, const size_t kmax,
                             const gsl_vector *alpha, gsl_matrix *L,
                             gsl_multifit_linear_workspace *work);

int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double * chisq,
                      gsl_multifit_linear_workspace * work);

int
gsl_multifit_wlinear_tsvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           const double tol,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double * chisq,
                           size_t * rank,
                           gsl_multifit_linear_workspace * work);

int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq,
                          gsl_multifit_linear_workspace * work);

int
gsl_multifit_wlinear_usvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           double tol,
                           size_t * rank,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double *chisq,
                           gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_est (const gsl_vector * x,
                         const gsl_vector * c,
                         const gsl_matrix * cov, double *y, double *y_err);

double
gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * w);

int
gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
                               const gsl_vector *c, gsl_vector *r);

/* gcv.c */
int
gsl_multifit_linear_gcv_init(const gsl_vector * y,
                             gsl_vector * reg_param,
                             gsl_vector * UTy,
                             double * delta0,
                             gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_gcv_curve(const gsl_vector * reg_param,
                              const gsl_vector * UTy,
                              const double delta0,
                              gsl_vector * G,
                              gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_gcv_min(const gsl_vector * reg_param,
                            const gsl_vector * UTy,
                            const gsl_vector * G,
                            const double delta0,
                            double * lambda,
                            gsl_multifit_linear_workspace * work);

double
gsl_multifit_linear_gcv_calc(const double lambda,
                             const gsl_vector * UTy,
                             const double delta0,
                             gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_gcv(const gsl_vector * y,
                        gsl_vector * reg_param,
                        gsl_vector * G,
                        double * lambda,
                        double * G_lambda,
                        gsl_multifit_linear_workspace * work);

typedef struct
{
  const char * name;     /* method name */
  int (*wfun)(const gsl_vector *r, gsl_vector *w);
  int (*psi_deriv)(const gsl_vector *r, gsl_vector *dpsi);
  double tuning_default; /* default tuning constant */
} gsl_multifit_robust_type;

typedef struct
{
  double sigma_ols;    /* OLS estimate of sigma */
  double sigma_mad;    /* MAD estimate of sigma */
  double sigma_rob;    /* robust estimate of sigma */
  double sigma;        /* final estimate of sigma */
  double Rsq;          /* R^2 coefficient of determination */
  double adj_Rsq;      /* degree of freedom adjusted R^2 */
  double rmse;         /* root mean squared error */
  double sse;          /* residual sum of squares */
  size_t dof;          /* degrees of freedom */
  size_t numit;        /* number of iterations */
  gsl_vector *weights; /* final weights */
  gsl_vector *r;       /* final residuals y - X c */
} gsl_multifit_robust_stats;

typedef struct
{
  size_t n;            /* number of observations */
  size_t p;            /* number of parameters */
  size_t numit;        /* number of iterations */
  size_t maxiter;      /* maximum iterations */
  const gsl_multifit_robust_type *type;
  double tune;         /* tuning parameter */

  gsl_vector *r;       /* residuals at current iteration */
  gsl_vector *weights; /* weights at current iteration */
  gsl_vector *c_prev;  /* coefficients from previous iteration */
  gsl_vector *resfac;  /* multiplicative factors for residuals */

  gsl_vector *psi;     /* psi(r) */
  gsl_vector *dpsi;    /* psi'(r) */

  gsl_matrix *QSI;     /* Q S^{-1} of original matrix X */
  gsl_vector *D;       /* balancing parameters of original matrix X */

  gsl_vector *workn;   /* workspace of length n */

  gsl_multifit_robust_stats stats; /* various statistics */

  gsl_multifit_linear_workspace *multifit_p;
} gsl_multifit_robust_workspace;

/* available types */
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_default;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_bisquare;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_cauchy;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_fair;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_huber;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_ols;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_welsch;

gsl_multifit_robust_workspace *gsl_multifit_robust_alloc(const gsl_multifit_robust_type *T,
                                                         const size_t n, const size_t p);
void gsl_multifit_robust_free(gsl_multifit_robust_workspace *w);
int gsl_multifit_robust_tune(const double tune,
                             gsl_multifit_robust_workspace *w);
int gsl_multifit_robust_maxiter(const size_t maxiter,
                                gsl_multifit_robust_workspace *w);
const char *gsl_multifit_robust_name(const gsl_multifit_robust_workspace *w);
gsl_multifit_robust_stats gsl_multifit_robust_statistics(const gsl_multifit_robust_workspace *w);
int gsl_multifit_robust_weights(const gsl_vector *r, gsl_vector *wts,
                                gsl_multifit_robust_workspace *w);
int gsl_multifit_robust(const gsl_matrix * X, const gsl_vector * y,
                        gsl_vector * c, gsl_matrix *cov,
                        gsl_multifit_robust_workspace *w);
int gsl_multifit_robust_est(const gsl_vector * x, const gsl_vector * c,
                            const gsl_matrix * cov, double *y, double *y_err);
int gsl_multifit_robust_residuals(const gsl_matrix * X,
                                  const gsl_vector * y,
                                  const gsl_vector * c, gsl_vector * r,
                                  gsl_multifit_robust_workspace * w);

__END_DECLS

#endif /* __GSL_MULTIFIT_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include "Watershed.h"
#include "math_functions.h"
#include "Globals.h"
#include "constitutive_equations.h"
#include "dgshed_numerical_flux.h"
#include "manufactured_solution.h"

/*****************************************************************************************//**
 * @file computeL.c
 *
 * This file contains code to evaluate the right hand side of the following discrete ODE obtained
 * from the DG discretization of the 2-D shallow water equations:
 *
 * \f$\frac{\partial \hat{\mathbf{w}}}{\partial t} = M^{-1}L(\hat{\mathbf{w}},t)\f$
 *
 * *******************************************************************************************/

/***************************************************************************/
/************ Compute the spatial discretization of the equation ***********/
/***************************************************************************/
extern void GetLGLWeights(int N, double *w);
extern double receive_floodwater(struct TwoDRegion *floodplain, int edg);

void computeChanL(struct channel* Chan, double time, double *dt, int channelNumber, int stage, double* RHSA, double* RHSQ)
{
	int fluxType = 1;
	int NumEdges = Chan->NumEdges;
	int NumEl = Chan->NumEl;
	int NumNodes = Chan->NumNodes;
	int Np = Chan->Np;
	int P = Chan->P;

#ifdef WDON
	/**** Compute the mass over an element to use later for ensuring ********/
	/**** positive mass with the flux in a wetting and drying treatment *************/
	double* mass = xcalloc(NumEl, sizeof(double));
	double LGLWeight[Np];
	GetLGLWeights(P, LGLWeight);

	for (int i = 0; i < NumEl; ++i)
	{
		double avgArea = 0;
		int begNode = i*Np + 1;
		for (int j =0; j < Np; j++)
		{
			avgArea += LGLWeight[j]*Chan->A[begNode+j];
		}
		mass[i] = avgArea;
	}
#endif

	/************** Compute the numerical flux at the faces *****************/
	double* Fhat1L = xcalloc(NumEdges, sizeof(double));
	double* Fhat2L = xcalloc(NumEdges, sizeof(double));
	double* Fhat1R = xcalloc(NumEdges, sizeof(double));
	double* Fhat2R = xcalloc(NumEdges, sizeof(double));

	for (int i=0; i < NumEdges; ++i)
	{
		double A_L, Q_L, A_R, Q_R, m1val, m2val;

		int leftNode = i*Np;
		int rightNode = leftNode+1;

		m1val = Chan->m1[i];
		m2val = Chan->m2[i];
		double bval = Chan->b[i];

		double tmpF[2];

#ifdef WDON
		// Check to see if the elements separated by this boundary are both dry
		if (i > 0 && i < NumEdges-1)
		{
			if (Chan->WD[i-1] == 0 && Chan->WD[i] == 0)
			{
				// Reflection flux for the left element
				A_L = Chan->A[leftNode];
				A_R = A_L;
				Q_L = Chan->Q[leftNode];
				Q_R = -Q_L;

				if (fluxType == 1)
					RoeFluxChan(tmpF, A_L, A_R, Q_L, Q_R,bval, m1val, m2val,0);
				else if (fluxType == 2)
					LFChan(tmpF, A_L, A_R, Q_L, Q_R,bval, m1val, m2val, 0);

				else
				{
					printf("Unknown flux type. Exiting now \n");
					exit(1);
				}

				Fhat1L[i] = tmpF[0];
				Fhat2L[i] = tmpF[1];

				// Reflection flux for the right element
				A_R = Chan->A[rightNode];
				A_L = A_R;
				Q_R = Chan->Q[rightNode];
				Q_L = -Q_R;

				if (fluxType == 1)
					RoeFluxChan(tmpF, A_L, A_R, Q_L, Q_R, bval, m1val, m2val, 0);
				else if (fluxType == 2)
					LFChan(tmpF, A_L, A_R, Q_L, Q_R, bval, m1val, m2val, 0);
				else
				{
					printf("Unknown flux type. Exiting now \n");
					exit(1);
				}
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];

				if (isnan(tmpF[0]) || isnan(tmpF[1]))
				{
					printf("both elements dry flux not a number, edge %d \n",i);
					exit(EXIT_FAILURE);
				}

			}

		}

		// if the elements are not both dry
		if ((i == 0) || (i == NumEl) || Chan->WD[i-1] == 1 || Chan->WD[i] == 1)
		{
			A_L = Chan->A[leftNode];
			A_R = Chan->A[rightNode];

			Q_L = Chan->Q[leftNode];
			Q_R = Chan->Q[rightNode];

			if (fluxType == 1)
				RoeFluxChan(tmpF, A_L, A_R, Q_L,Q_R,bval, m1val, m2val, g);
			else if (fluxType == 2)
				LFChan(tmpF, A_L, A_R, Q_L,Q_R,bval, m1val, m2val, g);

			else
			{
				printf("Unknown flux type. Exiting now \n");
				exit(1);
			}

			if (isnan(tmpF[0]) || isnan(tmpF[1]))
			{
				printf("A_L = %lf A_R = %lf Q_L = %lf Q_R = %lf\n", A_L, A_R, Q_L, Q_R);
				printf(" both elements wet flux not a number, edge %d, channel %d , time = %lf\n",i, channelNumber, time);
				exit(EXIT_FAILURE);
			}

			if (i==0 || Chan->WD[i-1] == 1)
			{
				Fhat1L[i] = tmpF[0];
				Fhat2L[i] = tmpF[1];

			}
			else
			{
				double newtmpF[2];
				if (fluxType == 1)
					RoeFluxChan(newtmpF, A_L, A_R, Q_L, Q_R, bval, m1val, m2val, 0);
				else if (fluxType == 2)
					LFChan(newtmpF, A_L, A_R, Q_L, Q_R, bval, m1val, m2val, 0);
				else
				{
					printf("Unknown flux type. Exiting now \n");
					exit(1);
				}
				Fhat1L[i] = newtmpF[0];
				Fhat2L[i] = newtmpF[1];
				if (isnan(newtmpF[0]) || isnan(newtmpF[1]))
				{
					printf("left element dry flux not a number, edge %d \n",i);
					exit(EXIT_FAILURE);
				}

			}
			if (i == NumEl || Chan->WD[i]==1)
			{
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];
			}
			else
			{
				if (fluxType == 1)
					RoeFluxChan(tmpF, A_L,A_R,Q_L,Q_R,bval, m1val, m2val, 0);
				else if (fluxType == 2)
					LFChan(tmpF, A_L,A_R,Q_L,Q_R,bval, m1val, m2val, 0);
				else
				{
					printf("Unknown flux type. Exiting now \n");
					exit(1);
				}
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];
				if (isnan(tmpF[0]) || isnan(tmpF[1]))
				{
					printf("right element dry flux not a number, edge %d \n",i);
					printf("A_L = %lf A_R = %lf Q_L = %lf Q_R = %lf bval = %lf m1val = %lf m2val = %lf\n", A_L, A_R, Q_L, Q_R, bval, m1val, m2val);
					printf("F = %lf\n", tmpF[0]);
					exit(EXIT_FAILURE);
				}

			}
		}

#else
		A_L = Chan->A[leftNode];
		A_R = Chan->A[rightNode];

		Q_L = Chan->Q[leftNode];
		Q_R = Chan->Q[rightNode];

		if (fluxType == 1)
			RoeFluxChan(tmpF, A_L, A_R, Q_L,Q_R,bval,m1val,m2val,g);
		else if (fluxType == 2)
			LFChan(tmpF, A_L, A_R, Q_L,Q_R,bval,m1val,m2val,g);
		else
		{
			printf("Unknown flux type. Exiting now \n");
			exit(1);
		}

		if (isnan(tmpF[0]) || isnan(tmpF[1]))
		{
			printf("flux not a number, edge %d \n",i);

			printf("A_L = %3.16f A_R = %3.16f Q_L = %3.16f Q_R = %31.6f, b = %lf, m1 = %lf, m2 = %lf \n", A_L, A_R, Q_L, Q_R, Chan->b[i], m1val, m2val);

			printf("Fhat = %lf Ghat = %lf\n", tmpF[0], tmpF[1]);
			exit(EXIT_FAILURE);
		}

		Fhat1L[i] = tmpF[0];
		Fhat2L[i] = tmpF[1];

		Fhat1R[i] = tmpF[0];
		Fhat2R[i] = tmpF[1];

#endif

		/**************************************************************************************************/
		double u_L = Q_L/A_L;
		double u_R = Q_R/A_R;
		//approximate c_L (exact for rectangular channels, but not for trapezoidal
		double c_L = sqrt(g*A_L/Chan->b[i]);
		double c_R = sqrt(g*A_R/Chan->b[i]);

		if (i==0)
			Chan->max_lambda = fmax((fabs(u_L)+c_L), (fabs(u_R)+c_R));
		// Compute maximum eigenvalue for the next time step
		double current_max =fmax((fabs(u_L) + c_L), (fabs(u_R) + c_R));
		Chan->max_lambda = max(Chan->max_lambda,current_max);
	}

#ifdef WDON
	for (int i =1; i < NumEdges-1; ++i)
	{
		int leftNode = i*Np;
		int rightNode = leftNode+1;
		double maxBetaOverAlpha = 2;
		// Check to see if this flux might possibly result in negative mass
		// If the mass will be negative on the left side
		//int cont = 1;
		//while(cont)
		//{
			if (Fhat1L[i]*maxBetaOverAlpha*(*dt) > mass[i-1] || -Fhat1R[i]*maxBetaOverAlpha*(*dt) > mass[i] )
			{
				Fhat1L[i] = 0.0;
				Fhat2L[i] = 0.0;
				Fhat1R[i] = 0.0;
				Fhat2R[i] = 0.0;

				//double A_L = Chan->A[leftNode];
				//double A_R = A_L;
				//double Q_L = Chan->Q[leftNode];
				//double Q_R = -Q_L;
				//double bval = Chan->NodalB[leftNode];
				//double m1val = Chan->Nodalm1[leftNode];
				//double m2val = Chan->Nodalm2[leftNode];

				//double tmpF[2];
				//double localG = g;

				//if (Chan->WD[i] == 0)
				//	localG = 0;

				//if (fluxType == 1)
				//	RoeFluxChan(tmpF, A_L, A_R, Q_L, Q_R, bval,m1val,m2val,localG);
				//else if (fluxType == 2)
				//	LFChan(tmpF, A_L, A_R, Q_L, Q_R, bval,m1val,m2val,localG);
				//else
				//{
				//	printf("Unknown flux type. Exiting now \n");
				//	exit(1);
				//}
				//Fhat1L[i] = tmpF[0];
				//Fhat2L[i] = tmpF[1];
				//printf("element %d of Channel %d will be dry. Reducing dt now.\n", i, channelNumber);
				//(*dt) = 0.5*(*dt);
			}
			//else
			//	cont = 0;
		//}

		//cont = 1;
		//while(cont)
		//{
		//	// If the mass will be negative on the right side
		//	if (-Fhat1R[i]*maxBetaOverAlpha*(*dt) > mass[i])
		//	{
		//		double A_R = Chan->A[rightNode];
		//		double A_L = A_R;
		//		double Q_R = Chan->Q[rightNode];
		//		double Q_L = -Q_R;
		//		double bval = Chan->NodalB[rightNode];
		//		double m1val = Chan->Nodalm1[rightNode];
		//		double m2val = Chan->Nodalm2[rightNode];

		//		double tmpF[2];
		//		double localG = g;
		//		if (Chan->WD[i] == 0)
		//			localG = 0;

		//		if (fluxType == 1)
		//			RoeFluxChan(tmpF, A_L, A_R, Q_L, Q_R, bval, m1val, m2val,localG);
		//		else if (fluxType == 2)
		//			LFChan(tmpF, A_L, A_R, Q_L, Q_R, bval,m1val,m2val,localG);
		//		else
		//		{
		//			printf("Unknown flux type. Exiting now \n");
		//			exit(1);
		//		}

		//		Fhat1R[i] = tmpF[0];
		//		Fhat2R[i] = tmpF[1];
		//		printf("element %d of channel %d will be dry. Decreasing dt now.\n",i, channelNumber);

		//		*dt = 0.5*(*dt);

		//	}
		//	else
		//		cont = 0;
		//}
	}
	free(mass);
#endif

	for (int k=0; k < NumEl; ++k)
	{

		double h = Chan->dh[k];
		gsl_vector *F1 = gsl_vector_alloc(Np);
		gsl_vector *F2 = gsl_vector_alloc(Np);
		gsl_vector *ST21 = gsl_vector_alloc(Np);
		gsl_vector *ST22 = gsl_vector_alloc(Np);
		gsl_vector *ST23 = gsl_vector_alloc(Np);

		//ST11 contains the qL term that comes from overland flow
		gsl_vector *ST11 = gsl_vector_calloc(Np);

		// ST12 contains term that takes flooding into account
		gsl_vector *ST12 = gsl_vector_calloc(Np);

		int begNode = k*Np;
		for (int i = 0; i < Np; i++)
		{
			double Aval = Chan->A[begNode+i+1];
			double Qval = Chan->Q[begNode+i+1];
			double qL = Chan->qL[begNode+i];
			double qM = Chan->qM[begNode+i];
			double bval = Chan->NodalB[begNode+i];
			double S0 = Chan->dz[begNode+i];
			double m1val = Chan->Nodalm1[begNode+i];
			double m2val = Chan->Nodalm2[begNode+i];
			double dm1val = Chan->dm1[begNode+i];
			double dm2val = Chan->dm2[begNode+i];
			double dbval = Chan->db[begNode+i];
			double nval = Chan->NodalnFriction[begNode+i];
			double I1val = getI1(Aval, bval, m1val, m2val);
			double I2val = getI2(Aval, bval, dbval, m1val, dm1val, m2val, dm2val);
			double Sfval = getS_f(Aval, Qval, bval, m1val, m2val, nval);
			double beta = Chan->beta[begNode+i];
			double localG;
#ifdef WDON
			if(Chan->WD[k]==1)
				localG = g;
			else
				localG = 0;
#else
			localG = g;
#endif

			gsl_vector_set(F1, i, Qval);
			gsl_vector_set(F2, i, beta*Qval*Qval/Aval + localG*I1val);
			gsl_vector_set(ST21, i, localG*I2val);
			gsl_vector_set(ST22, i, localG*Aval*S0);
			gsl_vector_set(ST23, i, localG*Aval*Sfval);
			gsl_vector_set(ST11, i, qL);
			gsl_vector_set(ST12, i, qM);

		}

		//gsl_blas_dgemv(CblasNoTrans, 2.0/h, InvM, ST, 1.0, ST11);

		gsl_vector *localFhat1 = gsl_vector_alloc(2);
		gsl_vector_set(localFhat1, 0, Fhat1R[k]);
		gsl_vector_set(localFhat1, 1, Fhat1L[k+1]);

		gsl_vector *localFhat2 = gsl_vector_alloc(2);
		gsl_vector_set(localFhat2, 0, Fhat2R[k]);
		gsl_vector_set(localFhat2, 1, Fhat2L[k+1]);

		// cacluate the volume integral of the flux
		gsl_vector *localRHS1 = gsl_vector_calloc(Np);
		gsl_vector *localRHS2 = gsl_vector_calloc(Np);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, VolMat, F1, 1.0, localRHS1);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, VolMat, F2, 1.0, localRHS2);

		// calculate the surface integral of the flux
		gsl_vector *SurfPart1 = gsl_vector_calloc(Np);
		gsl_vector *SurfPart2 = gsl_vector_calloc(Np);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, LIFT, localFhat1, 1.0, SurfPart1);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, LIFT, localFhat2, 1.0, SurfPart2);

		// calculate the RHS
		gsl_vector_add(localRHS1, SurfPart1);
		gsl_vector_add(localRHS1, ST11);
		gsl_vector_add(localRHS1, ST12);
		gsl_vector_add(localRHS2, SurfPart2);
		gsl_vector_add(localRHS2, ST21);
		gsl_vector_add(localRHS2, ST22);
		gsl_vector_sub(localRHS2, ST23);

		for(int i = 0; i < Np; i++)
		{
			RHSA[begNode+i+1] = gsl_vector_get(localRHS1, i);
			RHSQ[begNode+i+1] = gsl_vector_get(localRHS2, i);
		}

		gsl_vector_free(F1);
		gsl_vector_free(F2);
		gsl_vector_free(localFhat1);
		gsl_vector_free(localFhat2);
		gsl_vector_free(SurfPart1);
		gsl_vector_free(SurfPart2);
		gsl_vector_free(ST11);
		gsl_vector_free(ST12);
		gsl_vector_free(ST21);
		gsl_vector_free(ST22);
		gsl_vector_free(ST23);
		gsl_vector_free(localRHS1);
		gsl_vector_free(localRHS2);

	}

	free(Fhat1L);
	free(Fhat1R);
	free(Fhat2L);
	free(Fhat2R);

}

void computeLKinematicEls(double time, double* RHS, int fp)
{
	int NumEl = FloodplainList[fp]->NumEl;
	// calculate and store the flux (Q) at the downstream edge of each kinematic element
	// these fluxes will be used to calculate the upwinded flux values
	double *Fhat = xcalloc(NumEl, sizeof(double));
	for (int i = 0; i < NumEl; ++i)
	{
		struct kinematicEl *kinEl = KinematicElList[i];
		if (kinEl->isActive)
		{
			int Np = kinEl->Np;
			double myA = kinEl->A[Np-1];
			double myS0 = kinEl->dz[Np-1];
			double myNf = kinEl->NodalnFriction[Np-1];
			double myWeq = kinEl->weq;
			double myH = myA/myWeq;

			// using manning's N.
			Fhat[i] = myWeq*sqrt(myS0)*pow(myH, 5.0/3)/myNf;

			// using chezy's relationship. nf stands for C, chezy's coefficient
			//Fhat[i] = myWeq*sqrt(myS0*myH)*myH*myNf;
		}
	}

	// now calculate the right hand side
	for (int i = 0; i < NumEl; ++i)
	{
		struct kinematicEl *kinEl = KinematicElList[i];
		if (kinEl->isActive)
		{
			double h_el = kinEl->dh;
			int Np = kinEl->Np;

			gsl_vector *F = gsl_vector_calloc(Np);

			gsl_vector *ST = gsl_vector_calloc(Np);			 // will eventually store R-I (Rainfall - Infiltration)

			double weq = kinEl->weq;

			for (int j = 0; j < Np; ++j)
			{
				double A = kinEl->A[j];
				double nf = kinEl->NodalnFriction[j];
				double S0 = kinEl->dz[j];
				double H = A/weq;

				// using manning's N
				double myF = weq*sqrt(S0)*pow(H,5.0/3)/nf;

				// using Chezy's relationship (nf stands for C)
				//double myF = weq*sqrt(S0*H)*H*nf;

				if (S0 < 0 )
					printf("z1 = %lf z2 = %lf\n", kinEl->z1, kinEl->z2);
				if (isnan(myF))
				{
					printf("myF is nan; A = %lf, nf = %lf, S0 = %lf, H = %lf weq = %lf, kinElNum = %d\n", A, nf, S0, H, weq, i);
					exit(1);
				}
				gsl_vector_set(F, j, myF);

				double rainfall;
				if (time < 1800)
					rainfall = 1.0/12/3600;

				gsl_vector_set(ST, j, weq*rainfall);
				// second parking lot case
				//if (time  <= 180)
				//	gsl_vector_set(ST, j, weq*2.0/12/3600);
				//else if (time <= 360)
				//	gsl_vector_set(ST, j, weq*4.0/12/3600);

			}
			// Now calculate total upstream flux.
			// If the element is at the boundary of the watershed, flux coming in will be 0.
			int numUpstreamEls = kinEl->numUpstreamEls;
			double Fhat1 = 0.0;
			//if (numUpstreamEls > 1)
			//	printf(" i = %d UpEl1 = %d UpEl2 = %d\n", i, kinEl->upstreamEls[0], kinEl->upstreamEls[1]);
			for (int j = 0; j < numUpstreamEls; ++j)
			{
				int upstreamEl = kinEl->upstreamEls[j];
				//printf("upstream el = %d\t", upstreamEl);
				Fhat1 += Fhat[upstreamEl];
			}
			//printf("downstream el = %d\n", kinEl->el2);
			double Fhat2 = Fhat[i];
			gsl_vector *localFhat = gsl_vector_alloc(2);
			gsl_vector_set(localFhat, 0, Fhat1);
			gsl_vector_set(localFhat, 1, Fhat2);

			gsl_vector *localRHS = gsl_vector_calloc(Np);

			// calculate the volume integral
			gsl_blas_dgemv(CblasNoTrans, 2.0/h_el, VolMat, F, 1.0, localRHS);

			// calculate the surface integral
			gsl_vector *SurfPart = gsl_vector_calloc(Np);
			gsl_blas_dgemv(CblasNoTrans, 2.0/h_el, LIFT, localFhat, 1.0, SurfPart);

			// subtract the surface integral from the volume integral
			gsl_vector_add(localRHS, SurfPart);

			// add source term to the right hand side

			gsl_vector_add(localRHS, ST);

			for (int j = 0; j < Np; j++)
			{
				RHS[i*Np+j] = gsl_vector_get(localRHS, j);
				if (isnan(RHS[i*Np+j]))
				{
					printf("SurfPart = %lf \n", gsl_vector_get(SurfPart,j));
				}
			}
			gsl_vector_free(F);
			gsl_vector_free(localFhat);
			gsl_vector_free(SurfPart);
			gsl_vector_free(localRHS);
			gsl_vector_free(ST);

		}

	}
	free(Fhat);

}

/***********************************************************************************************//**
 * Function for evaluating the right hand side of the discrete ODE obtained from the DG discretization
 * of the 2-D shallow water equations
 * @param [in] currRegion a pointer to the junction structure corresponding to the junction that is
 * currently being worked on
 * @param [in] time  a double representing the current time of the simulation
 * @param [out] RHSZeta a pointer to an array of size NumEl x 3 in which the right hand side for
 * the water surface elevation will be stored
 * @param [out] RHSQx a pointer to an array of size NumEl x 3 in which the right hand side for the
 * momentum in the x-direction is stored
 * @param [out] RHSQy a pointer to an array of size NumEl x 3 in which the right hand side for the
 * momentum in the y-direction is stored
 *
 ***************************************************************************************************/

void compute2DL(struct TwoDRegion *currRegion, double time, double *RHSZeta, double *RHSQx, double *RHSQy, double dt)
{

	int NumEdges = currRegion->TotalNumEdges;
	int NumEl = currRegion->NumEl;
	int Nfp = currRegion->P + 1;
	int Np = currRegion->Np;

#ifdef WDON
	/**** Compute the mass over an element to use later for ensuring ********/
	/**** positive mass with the flux in a wetting and drying treatment *************/
	double* mass = xcalloc(NumEl, sizeof(double));
	gsl_vector *height = gsl_vector_calloc(Np);
	gsl_vector *tmp = gsl_vector_alloc(Np);
	gsl_vector *jac = gsl_vector_alloc(Np);

	for (int i = 0; i < NumEl; ++i)
	{
		for (int j = 0; j < Np; j++)
		{
			gsl_vector_set(height, j, currRegion->zeta[i][j] + currRegion->NodalZ[i][j]);
			gsl_vector_set(jac, j, currRegion->jac[i]);
		}

		gsl_blas_dgemv(CblasNoTrans, 1.0, MassMatrix2D, height, 0.0, tmp);
		gsl_blas_ddot(jac, tmp, &mass[i]);

	}

	gsl_vector_free(height);
	gsl_vector_free(tmp);
	gsl_vector_free(jac);
#endif

	/********************************* Compute Roe's flux at the face **********************************************/
	double **Fhat1dotn, **Fhat2dotn, **Fhat3dotn;
	Fhat1dotn = malloc(NumEl*sizeof(double*));
	Fhat2dotn = malloc(NumEl*sizeof(double*));
	Fhat3dotn = malloc(NumEl*sizeof(double*));
	for (int i = 0; i < NumEl; i++)
	{
		Fhat1dotn[i] = xcalloc(3*Nfp, sizeof(double));
		Fhat2dotn[i] = xcalloc(3*Nfp, sizeof(double));
		Fhat3dotn[i] = xcalloc(3*Nfp, sizeof(double));
	}

	for (int i=0; i<NumEdges; ++i)
	{
		int el1 = currRegion->EdgtoEls[i*2];
		int el2 = currRegion->EdgtoEls[i*2+1];

		int ledg1 = currRegion->GlobaltoLocalEdg[i*2];
		int ledg2 = currRegion->GlobaltoLocalEdg[i*2+1];

		// take the normal from element 1 side
		double nx = currRegion->nx[el1*3+ledg1];
		double ny = currRegion->ny[el1*3+ledg1];

		double tx = -ny;
		double ty = nx;

		for (int j = 0; j < Nfp; j++)
		{
			double zeta_in, Qx_in, Qy_in;
			double zeta_ex, Qx_ex, Qy_ex;
			double z_in, z_ex;			// should be the same as z_ex

			int lv1 = currRegion->GlobalEdgPosNegNodes[i][2*j];
			int lv2 = currRegion->GlobalEdgPosNegNodes[i][2*j+1];

			int pos1 = currRegion->PosInFVec[i][2*j];
			int pos2 = currRegion->PosInFVec[i][2*j+1];

			double Fn_in[3], Fn_ex[3];

			zeta_in = currRegion->zeta[el1][lv1];
			Qx_in = currRegion->Qx[el1][lv1];
			Qy_in = currRegion->Qy[el1][lv1];
			z_in = currRegion->NodalZ[el1][lv1];

			double Q_T_in = Qx_in*tx + Qy_in*ty;

			// need to change this later so that the boundary condition is given per node on the edge
			// if the edge is an exterior edge connected to a channel
			int bdrypres = currRegion->BdryPrescribed[i];

		//if (el1 == 531 || el2 == 531)
		//{
		//	printf("el1 = %d, el2 = %d\n", el1, el2);
		//	printf("WD1 = %d, WD2 = %d\n", currRegion->WD[el1], currRegion->WD[el2]);
		//	printf("bdrypres = %d\n", bdrypres);
		//}

			// manufactured solution
			if (bdrypres == 555)
				//if (bdrypres == 555 || (el1 == el2))
			{
				double xval = currRegion->NodalX[el1][lv1];
				double yval = currRegion->NodalY[el1][lv1];
				zeta_ex = getmanH(xval, yval, time);
				Qx_ex = getQx(xval, yval, time);
				Qy_ex = getQy(xval, yval, time);
			}
			// if inflow or outflow boundary, i.e. if bdrypres = 1 or 2
			//for floodplains
			else if (bdrypres == 111)
			{
				zeta_ex = currRegion->bzeta[i];
				Qx_ex = Qx_in;
				Qy_ex = Qy_in;
			}

			else if (bdrypres == 222)
			{
				zeta_ex = zeta_in;
				double Q_N_ex = currRegion->bQn[i];
				double Q_T_ex = Q_T_in;
				double denom = 1./(nx*ty-ny*tx);
				Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;
			}

			else if (bdrypres == 333)
			{
				zeta_ex = currRegion->bQn[i];
				double Q_N_ex = currRegion->bQn[i];
				double Q_T_ex = Q_T_in;
				double denom = 1./(nx*ty-ny*tx);
				Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;

			}
			// for junctions
			else if (bdrypres == 1 || bdrypres == 2)
			{
				zeta_ex = currRegion->bzeta[i];
				double Q_N_ex = currRegion->bQn[i];
				// for the first iteration
				if (zeta_ex == -10000.0 && Q_N_ex == -10000.0)
				{
					zeta_ex = zeta_in;
					Qx_ex = Qx_in;
					Qy_ex = Qy_in;
					z_ex = z_in;
				}
				else
				{
					double Q_T_ex = Q_T_in;
					double denom = 1./(nx*ty-ny*tx);
					Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
					Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;
					z_ex = z_in;
				}
			}

			// for floodplains
			else if (bdrypres == 3)
			{
				double Q_N_ex = receive_floodwater(currRegion, i);

				// element 1
				Qx_in = currRegion->Qx[el1][lv1];
				Qy_in = currRegion->Qy[el1][lv1];
				zeta_in = currRegion->zeta[el1][lv1];
				z_in = currRegion->NodalZ[el1][lv1];
				zeta_ex = zeta_in;
				Q_T_in = Qx_in*tx + Qy_in*ty;
				double Q_T_ex = Q_T_in;
				double denom = 1./(nx*ty-ny*tx);
				Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;

				double localG = g;

#ifdef WDON
				if (currRegion->WD[el1] == 0)
					localG = 0;
				else
					localG = g;
#endif

				double current_max_lam= RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex, z_in, nx, ny, localG,Fn_in);

				//if (el1 == 170 && fabs(Q_N_ex) > 0)
				//{
				//	printf("***************\n");
				//	printf("time = %lf\n", time);
				//	printf("j = %d el1 = %d, el2 = %d\n", j, el1, el2);
				//	printf("WD1 = %d, WD2 = %d\n", currRegion->WD[el1], currRegion->WD[el2]);
				//	printf("Q_N_ex = %lf, Qx_in = %lf, Qy_in = %lf, \n", Q_N_ex, Qx_in, Qy_in);
				//	printf("Fn_in = %lf, %lf, %lf\n", Fn_in[0], Fn_in[1], Fn_in[2]);
				//}

				// element 2
				nx = -nx;
				ny = -ny;
				tx = -ny;
				ty = nx;
				denom = 1./(nx*ty-ny*tx);
				zeta_in = currRegion->zeta[el2][lv2];
				z_in = currRegion->NodalZ[el2][lv2];
				zeta_ex = zeta_in;
				Qx_in = currRegion->Qx[el2][lv2];
				Qy_in = currRegion->Qy[el2][lv2];
				Q_T_in = Qx_in*tx + Qy_in*ty;
				Q_T_ex = Q_T_in;
				denom = 1./(nx*ty-ny*tx);
				Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;
				z_ex = z_in;

#ifdef WDON
				if (currRegion->WD[el2] == 0)
					localG = 0;
				else
					localG = g;
#endif

				current_max_lam= RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex, z_in, nx, ny, localG, Fn_ex);

				if (el2 == 170 && fabs(Q_N_ex)>0)
				{
					printf("***************\n");
					printf("el1 = %d, el2 = %d\n", el1, el2);

					printf("Q_N_ex = %lf\n", Q_N_ex);
					printf("Fn_in = %lf, %lf, %lf\n", Fn_in[0], Fn_in[1], Fn_ex[2]);
				}
			}

			// if edge i is a boundary edge but not connected to a channel, implement no flux boundary condition
			else if (el1 == el2)
			{
				zeta_ex  = currRegion->zeta[el2][lv2];
				// Compute the velocity in the normal direction
				double Q_N_in = Qx_in*nx + Qy_in*ny;
				double Q_T_in = Qx_in*tx + Qy_in*ty;

				// Reflect the velocity in the normal direction
				double Q_N_ex = -Q_N_in;
				double Q_T_ex = Q_T_in;

				// Compute the x and y components of the external state flow
				double denom = 1./(nx*ty - ny*tx);
				Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;

				z_ex = z_in;

			}

			// if the edge is not a boundary edge
			else
			{
				zeta_ex = currRegion->zeta[el2][lv2];
				Qx_ex = currRegion->Qx[el2][lv2];
				Qy_ex = currRegion->Qy[el2][lv2];
				z_ex = currRegion->NodalZ[el2][lv2];
			}

#ifdef WDON
			// Check to see if both of the elements separated by this edge are dry
			if (currRegion->WD[el1] == 0 && currRegion->WD[el2] == 0 && (bdrypres != 3) && (bdrypres != 222) && (bdrypres != 2) && (bdrypres !=1))
			{
				Fn_in[0] = 0.0;
				Fn_in[1] = 0.0;
				Fn_in[2] = 0.0;
				Fn_ex[0] = 0.0;
				Fn_ex[1] = 0.0;
				Fn_ex[2] = 0.0;

				//// Reflection flux for the interior element
				//double zeta_ex_ref = zeta_in;
				//double Q_N_in = Qx_in*nx + Qy_in*ny;
				//double Q_T_in = Qx_in*tx + Qy_in*ty;

				//double Q_N_ex = -Q_N_in;
				//double Q_T_ex = Q_T_in;
				//
				//double denom = 1./(nx*ty - ny*tx);
				//double Qx_ex_ref = (ty*Q_N_ex - ny*Q_T_ex)*denom;
				//double Qy_ex_ref = (-tx*Q_N_ex + nx*Q_T_ex)*denom;

				//double max_lam_in = RoeFluxJunc(zeta_in, zeta_ex_ref, Qx_in,Qx_ex_ref,
				//	Qy_in, Qy_ex_ref, z_in, nx, ny, 0, Fn_in);

				//// Reflection flux for the exterior element
				//double zeta_in_ref = zeta_ex;
				//Q_N_ex = Qx_ex*nx + Qy_ex*ny;
				//Q_T_ex = Qx_ex*tx + Qy_ex*ty;

				//Q_N_in = -Q_N_ex;
				//Q_T_in = Q_T_ex;
				//
				//double Qx_in_ref = (ty*Q_N_in - ny*Q_T_in)*denom;
				//double Qy_in_ref = (-tx*Q_N_in + nx*Q_T_ex)*denom;
				//double max_lam_ex = RoeFluxJunc(zeta_in_ref, zeta_ex, Qx_in_ref, Qx_ex, Qy_in_ref, Qy_ex, z_in, nx, ny, 0, Fn_ex);

				//if (i ==0)
				//	currRegion->max_lambda = max(max_lam_in, max_lam_ex);
				//else
				//{
				//	currRegion->max_lambda = fmax(currRegion->max_lambda, max_lam_in);
				//	currRegion->max_lambda = fmax(currRegion->max_lambda, max_lam_ex);
				//}

			}
			else if (bdrypres != 3 && bdrypres != 222)// if the elements aren't both dry
			{
				double Fhatdotn[3];
				double current_max_lam= RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex, z_in, nx, ny, g,Fhatdotn);

				if (isnan(Fhatdotn[0]) || isnan(Fhatdotn[1]) || isnan(Fhatdotn[2]))
				{
					printf("2D both elements wet flux not a number, edge %d, time %e, for domain type %d \n",i,time, currRegion->type);
					printf("ze_ex = %lf, ze_in = %lf, z_in = %lf, Qx_ex = %lf, Qx_in = %lf, Qy_ex = %lf, Qy_in = %lf, el_ex = %d, el_in = %d, bdrypres = %d\n", zeta_ex, zeta_in, z_in, Qx_ex, Qx_in, Qy_ex, Qy_in, el2, el1, bdrypres);
					printf("wdfalg1 = %d , wdflag2 = %d\n", currRegion->WD[el1], currRegion->WD[el2]);
					exit(EXIT_FAILURE);
				}

				// Check to see if the flux is large enough to dry up the elements
				// Calculate the sum of the length of the three edges of the interior and exterior element
				double edg_len_el1 = 0;
				for (int v = 0; v < 3; v++)
				{
					int v1 = currRegion->EltoVert[el1*3 + v];
					int v2 = currRegion->EltoVert[el1*3 + (v+1)%3];
					edg_len_el1 += sqrt(pow(currRegion->Vx[v1] - currRegion->Vx[v2], 2) + pow(currRegion->Vy[v1] - currRegion->Vy[v2],2));
				}

				double edg_len_el2 = 0;
				for (int v = 0; v < 3; v++)
				{
					int v1 = currRegion->EltoVert[el2*3 + v];
					int v2 = currRegion->EltoVert[el2*3 + (v+1)%3];
					edg_len_el2 += sqrt(pow(currRegion->Vx[v1] - currRegion->Vx[v2], 2) + pow(currRegion->Vy[v1] - currRegion->Vy[v2],2));
				}

				double maxBetaOverAlpha = 2;

				if ((1.01*Fhatdotn[0]*edg_len_el1*maxBetaOverAlpha*dt >= mass[el1]) || (-1.01*Fhatdotn[0]*edg_len_el2*maxBetaOverAlpha*dt >= mass[el2]))
				{
					//printf("either el %d or el %d will be dry \n", el1, el2);
					Fn_in[0] = 0.0;
					Fn_in[1] = 0.0;
					Fn_in[2] = 0.0;
					Fn_ex[0] = 0.0;
					Fn_ex[1] = 0.0;
					Fn_ex[2] = 0.0;

				}

				else // check to make sure the flux is not coming from a dry element
				{
					if (currRegion->WD[el1] == 1 || (el1 == el2))
					{
						for (int j =0; j < 3; ++j)
						{
							Fn_in[j] = Fhatdotn[j];
						}
					}
					else if (Fhatdotn[0] > 0)
					{
						//printf("initially, el1 = %d, j = %d, Fn_ex = %lf\n", el1, j, Fhatdotn[0]);
						zeta_ex = currRegion->zeta[el2][lv2];
						zeta_in = zeta_ex;
						Qx_ex = currRegion->Qx[el2][lv2];
						Qy_ex = currRegion->Qy[el2][lv2];
						double Q_N_ex = Qx_ex*nx + Qy_ex*ny;
						double Q_N_in = -Q_N_ex;
						double Q_T_ex = Qx_ex*tx + Qy_ex*ty;
						Q_T_in = Q_T_ex;
						double denom = 1./(nx*ty-ny*tx);
						Qx_in = (ty*Q_N_in - ny*Q_T_in)*denom;
						Qy_in = (-tx*Q_N_in + nx*Q_T_in)*denom;

						double max_lam_in = RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex,
								z_ex, nx, ny, g, Fn_in);
						if (isnan(Fn_in[0]) || isnan(Fn_in[1]) || isnan(Fn_in[2]))
						{
							printf("2D interior element dry flux not a number, edge %d, el_in = %d, el_ex = %d, domain type = %d\n",i, el1, el2, currRegion->type);
							printf("ze_ex = %lf, ze_in = %lf, z_in = %lf, z_ex = %lf, Qx_ex = %lf, Qx_in = %lf, Qy_ex = %lf, Qy_in = %lf, el_ex = %d, el_in = %d, bdrypres = %d\n", zeta_ex, zeta_in, z_in, z_ex, Qx_ex, Qx_in, Qy_ex, Qy_in, el2, el1, bdrypres);
							printf("wdfalg1 = %d , wdflag2 = %d\n", currRegion->WD[el1], currRegion->WD[el2]);
							exit(EXIT_FAILURE);
						}
						//printf("after reflection, el1 = %d, Fn_ex = %lf\n", el1, Fn_in[0]);
						current_max_lam = max(current_max_lam, max_lam_in);
					}

					else
					{
						double max_lam_in = RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex,
								z_in, nx, ny, 0, Fn_in);

					}

					if (currRegion->WD[el2] == 1 && (el1 != el2))
					{
						for (int j = 0; j < 3; ++j)
						{
							Fn_ex[j] = Fhatdotn[j];
						}
					}
					else if ((el1 != el2) && (Fhatdotn[0] < 0))
					{
					//	printf("initially, el2 = %d, Fn_ex = %lf\n", el2, -Fhatdotn[0]);
						Qx_in = currRegion->Qx[el1][lv1];
						Qy_in = currRegion->Qy[el1][lv1];
						zeta_in = currRegion->zeta[el1][lv1];
						double Q_N_in = Qx_in*nx + Qy_in*ny;
						double Q_T_in = Qx_in*tx + Qy_in*ty;
						double Q_N_ex = -Q_N_in;
						double Q_T_ex = Q_T_in;
						double denom = 1.0/(nx*ty - ny*tx);
						Qx_ex = (ty*Q_N_ex - ny*Q_T_ex)*denom;
						Qy_ex = (-tx*Q_N_ex + nx*Q_T_ex)*denom;
						zeta_ex = zeta_in;

						double max_lam_ex = RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex,
								Qy_in, Qy_ex, z_in, nx, ny, g, Fn_ex);
						current_max_lam = max(current_max_lam, max_lam_ex);

						if (isnan(Fn_ex[0]) || isnan(Fn_ex[1]) || isnan(Fn_ex[2]))
						{
							printf("2D exterior element dry flux not a number, edge %d \n",i);
							exit(EXIT_FAILURE);
						}
						//printf("el2 = %d, Fn_ex = %lf\n", el2, -Fn_ex[0]);

					//	printf("after reflection, el2 = %d, Fn_ex = %lf\n", el2, -Fn_ex[0]);
					}
					else if (el1 != el2 )
					{
						double max_lam_ex = RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex,
								Qy_in, Qy_ex, z_in, nx, ny, 0, Fn_ex);
						current_max_lam = max(current_max_lam, max_lam_ex);

					}

					if (i == 0)
						currRegion->max_lambda = current_max_lam;
					else
						currRegion->max_lambda = fmax(currRegion->max_lambda, current_max_lam);
				}

			}

#else

			double Fhatdotn[3];
			double current_max_lam= RoeFluxJunc(zeta_in, zeta_ex, Qx_in, Qx_ex, Qy_in, Qy_ex, z_in, nx, ny, g, Fhatdotn);

			if (isnan(Fhatdotn[0]) || isnan(Fhatdotn[1]) || isnan(Fhatdotn[2]))
			{
				printf("2D numerical flux not a number, edge %d \n",i);
				printf("zeta_in = %e, zeta_ex = %e, z_edge = %e, BdryPrescribed = %d \n", zeta_in, zeta_ex, z_in, bdrypres);
				exit(EXIT_FAILURE);
			}

			for (int j = 0; j < 3; ++j)
			{
				Fn_in[j] = Fhatdotn[j];
				Fn_ex[j] = Fhatdotn[j];
			}

			if (i == 0)
				currRegion->max_lambda = current_max_lam;
			else
				currRegion->max_lambda = max(currRegion->max_lambda, current_max_lam);

#endif

			// store Fhatdotn for the two elements connected by this edge
			Fhat1dotn[el1][pos1] = Fn_in[0];
			Fhat2dotn[el1][pos1] = Fn_in[1];
			Fhat3dotn[el1][pos1] = Fn_in[2];

			// store the value for the exterior element, only if the edge is not a boundary edge
			if (el2 != el1 && (bdrypres != 3))
			{
				Fhat1dotn[el2][pos2] = -Fn_ex[0];
				Fhat2dotn[el2][pos2] = -Fn_ex[1];
				Fhat3dotn[el2][pos2] = -Fn_ex[2];
			}
			else if (el2 != el1)  // for bdyrpes = 3, we have already taken sign of normals into account during flux calculation
			{
				Fhat1dotn[el2][pos2] = Fn_ex[0];
				Fhat2dotn[el2][pos2] = Fn_ex[1];
				Fhat3dotn[el2][pos2] = Fn_ex[2];
			}

		} // end node loop

	}	// end edge loop

	for (int k=0; k<NumEl; ++k)
	{
		gsl_vector *F1x = gsl_vector_alloc(Np);
		gsl_vector *F2x = gsl_vector_alloc(Np);
		gsl_vector *F3x = gsl_vector_alloc(Np);
		gsl_vector *F1y = gsl_vector_alloc(Np);
		gsl_vector *F2y = gsl_vector_alloc(Np);
		gsl_vector *F3y = gsl_vector_alloc(Np);
		gsl_vector *ST21 = gsl_vector_calloc(Np);
		gsl_vector *ST22 = gsl_vector_calloc(Np);
		gsl_vector *ST31 = gsl_vector_calloc(Np);
		gsl_vector *ST32 = gsl_vector_calloc(Np);
		//gsl_vector *ST11 = gsl_vector_calloc(Np);

		double localG = g;
#ifdef WDON
		if(currRegion->WD[k] == 0)
			localG = 0;
#endif

		for (int i = 0; i < Np; i++)
		{
			double zeta = currRegion->zeta[k][i];
			double Qx = currRegion->Qx[k][i];
			double Qy = currRegion->Qy[k][i];
			double z = currRegion->NodalZ[k][i];
			double H = zeta + z;
			double f2x = Qx*Qx/H + 0.5*localG*(H*H - z*z);
			double f2y = Qx*Qy/H;
			double f3x = f2y;
			double f3y = Qy*Qy/H + 0.5*localG*(H*H - z*z);
			gsl_vector_set(F1x, i, Qx);
			gsl_vector_set(F1y, i, Qy);
			gsl_vector_set(F2x, i, f2x);
			gsl_vector_set(F2y, i, f2y);
			gsl_vector_set(F3x, i, f3x);
			gsl_vector_set(F3y, i, f3y);

			double u = Qx/H;
			double v = Qy/H;
			double nfric = currRegion->NodalnFriction[k][i];
			double tau = localG*nfric*nfric*sqrt(u*u + v*v)/pow(H,4.0/3);
			double dzx = currRegion->Nodaldzx[k][i];
			double dzy = currRegion->Nodaldzy[k][i];

			gsl_vector_set(ST21, i, localG*zeta*dzx);
			gsl_vector_set(ST22, i, -tau*Qx);
			gsl_vector_set(ST31, i, localG*zeta*dzy);
			gsl_vector_set(ST32, i, -tau*Qy);

			// Source terms for manufactured solution
			//double xval = currRegion->NodalX[k][i];
			//double yval = currRegion->NodalY[k][i];

			//double st2 = getS2(xval, yval, time);
			//double st3 = getS3(xval, yval, time);
			//gsl_vector_set(ST21, i, st2);
			//gsl_vector_set(ST31, i, st3);

#ifdef WDON
			if (currRegion->WD[k] == 0)
			{
				gsl_vector_set(F1x, i, 0);
				gsl_vector_set(F1y, i, 0);
				gsl_vector_set(F2x, i, 0);
				gsl_vector_set(F2y, i, 0);
				gsl_vector_set(F3x, i, 0);
				gsl_vector_set(F3y, i, 0);
				gsl_vector_set(ST21, i, 0);
				gsl_vector_set(ST22, i, 0);
				gsl_vector_set(ST31, i, 0);
				gsl_vector_set(ST32, i, 0);
			}
#endif

		// rainfall
		//if (currRegion->type == 2)
		//	gsl_vector_set(ST11, i, 0.5*12/3600);

		}

		// Calculate the volume integral of the flux
		double rx = currRegion->rx[k];
		double ry = currRegion->ry[k];
		double sx = currRegion->sx[k];
		double sy = currRegion->sy[k];

		gsl_vector *localRHS1 = gsl_vector_calloc(Np);
		gsl_vector *localRHS2 = gsl_vector_calloc(Np);
		gsl_vector *localRHS3 = gsl_vector_calloc(Np);

		gsl_vector *lR12 = gsl_vector_calloc(Np);
		gsl_vector *lR13 = gsl_vector_calloc(Np);
		gsl_vector *lR14 = gsl_vector_calloc(Np);
		gsl_vector *lR22 = gsl_vector_calloc(Np);
		gsl_vector *lR23 = gsl_vector_calloc(Np);
		gsl_vector *lR24 = gsl_vector_calloc(Np);
		gsl_vector *lR32 = gsl_vector_calloc(Np);
		gsl_vector *lR33 = gsl_vector_calloc(Np);
		gsl_vector *lR34 = gsl_vector_calloc(Np);

		gsl_blas_dgemv(CblasNoTrans, rx, Drw, F1x, 0.0, localRHS1);
		gsl_blas_dgemv(CblasNoTrans, sx, Dsw, F1x, 0.0, lR12);
		gsl_blas_dgemv(CblasNoTrans, ry, Drw, F1y, 0.0, lR13);
		gsl_blas_dgemv(CblasNoTrans, sy, Dsw, F1y, 0.0, lR14);

		gsl_vector_add(localRHS1, lR12);
		gsl_vector_add(localRHS1, lR13);
		gsl_vector_add(localRHS1, lR14);

		gsl_blas_dgemv(CblasNoTrans, rx, Drw, F2x, 0.0, localRHS2);
		gsl_blas_dgemv(CblasNoTrans, sx, Dsw, F2x, 0.0, lR22);
		gsl_blas_dgemv(CblasNoTrans, ry, Drw, F2y, 0.0, lR23);
		gsl_blas_dgemv(CblasNoTrans, sy, Dsw, F2y, 0.0, lR24);

		gsl_vector_add(localRHS2, lR22);
		gsl_vector_add(localRHS2, lR23);
		gsl_vector_add(localRHS2, lR24);

		gsl_blas_dgemv(CblasNoTrans, rx, Drw, F3x, 0.0, localRHS3);
		gsl_blas_dgemv(CblasNoTrans, sx, Dsw, F3x, 0.0, lR32);
		gsl_blas_dgemv(CblasNoTrans, ry, Drw, F3y, 0.0, lR33);
		gsl_blas_dgemv(CblasNoTrans, sy, Dsw, F3y, 0.0, lR34);

		gsl_vector_add(localRHS3, lR32);
		gsl_vector_add(localRHS3, lR33);
		gsl_vector_add(localRHS3, lR34);

		// calculate the surface integral of the flux
		gsl_vector *localFhat1dotn = gsl_vector_alloc(3*Nfp);
		gsl_vector *localFhat2dotn = gsl_vector_alloc(3*Nfp);
		gsl_vector *localFhat3dotn = gsl_vector_alloc(3*Nfp);

		for (int i = 0; i < 3; i++)
		{
			double edgJac = currRegion->edgJac[k*3+i];
			for (int j = 0; j < Nfp; j++)
			{
				int index = i*Nfp+j;
				gsl_vector_set(localFhat1dotn, index, edgJac*Fhat1dotn[k][index]);
				gsl_vector_set(localFhat2dotn, index, edgJac*Fhat2dotn[k][index]);
				//printf("edgJac = %lf Fhat = %lf\n", edgJac, Fhat2dotn[k][index]);
				gsl_vector_set(localFhat3dotn, index, edgJac*Fhat3dotn[k][index]);

			}
		}

		double jac = currRegion->jac[k];

		gsl_vector *SurfPart1 = gsl_vector_calloc(Np);
		gsl_vector *SurfPart2 = gsl_vector_calloc(Np);
		gsl_vector *SurfPart3 = gsl_vector_calloc(Np);

		double fac = 1.0/jac;
		gsl_blas_dgemv(CblasNoTrans, fac, LIFT2D, localFhat1dotn, 0.0, SurfPart1);
		gsl_blas_dgemv(CblasNoTrans, fac, LIFT2D, localFhat2dotn, 0.0, SurfPart2);
		gsl_blas_dgemv(CblasNoTrans, fac, LIFT2D, localFhat3dotn, 0.0, SurfPart3);

		// subtract this from the RHS
		gsl_vector_sub(localRHS1, SurfPart1);
		gsl_vector_sub(localRHS2, SurfPart2);
		gsl_vector_sub(localRHS3, SurfPart3);

		// add source and sink terms
		//gsl_vector_add(localRHS1, ST11);
		//
		//gsl_vector_add(localRHS1, ST11);
		gsl_vector_add(localRHS2, ST21);
		gsl_vector_add(localRHS2, ST22);
		gsl_vector_add(localRHS3, ST31);
		gsl_vector_add(localRHS3, ST32);

		int begNode = k*Np;
		for(int i = 0; i < Np; i++)
		{
			RHSZeta[begNode+i] = gsl_vector_get(localRHS1,i);
			RHSQx[begNode+i] = gsl_vector_get(localRHS2,i);
			RHSQy[begNode+i] = gsl_vector_get(localRHS3, i);

		}

		// free all gsl vectors
		gsl_vector_free(localRHS1);
		gsl_vector_free(localRHS2);
		gsl_vector_free(localRHS3);
		gsl_vector_free(localFhat1dotn);
		gsl_vector_free(localFhat2dotn);
		gsl_vector_free(localFhat3dotn);
		gsl_vector_free(SurfPart1);
		gsl_vector_free(SurfPart2);
		gsl_vector_free(SurfPart3);
		gsl_vector_free(lR12);
		gsl_vector_free(lR13);
		gsl_vector_free(lR14);
		gsl_vector_free(lR22);
		gsl_vector_free(lR23);
		gsl_vector_free(lR24);
		gsl_vector_free(lR32);
		gsl_vector_free(lR33);
		gsl_vector_free(lR34);
		gsl_vector_free(F1x);
		gsl_vector_free(F2x);
		gsl_vector_free(F3x);
		gsl_vector_free(F1y);
		gsl_vector_free(F2y);
		gsl_vector_free(F3y);
		//gsl_vector_free(ST11);
		gsl_vector_free(ST21);
		gsl_vector_free(ST22);
		gsl_vector_free(ST31);
		gsl_vector_free(ST32);

	}

	// free allocated space
	for (int i = 0; i < NumEl; i++)
	{
		free(Fhat1dotn[i]);
		free(Fhat2dotn[i]);
		free(Fhat3dotn[i]);
	}
	free(Fhat1dotn);
	free(Fhat2dotn);
	free(Fhat3dotn);

#ifdef WDON
	free(mass);
#endif

}

ext: c

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
/**
 * author: Jochen K"upper
 * created: Jan 2002
 * file: pygsl/src/statistics/longmodule.c
 * $Id: ucharmodule.c,v 1.5 2004/03/24 08:40:45 schnizer Exp $
 *
 *"
 */

#include <Python.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_statistics.h>

/* include real functions for different data-types */

#define STATMOD_APPEND_PY_TYPE(X) X ## Int
#define STATMOD_APPEND_PYC_TYPE(X) X ## UBYTE
#define STATMOD_FUNC_EXT(X, Y) X ## _uchar ## Y
#define STATMOD_PY_AS_C PyInt_AsLong
#define STATMOD_C_TYPE unsigned char
#include "functions.c"

PyGSL_STATISTICS_INIT(uchar, "uchar")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "Stroustrup"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/*
* Copyright (c) 2014, Newcastle University, UK.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
*    this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
*    this list of conditions and the following disclaimer in the documentation
*    and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/

// Linear Regression
// Dan Jackson, 2014 (but see below)

#include <stdlib.h>
//#include <math.h>

#include "linearregression.h"

// Linear regression with one independent variable
double *LinearModelFitOneIndependent(int n, double *y, double *x)
{
	static double coef[3];		// offset(intersect), scale(gradient), spare(to be compatible with two-variable version)
	int i;

	// sum(Xi * Yi)
	// sum(Xi^2)
	double sumx = 0, sumy = 0;
	double sumxy = 0, sumxsq = 0;
	for (i = 0; i < n; i++)
	{
		sumx += x[i];
		sumy += y[i];
		sumxy += x[i] * y[i];
		sumxsq += x[i] * x[i];
	}

	// Calculate mean of x and y
	double xmean = n > 0 ? sumx / n : 0;
	double ymean = n > 0 ? sumy / n : 0;

	//                 sum(x*y) - (n * xm * ym)
	// Slope (b):  b = ------------------------
	//                  sum(x^2) - (n * xm^2)
	double bNumerator = sumxy - (n * xmean * ymean);
	double bDenominator = sumxsq - (n * xmean * xmean);
	double b = bDenominator != 0.0 ? bNumerator / bDenominator : 0.0;

	// Intercept (a):  a = ym - b * xm
	double a = ymean - (b * xmean);

	// Line of best fit: y = a + b * X
	coef[0] = a;
	coef[1] = b;
	coef[2] = 0.0;		// just to be compatible with two-variable version

	return &coef[0];
}

// Linear regression with two independent variables
double *LinearModelFitTwoIndependent(int n, double *y, double *x1, double *x2)
{
	// Implemented from information from: http://faculty.cas.usf.edu/mbrannick/regression/Reg2IV.html
	static double coef[3];		// offset (intersect), scale 1 (gradient 1), scale 2 (gradient 2)
	int i;

	double sumx1 = 0, sumx2 = 0, sumy = 0;
	double sumx1y = 0, sumx2y = 0;
	double sumx1x2 = 0;
	double sumx1sq = 0, sumx2sq = 0;
	for (i = 0; i < n; i++)
	{
		sumx1 += x1[i];
		sumx2 += x2[i];
		sumy += y[i];
		sumx1y += x1[i] * y[i];
		sumx2y += x2[i] * y[i];
		sumx1x2 += x1[i] * x2[i];
		sumx1sq += x1[i] * x1[i];
		sumx2sq += x2[i] * x2[i];
	}

	// Calculate mean of x1, x2, and y
	double x1mean = n > 0 ? sumx1 / n : 0;
	double x2mean = n > 0 ? sumx2 / n : 0;
	double ymean = n > 0 ? sumy / n : 0;

	//                   sum(x2^2).sum(x1*y) - sum(x1*x2).sum(x2*y)
	// Slope (b1):  b1 = ------------------------------------------
	//                       sum(x1^2).sum(x2^2) - sum(x1*x2)^2
	double b1Numerator = (sumx2sq * sumx1y) - (sumx1x2 * sumx2y);
	double b1Denominator = (sumx1sq * sumx2sq) - (sumx1x2 * sumx1x2);
	double b1 = b1Numerator != 0.0 ? b1Numerator / b1Denominator : 0.0;

	//                   sum(x1^2).sum(x2*y) - sum(x1*x2).sum(x1*y)
	// Slope (b2):  b2 = ------------------------------------------
	//                       sum(x1^2).sum(x2^2) - sum(x1*x2)^2
	double b2Numerator = (sumx1sq * sumx2y) - (sumx1x2 * sumx1y);
	double b2Denominator = b1Denominator;		// Same as denominator of b1
	double b2 = b2Denominator != 0.0 ? b2Numerator / b2Denominator : 0.0;

	// Intercept (a):  a = ym - b1 * x1m - b2 * x2m
	double a = ymean - (b1 * x1mean) - (b2 * x2mean);

	// Line of best fit: y = a + b1 * x1 + b2 * x2
	coef[0] = a;
	coef[1] = b1;
	coef[2] = b2;

	return &coef[0];
}

#ifdef ENABLE_APPROXIMATE
// Linear regression with two independent variables, weighted
// NOTE: This is not really correct...
double *LinearModelFitTwoIndependentWeightedApproximately(int n, double *y, double *x1, double *x2, double *weights)
{
	double *weightedY = (double *)malloc(sizeof(double) * n);
	double *weightedX1 = (double *)malloc(sizeof(double) * n);
	double *weightedX2 = (double *)malloc(sizeof(double) * n);
	int i;

	// HACK: To reduce influence on intercept
	double influenceSum = 0;
	for (i = 0; i < n; i++)
	{
		double sw = sqrt(weights[i]);
		weightedY[i] = y[i] * sw;
		weightedX1[i] = x1[i] * sw;
		weightedX2[i] = x2[i] * sw;
		influenceSum += sw;
	}

	double *coef = LinearModelFitTwoIndependent(n, weightedY, weightedX1, weightedX2);
	double influence = n != 0 ? influenceSum / n : 0;
	if (influence != 0.0) { coef[0] /= influence; }

	free(weightedY);
	free(weightedX1);
	free(weightedX2);
	return coef;
}
#endif

#ifdef ENABLE_GSL

#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit.h>
//#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

/* General weighted case in gsl/multifit/multilinear.c */

double *LinearModelFitTwoIndependentWeighted(int n, double *y, double *x1, double *x2, double *weights)
{
	#define NPARAMS 2

	// Allocate
	gsl_matrix *matrixX = gsl_matrix_alloc(n, NPARAMS);
	gsl_vector *vectorW = gsl_vector_alloc(n);
	gsl_vector *vectorY = gsl_vector_alloc(n);
	gsl_vector *vectorC = gsl_vector_alloc(1 + NPARAMS);

	// Copy data
	int i;
	for (i = 0; i < n; i++)
	{
		// NPARAMS
		gsl_matrix_set(matrixX, i, 0, x1[i]);
		gsl_matrix_set(matrixX, i, 1, x2[i]);
		gsl_vector_set(vectorW, i, weights[i]);
		gsl_vector_set(vectorY, i, y[i]);
	}

	// Compute
	gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n, NPARAMS);
	int ret = gsl_multifit_wlinear(matrixX, vectorW, vectorY, vectorC, /*gsl_matrix * cov*/ NULL, /*double * chisq*/ NULL, work);
	gsl_multifit_linear_free(work);

	static double coef[NPARAMS + 1];
	// NPARAMS+1
	coef[0] = gsl_vector_get(vectorC, 0);
	coef[1] = gsl_vector_get(vectorC, 1);
	coef[2] = gsl_vector_get(vectorC, 2);

	// Free
	gsl_matrix_free(matrixX);
	gsl_vector_free(vectorW);
	gsl_vector_free(vectorY);
	gsl_vector_free(vectorC);

	return coef;
}
#endif

max_stars_repo_name: digitalinteraction/openmovement

max_issues_repo_name: digitalinteraction/openmovement

max_forks_repo_name: digitalinteraction/openmovement


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "nrsrc/nrutil.h"

#include "prototypes.h"
#include "globvars.h"

double comp_DF_halo_exact(double E, double L);
double comp_DF_bulge_exact(double E, double L);

double comp_DF_H_iso(double E);
double comp_DF_H_ani(double Q);

double comp_DF_halo_Eddington(double Q);
double comp_DF_bulge_Eddington(double Q);

void comp_DF_init(void);
void compute_d2rhodpsi2_halo(void);
void compute_d2rhodpsi2_bulge(void);

double eddington_integrand(double t, void *params);

#define WORKSIZE 100000
gsl_integration_workspace *Workspace;

double *xi, *yi, *eddint;
double *list_radius;
double *list_E;                            /* tabulated energies */
double *DistFunc_halo, *DistFunc_bulge;    /* distribution function(s) */
double *psi_R;
double *rho_R_halo,*rho_R_bulge;
double *drhodpsi_halo,*d2rhodpsi2_halo;
double *drhodpsi_bulge,*d2rhodpsi2_bulge;

/* size of DF look-up table */
int DFSIZE= RSIZE * 2;

void compute_DF_lookuptable(void)
{
  FILE *fd;
  char dffile[100]="";
  int i;

  printf("Start computing distribution function.\n"); fflush(stdout);

#ifdef MAXWELLIAN
  AnisotropyRadius= -1;
  printf("MAXWELLIAN is turned on\n");
#else
  printf("MAXWELLIAN is turned off\n");
#endif

#ifdef DF_H_MODEL
  printf("DF_H_MODEL is turned on\n");
#else
  printf("DF_H_MODEL is turned off\n");
#endif

#ifdef DF_EDDINGTON
  printf("DF_EDDINGTON is turned on\n");
#else
  printf("DF_EDDINGTON is turned off\n");
#endif

  printf("AnisotropyRadius= %g\n", AnisotropyRadius);
  printf("ra= %g\n", AnisotropyRadius * RH);

  comp_DF_init();

  for(i = 0; i <= DFSIZE; i++)
    {
	printf("Computing DF table bin %d/%d (e= %g)\n", i, DFSIZE, -1.0*psi_R[i]); fflush(stdout);

	list_E[i]= -1.0 * psi_R[i];

	DistFunc_halo[i] = comp_DF_halo_exact(list_E[i], 0.0);
	DistFunc_bulge[i] = comp_DF_bulge_exact(list_E[i], 0.0);
    }
  printf("\n");

  /* write this to a file so we have a record of it */
  if(strstr(OutputFile,".hdf5"))
        strncpy(dffile, OutputFile, strlen(OutputFile)-5);
  if(strstr(OutputFile,".dat"))
        strncpy(dffile, OutputFile, strlen(OutputFile)-4);
  strcat(dffile, ".df");
  printf("DF look-up table saved in file: %s\n",dffile);
  if((fd = fopen(dffile,"w")))
    {
	fprintf(fd,"   Energy       f_halo(E)    f_bulge(E)\n");
	for(i= 0; i<= DFSIZE; i++)
	  fprintf(fd," %10.5e   %10.5e   %10.5e\n",list_E[i],DistFunc_halo[i],DistFunc_bulge[i]);
	fclose(fd);
    }

}

double comp_DF_halo_exact(double E, double L)
{
  double Q, f_e= 0.0;

  if((N_HALO == 0) || (M_HALO == 0))
    return 0;

  if (AnisotropyRadius > 0.0)
	Q= - E - (L*L/2./(AnisotropyRadius * RH)/(AnisotropyRadius * RH));
  else
	Q= - E;

#ifdef DF_H_MODEL
	/* calculate f_e exactly from the H90 analytic formula */
      if (AnisotropyRadius > 0.0)
        f_e= comp_DF_H_ani(Q);
      else
        f_e= comp_DF_H_iso(E);
#endif

#ifdef DF_EDDINGTON
	/* calculate f_e using Eddington's formula */
	f_e= comp_DF_halo_Eddington(Q);
	printf("f_e_edd(Q)/f_H_iso(E)= %g\n",f_e/comp_DF_H_iso(E));
	printf("f_e_edd(Q)/f_H_ani(Q)= %g\n",f_e/comp_DF_H_ani(Q)); fflush(stdout);
	printf("------\n");
#endif

  return f_e;
}

double comp_DF_bulge_exact(double E, double L)
{
  double Q, f_e= 0.0;

  if((N_BULGE == 0) || (M_BULGE == 0))
    return 0;

/* for the moment, we only assume isotropic DF and use Eddington's formula */

/*
  if (AnisotropyRadius > 0.0)
        Q= - E - (L*L/2./(AnisotropyRadius * RH)/(AnisotropyRadius * RH));
  else
*/
        Q= - E;

  f_e= comp_DF_bulge_Eddington(Q);

  return f_e;
}

/*
   interpolate the halo f_e from the pre-computed table
*/
double comp_DF_halo(double E, double L)
{
  double f_e= 0.0, Q, ee, ue;
  int ie;

  ee= E;

  if(AnisotropyRadius > 0.0)
    {
      Q= - E - (L*L/2./(AnisotropyRadius * RH)/(AnisotropyRadius * RH));

      /* in this case, table is stored in terms of -Q */
      ee= -1.0 * Q;
    }

  /* if e>0 or q<0 then dist. func. is zero */
  if(ee>0.0)
    return 0.0;

  ie= find_idx(ee, list_E, DFSIZE);

/*
  printf("\nE= %g  Q= %g ee= %g\n",E, Q, ee);
  printf("ie= %d  list_E[0]= %g   list_E[DFSIZE]= %g \n",ie, list_E[0], list_E[DFSIZE]);
  printf("list_E[%d]= %g\n", ie-1, list_E[ie-1]);
  printf("list_E[%d]= %g\n", ie, list_E[ie]);
  printf("list_E[%d]= %g\n", ie+1, list_E[ie+1]);
  fflush(stdout);
*/

  if(ie < 0 || ie > DFSIZE)
    ie= 0;

  ue = (ee - list_E[ie]) / (list_E[ie + 1] - list_E[ie]);

  if(ee>list_E[0])
    {
      /*
      printf("\nparticle has energy very close to zero\n");
      printf("ee=  %g, table min list_E[0]= %g\n",ee,list_E[0]);
      printf("return 0.\n");
      */
      ue= 0.0;
    }

  if(ie < 0 || ie >= DFSIZE || ue < 0 || ue > 1)
    {
      printf("\nee= %g   (in bounds? list_E[0]= %g  list_E[DFSIZE]= %g\n",ee, list_E[0], list_E[DFSIZE]);
      printf("in: comp_DF_fromtable\nerror: ie=%d out of range ue=%g\nstopping\n", ie, ue);
      exit(0);
    }

  f_e= DistFunc_halo[ie] * (1-ue) + DistFunc_halo[ie+1] * ue;

  return f_e;
}

/*
   interpolate the bulge f_e from the pre-computed table
*/
double comp_DF_bulge(double E, double L)
{
  double f_e= 0.0, Q, ee, ue;
  int ie;

  ee= E;

  /* for the moment, we'll assume the bulge is isotropic */
/*
  if(AnisotropyRadius > 0.0)
    {
      Q= - E - (L*L/2./(AnisotropyRadius * RH)/(AnisotropyRadius * RH));
      ee= -1.0 * Q;
    }
*/

  /* if e>0 or q<0 then dist. func. is zero */
  if(ee>0.0)
    return 0.0;

  ie= find_idx(ee, list_E, DFSIZE);

  if(ie < 0 || ie > DFSIZE)
    ie= 0;

  ue = (ee - list_E[ie]) / (list_E[ie + 1] - list_E[ie]);

  if(ee>list_E[0])
    {
      /*
      printf("\nparticle has energy very close to zero\n");
      printf("ee=  %g, table min list_E[0]= %g\n",ee,list_E[0]);
      printf("return 0.\n");
      */
      ue= 0.0;
    }

  if(ie < 0 || ie >= DFSIZE || ue < 0 || ue > 1)
    {
      printf("\nee= %g   (in bounds? list_E[0]= %g  list_E[DFSIZE]= %g\n",ee, list_E[0], list_E[DFSIZE]);
      printf("in: comp_DF_fromtable\nerror: ie=%d out of range ue=%g\nstopping\n", ie, ue);
      exit(0);
    }

  f_e= DistFunc_bulge[ie] * (1-ue) + DistFunc_bulge[ie+1] * ue;

  return f_e;
}

/* ----------------------------------------------------------------- */

/* the following are analytic forms of the distribution
   function, currently, it's only Hernquist 1990 */

/* Isotropic Hernquist Model */

double comp_DF_H_iso(double E)
{
  double vg, q, prefac, f_e;

  if (E >= 0.0)
	return 0.0;

  vg= sqrt(G * M_HALO / RH);
  q= sqrt(-1.0 * E * RH / G / M_HALO);

  if (q >= 1.0)
	return (1.0e+30);

  prefac= M_HALO / (8.0 * sqrt(2.0) * PI * PI * PI * RH * RH * RH * vg * vg * vg);

  f_e= pow((1-q*q), -5./2.) * ( 3.*asin(q) + q*pow((1.-q*q),0.5)*(1.-2.*q*q) * (8.*q*q*q*q - 8.*q*q - 3.));
/*
  printf("               E= %g  vg= %g  q= %g  prefac= %g f_e %g\n", E, vg, q, prefac, f_e); fflush(stdout);
*/

  return prefac * f_e;
}

/* Anisotropic Hernquist Model */

double comp_DF_H_ani(double Q)
{

  if (AnisotropyRadius <= 0.0)
	return 0;

  double vg, qbar, prefac, f_e, f_iso;
  double ra;

  ra= AnisotropyRadius * RH;

  if (Q <= 0.0)
	return 0;

  vg= sqrt(G * M_HALO / RH);
  qbar= sqrt(RH * Q / G / M_HALO);

  if (qbar >= 1.0)
	return 1.0e+30;

  prefac= M_HALO / (sqrt(2.0) * PI * PI * PI * RH * RH * RH * vg * vg * vg);

  f_e= (RH * RH / ra / ra) * qbar * (1. - 2.*qbar*qbar);

  f_iso= comp_DF_H_iso(-1.0*Q);

/*
printf("DF_ani  Q= %g   qbar= %g  f_iso= %g  prefac= %g  f_e= %g\n",Q,qbar,f_iso,prefac,f_e); fflush(stdout);
*/

  return f_iso + prefac * f_e;
}

/* -----------------------------------------------------------------

     Now, calculate DF using the Eddington formula.

   ----------------------------------------------------------------- */

gsl_spline * d2rhodpsi2_spline_halo;
gsl_interp_accel * d2rhodpsi2_spline_acc_halo;

gsl_spline * d2rhodpsi2_spline_bulge;
gsl_interp_accel * d2rhodpsi2_spline_acc_bulge;

// Parameter struct for the Eddington integration function.
typedef struct {
  const gsl_spline* s;
  gsl_interp_accel* a;
  double Q;
} spline_params;

/* The integrand for the Eddington integration. Now uses
   the GSL spline functionality rather than the NR, which is evil. */
double eddington_integrand(double t, void *params)
{
  spline_params* p = params;

//printf("t= %g\t", t);
//printf("Q= %g\t",p->Q);
//printf("sqrt(p->Q-t)= %g   \t", sqrt(p->Q-t));
//printf("d2rhodpsi2= %g  \t", gsl_spline_eval(p->s, t, p->a));
  double v2 = gsl_spline_eval(p->s, t, p->a) / sqrt(p->Q - t);
//printf("v2= %g  \n",v2); fflush(stdout);

  return v2;
}

double comp_DF_halo_Eddington(double Q)
{
  double f_e;
  int i;
  double result, abserr;

  /* first, we'll generate the spline for d2rho/dpsi2 term in integrand */

  if(Q<=0)
    return 0.0;

  /* load spline for d2rho/dpsi2 term into function parameters */
  spline_params p;
  p.s = d2rhodpsi2_spline_halo;
  p.a = d2rhodpsi2_spline_acc_halo;
  p.Q = Q;

  /* error checking: write entire DF integrand and the spline */
/*
  FILE *fd;
  char tempfile[100]="";
  sprintf(tempfile, "df_integrands/%g.txt", Q);
  if((fd = fopen(tempfile,"w")))
    {
        fprintf(fd,"Q= %g \n",Q);
        fprintf(fd,"       xi            yi          d2rho_dpsi2            Q           Q-psiR     sqrt(Q-psiR) \n");
        for(i = 0; i <= DFSIZE; i++)
          fprintf(fd," %8.5e     %8.5e     %8.5e    %8.5e    %8.5e      %8.5e     \n",xi[i],yi[i],d2rho_dpsi2[i], Q, Q-psi_R[i], sqrt(Q-psi_R[i]));
        fclose(fd);
    }
*/

  gsl_function F;
  F.function = &eddington_integrand;
  F.params = &p;

  double epsabs=0;
  double epsrel=1e-4;
  int intstatus;

  /* do our own error handling, so that we can adjust the accuracy if needed (i.e., at small r) */
  gsl_error_handler_t * old_handler = gsl_set_error_handler_off(); /* turn off the error-handler */

  do
  {
    intstatus= gsl_integration_qags(&F, 0.0, Q, epsabs, epsrel, WORKSIZE, Workspace, &result, &abserr);
    printf("epsrel = %g     result = %g +/- %g   (%g %)\n", epsrel, result, abserr, 100.0*abserr/result);
    epsrel *= 2.0;
  } while(intstatus == GSL_EROUND);

  //printf("Workspace.intervals = %d\n", Workspace->size);
  //printf("Workspace.maximum_level = %d\n", Workspace->maximum_level);

  gsl_set_error_handler(old_handler);          /*  turn it on again */

  f_e = result / sqrt(8.0) / PI / PI;

  if(f_e<0)
    f_e= 0.0;

  return f_e;
}

double comp_DF_bulge_Eddington(double Q)
{
  double f_e;
  int i;
  double result, abserr;

  if(Q<=0)
    return 0.0;

  spline_params p;
  p.s = d2rhodpsi2_spline_bulge;
  p.a = d2rhodpsi2_spline_acc_bulge;
  p.Q = Q;

  gsl_function F;
  F.function = &eddington_integrand;
  F.params = &p;

  double epsabs=0;
  double epsrel=1e-4;
  int intstatus;

  gsl_error_handler_t * old_handler = gsl_set_error_handler_off(); /* turn off the error-handler */

  do
  {
    intstatus= gsl_integration_qags(&F, 0.0, Q, epsabs, epsrel, WORKSIZE, Workspace, &result, &abserr);
    //printf("epsrel = %g     result = %g +/- %g   (%g %)\n", epsrel, result, abserr, 100.0*abserr/result);
    //printf("error = %g\n", abserr);
    epsrel *= 2.0;
  } while(intstatus == GSL_EROUND);

  gsl_set_error_handler(old_handler);          /*  turn it on again */

  f_e = result / sqrt(8.0) / PI / PI;

  if(f_e<0)
    f_e= 0.0;

  return f_e;
}

/*
   Computes d2_rho / dpsi2, which is the primary term
   in the integral when using the Eddington formulation
   to calculate the distribution function.

   We do this for the halo and bulge components separately.

*/
void compute_d2rhodpsi2_halo()
{
  FILE *fd;
  char drhodpsifile[100]="";
  int i;

  double R, RplusdR, RminusdR, dR;
  double rho_RplusdR, rho_RminusdR;
  double psi_RplusdR, psi_RminusdR;
  double slope1, slope2, ra;

  // normalize just in case we get really small densities
  double rho_norm, psi_norm;

  for(i = 0; i <= DFSIZE; i++)
    {
	if(i < RSIZE)
	    /* continues same scaling as list_R
            R= list_R[RSIZE] * exp((RSIZE-i) * (log(LL) - log(Baselen)) / (RSIZE - 1));   */
	    /* more modest scaling to zero energy (ie larger radii) */
            R= list_R[RSIZE] * pow(10., 5.0 * (RSIZE-i) / RSIZE);
            //R= list_R[RSIZE] * pow(10., 7.0 * (RSIZE-i) / RSIZE);
	else
            R= list_R[DFSIZE-i];

        if(i == DFSIZE)
            R= list_R[1] * 0.5;

	/* interval of integration */
	dR= R * exp(1.0e-2 * (log(LL) - log(Baselen)) / (RSIZE - 1)) - R;
	if(dR>1.5)
	  dR= 1.5;
	if((dR/R)<1.0e-4)
	  dR=1.0e-4*R;
        RplusdR= R + dR;
        RminusdR= R - dR;

	list_radius[i]= R;

        /* compute psi (=-phi, i.e., negative the potential) along the z axis, although
           the full DF assumes spherical symmetry */
        psi_R[i]= -1.0 * comp_phi(0, R);;
        psi_RplusdR= -1.0 * comp_phi(0, RplusdR);
        psi_RminusdR= -1.0 * comp_phi(0, RminusdR);

        /* -------------------------------
	    DF for halo component  */
        rho_R_halo[i]= comp_rho_halo(R,0);
        rho_RplusdR= comp_rho_halo(RplusdR,0);
        rho_RminusdR= comp_rho_halo(RminusdR,0);

	if(AnisotropyRadius>0)
	  {
	    ra= AnisotropyRadius * RH;
            rho_R_halo[i] *= (1.0 + (R * R / ra / ra ));
            rho_RplusdR *= (1.0 + (RplusdR * RplusdR / ra / ra));
            rho_RminusdR *= (1.0 + (RminusdR * RminusdR / ra / ra));
	  }

	// normalize
	rho_norm= rho_R_halo[i];
	rho_R_halo[i] /= rho_norm;
	rho_RplusdR /= rho_norm;
	rho_RminusdR /= rho_norm;
	psi_norm= psi_R[i];
	psi_R[i] /= psi_norm;
	psi_RplusdR /= psi_norm;
	psi_RminusdR /= psi_norm;

        slope1= (rho_RplusdR - rho_R_halo[i]) / (psi_RplusdR - psi_R[i]);
        slope2= (rho_R_halo[i] - rho_RminusdR) / (psi_R[i] - psi_RminusdR);

        drhodpsi_halo[i]= 0.5*(slope1+slope2);

        d2rhodpsi2_halo[i]= (rho_RplusdR+rho_RminusdR-2.*rho_R_halo[i])/(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR);

	// now put units back in
	psi_R[i] *= psi_norm;
	rho_R_halo[i] *= rho_norm;
	drhodpsi_halo[i] *= rho_norm / psi_norm;
	d2rhodpsi2_halo[i] *= rho_norm / psi_norm / psi_norm;

//printf("i= %d, R= %g,  second derivative denominators: %g, %g    diff= %g\n",i,R,(psi_RplusdR-psi_R[i]),(psi_R[i]-psi_RminusdR),(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR));
//printf("i= %d, R= %g,  rho= %g   drhodpsi_halo= %g    d2rhodpsi2_halo= %g     second derivative denominators: %g, %g    diff= %g\n",i,R,rho_R_halo[i],drhodpsi_halo[i],d2rhodpsi2_halo[i],(psi_RplusdR-psi_R[i]),(psi_R[i]-psi_RminusdR),(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR));
//printf("i= %d, R= %g,    rho_0= %g   ani_factor= %g  rho_R_halo= %g\n",i,R,comp_rho_halo(R,0),(1.0 + (R * R / ra / ra )),rho_R_halo[i]);

    }

  xi[0] = 0;
  yi[0] = 0;
  /* generate spline for the d2rho / dpsi2 term */
  for(i = 0; i <= DFSIZE; i++)
    {
	xi[i + 1] = psi_R[i];
	yi[i + 1] = d2rhodpsi2_halo[i];
	if( (xi[i+1] - xi[i]) < 0)  xi[i+1] =  xi[i] + 1e-5;

	printf("xi[i] = %f  yi[i] = %f delta_x = %f\n",xi[i + 1], yi[i + 1], xi[i+1] - xi[i]);
    }

  printf("check_abcdefg\n");

  gsl_spline_init(d2rhodpsi2_spline_halo, xi, yi, DFSIZE+1);

  printf("check_abcdefg done\n");

  /* write this to a file so we have a record of it */
  strcpy(drhodpsifile,"");
  if(strstr(OutputFile,".hdf5"))
        strncpy(drhodpsifile, OutputFile, strlen(OutputFile)-5);
  if(strstr(OutputFile,".dat"))
        strncpy(drhodpsifile, OutputFile, strlen(OutputFile)-4);
  strcat(drhodpsifile, ".drhodpsi_halo");
  if((fd = fopen(drhodpsifile,"w")))
    {
        fprintf(fd,"# drhodpsi file, halo component\n");
        fprintf(fd,"# n= %d \n",DFSIZE);
        fprintf(fd,"# \n");
        fprintf(fd,"#   R (kpc)         psi               rho         drho/dpsi      d2rho/dpsi2   spline(psi)   \n");
        fprintf(fd,"# \n");
        for(i = 0; i <= DFSIZE; i++)
          fprintf(fd," %8.5e     %8.5e      %8.5e     %8.5e     %8.5e    %8.5e  \n",list_radius[i],psi_R[i],rho_R_halo[i],drhodpsi_halo[i],d2rhodpsi2_halo[i],gsl_spline_eval(d2rhodpsi2_spline_halo, 0.999*psi_R[i], d2rhodpsi2_spline_acc_halo));

        fclose(fd);
    }

}

void compute_d2rhodpsi2_bulge()
{
  FILE *fd;
  char drhodpsifile[100]="";
  int i;

  double R, RplusdR, RminusdR, dR;
  double rho_RplusdR, rho_RminusdR;
  double psi_RplusdR, psi_RminusdR;
  double slope1, slope2, ra;

  // normalize just in case we get really small densities
  double rho_norm, psi_norm;

  if(M_BULGE == 0)
    return;

  for(i = 0; i <= DFSIZE; i++)
    {
	if(i < RSIZE)
	    /* continues same scaling as list_R
            R= list_R[RSIZE] * exp((RSIZE-i) * (log(LL) - log(Baselen)) / (RSIZE - 1));   */
	    /* more modest scaling to zero energy (ie larger radii) */
            R= list_R[RSIZE] * pow(10., 5.0 * (RSIZE-i) / RSIZE);
            //R= list_R[RSIZE] * pow(10., 7.0 * (RSIZE-i) / RSIZE);
	else
            R= list_R[DFSIZE-i];

        if(i == DFSIZE)
            R= list_R[1] * 0.5;

	/* interval of integration */
	dR= R * exp(1.0e-2 * (log(LL) - log(Baselen)) / (RSIZE - 1)) - R;
	if(dR>1.5)
	  dR= 1.5;
	if((dR/R)<1.0e-4)
	  dR=1.0e-4*R;
        RplusdR= R + dR;
        RminusdR= R - dR;

	list_radius[i]= R;

        /* compute psi (=-phi, i.e., negative the potential) along the z axis, although
           the full DF assumes spherical symmetry */
        psi_R[i]= -1.0 * comp_phi(0, R);;
        psi_RplusdR= -1.0 * comp_phi(0, RplusdR);
        psi_RminusdR= -1.0 * comp_phi(0, RminusdR);

        /* -------------------------------
	    DF for bulge component  */
        rho_R_bulge[i]= comp_rho_bulge(R,0);
        rho_RplusdR= comp_rho_bulge(RplusdR,0);
        rho_RminusdR= comp_rho_bulge(RminusdR,0);

	if(AnisotropyRadius>0)
	  {
	    ra= AnisotropyRadius * RH;
            rho_R_bulge[i] *= (1.0 + (R * R / ra / ra ));
            rho_RplusdR *= (1.0 + (RplusdR * RplusdR / ra / ra));
            rho_RminusdR *= (1.0 + (RminusdR * RminusdR / ra / ra));
	  }

	// normalize
	rho_norm= rho_R_bulge[i];
	rho_R_bulge[i] /= rho_norm;
	rho_RplusdR /= rho_norm;
	rho_RminusdR /= rho_norm;
	psi_norm= psi_R[i];
	psi_R[i] /= psi_norm;
	psi_RplusdR /= psi_norm;
	psi_RminusdR /= psi_norm;

        slope1= (rho_RplusdR - rho_R_bulge[i]) / (psi_RplusdR - psi_R[i]);
        slope2= (rho_R_bulge[i] - rho_RminusdR) / (psi_R[i] - psi_RminusdR);

        drhodpsi_bulge[i]= 0.5*(slope1+slope2);

        d2rhodpsi2_bulge[i]= (rho_RplusdR+rho_RminusdR-2.*rho_R_bulge[i])/(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR);

	// now put units back in
	psi_R[i] *= psi_norm;
	rho_R_bulge[i] *= rho_norm;
	drhodpsi_bulge[i] *= rho_norm / psi_norm;
	d2rhodpsi2_bulge[i] *= rho_norm / psi_norm / psi_norm;

//printf("i= %d, R= %g,  second derivative denominators: %g, %g    diff= %g\n",i,R,(psi_RplusdR-psi_R[i]),(psi_R[i]-psi_RminusdR),(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR));
//printf("i= %d, R= %g,  rho= %g   drhodpsi_bulge= %g    d2rhodpsi2_bulge= %g     second derivative denominators: %g, %g    diff= %g\n",i,R,rho_R_bulge[i],drhodpsi_bulge[i],d2rhodpsi2_bulge[i],(psi_RplusdR-psi_R[i]),(psi_R[i]-psi_RminusdR),(psi_RplusdR-psi_R[i])/(psi_R[i]-psi_RminusdR));
//printf("i= %d, R= %g,    rho_0= %g   ani_factor= %g  rho_R_bulge= %g\n",i,R,comp_rho_bulge(R,0),(1.0 + (R * R / ra / ra )),rho_R_bulge[i]);

    }

  /* generate spline for the d2rho / dpsi2 term */
  for(i = 0; i <= DFSIZE; i++)
    {
	xi[i + 1] = psi_R[i];
	yi[i + 1] = d2rhodpsi2_bulge[i];
        if( (xi[i+1] - xi[i]) < 0)  xi[i+1] =  xi[i] + 1e-5;

    }

  gsl_spline_init(d2rhodpsi2_spline_bulge, xi, yi, DFSIZE+1);

  /* write this to a file so we have a record of it */
  strcpy(drhodpsifile,"");
  if(strstr(OutputFile,".hdf5"))
        strncpy(drhodpsifile, OutputFile, strlen(OutputFile)-5);
  if(strstr(OutputFile,".dat"))
        strncpy(drhodpsifile, OutputFile, strlen(OutputFile)-4);
  strcat(drhodpsifile, ".drhodpsi_bulge");
  if((fd = fopen(drhodpsifile,"w")))
    {
        fprintf(fd,"# drhodpsi file, bulge component\n");
        fprintf(fd,"# n= %d \n",DFSIZE);
        fprintf(fd,"# \n");
        fprintf(fd,"#   R (kpc)         psi               rho         drho/dpsi      d2rho/dpsi2   spline(psi)   \n");
        fprintf(fd,"# \n");
        for(i = 0; i <= DFSIZE; i++)
          fprintf(fd," %8.5e     %8.5e      %8.5e     %8.5e     %8.5e    %8.5e  \n",list_radius[i],psi_R[i],rho_R_bulge[i],drhodpsi_bulge[i],d2rhodpsi2_bulge[i],gsl_spline_eval(d2rhodpsi2_spline_bulge, 0.999*psi_R[i], d2rhodpsi2_spline_acc_bulge));

        fclose(fd);
    }

}

void comp_DF_init(void)
{
  /* stores spline on the Eddington integrand */
  xi= vector(1, DFSIZE+1);
  yi= vector(1, DFSIZE+1);
  eddint= vector(1, DFSIZE+1);

  /* arrays for energies, radii, d2rho/dpsi2 term in Eddington integrand, and pre-computed DF  */
  list_radius= vector(0, DFSIZE);
  list_E= vector(0, DFSIZE);
  psi_R= vector(0, DFSIZE);

  /* halo specific */
  rho_R_halo= vector(0, DFSIZE);
  drhodpsi_halo= vector(0, DFSIZE);
  d2rhodpsi2_halo= vector(0, DFSIZE);
  DistFunc_halo= vector(0, DFSIZE);

  /* bulge specific */
  rho_R_bulge= vector(0, DFSIZE);
  drhodpsi_bulge= vector(0, DFSIZE);
  d2rhodpsi2_bulge= vector(0, DFSIZE);
  DistFunc_bulge= vector(0, DFSIZE);

  /* setup the d2rho_dpsi2 array */

  /* for the halo */
  printf("allocating for halo\n");
  d2rhodpsi2_spline_halo = gsl_spline_alloc(gsl_interp_cspline, DFSIZE+1);
  d2rhodpsi2_spline_acc_halo = gsl_interp_accel_alloc();
  printf("evaluating for halo\n");
  compute_d2rhodpsi2_halo();

  /* for the bulge */
  printf("allocating for bulge\n");
  d2rhodpsi2_spline_bulge = gsl_spline_alloc(gsl_interp_cspline, DFSIZE+1);
  d2rhodpsi2_spline_acc_bulge = gsl_interp_accel_alloc();
  compute_d2rhodpsi2_bulge();

  printf("some workplace alloc\n");
  Workspace = gsl_integration_workspace_alloc(WORKSIZE);
}

void comp_DF_Eddington_close(void)
{
  free_vector(xi, 1, DFSIZE+1);
  free_vector(yi, 1, DFSIZE+1);
  free_vector(eddint, 1, DFSIZE+1);
  free_vector(psi_R, 1, DFSIZE);
  free_vector(list_radius, 1, DFSIZE);
  free_vector(list_E, 1, DFSIZE);
  free_vector(rho_R_halo, 1, DFSIZE);
  free_vector(rho_R_bulge, 1, DFSIZE);
  free_vector(drhodpsi_halo, 1, DFSIZE);
  free_vector(drhodpsi_bulge, 1, DFSIZE);
  free_vector(d2rhodpsi2_halo, 1, DFSIZE);
  free_vector(d2rhodpsi2_bulge, 1, DFSIZE);
  free_vector(DistFunc_halo, 1, DFSIZE);
  free_vector(DistFunc_bulge, 1, DFSIZE);

  gsl_integration_workspace_free(Workspace);

  gsl_spline_free(d2rhodpsi2_spline_halo);
  printf("check_1234\n");
  gsl_interp_accel_free(d2rhodpsi2_spline_acc_halo);

  gsl_spline_free(d2rhodpsi2_spline_bulge);
  printf("check_abcd\n");
  gsl_interp_accel_free(d2rhodpsi2_spline_acc_bulge);
}

/* ----------------------------------------------------------------- */

double compute_ani_beta(double r)
{
  double beta;

  /* isotropic model */
  beta = 0.0;

  if (AnisotropyRadius > 0)
    {
        double ra;

        ra= AnisotropyRadius * RH;

        /* Osipkov & Merritt anisotropy dependency */
        beta = r * r / (r * r + ra * ra);
    }

  return beta;
}


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

int
func (double t, const double y[], double f[],
      void *params)
{
  double mu = *(double *)params;
  f[0] = y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

int
main (void)
{
  const gsl_odeiv_step_type * T
    = gsl_odeiv_step_rk8pd;

  gsl_odeiv_step * s
    = gsl_odeiv_step_alloc (T, 2);
  gsl_odeiv_control * c
    = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e
    = gsl_odeiv_evolve_alloc (2);

  double mu = 10;
  gsl_odeiv_system sys = {func, jac, 2, &mu};

  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y[2] = { 1.0, 0.0 };

  while (t < t1)
    {
      int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys,
                                           &t, t1,
                                           &h, y);

      if (status != GSL_SUCCESS)
          break;

      printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
//
//  mcmc_lotkavolterra.c
//
//
//  Created by Juan Pablo Molano on 15/11/15.
//
//

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#define USAGE "./mcmc_lotkavolterra.x n_burn n_steps"

//------------------------------------------------------------------------------------------------------------

void load_matrix(char *filename, float *time, float *pre, float *dep , int n);
void print_data(float *array, float *array2 , float *array3 , float *array4 , float *array5 , float *array6, float *array7, float *array8 , float *array9 , float *array10,  int n_puntos);
float *reserva(int n);
float likelihood(float *y_obs, float *y_model, int n_row);
void ec_dep(float *array, float x, float y, float gamma, float delta);
void ec_pre(float *array, float x, float y, float alpha, float beta);
int min_likelihood( float *array , int n_row);

//------------------------------------------------------------------------------------------------------------

int main(int argc, char **argv){

    float *tim, *pre, *dep;
    float  m_x , m_y , k1_x , k1_y , k2_x , k2_y ,  k3_x , k3_y ,  k4_x , k4_y;
    float x1 , x2 , x3, x4, y1, y2, y3, y4 , h;
    float alpha_prime, beta_prime, gamma_prime, delta_prime, l_prime, l_init, gam, alph, bet , best_alpha , best_beta, best_gamma, best_delta;
    float *fit_pre, *fit_dep , *dep_init, *dep_prime , *pre_init , *pre_prime;
    int i, j , n_row = 96;
    float *alpha_a, *beta_a , *gamma_a , *delta_a , *walk_a, *best_dep, *best_pre;
    srand48(time(NULL));

    alpha_a = reserva( atof(argv[1]) + atof(argv[2]));
    beta_a = reserva( atof(argv[1]) + atof(argv[2]));
    gamma_a= reserva( atof(argv[1]) + atof(argv[2]));
    delta_a = reserva( atof(argv[1]) + atof(argv[2]));
    walk_a = reserva( atof(argv[1]) + atof(argv[2]));

    alpha_a[0] = drand48();
    beta_a[0] = drand48();
    gamma_a[0] = drand48();
    delta_a[0] = drand48();

    tim = reserva(96);
    pre = reserva(96);
    dep = reserva(96);
    fit_dep = reserva(96);
    fit_pre = reserva(96);
    dep_init = reserva(96);
    dep_prime = reserva(96);
    pre_init = reserva(96);
    pre_prime = reserva(96);
    best_dep = reserva(96);
    best_pre = reserva(96);

    load_matrix( "new_data.txt" , tim , pre , dep, 96);

    h = tim[1] - tim[0];

    fit_pre[0] = pre[0];
    fit_dep[0] = dep[0];

    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    for(i=1;i<96;i++){

        ec_pre( &k1_x , fit_pre[i-1] , fit_dep[i-1] , alpha_a[0] , beta_a[0]);
        ec_dep( &k1_y , fit_pre[i-1] , fit_dep[i-1] , gamma_a[0] , delta_a[0]);

        x1 = fit_pre[i-1] + (h/2.0) * k1_x;
        y1 = fit_dep[i-1] + (h/2.0) * k1_y;

        ec_pre( &k2_x , fit_pre[i-1] , fit_dep[i-1] , alpha_a[0] , beta_a[0]);
        ec_dep( &k2_y , fit_pre[i-1] , fit_dep[i-1] , gamma_a[0] , delta_a[0]);

        x2 = fit_pre[i-1] + (h/2.0) * k2_x;
        y2 = fit_dep[i-1] + (h/2.0) * k2_y;

        ec_pre( &k3_x , fit_pre[i-1] , fit_dep[i-1] , alpha_a[0] , beta_a[0]);
        ec_dep( &k3_y , fit_pre[i-1] , fit_dep[i-1] , gamma_a[0] , delta_a[0]);

        x3 = fit_pre[i-1] + h * k3_x;
        y3 = fit_dep[i-1] + h * k3_y;

        ec_pre( &k4_x , fit_pre[i-1] , fit_dep[i-1] , alpha_a[0] , beta_a[0]);
        ec_dep( &k4_y , fit_pre[i-1] , fit_dep[i-1] , gamma_a[0] , delta_a[0]);

        m_x = (1.0/6.0)*(k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
        m_y = (1.0/6.0)*(k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);

        fit_pre[i] = fit_pre[i-1] + h * m_x;
        fit_dep[i] = fit_dep[i-1] + h * m_y;
    }

    for (i=0; i < atof(argv[1]) + atof(argv[2]); i++) {

        alpha_prime = gsl_ran_gaussian(r, 0.1) + alpha_a[i];
        beta_prime = gsl_ran_gaussian(r, 0.1) + beta_a[i];
        gamma_prime = gsl_ran_gaussian(r, 0.1) + gamma_a[i];
        delta_prime = gsl_ran_gaussian(r, 0.1) + delta_a[i];

        for(j=1;j<96;j++){

            dep_init[0] = dep[0];
            pre_init[0] = pre[0];

            ec_pre( &k1_x , pre_init[j-1] , dep_init[j-1] , alpha_a[i] , beta_a[i]);
            ec_dep( &k1_y , pre_init[j-1] , dep_init[j-1] , gamma_a[i] , delta_a[i]);

            x1 = pre_init[j-1] + (h/2.0) * k1_x;
            y1 = dep_init[j-1] + (h/2.0) * k1_y;

            ec_pre( &k2_x , pre_init[j-1] , dep_init[j-1] , alpha_a[i] , beta_a[i]);
            ec_dep( &k2_y , pre_init[j-1] , dep_init[j-1] , gamma_a[i] , delta_a[i]);

            x2 = pre_init[j-1] + (h/2.0) * k2_x;
            y2 = dep_init[j-1] + (h/2.0) * k2_y;

            ec_pre( &k3_x , pre_init[j-1] , dep_init[j-1] , alpha_a[i] , beta_a[i]);
            ec_dep( &k3_y , pre_init[j-1] , dep_init[j-1] , gamma_a[i] , delta_a[i]);

            x3 = pre_init[j-1] + h * k3_x;
            y3 = dep_init[j-1] + h * k3_y;

            ec_pre( &k4_x , pre_init[j-1] , dep_init[j-1] , alpha_a[i] , beta_a[i]);
            ec_dep( &k4_y , pre_init[j-1] , dep_init[j-1] , gamma_a[i] , delta_a[i]);

            m_x = (1.0/6.0)*(k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
            m_y = (1.0/6.0)*(k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);

            pre_init[j] = pre_init[i-1] + h * m_x;
            dep_init[j] = dep_init[i-1] + h * m_y;
        }

        for(j=1;j<96;j++){

            dep_prime[0] = dep[0];
            pre_prime[0] = pre[0];

            ec_pre( &k1_x , pre_prime[j-1] , dep_prime[j-1] , alpha_prime , beta_prime);
            ec_dep( &k1_y , pre_prime[j-1] , dep_prime[j-1] , gamma_prime , delta_prime);

            x1 = pre_prime[j-1] + (h/2.0) * k1_x;
            y1 = dep_prime[j-1] + (h/2.0) * k1_y;

            ec_pre( &k2_x , pre_prime[j-1] , dep_prime[j-1] ,  alpha_prime , beta_prime);
            ec_dep( &k2_y , pre_prime[j-1] , dep_prime[j-1] , gamma_prime , delta_prime);

            x2 = pre_prime[j-1] + (h/2.0) * k2_x;
            y2 = dep_prime[j-1] + (h/2.0) * k2_y;

            ec_pre( &k3_x , pre_prime[j-1] , dep_prime[j-1] ,  alpha_prime , beta_prime);
            ec_dep( &k3_y , pre_prime[j-1] , dep_prime[j-1] , gamma_prime , delta_prime);

            x3 = pre_prime[j-1] + h * k3_x;
            y3 = dep_prime[j-1] + h * k3_y;

            ec_pre( &k4_x , pre_prime[j-1] , dep_prime[j-1] ,  alpha_prime , beta_prime);
            ec_dep( &k4_y , pre_prime[j-1] , dep_prime[j-1] , gamma_prime , delta_prime);

            m_x = (1.0/6.0)*(k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
            m_y = (1.0/6.0)*(k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);

            pre_prime[j] = pre_prime[j-1] + h * m_x;
            dep_prime[j] = dep_prime[j-1] + h * m_y;
        }

        l_prime = likelihood( dep, dep_prime , n_row) + likelihood( pre , pre_prime , n_row);
        l_init = likelihood( dep , dep_init ,  n_row) + likelihood( pre , pre_init , n_row);

        gam = l_init - l_prime;

        if (gam>=0.0) {
            alpha_a[i+1] = alpha_prime;
            beta_a[i+1] = beta_prime;
            gamma_a[i+1] = gamma_prime;
            delta_a[i+1] = delta_prime;
            walk_a[i+1] = l_prime;
        }
        else{

            bet = drand48();
            alph = exp(gam);

            if (bet <= alph) {

                alpha_a[i+1] = alpha_prime;
                beta_a[i+1] = beta_prime;
                gamma_a[i+1] = gamma_prime;
                delta_a[i+1] = delta_prime;
                walk_a[i+1] = l_prime;

            }
            else{

                alpha_a[i+1] = alpha_a[i];
                beta_a[i+1] = beta_a[i];
                gamma_a[i+1] = gamma_a[i];
                delta_a[i+1] = delta_a[i];
                walk_a[i+1] = l_prime;

            }

        }
    }

    j = min_likelihood( walk_a , n_row);
    best_alpha = alpha_a[j];
    best_beta = beta_a[j];
    best_gamma = gamma_a[j];
    best_delta = delta_a[j];

    for(i=1;i<96;i++){

        best_dep[0] = dep[0];
        best_pre[0] = pre[0];

        ec_pre( &k1_x , best_pre[i-1] , best_dep[i-1] , best_alpha , best_beta);
        ec_dep( &k1_y ,  best_pre[i-1] , best_dep[i-1], best_gamma , best_delta);

        x1 = fit_pre[i-1] + (h/2.0) * k1_x;
        y1 = fit_dep[i-1] + (h/2.0) * k1_y;

        ec_pre( &k2_x ,  best_pre[i-1] , best_dep[i-1] ,  best_alpha , best_beta);
        ec_dep( &k2_y ,  best_pre[i-1] , best_dep[i-1] , best_gamma , best_delta);

        x2 = fit_pre[i-1] + (h/2.0) * k2_x;
        y2 = fit_dep[i-1] + (h/2.0) * k2_y;

        ec_pre( &k3_x ,  best_pre[i-1] , best_dep[i-1],  best_alpha , best_beta);
        ec_dep( &k3_y ,  best_pre[i-1] , best_dep[i-1], best_gamma , best_delta);

        x3 = fit_pre[i-1] + h * k3_x;
        y3 = fit_dep[i-1] + h * k3_y;

        ec_pre( &k4_x ,  best_pre[i-1] , best_dep[i-1],  best_alpha , best_beta);
        ec_dep( &k4_y ,  best_pre[i-1] , best_dep[i-1] , best_gamma , best_delta);

        m_x = (1.0/6.0)*(k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
        m_y = (1.0/6.0)*(k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);

        best_pre[i] = best_pre[i-1] + h * m_x;
        best_dep[i-1] = best_dep[i-1] + h * m_y;
    }

    print_data( tim , pre, dep, best_pre , best_pre, walk_a , alpha_a , beta_a , gamma_a , delta_a , n_row);

    gsl_rng_free (r);

    return 0;

}

//------------------------------------------------------------------------------------------------------------

void load_matrix(char *filename, float *time, float *pre, float *dep , int n){

    FILE *in;
    int i;
    in=fopen(filename,"r");
    if(!in){
        printf("problems opening the file %s\n", filename);
        exit(1);
    }
    for(i=0;i<(n);i++){

        fscanf(in, "%f %f %f\n", &time[i], &pre[i], &dep[i]);

    }
    fclose(in);
}

//------------------------------------------------------------------------------------------------------------

void print_data(float *array, float *array2 , float *array3 , float *array4 , float *array5 , float *array6, float *array7, float *array8 , float *array9 ,  float *array10 ,   int n_puntos){
    int i;
    for(i=0;i<n_puntos;i++){
        printf("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n", array[i] , array2[i] , array3[i] , array4[i] , array5[i] , array6[i], array7[i] , array8[i] , array9[i] , array10[i]);
    }
}

//------------------------------------------------------------------------------------------------------------

float *reserva(int n){
    float *array;
    int i;
    if(!(array = malloc(n * sizeof(float)))){
        printf("Problema en reserva\n");
        exit(1);
    }
    for(i=0;i<n ;i++){
        array[i] = 0.0;
    }
    return array;
}

//------------------------------------------------------------------------------------------------------------

int min_likelihood( float *array , int n_row){

    int i;
    int ind = 0;

    for(i = 1; i < n_row; i ++){
        if ( array[i] < ind) {
            ind = i;
        }

    }

    return ind;
}

//------------------------------------------------------------------------------------------------------------

void ec_pre(float *array, float x, float y, float alpha, float beta){

    *array = x*(alpha - beta*y);

}

//------------------------------------------------------------------------------------------------------------

void ec_dep(float *array, float x, float y, float gamma, float delta){

    *array = -y*(gamma - delta*x);

}

//------------------------------------------------------------------------------------------------------------

float likelihood(float *y_obs, float *y_model, int n_row){

    int i;
    float c=0;

    for (i=0; i<n_row; i++) {

        c = c + ( y_obs[i] - y_model[i]);
    }
    c = (1.0/2.0)*pow(c,2);

    return c;

}