//		ret.terms[i]=complex(0.0,0.0);
//
//	C++ Routines to manipulate Taylor series for solving a differential
//  equation.
//
//  By Dennis J. Darland 
//  4/3/91
//  Public Domain
//  Version 1.6
//  Many enhancements expected.
//  PROVIDED WITH NO WARRANTY.
//  
#include "ats.h"
t_s::t_s()
	{
	alloc_terms=MAX_TERMS;
	}
t_s::~t_s()
	{
	}
t_s::t_s(t_s &a)
	{
	alloc_terms=MAX_TERMS;
	set_first(a.first_term);
	set_last(a.last_term);
	set_x0(a.x0);
	set_h(a.h);
	set_terms(*this,a.terms);
	}
void t_s::set_first(int i)
	{
	first_term=i;
	}
void t_s::set_last(int i)
	{
	last_term=i;
	}
void t_s::set_h(complex h_in)
	{
	h=h_in;
	}
void t_s::set_x0(complex x0_in)
	{
	x0=x0_in;
	}
void set_data(t_s &a,t_s &b)
	{
	a.set_first(b.first_term);
	a.set_last(b.last_term);
	a.set_x0(b.x0);
	a.set_h(b.h);
	a.alloc_terms=MAX_TERMS;
	}
void set_terms(t_s &ret,complex *values)
	{
	register int i;
	
	for (i=0;i<MAX_TERMS;i++)
		ret.terms[i]=values[i];
	}
complex ats(int m,t_s &a,t_s &b,int j)
	{
	// function to implement Liebniz rule
	complex cret;
	register int i;
	
	cret=complex(0.0,0.0);
	
	for (i=j;i<=m;i++)
		{
		cret += a.terms[i]*b.terms[m-i];
		}
	return(cret);
	}

complex att(int m,t_s &a,t_s &b,int j)
	{
	// function for modified Liebniz rule
	complex ret;
	register int i;
	
	ret=complex(0.0,0.0);
	if (j > m)
		return(ret);	
	for (i=j;i<=m;i++)
		{
		ret += a.terms[i]*b.terms[m-i+1]*complex(double(m-i+1),0.0);
		}
	ret /= complex(double(m+1),0.0);
	return(ret);
	}
t_s & t_s::operator+=(t_s &b)
	{
	// Add taylor series
	register int i;
	
	if ((first_term != b.first_term)
	|| (last_term != b.last_term))
		report_error(0);
	set_data(*this,b);
	for (i=first_term; i < last_term;i++)
		terms[i]+=b.terms[i];
	return (*this);
	}
t_s & t_s::operator-=(t_s &b)
	{
	// Subtract taylor series
	register int i;
	
	if ((first_term != b.first_term)
	|| (last_term != b.last_term))
		report_error(0);
	set_data(*this,b);
	for (i=first_term; i < last_term;i++)
		terms[i]-=b.terms[i];
	return (*this);
	}
t_s& t_s::operator=(t_s &b)
	{
	// assign taylor series
	register int i;
	set_data(*this,b);
	for (i=0; i < b.last_term;i++)
		this->terms[i]=b.terms[i];
	return(*this);
	}
void display(t_s &a,int plot)
	{
	static pt=0;
	FILE *plot_cmds;
	char fname[64];
	char cmd[81];
	
	// display taylor series
	register int i;
	
	sprintf(fname,"ram:mapleplot_%d",pt);
	pt++;
	plot_cmds = fopen("ram:plottmp","w");
	fprintf(plot_cmds,"plotdevice := amiga;\n");
	fprintf(plot_cmds,"plot([");
	cout << "alloc_terms=" << a.alloc_terms << "\n"; 
	cout << "first_term =" << a.first_term << "\n"; 
	cout << "last_term=" << a.last_term << "\n";
	cout << "h=" << a.h << "\n";
	cout << "x0=" << a.x0 << "\n";
	for (i=0;i<a.last_term;i++)
		{
		cout << "terms[" << i << "]= " << real(a.terms[i]) << " + "
		<< imag(a.terms[i]) << " I\n";
		if (i != 0)
			fprintf(plot_cmds,",");
		fprintf(plot_cmds,"%d,sign(%g)*(ln(1+abs(%g)))\n",
		i,real(a.terms[i]),real(a.terms[i]));
		}
	fprintf(plot_cmds,"]);\n");
	fprintf(plot_cmds,"quit\n");
	fclose(plot_cmds);
	sprintf(cmd,"iconz <ram:plottmp >%s pltcvt",fname);
	system(cmd);
	if (plot>1)
		{
		sprintf(cmd,"maple <%s",fname);
		system(cmd);
		}
	}
