/* Steven Andrews, 3/02 */
/* List of basis functions to be used with BasisFn.c */
/* See documentation called BasisFns doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <math.h>
#include <string.h>
#include "BasisFns.h"
#include "math2.h"
#include "Rn.h"
#include "RnSort.h"
#include "Set.h"
#include "Spectra.h"
#include "string2.h"

#define MAXPARAM 20

basisptr DeclareBasis(set bset,int n,char *name,char *desc,float (*addr)(float,float *,sptr,float *),float *param,char **pname);

basisptr DeclareBasis(set bset,int n,char *name,char *desc,float (*addr)(float,float *,sptr,float *),float *param,char **pname) {
	basisptr b;
	int i,ok;
	
	if(n<0||n>=MAXPARAM) return NULL;
	b=BasisAlloc(n);
	if(!b) return NULL;
	strncpy(b->name,name,STRCHAR);
	strncpy(b->proc,name,STRCHAR);
	strncpy(b->desc,desc,STRCHAR);
	b->addr=addr;
	for(i=0;i<n;i++) {
		b->param[i]=b->pold[i]=param[i];
		strncpy(b->pname[i],pname[i],STRCHAR); }
	ok=SetInsert((void *)b->name,(void *)b,bset);
	if(ok!=1) {BasisFree(b);return NULL;}
	return b; }


int getbasis(set bset) {
	basisptr b;
	int n;
	char *name;												/* name and procedure							*/
	char *desc;												/* short text description					*/
	float (*addr)(float,float *,sptr,float *);	/* pointer to function	*/
	float param[MAXPARAM];						/* initial parameter values				*/
	char *pname[MAXPARAM];						/* parameter names 								*/
	
	if(!bset) return 6;

	/* constant */
	name="constant";
	desc="y=offset";
	addr=&constbasis;
	n=1;
	pname[0]="offset";
	param[0]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* spectrum */
	name="spectrum";
	desc="y=value of spectrum";
	addr=&spectbasis;
	n=1;
	pname[0]="weight";
	param[0]=0.01;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->isspec=1;

	/* line */
	name="line";
	desc="y=slope*(x-x0)";
	addr=&linebasis;
	n=2;
	pname[0]="slope";
	param[0]=0.001;
	pname[1]="x0";
	param[1]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* exp */
	name="exp";
	desc="y=factor*exp(slope*x)";
	addr=&expbasis;
	n=2;
	pname[0]="factor";
	param[0]=1;
	pname[1]="slope";
	param[1]=0.1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* log */
	name="log";
	desc="y=weight*ln(slope*x+intercept)";
	addr=&logbasis;
	n=3;
	pname[0]="weight";
	param[0]=0.0001;
	pname[1]="slope";
	param[1]=1;
	pname[2]="intercept";
	param[2]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* quad */
	name="quad";
	desc="y=curve*(x-x0)^2+slope*(x-x0)+intercept";
	addr=&quadbasis;
	n=4;
	pname[0]="curve";
	param[0]=0.000001;
	pname[1]="slope";
	param[1]=0.00001;
	pname[2]="intercept";
	param[2]=0.01;
	pname[3]="x0";
	param[3]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	b->freeze[3]=1;
	if(!b) return 1;

	/* asinh */
	name="asinh";
	desc="y=weight*asinh(slope*x+intercept)";
	addr=&asinhbasis;
	n=3;
	pname[0]="weight";
	param[0]=0.0001;
	pname[1]="slope";
	param[1]=1;
	pname[2]="intercept";
	param[2]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* gaussian */
	name="gaussian";
	desc="y=area/(std_dev*Ã2¹)*exp(-(x-mean)^2/(2*std_dev^2))";
	addr=&gaussbasis;
	n=3;
	pname[0]="area";
	param[0]=1;
	pname[1]="mean";
	param[1]=0;
	pname[2]="std_dev";
	param[2]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;

	/* xgauss */
	name="xgauss";
	desc="y=x*gaussian(x)";
	addr=&xgaussbasis;
	n=3;
	pname[0]="area";
	param[0]=0.0002;
	pname[1]="mean";
	param[1]=1945;
	pname[2]="std_dev";
	param[2]=4;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;
	
	/* sine */
	name="sine";
	desc="y=amplitude*sin(frequency*x+shift)";
	addr=&sinbasis;
	n=3;
	pname[0]="amp";
	param[0]=1;
	pname[1]="freq";
	param[1]=1;
	pname[2]="shift";
	param[2]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* lorentz */
	name="lorentz";
	desc="y=max/(1+((x-mean)/(fwhm/2))^2)";
	addr=&lorentzbasis;
	n=3;
	pname[0]="max";
	param[0]=1;
	pname[1]="mean";
	param[1]=0;
	pname[2]="fwhm";
	param[2]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;

	/* peak */
	name="peak";
	desc="y=(1-shape)*xgaussian(x)+shape*xlorentzian(x)";
	addr=&peakbasis;
	n=4;
	pname[0]="height";
	param[0]=1;
	pname[1]="position";
	param[1]=2250;
	pname[2]="fwhm";
	param[2]=10;
	pname[3]="shape";
	param[3]=0.5;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;
	b->paramlo[3]=0;
	b->paramhi[3]=1;

	/* peakd1 */
	name="peakd1";
	desc="y=x weighted derivative of a peak";
	addr=&peak1basis;
	n=4;
	pname[0]="height";
	param[0]=1;
	pname[1]="position";
	param[1]=2250;
	pname[2]="fwhm";
	param[2]=10;
	pname[3]="shape";
	param[3]=0.5;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;
	b->paramlo[3]=0;
	b->paramhi[3]=1;

	/* peakd2 */
	name="peakd2";
	desc="y=x weighted second derivative of a peak";
	addr=&peak2basis;
	n=4;
	pname[0]="height";
	param[0]=1;
	pname[1]="position";
	param[1]=2250;
	pname[2]="fwhm";
	param[2]=10;
	pname[3]="shape";
	param[3]=0.5;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[2]=0;
	b->paramlo[3]=0;
	b->paramhi[3]=1;

	/* peakz */
	name="peakz";
	desc="y=z0*peak(x)+z1*peakd1(x)+z2*peakd2(x)";
	addr=&peakzbasis;
	n=7;
	pname[0]="z0";
	param[0]=0.0001;
	pname[1]="z1";
	param[1]=0.001;
	pname[2]="z2";
	param[2]=0.01;
	pname[3]="height";
	param[3]=1;
	pname[4]="position";
	param[4]=2250;
	pname[5]="fwhm";
	param[5]=10;
	pname[6]="shape";
	param[6]=0.5;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[0]=0;
	b->paramlo[2]=0;
	b->freeze[3]=b->freeze[4]=b->freeze[5]=b->freeze[6]=1;

	/* diffuse */
	name="diffuse";
	desc="y=spect(x)*area*exp(-dt*x^2)";
	addr=&diffusebasis;
	n=2;
	pname[0]="area";
	param[0]=1;
	pname[1]="dt";
	param[1]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->isspec=1;

	/* diffuse2 */
	name="diffuse2";
	desc="y=c0*erf(Ã(td/(x-t0))/4)^2";
	addr=&diffuse2basis;
	n=3;
	pname[0]="c0";
	param[0]=100;
	pname[1]="td";
	param[1]=1;
	pname[2],"t0";
	param[2]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* decay */
	name="decay";
	desc="y=convolution of gaussian with truncated exponential";
	addr=&convexpbasis;
	n=4;
	pname[0]="height";
	param[0]=1;
	pname[1]="fwhm";
	param[1]=0.16;
	pname[2]="tau";
	param[2]=1;
	pname[3]="shift";
	param[3]=10;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;
	b->paramlo[1]=0;
	b->paramlo[2]=0;
	b->freeze[1]=1;

	/* rational */
	name="rational";
	desc="y=(n0+n1*x+n2*x^2+n3*x^3+n4*x^4+n5*x^5)/(d0+d1*x+d2*x^2+d3*x^3+d4*x^4+d5*x^5)";
	addr=&rationbasis;
	n=12;
	pname[0]="n0";
	pname[1]="n1";
	pname[2]="n2";
	pname[3]="n3";
	pname[4]="n4";
	pname[5]="n5";
	pname[6]="d0";
	pname[7]="d1";
	pname[8]="d2";
	pname[9]="d3";
	pname[10]="d4";
	pname[11]="d5";
	param[0]=param[6]=1;
	param[1]=param[2]=param[3]=param[4]=param[5]=0;
	param[7]=param[8]=param[9]=param[10]=param[11]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	b->freeze[6]=b->freeze[8]=b->freeze[9]=b->freeze[10]=b->freeze[11]=1;
	if(!b) return 1;

	/* sigmoid */
	name="sigmoid";
	desc="y=min+(max-min)/(1+10^(slope*(ec50-x)))";
	addr=&sigmoidbasis;
	n=4;
	pname[0]="min";
	param[0]=0;
	pname[1]="max";
	param[1]=1;
	pname[2]="ec50";
	param[2]=-1;
	pname[3]="slope";
	param[3]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	/* hill */
	name="hill";
	desc="y=amp*x^n/(x^n+ec50^n)";
	addr=&hillbasis;
	n=3;
	pname[0]="amp";
	param[0]=1;
	pname[1]="ec50";
	param[1]=1;
	pname[2]="n";
	param[2]=1;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	b->paramlo[1]=0;
	b->paramlo[2]=0;
	if(!b) return 1;

	/* lines */
	name="lines";
	desc="y= slope a0 to x1,y1, then to x2,y2, then to x3,y3, then slope a3";
	addr=&linesbasis;
	n=8;
	pname[0]="a0";
	param[0]=1;
	pname[1]="x1";
	param[1]=1;
	pname[2]="y1";
	param[2]=1;
	pname[3]="x2";
	param[3]=2;
	pname[4]="y2";
	param[4]=1;
	pname[5]="x3";
	param[5]=3;
	pname[6]="y3";
	param[6]=2;
	pname[7]="a3";
	param[7]=0;
	b=DeclareBasis(bset,n,name,desc,addr,param,pname);
	if(!b) return 1;

	return 0; }


