/* Steven Andrews, 8/4/1992; modified 8/7/96, 9/1/98, 10/18/01 */
/* See documentation called math2 doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <math.h>
#include <stdlib.h>
#include <float.h>
#include "math2.h"


int sign(float x)	{
	return((x==0)?0:(x>0)?1:-1);	}


int signD(double x)	{
	return((x==0)?0:(x>0)?1:-1);	}


float sinc(float x)	{
	return (x==0)?1:sin(x)/x;	}


float box(float x)	{
	return (fabs(x)>1)?0:1;	}


float bessj0(float x) { /* Copied from Numerical Recipies */
	double ax,z;
	double xx,y,ans,ans1,ans2;
	
	if((ax=fabs(x))<8.0) {
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+y*(9494680.718+y*(59272.64853+y*(267.8532712+y*1.0))));
	ans=ans1/ans2; }
	else {
		z=8.0/ax;
		y=z*z;
		xx=ax-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2=-0.1562499995e-1+y*(0.1430488765e-3+y*(-0.6911147651e-5+y*(0.7621095161e-6-y*0.934935152e-7)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); }
	return ans; }


float bessj1(float x)	{	/* Copied exactly from Numerical Recipies */
	double ax,z;
	double xx,y,ans,ans1,ans2;

	if((ax=fabs(x))<8.0) {
		y=x*x;
		ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
		ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
		ans=ans1/ans2;	}
	else	{
		z=8.0/ax;
		y=z*z;
		xx=ax-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
		if(x<0.0) ans=-ans;	}
	return ans;	}


float gauss(float x,float mean,float sd)	{
	return 1.0/(sd*SQRT2PI)*exp(-1.0/2*((x-mean)/sd*(x-mean)/sd));	}


float gammaln(float x)	{	/* Returns natural log of gamma function, partly from Numerical Recipies and part from Math CRC */
	double sum,t;
	int j;
	static double c[6]={76.18009173,-86.50532033,24.01409822,-1.231739516,0.120858003e-2,-0.536382e-5};
	
	if(x==floor(x)&&x<=0) sum=DBL_MAX;					// 0 or negative integer
	else if(x==floor(x))	{											// positive integer
		sum=0;
		for(t=2;t<x-0.1;t+=1.0)	sum+=log(t);	}
	else if(x==0.5)	sum=0.572364942;						// 1/2
	else if(2*x==floor(2*x)&&x>0)	{							// other positive half integer
		sum=0.572364942;
		for(t=0.5;t<x-0.1;t+=1)	sum+=log(t);	}
	else if(2*x==floor(2*x))	{									// negative half integer
		sum=0.572364942;
		for(t=0.5;t<-x+0.1;t+=1)	sum-=log(t);	}
	else if(x<0)																// other negative
		sum=gammaln(x+1)-log(-x);
	else	{																			// other positive
		x-=1.0;
		t=x+5.5;
		t-=(x+0.5)*log(t);
		sum=1.0;
		for(j=0;j<=5;j++)	{
			x+=1.0;
			sum+=c[j]/x;	}
		sum=-t+log(2.50662827465*sum);	}
	return(sum);	}


float gammp(float a,float x)	{ // Returns incomplete gamma function, partly from Numerical Recipes
	double sum,del,ap,eps;
	double gold=0,g=1,fac=1,b1=1,b0=0,anf,ana,an=0,a1,a0=1;

	eps=3e-7;
	if(x<0||a<=0) return -1;			// out of bounds
	else if(x==0) return 0;
	else if(x<a+1)	{
		ap=a;
		del=sum=1.0/a;
		while(fabs(del)>fabs(sum)*eps&&ap-a<100)	{
			ap+=1;
			del*=x/ap;
			sum+=del;	}
		return sum*exp(-x+a*log(x)-gammaln(a));	}
	else {
		a1=x;
		for(an=1;an<100;an++)	{
			ana=an-a;
			a0=(a1+a0*ana)*fac;
			b0=(b1+b0*ana)*fac;
			anf=an*fac;
			a1=x*a0+anf*a1;
			b1=x*b0+anf*b1;
			if(a1)	{
				fac=1.0/a1;
				g=b1*fac;
				if(fabs((g-gold)/g)<eps)
					return 1.0-exp(-x+a*log(x)-gammaln(a))*g;
				gold=g;	}}}
	return -1;	}							// shouldn't ever get here


float erfn(float x)	{				// Numerical Recipies
	return (x<0?-gammp(0.5,x*x):gammp(0.5,x*x));	}


float erfnc(float x)	{			// Numerical Recipies
	return 1.0-(x<0?-gammp(0.5,x*x):gammp(0.5,x*x));	}


float betaln(float x1,float x2)	{				// Math CRC
	return(gammaln(x1)+gammaln(x2)-gammaln(x1+x2));	}


int iseven(int x)	{
	return(1-abs(x)%2);	}


int is2ton(int x)	{
	return (x&-x)==x; }


int isinteger(float x) {
	return(floor(x)==x); }


int next2ton(int x)	{
	int i;
	
	if(x<0) return 0;
	if(!x) return 1;
	for(i=0;x!=1;i++)	x=x>>1;
	return x<<(i+1);	}


float factorial(int n)	{
	double y;
	
	y=1;
	for(;n>1;n--)	y*=n;
	return y;	}


float choose(int n,int m)	{
	double y;
	
	if(m>n/2) m=n-m;
	y=1;
	while(m>0)	y*=n--/m--;
	return y;	}


int minus1to(int x)	{
	return((abs(x)%2)?(-1):1);	}