void t_s::report_error(int no)
	{
	cout << "ats error = " << no;
	}
t_s operator+(t_s &a,t_s &b)
	{
	// Add two taylor series
	register int i;
	t_s ret;
	complex ac,bc;
	
	set_data(ret,a);
	for (i=a.first_term;i<a.last_term;i++)
		{
		ac=a.terms[i];
		bc=b.terms[i];
		ret.terms[i]=ac+bc;
		}
	return(ret);
	}

t_s operator-(t_s &a,t_s &b)
	{
	// Subtract two taylor series
	register int i;
	t_s ret;
	complex ac,bc;
		
	set_data(ret,a);
	
	for (i=a.first_term;i<a.last_term;i++)
		{
		ac=a.terms[i];
		bc=b.terms[i];
		ret.terms[i]=ac-bc;
		}
	return(ret);
	}

t_s operator*(t_s &a,t_s &b)
	{
	// Multiply two taylor series
	register int i;
	t_s ret;
	
	set_data(ret,a);
	ret.terms[0]=a.terms[0]*b.terms[0];
	for (i=1;i<a.last_term;i++)
		ret.terms[i]=ats(i,a,b,0);
	return(ret);
	}


t_s operator/(t_s &a,t_s &b)
	{
	// divide two taylor series
	register int i;
	t_s ret;
	
	set_data(ret,a);
	ret.terms[0]=a.terms[0]/b.terms[0];
	for (i=1;i<a.last_term;i++)
		{
		ret.terms[i]=(a.terms[i]-ats(i,b,ret,1))/b.terms[0];
		}
	return(ret);
	}

t_s x(t_s & a)
	{
	// independant variable - x
	register int i;
	t_s ret;
	
	set_data(ret,a);
	ret.terms[0]=a.x0;
	ret.terms[1]=complex(1.0,0.0)*ret.h;
	for (i=2;i<a.last_term;i++)
		ret.terms[i]=complex(0.0,0.0);
	return(ret);
	}

t_s constant(t_s & a,complex b)
	{
	// constant - b
	register int i;
	t_s ret;
	
	set_data(ret,a);
	ret.terms[0]=b;
	for (i=1;i<a.last_term;i++)
		ret.terms[i]=complex(0.0,0.0);
	return(ret);
	}

t_s expt(t_s & y,t_s & y1)
	{
	// y to power y1
	register int k,ka;
	t_s ret,a1,a2;
	
	set_data(ret,y);
	set_data(a1,y);
	set_data(a2,y);
	ret.terms[0]=expt(y.terms[0],y1.terms[0]);
	a1.terms[0]=lnn(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		a1.terms[k]=(y.terms[k] - att(ka,y,a1,1))/y.terms[0];
		a2.terms[ka]=ats(k,y,a1,0)*complex(double(k),0.0)/y.h;
		ret.terms[k]=ats(ka,ret,a2,0)*y.h/complex(double(k),0.0);
		}
	return(ret);
	}

t_s sqrt(t_s & y)
	{
	// square root
	register int k;
	t_s ret;
	
	set_data(ret,y);
	ret.terms[0]=my_sqrt(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ret.terms[k]=complex(0.0,0.0);
		ret.terms[k]=(y.terms[k]-ats(k,ret,ret,1))/(ret.terms[0]*
		complex(2.0,0.0));
		}
	return(ret);
	}
	
t_s exp(t_s & y)
	{
	// e to y power
	register int k,ka;
	t_s ret;
	
	set_data(ret,y);
	ret.terms[0]=exp(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		ret.terms[k]=att(ka,ret,y,0);
		}
	return(ret);
	}
t_s log(t_s & y)
	{
	// ln of y
	register int k,ka;
	t_s ret;
	
	set_data(ret,y);
	ret.terms[0]=lnn(y.terms[0]);
	ret.terms[1]=y.terms[1]/y.terms[0];
	for (k=2;k<y.last_term;k++)
		{
		ka = k - 1;
		ret.terms[k]=(y.terms[k] - att(ka,y,ret,1))/y.terms[0];
		}
	return(ret);
	}