float constbasis(float x,float *param,sptr spec,float *deriv) {
	/* constant(x)=offset */
	/* param[0] = offset	*/
	if(deriv) deriv[0]=1;
 	return param[0]; }

float spectbasis(float x,float *param,sptr spec,float *deriv) {
	/* spectrum(x)=value of spectrum */
	/* param[0] = weight */
	int j=-2;
	float y;

	y=interpolate1(spec->x,spec->y,spec->n,&j,x);
	if(deriv) deriv[0]=y;
	return(param[0]*y); }

float linebasis(float x,float *param,sptr spec,float *deriv) {
	/* line(x)=slope*(x-x0) */
	/* param[0] = slope */
	/* param[1] = x0 */
	if(deriv) {
		deriv[0]=x-param[1];
		deriv[1]=-param[0]; }
	return param[0]*(x-param[1]); }

float expbasis(float x,float *param,sptr spec,float *deriv) {
	/* exp(x)=factor*exp(slope*x) */
	/* param[0] = factor */
	/* param[1] = slope */
	if(deriv) {
		deriv[0]=exp(param[1]*x);
		deriv[1]=param[0]*x*exp(param[1]*x); }
	return param[0]*exp(param[1]*x); }

float logbasis(float x,float *param,sptr spec,float *deriv) {
	/* log(x)=weight*ln(slope*x+intercept) */
	/* param[0] = weight */
	/* param[1] = slope */
	/* param[2] = intercept */
	float f;
	
	f=param[1]*x+param[2];
	if(f<=0)	{
		if(deriv) deriv[0]=deriv[1]=deriv[2]=0;
		return 0; }
	if(deriv)	{
		deriv[0]=log(f);
		deriv[1]=param[0]*x/f;
		deriv[2]=param[0]/f; }
	return param[0]*log(f); }

