//		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"
complex radius_cp(t_s & y,complex & order)
	{
	register int m;
	complex rm0,rm1,rm2,rm3,rm4,rcs,c;
	complex t1,t2,t3,t4;
		
	m = y.last_term - 1;
	rm0 = y.terms[m]/y.terms[m-1];
	cout << "rm0 = " << rm0 << "\n";
	rm1 = y.terms[m-1]/y.terms[m-2];
	cout << "rm1 = " << rm1 << "\n";
	rm2 = y.terms[m-2]/y.terms[m-3];
	cout << "rm2 = " << rm2 << "\n";
	rm3 = y.terms[m-3]/y.terms[m-4];
	cout << "rm3 = " << rm3 << "\n";
	rm4 = y.terms[m-4]/y.terms[m-5];
	cout << "rm4 = " << rm4 << "\n";
	t1 = complex(12.0,0.0)*rm2*rm1*rm1*rm4;
	t1 = t1 -complex(6.0,0.0)*rm3*rm3*rm2*rm4;
	t1 = t1 -complex(10.0,0.0)*rm1*rm1*rm2*rm3;
	t1 = t1 +complex(8.0,0.0)*rm2*rm2*rm3*rm4;
	t1 = t1 +complex(10.0,0.0)*rm2*rm2*rm1*rm3;
	t1 = t1 +complex(18.0,0.0)*rm1*rm3*rm3*rm4;
	t1 = t1 +complex(25.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
	t1 = t1 -complex(5.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm3;
	t1 = t1 +complex(12.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
	t1 = t1 -complex(7.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm1*rm3;
	t1 = t1 -complex(32.0,0.0)*rm2*rm1*rm3*rm4;
	t1 = t1 -complex(8.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
	t1 = t1 +complex(5.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm2*rm4;
	t1 = t1 -complex(12.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
	t1 = t1 -complex(15.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
	t1 = t1 +complex(8.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
	t1 = t1 -complex(3.0,0.0)*complex(double(m),0.0)*rm0*rm1*rm3*rm4;
	t1 = t1 +rm1*rm3*rm4*rm0*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 +rm2*rm1*rm3*rm0*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 -complex(2.0,0.0)*rm2*rm1*rm4*rm0*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 -complex(2.0,0.0)*rm1*rm1*rm2*rm3*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 +complex(3.0,0.0)*rm2*rm1*rm1*rm4*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 +rm2*rm2*rm1*rm3*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 -complex(5.0,0.0)*rm2*rm1*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 +complex(2.0,0.0)*rm1*rm3*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 +complex(2.0,0.0)*rm2*rm2*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
	t1 = t1 -rm3*rm3*rm2*rm4*complex(double(m),0.0)*complex(double(m),0.0);
	t2 = -complex(5.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
	t2 = t2 +complex(3.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
	t2 = t2 +complex(10.0,0.0)*rm2*rm1*rm3*rm4;
	t2 = t2 -complex(3.0,0.0)*rm2*rm1*rm1*rm4;
	t2 = t2 +complex(2.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
	t2 = t2 -complex(4.0,0.0)*rm2*rm2*rm3*rm4;
	t2 = t2 +complex(2.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
	t2 = t2 -complex(double(m),0.0)*rm3*rm3*rm2*rm4;
	t2 = t2 -complex(6.0,0.0)*rm1*rm3*rm3*rm4;
	t2 = t2 +complex(3.0,0.0)*rm3*rm3*rm2*rm4;
	t2 = t2 -complex(2.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
	t2 = t2 +complex(double(m),0.0)*rm0*rm1*rm3*rm4;
	t2 = t2 +complex(double(m),0.0)*rm0*rm2*rm1*rm3;
	t2 = t2 -complex(2.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
	t2 = t2 +complex(2.0,0.0)*rm1*rm1*rm2*rm3;
	t2 = t2 +complex(double(m),0.0)*rm2*rm2*rm1*rm3;
	t2 = t2 -complex(2.0,0.0)*rm2*rm2*rm1*rm3;
	order = -t1/(complex(2.0,0.0)*t2);
	t3 = rm3*rm4;
	t3 = t3 -complex(4.0,0.0)*rm2*rm4;
	t3 = t3 +rm2*rm1;
	t3 = t3 -complex(4.0,0.0)*rm1*rm3;
	t3 = t3 +complex(3.0,0.0)*rm1*rm4;
	t3 = t3 +complex(3.0,0.0)*rm2*rm3;
	t4 = -complex(5.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
	t4 = t4 +complex(3.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
	t4 = t4 +complex(10.0,0.0)*rm2*rm1*rm3*rm4;
	t4 = t4 -complex(3.0,0.0)*rm2*rm1*rm1*rm4;
	t4 = t4 +complex(2.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
	t4 = t4 -complex(4.0,0.0)*rm2*rm2*rm3*rm4;
	t4 = t4 +complex(2.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
	t4 = t4 -complex(double(m),0.0)*rm3*rm3*rm2*rm4;
	t4 = t4 -complex(6.0,0.0)*rm1*rm3*rm3*rm4;
	t4 = t4 +complex(3.0,0.0)*rm3*rm3*rm2*rm4;
	t4 = t4 -complex(2.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
	t4 = t4 +complex(double(m),0.0)*rm0*rm1*rm3*rm4;
	t4 = t4 +complex(double(m),0.0)*rm0*rm2*rm1*rm3;
	t4 = t4 -complex(2.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
	t4 = t4 +complex(2.0,0.0)*rm1*rm1*rm2*rm3;
	t4 = t4 +complex(double(m),0.0)*rm2*rm2*rm1*rm3;
	t4 = t4 -complex(2.0,0.0)*rm2*rm2*rm1*rm3;
	rcs  = -t3/t4;
	c = my_sqrt(rcs)*y.h;
	return(c);
	}
complex radius_xx(t_s & y,complex & order)
	{
	register int m;
	complex rm0,rm1,hdrc,rcs;
		
	m = y.last_term - 1;
	rm0 = y.terms[m]/y.terms[m-1];
	cout << "rm0 = " << rm0 << "\n";
	rm1 = y.terms[m-1]/y.terms[m-2];
	cout << "rm1 = " << rm1 << "\n";
	hdrc = complex(double(m),0.0)*rm0-complex(double(m-1),0.0)*rm1;
	cout << "hdrc = " << hdrc << "\n";
	rcs = y.h/hdrc;
	cout << "rcs = " << rcs << "\n";
	order = complex(double(m),0.0)*rm0*rcs/y.h-(complex(double(m-1),0.0));
	cout << "order = " << order << "\n";
	return(rcs);
	}