t_s sin(t_s & y)
	{
	// sin of y
	register int k,ka;
	t_s f,g;
	
	set_data(f,y);
	set_data(g,y);
	f.terms[0]=sin(y.terms[0]);
	g.terms[0]=cos(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=att(ka,g,y,0);
		g.terms[k]=complex(-1.0,0.0)*att(ka,f,y,0);
		}
	return(f);
	}
	
t_s cos(t_s & y)
	{
	// cos of y
	register int k,ka;
	t_s f,g;
	
	set_data(f,y);
	set_data(g,y);
	f.terms[0]=sin(y.terms[0]);
	g.terms[0]=cos(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=att(ka,g,y,0);
		g.terms[k]=complex(-1.0,0.0)*att(ka,f,y,0);
		}
	return(g);
	}
	
t_s tan(t_s & y)
	{
	// cos of y
	register int k,ka;
	t_s f,a1,a2;
	
	set_data(f,y);
	set_data(a1,y);
	set_data(a2,y);
	a1.terms[0]=sin(y.terms[0]);
	a2.terms[0]=cos(y.terms[0]);
	f.terms[0]=a1.terms[0]/a2.terms[0];
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		a1.terms[k]=att(ka,a2,y,0);
		a2.terms[k]=complex(-1.0,0.0)*att(ka,a1,y,0);
		f.terms[k]=(a1.terms[k] - ats(k,a2,f,1))/a2.terms[0];
		}
	return(f);
	}
	
t_s sinh(t_s & y)
	{
	// sinh of y
	register int k,ka;
	t_s f,g;
	
	set_data(f,y);
	set_data(g,y);
	f.terms[0]=sinh(y.terms[0]);
	g.terms[0]=cosh(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=att(ka,g,y,0);
		g.terms[k]=att(ka,f,y,0);
		}
	return(f);
	}
	
t_s cosh(t_s & y)
	{
	// cosh of y
	register int k,ka;
	t_s f,g;
	
	set_data(f,y);
	set_data(g,y);
	f.terms[0]=sinh(y.terms[0]);
	g.terms[0]=cosh(y.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=att(ka,g,y,0);
		g.terms[k]=att(ka,f,y,0);
		}
	return(g);
	}
	
t_s tanh(t_s & y)
	{
	// tanh of y
	register int k,ka;
	t_s f,a1,a2;
	
	set_data(f,y);
	set_data(a1,y);
	set_data(a2,y);
	a1.terms[0]=sinh(y.terms[0]);
	a2.terms[0]=cosh(y.terms[0]);
	f.terms[0]=a1.terms[0]/a2.terms[0];
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		a1.terms[k]=att(ka,a2,y,0);
		a2.terms[k]=att(ka,a1,y,0);
		f.terms[k]=(a1.terms[k] - ats(k,a2,f,1))/a2.terms[0];
		}
	return(f);
	}
	
t_s asin(t_s & y)
	{
	// asin of y
	register int k,ka;
	t_s f,a1;
	
	set_data(f,y);
	set_data(a1,y);
	f.terms[0]=asin(y.terms[0]);
	a1.terms[0]=cos(f.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=(y.terms[k] - att(ka,a1,f,1))/a1.terms[0];
		a1.terms[k]=complex(-1.0,0.0)*att(ka,y,f,0);
		}
	return(f);
	}
	
t_s acos(t_s & y)
	{
	// acos of y
	register int k,ka;
	t_s f,a1;
	
	set_data(f,y);
	set_data(a1,y);
	f.terms[0]=acos(y.terms[0]);
	a1.terms[0]=sin(f.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=complex(-1.0,0.0)*
		(y.terms[k] + att(ka,a1,f,1))/a1.terms[0];
		a1.terms[k]=att(ka,y,f,0);
		}
	return(f);
	}
	
t_s atan(t_s & y)
	{
	// atan of y
	register int k,ka;
	t_s f,a1,a2;
	
	set_data(f,y);
	set_data(a1,y);
	set_data(a2,y);
	f.terms[0]=atan(y.terms[0]);
	a1.terms[0]=sin(f.terms[0]);
	a2.terms[0]=cos(f.terms[0]);
	for (k=1;k<y.last_term;k++)
		{
		ka = k - 1;
		f.terms[k]=(ats(k,y,a2,1)-y.terms[0]*att(ka,a1,f,1) -
			att(ka,a2,f,1))/(a2.terms[0]+y.terms[0]*a1.terms[0]);
		a1.terms[k]=att(ka,a2,f,0);
		a2.terms[k]=complex(-1.0,0.0)*att(ka,a1,f,0);
		}
	return(f);
	}