float quadbasis(float x,float *param,sptr spec,float *deriv) {
	/* quad(x)=a*x^2+b*x+c */
	/* param[0] = a */
	/* param[1] = b */
	/* param[2] = c */
	/* param[3] = x0 */
	float xdif;

	xdif=x-param[3];
	if(deriv)	{
		deriv[0]=xdif*xdif;
		deriv[1]=xdif;
		deriv[2]=1;
		deriv[3]=-2*param[0]*xdif-param[1]; }
	return param[0]*xdif*xdif+param[1]*xdif+param[2]; }

float asinhbasis(float x,float *param,sptr spec,float *deriv) {
	/* asinh(x)=weight*asinh(slope*x+intercept) */
	/* param[0] = weight */
	/* param[1] = slope */
	/* param[2] = intercept */
	float f;
	
	f=param[1]*x+param[2];
	if(deriv)	{
		deriv[0]=asinh(f);
		deriv[1]=param[0]*x/sqrt(f*f+1);
		deriv[2]=param[0]/sqrt(f*f+1); }
	return param[0]*asinh(f); }

float gaussbasis(float x,float *param,sptr spec,float *deriv) {
	/* gauss(x)=area/(std_dev*Ã2¹)*exp(-(x-mean)^2/(2*std_dev^2)) */
	/* param[0] = area */
	/* param[1] = mean */
	/* param[2] = std_dev */
	float xdif,g;
	
	if(param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=0;
		return 0; }
	xdif=x-param[1];
 	g=1.0/(param[2]*sqrt(2*PI))*exp(-1.0/2*(xdif/param[2]*xdif/param[2]));
	if(deriv)	{
		deriv[0]=g;
		deriv[1]=param[0]*g*xdif/(param[2]*param[2]);
		deriv[2]=param[0]*g*(-1.0/param[2]+xdif/param[2]*xdif/param[2]/param[2]); }
	return param[0]*g; }

