//
//	angular function, pi, tau by using recursive call of function
//
//	compile:
//	c++ pi_tau.cpp
//
//	calc speed is significantly decreased beyond n=30
//

#include <iostream>
#include <cmath>

using std::cout;
using std::cerr;
using std::endl;

const double PI = 2.0*asin(1.0);

double pi(const int n, const double x)
//
//	n must be n>=1
//	but n==0 is necessary for the computation
//
{
	if(n == 0)
		return 0.0;
	if(n == 1)
		return 1.0;
	return (2.0*n-1.0)/(n-1.0) * x * pi(n-1,x) - n/(n-1.0) * pi(n-2,x);
}

double tau(const int n, const double x)
//
//	n must be n>=1
//
{
	return n*x*pi(n,x)-(n+1)*pi(n-1,x);
}


int main(int argc, char* argv[])
{
	const int Ntheta = 91;	//	div of theta

	const int MAX_N = 20;	//	maximum value of n

	cout << "pi_n(cos(theta))\n\n";
	cout << ",n=,";
	for(int n=1; n<=MAX_N; n++) {
		cout << n << ',';
	}

	const double dtheta = PI/2 / (Ntheta-1);	//	angular resolution

	cout << "\ntheta(deg)\n";
	for(double theta=0; theta<=PI/2+dtheta/1000; theta += dtheta) {

		const double x = cos(theta);

		cout << theta*180/PI << ",,";
		for(int n=1; n<=MAX_N; n++) {
			cout << pi(n, x) << ",";
		}
		cout << "\n";
	}
	cout << "\n\n";


	cout << "tau_n(cos(theta))\n\n";
	cout << ",n=,";
	for(int n=1; n<=MAX_N; n++) {
		cout << n << ',';
	}

	cout << "\ntheta(rad)\n";
	for(double theta=0; theta<=PI/2+dtheta/1000; theta += dtheta) {

		const double x = cos(theta);

		cout << theta*180/PI << ",,";
		for(int n=1; n<=MAX_N; n++) {
			cout << tau(n, x) << ",";
		}
		cout << "\n";
	}

	return 0;

}