int intpower(int n,int p)	{
	int y=1;

	if(p<0) return 0;
	for(;p;p--) y*=n;
	return y;	}


int gcomdiv(int m,int n) {
	int temp;
	
	if(m==0||n==0) return 1;
	if(m<0) m=-m;
	if(n<0) n=-n;
	while(m>0) {
		if(m<n) {temp=m;m=n;n=temp;}
		m%=n; }
	return n; }


float hermite(float x,int n)	{
	if(n==0)	return 1;
	else if(n==1)	return 2*x;
	else if(n>1)	return(2*x*hermite(x,n-1)-2*(n-1)*hermite(x,n-2));
	else return 0;	}


float constrain(float x,float lo,float hi) {
	if(x<lo) return lo;
	if(x>hi) return hi;
	return x; }


float reflect(float x,float lo,float hi) {
	if(x>=lo&&x<=hi) return x;
	if(x<lo) return reflect(2*lo-x,lo,hi);
	return reflect(2*hi-x,lo,hi); }


float inversefn(float (*fn)(float),float y,float x1,float x2,int n) {
	float dx,y2;

	if((*fn)(x1)<(*fn)(x2))	dx=x2-x1;
	else {
		dx=x1-x2;
		x1=x2; }
	for(;n>0;n--) {
		dx*=0.5;
		x2=x1+dx;
		y2=(*fn)(x2);
		if(y2<y) x1=x2; }
	return x1+0.5*dx; }


double quadpts2paramsD(double *x,double *y,double *abc) {
	double det;

	det=x[0]*x[0]*x[1]-x[0]*x[0]*x[2]-x[0]*x[1]*x[1]+x[0]*x[2]*x[2]+x[1]*x[1]*x[2]-x[1]*x[2]*x[2];
	if(det>0||det<0) {
		abc[0]=((x[1]-x[2])*y[0]+(-x[0]+x[2])*y[1]+(x[0]-x[1])*y[2])/det;
		abc[1]=((-x[1]*x[1]+x[2]*x[2])*y[0]+(x[0]*x[0]-x[2]*x[2])*y[1]+(-x[0]*x[0]+x[1]*x[1])*y[2])/det;
		abc[2]=((x[1]*x[1]*x[2]-x[1]*x[2]*x[2])*y[0]+(-x[0]*x[0]*x[2]+x[0]*x[2]*x[2])*y[1]+(x[0]*x[0]*x[1]-x[0]*x[1]*x[1])*y[2])/det; }
	return det; }


double fouriersumD(double *a,double *b,int n,double l,double x) {
	double sum;
	int j;

	sum=a[0]/2;
	for(j=1;j<n;j++)
		sum+=a[j]*cos(j*PI*x/l)+b[j]*sin(j*PI*x/l);
	return sum; }


void SetHillParam(float *hp,float a,float e,float n) {
	hp[0]=a;
	hp[1]=e;
	hp[2]=n;
	return; }


/* Hill function.  hp is parameters in 3-element vector: A,E,n.  x is
independent variable.  Function: H(x)=A*x^n/(E^n+x^n) */
float HillFn(float *hp,float x) {
	return hp[0]*pow(x,hp[2])/(pow(hp[1],hp[2])+pow(x,hp[2])); }


/* Composition of Hill functions 1 and 2: H_tot(x)=H2(H1(x)).  The composition
is not exactly a Hill function but is sufficiently close that H_tot is a "close"
Hill function.  Parameters for H1, H2, and H_tot are the same as for HillFn. */
void HillFnCompose(float *hp1,float *hp2,float *hp12) {
	float a1,a2,a12,e1,e2,e12,n1,n2,n12,e2p,a12p,e12p;

	a1=hp1[0];
	a2=hp2[0];
	e1=hp1[1];
	e2=hp2[1];
	n1=hp1[2];
	n2=hp2[2];
	e2p=e2/a1;
	a12=a2/(1.0+pow(e2p,n2));
	e12=e1/pow(1/e2p*pow(1.0+2*pow(e2p,n2),1.0/n2)-1.0,1.0/n1);
	a12p=a12/a2;
	e12p=e12/e1;
	n12=4.0*n1*n2*pow(e2p,n2)/(a12p*pow(e12p,n1))/pow(1.0+pow(e12p,-n1),n2+1)/pow(pow(1.0+pow(e12p,-n1),-n2)+pow(e2p,n2),2);
	hp12[0]=a12;
	hp12[1]=e12;
	hp12[2]=n12;
	return; }


/* Composition of Hill functions 1 and 2, including feedback from 2 to 1.  The
composition is not exactly a Hill function but is a "close" fit to one.  The
reaction mechanism is that the activated product of function 2 is a reactant for
the inactivation of reaction 1.  Either, both, or neither hpf1 and hpf2 are
returned; for one to be returned send in a non-NULL address.  It is assumed, but
not checked, that both input n values (hp1[2], hp2[2]) are equal to 1. */
void HillFnComposeNF1(float *hp1,float *hp2,float *hpf1,float *hpf12) {
	float a1,e1,a2,e2;

	a1=hp1[0];
	e1=hp1[1];
	a2=hp2[0];
	e2=hp2[1];
	if(hpf1) {
		hpf1[0]=a1;
		hpf1[1]=e1*a1/(2*e2+a1);
		hpf1[2]=(2*e2+a1)/(3*e2+a1); }
	if(hpf12) {
		hpf12[0]=a1*a2/(a1+e2);
		hpf12[1]=e1*e2*a1/(2.0*(e2+a1)*(e2+a1));
		hpf12[2]=2.0/3.0; }
	return; }