#define VERN "spectr2.cpp"

// spectr2.cpp   Jul 11  allow Poisson ratio as input
// spectr1.cpp  Jun 30  enumerate all frequencies
//              and write to a data file
// - beta2b.cpp  June 30   slight changes


#include <iostream.h>
#include <conio.h>
#include <graphics.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

double j(int n,double x);
double factorial(double arg);
double SeriesSphBess1j(int n,double x);
double SphBess1j(int n,double x);
double dubfactorial(double arg);
double evalan(int n,double alpha,double beta);
double evalbn(int n,double alpha,double beta);
double evalcn(int n,double alpha,double beta);
double evaldn(int n,double alpha,double beta);

void main(void)
{ // main
clrscr();
int ipx,ipy,i1;
double beta,f1,f2,sfx,alpha;
gotoxy(1,1);cout<<VERN;
gotoxy(1,2);cout<<"Calculation of spheroidal and toroidal frequencies of an elastic sphere";
gotoxy(1,3);cout<<"based on correction of equation (24) in Z. Ye's Chin. J. Phys 2000 paper.";
double nuPoisson,alphabybeta;

gotoxy(2,12);cout<<"nu must be from -1 to 0.5     Enter -2 to exit";
gotoxy(2,10);cout<<"What Poisson ratio? ";cin>>nuPoisson;
if (nuPoisson<=-2) {closegraph();exit(0);}
gotoxy(2,12);cout<<"                                              ";
alphabybeta=sqrt((1-2.0*nuPoisson)/(2-2.0*nuPoisson));
//gotoxy(2,10);cout<<"What is alpha/beta? ";cin>>alphabybeta;
if (alphabybeta<=0) return;
double db,func,func0,oldf1,oldf2; int countroots=0,nn;
db=0.001;
countroots=0;
char foutname[40];
gotoxy(2,12);
cout<<"Name of output file? (.txt extension; q to quit) ";cin>>foutname;
if (foutname[0]=='q') {closegraph();exit(0);}
FILE *fout;
fout=fopen(foutname,"wb");
fprintf(fout,"%s%s%s%c%c","// ",foutname," July 11, 2002   ",13,10);
fprintf(fout,"%s%f%s%f%s%c%c","// Poisson=",nuPoisson,"  alpha/beta=",alphabybeta,"   ",13,10);
//--------------------------------------------------
nn=0;
//  n=0 spheroidal
func0=0; oldf1=0;oldf2=0;
for (beta=db;beta<12;beta=beta+db)
{ // beta
alpha=alphabybeta*beta;
f1=tan(alpha)/alpha;
f2=1.0/(1.0-0.25*beta*beta);
func=f1-f2;
if ((func*func0<0)&&(fabs(f1-oldf1)<1)&&(fabs(f2-oldf2)<1))
{ // found a root
countroots++;
gotoxy(40+countroots/25,1+(countroots%25));cout<<" sph n=0 beta="<<beta<<"  ";
fprintf(fout,"%s%d%c%f%s%c%c","sph ",nn,' ',beta,"   ",13,10);
} // found a root
func0=func;
oldf1=f1;
oldf2=f2;
} // beta
//-----------------------------------------
//  n>=1  spheroidal
for (nn=1;nn<13;nn++)
{ // nn
func0=0;
for (beta=db;beta<12;beta=beta+db)
{ // beta
double ra,rb,rc,rd;
alpha=beta*alphabybeta;
ra=evalan(nn,alpha,beta);
rb=evalbn(nn,alpha,beta);
rc=evalcn(nn,alpha,beta);
rd=evaldn(nn,alpha,beta);
func=ra*rd-rb*rc;
if (func*func0<0)
{ // found a root
countroots++;
gotoxy(40+countroots/25,1+(countroots%25));cout<<" sph "<<nn<<" beta="<<beta<<"  ";
fprintf(fout,"%s%d%c%f%s%c%c","sph ",nn,' ',beta,"   ",13,10);
} // found a root
func0=func;
} // beta
} // nn
//---------------------------------
// toroidal   (n>=1 only;   n=0 are not physical modes)
for (nn=1;nn<13;nn++)
{ // n
gotoxy(67,2);cout<<"tor n="<<nn<<" ";
func0=0;
for (beta=db;beta<12;beta=beta+db)
{ //
double ra,rb,rc,rd,nd,det;  nd=nn;
det=(nd-1.0)*j(nn,beta)-beta*j(nn+1,beta);
if (det*func0<0)
{ // found a root
countroots++;
gotoxy(40+countroots/25,1+(countroots%25));cout<<" tor "<<nn<<" beta="<<beta<<" ";
fprintf(fout,"%s%d%c%f%s%c%c","tor ",nn,' ',beta,"   ",13,10);
} // found a root
func0=det;
} //
} // nn
//--------------------------------------

fprintf(fout,"%s%c%c","--------------------------",13,10);
fclose(fout);
gotoxy(1,25);cout<<"Done. Press any key to exit. ";
getch();
return;
} // main---------------------------------------------------------

