#include "iostream.h"
#include "math.h"
#include "stdlib.h"

double c2, factr;
int m, n;

int sgn(double f)
{
	if (f>0)
	{
		return 1;
	}
	if (f<0)
	{
		return -1;
	}
	return 0;
}

void lubksb(double a[], int n, int an,int indx[], double b[])
{
	int i,j,ll;
	double sum;
    int ii = 0;
    for (i = 1; i<=n; i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0 )
		{
            for (j = ii; j<=i-1; j++)
			{
                sum = sum - a[(i-1)*an+j] * b[j];
            }
		}
        else
		{
			if (sum != 0.0)
			{
				ii = i;
			}
        }
        b[i] = sum;
    }
    for (i = n; i>=1; i--)
	{
        sum = b[i];
        if (i < n)
		{
            for (j = i + 1; j<=n; j++)
			{
                sum = sum - a[(i-1)*an+j] * b[j];
            }
        }
        b[i] = sum / a[(i-1)*an+i];
    }
}

void ludcmp(double a[], int n, int an,int indx[], int& d)
{
	const int nmax=100;
	const double tiny=1e-20;
	double vv[100];
	double sum,dum,aamax;
	int i,j,k,imax;
	d=1;
	for (i=1; i<=n; i++)
	{
		aamax=0.0;
		for (j=1; j<=n;j++)
		{
			if (fabs(a[(i-1)*an+j])>aamax)
			{
				aamax=fabs(a[(i-1)*an+j]);
			}
		}
		if (aamax==0.0)
		{
			cout<<"singular matrix"<<endl;
		}
		vv[i]=1.0/aamax;
	}
	for (j=1; j<=n;j++)
	{
        if (j > 1)
		{
            for (i=1; i<=j-1;i++)
			{
                sum=a[(i-1)*an+j];
                if (i>1)
				{
                    for (int k=1; k<=i-1; k++)
					{
                        sum=sum-a[(i-1)*an+k]*a[(k-1)*an+j];
                    }
                    a[(i-1)*an+j]=sum;
                }
            }
        }
        aamax=0.0;
        for (i=j; i<=n; i++)
		{
            sum=a[(i-1)*an+j];
            if (j>1)
			{
                for (k=1; k<=j-1; k++)
				{
                    sum=sum-a[(i-1)*an+k]*a[(k-1)*an+j];
				}
                a[(i-1)*an+j]=sum;
            }
            dum=vv[i]*fabs(sum);
            if (dum>=aamax)
			{
                imax=i;
                aamax=dum;
            }
        }
        if (j!=imax)
		{
            for (k=1; k<=n; k++)
			{
                dum=a[(imax-1)*an+k];
                a[(imax-1)*an+k]=a[(j-1)*an+k];
                a[(j-1)*an+k]=dum;
            }
            d=-d;
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (j!=n)
		{
            if (a[(j-1)*n+j]==0.0)
			{
				a[(j-1)*n+j]=tiny;
			}
            dum=1.0/a[(j-1)*an+j];
            for (i=j+1; i<=n; i++)
			{
                a[(i-1)*an+j]=a[(i-1)*an+j]*dum;
            }
		}
    }
    if (a[(n-1)*an+n]==0.0)
	{
		a[(n-1)*an+n]=tiny;
	}
}
void load(double x1, double v[], double y[])
{
	double dx=0;
    y[3] = v[1];
    y[2] = -(y[3] - c2) * factr / 2.0 / (m + 1.0);
    y[1] = factr + y[2] * dx;
    y[6] = v[2];
    y[5] = -(y[6] + c2) * factr / 2.0 / (m + 1.0);
    y[4] = factr + y[5] * dx;
}

void score(double x2, double y[], double f[])
{
    if (((n - m) % 2) == 0)
	{
        f[1] = y[2];
        f[2] = y[5];
	}
    else
	{
        f[1] = y[1];
        f[2] = y[4];
	}
}

void derivs(double x, double y[], double dydx[])
{
    dydx[1] = y[2];
    dydx[3] = 0.0;
    dydx[2] = (2.0 * x * (m + 1.0) * y[2] - (y[3] - c2 * x * x) * y[1]) / (1.0 - x * x);
    dydx[4] = y[5];
    dydx[6] = 0.0;
    dydx[5] = (2.0 * x * (m + 1.0) * y[5] - (y[6] + c2 * x * x) * y[4]) / (1.0 - x * x);
}

void rk4(double y[], double dydx[], int n, double x, double h, double yout[])
{
	double yt[11], dyt[11], dym[11];
	double hh = h * 0.5;
	double h6 = h / 6.0;
	double xh = x + hh;
	int i;
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + hh * dydx[i];
	}
	derivs(xh, yt, dyt);
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + hh * dyt[i];
	}
	derivs(xh, yt, dym);
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + h * dym[i];
		dym[i] = dyt[i] + dym[i];
	}
	derivs(x + h, yt, dyt);
	for (i = 1; i<=n; i++)
	{
		yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
	}
	//erase dym, dyt, yt
}