float xgaussbasis(float x,float *param,sptr spec,float *deriv) {
	/* xgauss(x)=x*area/(std_dev*Ã2¹)*exp(-(x-mean)^2/(2*std_dev^2)) */
	/* param[0] = area */
	/* param[1] = mean */
	/* param[2] = std_dev */
	float xdif,g;
	
	if(param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=0;
		return 0; }
	xdif=x-param[1];
 	g=x/(param[2]*sqrt(2*PI))*exp(-1.0/2*(xdif/param[2]*xdif/param[2]));
	if(deriv)	{
		deriv[0]=g;
		deriv[1]=param[0]*g*xdif/(param[2]*param[2]);
		deriv[2]=param[0]*g*(-1.0/param[2]+xdif/param[2]*xdif/param[2]/param[2]); }
	return param[0]*g; }

float sinbasis(float x,float *param,sptr spec,float *deriv) {
	/* sine(x)=amplitude*sin(frequency*x+shift) */
	/* param[0] = amplitude */
	/* param[1] = frequency */
	/* param[2] = shift */
	float f;
	
	f=param[1]*x+param[2];
	if(deriv)	{
		deriv[0]=sin(f);
		deriv[1]=param[0]*x*cos(f);
		deriv[2]=param[0]*cos(f); }
	return param[0]*sin(f); }

float lorentzbasis(float x,float *param,sptr spec,float *deriv) {
	/* lorentz(x)=max/(1+((x-mean)/(fwhm/2))^2) */
	/* param[0] = max */
	/* param[1] = mean */
	/* param[2] = fwhm */
	float a,x2,l;
	
	if(param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=0;
		return 0; }
	a=2/param[2];
	x2=x-param[1];
	l=1.0/(1.0+a*x2*a*x2);
	if(deriv)	{
		deriv[0]=l;
		deriv[1]=2.0*param[0]*l*l*a*a*x2;
		deriv[2]=param[0]*l*l*x2*a*x2*a*a; }
	return param[0]*l; }

float peakbasis(float x,float *param,sptr spec,float *deriv) {
	/* peak(x)=(1-shape)*xgaussian(x)+shape*xlorentzian(x) */
	/* see documentation for exact equation */
	/* param[0] = height */
	/* param[1] = position */
	/* param[2] = fwhm */
	/* param[3] = shape, 0=gauss, 1=lorentz */
	float g,l,g1,l1,x2,w,a,b,m,h,s;
	
	if(param[1]==0||param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=deriv[3]=0;
		return 0; }
	h=param[0];
	m=param[1];
	w=param[2];
	s=param[3];
	a=h*(1.0-s);
	b=h*s;
	x2=x-m;
	g=exp(-x2/w*x2/w*4*0.69314718056);
	l=1.0/(1.0+4.0*x2/w*x2/w);
	if(deriv)	{
		g1=-8.0*0.69314718056/(w*w)*x2*g;
		l1=-8.0/(w*w)*l*l*x2;
		deriv[0]=x/m*((1-s)*g+s*l);
		deriv[1]=-x/m*(a*g/m+a*g1+b*l/m+b*l1);
		deriv[2]=-x*x2/(m*w)*(a*g1+b*l1);
		deriv[3]=h*x/m*(-g+l); }
	return x/m*(a*g+b*l); }