t_s deriv(t_s &a) 
	{
	// derivative of a taylor series
	register int i;
	t_s ret;
		
	set_data(ret,a);
	
	for (i=a.first_term;i<a.last_term-1;i++)
		{
		if (i > 0)
			ret.terms[i]=a.terms[i+1]*complex(double(i),0.0)/a.h;
		else 
			ret.terms[i]=a.terms[i+1]/a.h;
		}
	ret.terms[a.last_term]=complex(0.0,0.0); // unknown but small
	return(ret);
	}
void diff_eq(t_s &ret2,t_s &y)
	{
	// solve 1st order ordinary differential equation provided in "equation".
	register int k;
	t_s work;
	t_s work2;
	
	work=y;
	work2=y;
	for (k=work.last_term;k<MAX_TERMS;k++)
		{
		// get derivative t_s (work) from function t_s work2
		equation(work,work2); 
		set_term_from_deriv(k,work,work2);
		}
	ret2 = work2;
	}
void higher_diff_eq(t_s &ret2,t_s &y,int n)
	{
	// solve 1st order ordinary differential equation provided in "equation".
	register int k;
	t_s work;
	t_s work2;
	
	work=y;
	work2=y;
	for (k=work.last_term;k<MAX_TERMS;k++)
		{
		// get derivative t_s (work) from function t_s work2
		equation(work,work2); 
		higher_set_term_from_deriv(k,work,work2,n);
		}
	ret2 = work2;
	}
void set_term_from_deriv(int i,t_s &dy,t_s &y)
	{
	if (i > 0) 
		y.terms[i]=dy.terms[i-1]*y.h/complex(double(i),0.0);
	else 
		y.terms[i]=dy.terms[i-1]*y.h;
	y.last_term++;
	}
void higher_set_term_from_deriv(int i,t_s &dy,t_s &y,int n)
	{
	complex adj1,adj2;
	
	adj1 = complex(1.0,0.0);
	adj2 = complex(1.0,0.0);
	
	for (int a=0;a<n;a++)
		adj1 *= dy.h;
	for (int b=0;b<n;b++)
		{
		if (i+b > 0)
			{
			adj2 *= complex(double(i+b),0.0);
			}
		}
	y.terms[i]=dy.terms[i-n]*adj1/adj2;
	y.last_term++;
	}
void jump(t_s &y)
	{
	t_s y2;
	
	y2 	= y;
	y2.terms[0] = complex(0.0,0.0);
	for (int i = 0;i < y.last_term;i++)
		{
		y2.terms[0] += y.terms[i];
		}
	y2.x0 = y.x0 + y.h;
	y2.last_term = 1;
	y = y2;
	}

void leap(t_s &y,int n)
	{
	t_s y2;
	complex adj1,adj2;
		
	y2 	= y;

	adj2 = complex(1.0,0.0); 
	for (int j = 0;j < n; j++)
		{
		y2.terms[j] = complex(0.0,0.0);
		for (int i = 0;i < y.last_term-j;i++)
			{
			adj1 = complex(1.0,0.0);
			for (int k = 1 ; k <= j; k++)
				{
				adj1 *= complex(double(i+k),0.0);
				}
			y2.terms[j] += y.terms[i+j]*adj1;
			}
		if (j > 0)
			{
			adj2 *= j;
			y2.terms[j] /= adj2;
			}
		}
	y2.x0 = y.x0 + y.h;
	y2.last_term = n;
	y = y2;
	}

int get_last_term(t_s &y)
	{
	return(y.last_term);
	}
	
double get_x0(t_s &y)
	{
	return(real(y.x0));
	}
double get_y0(t_s &y)
	{
	return(real(y.terms[0]));
	}

void adjust_h(t_s &y,complex new_h)
	{
	t_s y2;
	complex t,adj;
	
	t = new_h / y.h;
	adj = t;
	y2 	= y;
	for (int i = 1;i < y.last_term;i++)
		{
		y2.terms[i] = y.terms[i]*adj;
		adj *= t;
		}
	y2.h = new_h;
	y = y2;
	}