void rkqc(double y[], double dydx[], int n, double& x, double htry, double eps,
		  double yscal[], double& hdid, double& hnext)
{
	double one = 1.0;
	double safety = 0.9;
	double errcon = 0.0006;
	double fcor = 0.066667;
	double ytemp[11], ysav[11], dysav[11];
	double pgrow = -0.2;
	double pshrnk = -0.25;
	double xsav = x;
	int i,flag;
	for (i = 1; i<=n; i++)
	{
		ysav[i] = y[i];
		dysav[i] = dydx[i];
	}
	double hh,h = htry;
	do
	{
		hh = 0.5 * h;
		rk4(ysav, dysav, n, xsav, hh, ytemp);
		x = xsav + hh;
		derivs(x, ytemp, dydx);
		rk4(ytemp, dydx, n, x, hh, y);
		x = xsav + h;
		if (x == xsav)
		{
			cout<<" stepsize not significant in rkqc."<<endl;
			return;
		}
		rk4(ysav, dysav, n, xsav, h, ytemp);
		double errmax = 0.0;
		for (i = 1; i<=n; i++)
		{
			ytemp[i] = y[i] - ytemp[i];
			if (errmax < fabs(ytemp[i] / yscal[i]))
			{
				errmax = fabs(ytemp[i] / yscal[i]);
			}
		}
		errmax = errmax / eps;
		if (errmax > one)
		{
			h = safety * h * pow(errmax , pshrnk);
			flag = 1;
		}
		else
		{
			hdid = h;
			if (errmax > errcon)
			{
				hnext = safety * h * pow(errmax , pgrow);
			}
			else
			{
				hnext = 4.0 * h;
			}
			flag = 0;
		}
	}while(flag == 1);
	for (i = 1; i<=n; i++)
	{
		y[i] = y[i] + ytemp[i] * fcor;
	}
	//erase dysav, ysav, ytemp
}

void odeint(double ystart[], int nvar, double x1, double x2, double eps,
			double h1, double hmin, int& nok, int& nbad, double xp[],
			double yp[], int ypn)
{
	int kmax=0;
    int maxstp = 10000;
    double two = 2.0;
    double zero = 0.0;
    double tiny = 1e-30;
	double xsav=0;
	double dxsav=0;
    double yscal[11], y[11], dydx[11];
    double x = x1;
    double h = fabs(h1) * sgn(x2 - x1);
    nok = 0;
    nbad = 0;
    int i,kount = 0;
    for (i = 1; i<=nvar; i++)
	{
        y[i] = ystart[i];
    }
    if (kmax > 0)
	{
		xsav = x - dxsav * two;
	}
	int nstp;
    for (nstp = 1; nstp<=maxstp; nstp++)
	{
        derivs(x, y, dydx);
        for (i = 1; i<=nvar; i++)
		{
            yscal[i] = fabs(y[i]) + fabs(h * dydx[i]) + tiny;
        }
        if (kmax > 0)
		{
            if (fabs(x - xsav) > fabs(dxsav))
			{
                if (kount < kmax - 1)
				{
                    kount = kount + 1;
                    xp[kount] = x;
                    for (i = 1; i<=nvar; i++)
					{
                        yp[(i-1)*ypn+kount] = y[i];
                    }
                    xsav = x;
                }
            }
        }
        if ((x + h - x2) * (x + h - x1) > zero)
		{
			h = x2 - x;
		}
		double hnext,hdid;
        rkqc(y, dydx, nvar, x, h, eps, yscal, hdid, hnext);
        if (hdid == h)
		{
            nok = nok + 1;
		}
        else
		{
            nbad = nbad + 1;
        }
        if ((x - x2) * (x2 - x1) >= zero)
		{
            for (i = 1; i<=nvar; i++)
			{
                ystart[i] = y[i];
            }
            if (kmax != 0)
			{
                kount = kount + 1;
                xp[kount] = x;
                for (i = 1; i<=nvar; i++)
				{
                    yp[(i-1)*ypn+kount] = y[i];
                }
            }
            //erase dydx, y, yscal
            return;
        }
        if (fabs(hnext) < hmin)
		{
            cout<<" stepsize smaller than minimum."<<endl;
            return;
        }
        h = hnext;
    }
    cout<<" too many steps."<<endl;
}