float peak1basis(float x,float *param,sptr spec,float *deriv) {
	/* peakd1(x)=x weighted derivative of a peak */
	/* see documentation for exact equation */
	/* param[0] = height */
	/* param[1] = position */
	/* param[2] = fwhm */
	/* param[3] = shape, 0=gauss, 1=lorentz */
	float g,l,g1,l1,g2,l2,x2,w,a,b,m,h,s;
	
	if(param[1]==0||param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=deriv[3]=0;
		return 0; }
	h=param[0];
	m=param[1];
	w=param[2];
	s=param[3];
	a=h*(1.0-s);
	b=h*s;
	x2=x-m;
	g=exp(-x2/w*x2/w*4*0.69314718056);
	l=1.0/(1.0+4*x2/w*x2/w);
	g1=-8.0*0.69314718056/(w*w)*x2*g;
	l1=-8.0/(w*w)*l*l*x2;
	if(deriv)	{
		g2=-8.0*0.69314718056/(w*w)*(g1*x2+g);
		l2=-8.0/(w*w)*(2*l*l1*x2+l*l);
		deriv[0]=x/m*((1-s)*g1+s*l1);
		deriv[1]=-x/m*(a*g1/m+a*g2+b*l1/m+b*l2);
		deriv[2]=-x/(m*w)*(x2*(a*g2+b*l2)+(a*g1+b*l1));
		deriv[3]=h*x/m*(-g1+l1); }
	return x/m*(a*g1+b*l1); }

float peak2basis(float x,float *param,sptr spec,float *deriv) {
	/* peakd2(x)=x weighted second derivative of a peak */
	/* see documentation for exact equation */
	/* param[0] = height */
	/* param[1] = position */
	/* param[2] = fwhm */
	/* param[3] = shape, 0=gauss, 1=lorentz */
	float g,l,g1,l1,g2,l2,g3,l3,x2,w,a,b,m,h,s;
	
	if(param[1]==0||param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=deriv[3]=0;
		return 0; }
	h=param[0];
	m=param[1];
	w=param[2];
	s=param[3];
	a=h*(1.0-s);
	b=h*s;
	x2=x-m;
	g=exp(-x2/w*x2/w*4*0.69314718056);
	l=1.0/(1.0+4.0*x2/w*x2/w);
	g1=-8.0*0.69314718056/(w*w)*x2*g;
	l1=-8.0/(w*w)*l*l*x2;
	g2=-8.0*0.69314718056/(w*w)*(g1*x2+g);
	l2=-8.0/(w*w)*(2.0*l*l1*x2+l*l);
	if(deriv)	{
		g3=-8.0*0.69314718056/(w*w)*(g2*x2+2*g1);
		l3=-8.0/(w*w)*(2*l1*l1*x2+2*l*l2*x2+4*l*l1);
		deriv[0]=x/m*((1.0-s)*g2+s*l2);
		deriv[1]=-x/m*(a*g2/m+a*g3+b*l2/m+b*l3);
		deriv[2]=-x/(m*w)*(x2*(a*g3+b*l3)+2*(a*g2+b*l2));
		deriv[3]=h*x/m*(-g2+l2); }
	return x/m*(a*g2+b*l2); }

float peakzbasis(float x,float *param,sptr spec,float *deriv) {
	/* peakz(x)=z0*peak(x)+z1*peakd1(x)+z2*peakd2(x)	*/
	/* see documentation for exact equation */
	/* param[0] = z0 */
	/* param[1] = z1 */
	/* param[2] = z2 */
	/* param[3] = height */
	/* param[4] = position */
	/* param[5] = fwhm */
	/* param[6] = shape, 0=gauss, 1=lorentz */
	float g,l,g1,l1,g2,l2,g3,l3,x2,w,a,b,m,h,s,z0,z1,z2,zg,zl,zgoff,zloff;
	
	if(param[4]==0||param[5]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=deriv[3]=deriv[4]=deriv[5]=deriv[6]=0;
		return 0; }
	z0=param[0];
	z1=param[1];
	z2=param[2];
	h=param[3];
	m=param[4];
	w=param[5];
	s=param[6];
	a=h*(1.0-s);
	b=h*s;
	x2=x-m;
	g=exp(-x2/w*x2/w*4*0.69314718056);
	l=1.0/(1.0+4.0*x2/w*x2/w);
	g1=-8.0*0.69314718056/(w*w)*x2*g;
	l1=-8.0/(w*w)*l*l*x2;
	g2=-8.0*0.69314718056/(w*w)*(g1*x2+g);
	l2=-8.0/(w*w)*(2*l*l1*x2+l*l);
	zg=z0*g+z1*g1+z2*g2;
	zl=z0*l+z1*l1+z2*l2;
	if(deriv)	{
		g3=-8.0*0.69314718/(w*w)*(g2*x2+2*g1);
		l3=-8.0/(w*w)*(2*l1*l1*x2+2*l*l2*x2+4*l*l1);
		zgoff=z0*g1+z1*g2+z2*g3;
		zloff=z0*l1+z1*l2+z2*l3;
		deriv[0]=x/m*(a*g+b*l);
		deriv[1]=x/m*(a*g1+b*l1);
		deriv[2]=x/m*(a*g2+b*l2);
		deriv[3]=x/m*((1-s)*zg+s*zl);
		deriv[4]=-x/m*(a*zg/m+a*zgoff+b*zl/m+b*zloff);
		deriv[5]=-x/(m*w)*(x2*(a*zgoff+b*zloff)+a*(z1*g1+2*z2*g2)+b*(z1*l1+2*z2*l2));
		deriv[6]=h*x/m*(-zg+zl); }
	return x/m*(a*zg+b*zl); }