double j(int n,double x)
{ // spherical bessel functions of first kind
// Using Ye's normalization as given in 2000 Chin. J. Phys. paper
double res,nd;nd=n;
res=SphBess1j(n,x);
// Fixup to try to match Ye's normalization
res=res/factorial(nd);
res=res/pow(2.0,nd);
res=res*factorial(2.0*nd+1.0);
res=res/dubfactorial(2.0*nd+1.0);
return res;
} // spherical bessel functions of first kind

double dubfactorial(double arg)
{ // double factor  n!! = n(n-2)(n-4)...
double ans,arg2,r2; int i1,i2;
ans=1.0;
i1 = floor(arg+0.5);
for (i2=i1;i2>1;i2=i2-2)
{ //
r2=i2;
ans=ans*r2;
} //
return ans;
} // double factor  n!! = n(n-2)(n-4)...

double factorial(double arg)
{ //
double ans,arg2,r2; int i1,i2;
ans=1.0;
i1 = floor(arg+0.5);
for (i2=2;i2<=i1;i2++)
{ //
r2=i2;
ans=ans*r2;
} //
return ans;
} //

double SeriesSphBess1j(int n,double x)
{ // infinite series spherical bessel function of 1st kind
// Note: summing to 15 gives optimal precision out to
//  about x=15, after which the series blow up, if "double"
//  is used.  If double precision is used, then it is possible
//  to sum more terms to extend the range of accurate convergence
//  without getting the blow up problem.
// (Double precision as used here is not really necessary.)
if (x>14) {cout<<"x ="<<x<<" too large for series sph bess func.";exit(0);}
double ans;
double ansd=0,nr,sr,xd;
xd=x;
nr=n;
int s;
for (s=0;s<15;s++)
{ //
sr=s;
ansd=ansd+pow(-1.0,sr)*factorial(sr+nr)*pow(xd,2.0*sr)/(factorial(sr)*factorial(2.0*s+2.0*n+1));
} //
ansd = ansd * pow(2.0,nr)*pow(xd,nr);
ans=ansd;
return ans;
} //

double SphBess1j(int n,double x)
{ //  SphBess1j
if (n<0) {cout<<"Negative n for Sph Bess 1";exit(0);}
if (n==0) return sin(x)/x;
if (n==1) return (sin(x)/(x*x))-(cos(x)/x);
if (n==2) return ( (3.0/(x*x*x)) - (1.0/x) )*sin(x) - (3.0/(x*x))*cos(x);
return SeriesSphBess1j(n,x); // use series expression (valid x<14)
} //  SphBess1j

double evalan(int n,double alpha,double beta)
{ //  evalan
double svan,nd; nd=n;
svan = pow(-1.0,nd)*
(
beta*beta*j(n,alpha)-2*(nd-1.0)*alpha*j(n-1,alpha)
);
//l//'/// (2.0*nd+1.0);
return svan;
} //  evalan

double evalbn(int n,double alpha,double beta)
{ //  evalan
double svbn,nd; nd=n;
svbn=pow(-1.0,nd+1.0)*
(beta*beta*j(n,alpha)-2.0*(nd+2.0)*j(n+1,alpha)*alpha);
///(2.0*nd+1.0);
return svbn;
} //  evalbn

double evalcn(int n,double alpha,double beta)
{ //  evalcn
double svcn,nd;nd=n;
svcn=(beta*beta*j(n,beta)-2.0*(nd-1.0)*beta*j(n-1,beta));
return svcn;
} //  evalan

double evaldn(int n,double alpha,double beta)
{ //  evalan
double svdn,nd; nd=n;
//svdn=pow((-1.0/beta),nd)*(nd/(nd+1.0))*
svdn=(nd/(nd+1.0))*
(
beta*beta*j(n,beta)-2.0*(nd+2.0)*beta*j(n+1,beta)
);
//alpha=kl*ObjRad;
//beta=ks*ObjRad;
return svdn;
} //  evalan
// Incompressibility approximation: cl infinity, so alpha=0.