#include <stdio.h>
#include <math.h>

#define		FAILURE			0
#define		SUCCESS			1

#define		GAMMA_PLUS		1
#define		GAMMA_MINUS		2
#define		LF_SHOCK		3
#define		RF_SHOCK		4
#define		CONTACT			5

struct _STATE
{
	double d,p,u;
};
typedef struct _STATE STATE;

struct _CENTERED_WAVE
{
	int wave_type;		/* simple wave or shock */
	double	speed_leading_edge,speed_trailing_edge;	/* for simple wave */
	double	speed_shock;				/* for shock */
	double	speed_contact;				/* for contact line */
};
typedef struct _CENTERED_WAVE CENTERED_WAVE;

struct _RIEMANN_SOLN
{
	STATE left_state;
	STATE right_state;
	STATE left_center_state;
	STATE right_center_state;
	CENTERED_WAVE left_wave;
	CENTERED_WAVE right_wave;
	CENTERED_WAVE contact;
};
typedef struct _RIEMANN_SOLN RIEMANN_SOLN;

struct _RIEMANN_INPUT
{
	STATE left_state;
	STATE right_state;
	double gamma;
};
typedef struct _RIEMANN_INPUT RIEMANN_INPUT;

static int BisectionFindRoot(double(*f)(double,void*),void*,double*,
			double,double,double,double);
static int SecantFindRoot(double(*f)(double,void*),void*,double*,
			double,double,double,double);
static int RiemannSolution(RIEMANN_INPUT,RIEMANN_SOLN*);
static double UlMinusUr(double,void*);
static double u_left_centerd_wave(STATE,double);
static double u_right_centerd_wave(STATE,double);

static int RiemannSolnAtXi(
	RIEMANN_SOLN riem_soln,
	STATE *soln_at_xi)
{
}	/* end RiemannSolnAtXi */

static int RiemannSolution(
	RIEMANN_INPUT input,
	RIEMANN_SOLN *riem_soln)
{
	void *params;
	double (*func)(double,void*);
	double epsilon,delta;
	double p_mid,p0,p1;

	params = (void*)&input;
	func = UlMinusUr;
	SecantFindRoot(func,params,&p_mid,p0,p1,epsilon,delta);
	/*
	BisectionFindRoot(func,params,&p_mid,p0,p1,epsilon,delta);
	*/
}	/* end RiemannSolution */

/*	bisection method to solve equation f(x) = 0        */

static int BisectionFindRoot(        	/* return SUCCESS and FAILURE */ 
	double (*f)(double,void*),	/* function f(x) */
	void *f_params,         	/* function parameters, unknown type */
	double *px,         	/* pointer to solution */
	double x0,		/* input of initial guess 0 */
	double x1,         	/* input of initial guess 1 */
	double epsilon,		/* tolerance for x */
	double delta) 		/* tolerance for f(x) */
{
}	/* end bisection_find_root */

/*	secant method to solve equation f(x) = 0	*/

static int SecantFindRoot(        	/* return SUCCESS and FAILURE */ 
	double (*f)(double,void*),	/* function f(x) */
	void *f_params,         	/* function parameters, unknown type */
	double *px,         	/* pointer to solution */
	double x0,		/* input of initial guess 1 */
	double x1,         	/* input of initial guess 0 */
	double epsilon,		/* tolerance for x */
	double delta) 		/* tolerance for f(x) */
{
}	/* end SecantFindRoot */

static double UlMinusUr(
	double p,
	void *params)
{
	double ul,ur;

	RIEMANN_INPUT *riem_input = (RIEMANN_INPUT*)params;
	ul = u_left_centerd_wave(riem_input->left_state,p);
	ur = u_right_centerd_wave(riem_input->right_state,p);
	return ul - ur;
}	/* end UlMinusUr */

static double u_left_centerd_wave(
	STATE lstate,
	double p_star)
{
	double u_across_wave;
	if (p_star < lstate.p) /* for rarefaction */
	{
	}
	else /* (p_star >= lstate.p), for shock */
	{
	}
	return ul_star;
}	/* end u_left_centerd_wave */

static double u_right_centerd_wave(
	STATE state,
	double p_star)
{
	double u_across_wave;
	if (p_star < state.p) /* for rarefaction */
	{
	}
	else /* (p_star >= state.p), for shock */
	{
	}
	return ur_star;
}	/* end u_right_centerd_wave */