void shoot(int nvar, double v[], double delv[], int n2, double x1, double x2,
		   double eps, double h1, double hmin, double f[], double dv[])
{
    double y[21], dfdv[401],xp[201], yp[2001];
	int ypn=100;
	int indx[21]; 
    load(x1, v, y);
	int i,iv,nok,nbad;
    odeint(y, nvar, x1, x2, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(x2, y, f);
	double sav;
    for (iv = 1; iv<=n2; iv++)
	{
        sav = v[iv];
        v[iv] = v[iv] + delv[iv];
        load(x1, v, y);
        odeint(y, nvar, x1, x2, eps, h1, hmin, nok, nbad, xp, yp,ypn);
        score(x2, y, dv);
        for (i = 1; i<=n2; i++)
		{
            dfdv[(i-1)*20+iv] = (dv[i] - f[i]) / delv[iv];
        }
        v[iv] = sav;
    }
    for (iv = 1; iv<=n2; iv++)
	{
        dv[iv] = -f[iv];
    }
	int det;
    ludcmp(dfdv, n2, 20,indx, det);
    lubksb(dfdv, n2, 20,indx, dv);
    for (iv = 1; iv<=n2; iv++)
	{
        v[iv] = v[iv] + dv[iv];
    }
}

void main()
{
    //program d15r1
    //driver for routine shoot
    //Solves for eigenvalues of spheroidal harmonics. Both
    //Prolate and Oblate case are handled simultaneously,
    //leading to six first-order equations. Unknown to
    //SHOOT, these are actually two independent sets of
    //three coupled equations, one set with c^2 positive
    //and the other with c^2 negative.
    int nvar = 6; int n2 = 2; double delta = 0.001; double eps = 0.000001;
    double v[3], delv[3], f[3], dv[3];
    double dx = 0.0001;
    do
	{
		cout<<"input m,n,c-squared (999 to end)"<<endl;
		m = 2;
		n = 2;
		c2 = 0.1;
		if (c2 == 999)
		{
			exit(1);
		}
	}while ((n < m) || (m < 0) || (n < 0));
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
	cout.width(1);
    cout<<m;
	cout.width(6);
	cout<<n;
	cout.width(12);
	cout<<c2<<endl;
    factr = 1;
	int q1;
    if (m != 0)
	{
        q1 = n;
        for (int i = 1; i<=m; i++)
		{
            factr = -0.5 * factr * (n + i) * ((double)q1 / (double)i);
            q1 = q1 - 1;
        }
    }
    v[1] = n * (n + 1) - m * (m + 1) + c2 / 2.0;
    v[2] = n * (n + 1) - m * (m + 1) - c2 / 2.0;
    delv[1] = delta * v[1];
    delv[2] = delv[1];
    double h1 = 0.1;
    double hmin = 0.0;
    double x1 = -1.0 + dx;
    double x2 = 0.0;
	cout<<endl;
    cout<<"            prolate                   oblate"<<endl;
    cout<<"      mu(m,n)    error est.    mu(m,n)      error est."<<endl;
    do
	{
		shoot(nvar, v, delv, n2, x1, x2, eps, h1, hmin, f, dv);
		cout.width(13);
		cout<<v[1];
		cout.width(13);
		cout<<dv[1];
		cout.width(13);	
		cout<<v[2];
		cout.width(13);
		cout<<dv[2]<<endl;
    }while (fabs(dv[1]) > fabs(eps * v[1]) || fabs(dv[2]) > fabs(eps * v[2]));
}