float diffusebasis(float x,float *param,sptr spec,float *deriv) {
	/* diffuse(x)=spect(x)*area*exp(-dt*x^2) */
	/* param[0] = area */
	/* param[1] = dt */
	int j=-2;
	float y,g;
	
	y=interpolate1(spec->x,spec->y,spec->n,&j,x);
 	g=y*exp(-param[1]*x*x);
	if(deriv)	{
		deriv[0]=g;
		deriv[1]=-param[0]*x*x*g; }
	return param[0]*g; }

float diffuse2basis(float x,float *param,sptr spec,float *deriv) {
	/* diffuse2(x)=c0*erf(Ã(td/(x-t0))/4)^2 */
	/* param[0] = c0 */
	/* param[1] = td */
	/* param[2] = t0 */
	float c0,d,t0,e;
	
	c0=param[0];
	d=param[1];
	t0=param[2];
	if(d*(x-t0)<=0)	{
		if(deriv)	{
			deriv[0]=1;
			deriv[1]=0;
			deriv[2]=0; }
		return c0; }
	else	{
		e=erfn(sqrt(d/(x-t0))/4);
		if(deriv)	{
			deriv[0]=e*e;
			deriv[1]=c0*e/2/sqrt(d*PI*(x-t0))*exp(-d/16/(x-t0));
			deriv[2]=c0*e/2*sqrt(d/(PI*(x-t0)*(x-t0)*(x-t0)))*exp(-d/16/(x-t0)); }
		return c0*e*e; }}

float convexpbasis(float x,float *param,sptr spec,float *deriv) {
	/* decay(x)=convolution of gaussian with truncated exponential */
	/* see documentation for exact equation */
	/* param[0] = height */
	/* param[1] = fwhm */
	/* param[2] = tau */
	/* param[3] = shift */
	float s,k,a,y,t;
	
	if(param[1]==0||param[2]==0) {
		if(deriv) deriv[0]=deriv[1]=deriv[2]=deriv[3]=0;
		return 0; }
	a=param[0];
	s=param[1]*0.424660900145;
	k=1/param[2];
	t=x-param[3];
	y=1.0/2*exp(-k*t+s*s*k*k/2)*(1+erfn((t-s*s*k)/(s*SQRT2)));
	if(deriv)	{
		deriv[0]=y;
		deriv[1]=a*0.424660900145*(-(k*s*s+t)/(s*SQRT2PI)*exp(-t*t/(2*s*s))+k*k*s*y);
		deriv[2]=a*k*k*(s/SQRT2PI*exp(-t*t/(2*s*s))-(k*s*s-t)*y);
		deriv[3]=a*(-1/(s*SQRT2PI)*exp(-t*t/(2*s*s))+k*y); }
	return a*y; }

float rationbasis(float x,float *param,sptr spec,float *deriv) {
	/* y=(n0+n1*x+n2*x^2+n3*x^3+n4*x^4+n5*x^5)/(d0+d1*x+d2*x^2+d3*x^3+d4*x^4+d5*x^5) */
	/* param[0] = n0 */
	/* param[1] = n1 */
	/* param[2] = n2 */
	/* param[3] = n3 */
	/* param[4] = n4 */
	/* param[5] = n5 */
	/* param[6] = d0 */
	/* param[7] = d1 */
	/* param[8] = d2 */
	/* param[9] = d3 */
	/* param[10] = d4 */
	/* param[11] = d5 */
	float num,den,*d,*p;

	p=param;
	d=deriv;
	num=p[0]+x*(p[1]+x*(p[2]+x*(p[3]+x*(p[4]+x*p[5]))));
	den=p[6]+x*(p[7]+x*(p[8]+x*(p[9]+x*(p[10]+x*p[11]))));
	if(den==0) {
		if(d) d[0]=d[1]=d[2]=d[3]=d[4]=d[5]=d[6]=d[7]=d[8]=d[9]=d[10]=d[11]=0;
		return 0; }
	if(deriv) {
		d[5]=x*(d[4]=x*(d[3]=x*(d[2]=x*(d[1]=x*(d[0]=1.0/den)))));
		d[11]=x*(d[10]=x*(d[9]=x*(d[8]=x*(d[7]=x*(d[6]=-num/den/den))))); }
	return num/den; }

float sigmoidbasis(float x,float *param,sptr spec,float *deriv) {
	/* y=b+(t-b)*f(x) with f(x)=1/(1+10^(h*(e-x))) */
	/* param[0] = b = min    */
	/* param[1] = t = max    */
	/* param[2] = e = ec50   */
	/* param[3] = h = slope  */
	float b,t,e,h,f;
	
	b=param[0];
	t=param[1];
	e=param[2];
	h=param[3];
	f=1.0/(1.0+exp(log(10)*h*(e-x)));
	if(deriv) {
		deriv[0]=1.0-f;
		deriv[1]=f;
		deriv[2]=h*(b-t)*f*(1.0-f)*log(10);
		deriv[3]=(e-x)*(b-t)*f*(1.0-f)*log(10); }
	return b+(t-b)*f; }
	
float hillbasis(float x,float *param,sptr spec,float *deriv) {
	/* y=a*x^n/(x^n+e^n)            */
	/* param[0] = a = amplitude     */
	/* param[1] = e = ec50          */
	/* param[2] = n = cooperativity */
	float a,e,n,f;

	a=param[0];
	e=param[1];
	n=param[2];
	if(x<=0) f=0;
	else f=1.0/(pow(e/x,n)+1.0);
	if(deriv) {
		deriv[0]=f;
		deriv[1]=-a*f*f*n*pow(e/x,n-1)/x;
		deriv[2]=-a*f*f*log(e/x)*pow(e/x,n); }
	return a*f; }

float linesbasis(float x,float *param,sptr spec,float *deriv) {
	/* y=slope a0 to x1,y1, then to x2,y2, then to x3,y3, then slope a3 */
	/* param[0] = a0 */
	/* param[1] = x1 */
	/* param[2] = y1 */
	/* param[3] = x2 */
	/* param[4] = y2 */
	/* param[5] = x3 */
	/* param[6] = y3 */
	/* param[7] = a3 */
	float a0,x1,y1,x2,y2,x3,y3,a3,y;

	a0=param[0];
	x1=param[1];
	y1=param[2];
	x2=param[3];
	y2=param[4];
	x3=param[5];
	y3=param[6];
	a3=param[7];
	if(deriv)
		deriv[0]=deriv[1]=deriv[2]=deriv[3]=deriv[4]=deriv[5]=deriv[6]=deriv[7]=0;

	if(x<x1) {
		y=a0*(x-x1)+y1;
		if(deriv) {
			deriv[0]=x-x1;
			deriv[1]=-a0;
			deriv[2]=1; }}
	else if(x<x2) {
		y=(y2-y1)/(x2-x1)*(x-x1)+y1;
		if(deriv) {
			deriv[1]=(y2-y1)*(x-x2)/(x2-x1)/(x2-x1);
			deriv[2]=1-(x-x1)/(x2-x1);
			deriv[3]=-(y2-y1)*(x-x1)/(x2-x1)/(x2-x1);
			deriv[4]=(x-x1)/(x2-x1); }}
	else if(x<x3) {
		y=(y3-y2)/(x3-x2)*(x-x2)+y2;
		if(deriv) {
			deriv[3]=(y3-y2)*(x-x3)/(x3-x2)/(x3-x2);
			deriv[4]=1-(x-x2)/(x3-x2);
			deriv[5]=-(y3-y2)*(x-x2)/(x3-x2)/(x3-x2);
			deriv[6]=(x-x2)/(x3-x2); }}
	else {
		y=a3*(x-x3)+y3;
		if(deriv) {
			deriv[5]=-a3;
			deriv[6]=1;
			deriv[7]=x-x3; }}
	return y; }