// wetdog4.cpp  July 5, 2002
// wetdog3 July 4  - generalize to order n
// wetdog2.cpp  July 3, 2002   Find beta for wet dog mode exactly,
//  find finding zeroes in the strain matrix.

#include <iostream.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <dos.h>

double LegenP(int n,double x);
double SphBess1j(int n,double x);
double SeriesSphBess1j(int n,double x);
double factorial(double arg);
void EvalWetDog_u(double x1,double y1,double z1);
double ang_delta,wdX,wdY,wdZ;
double divm[3][3];
double emat[3][3];
int nn; int showpartialresults;
double dtheta=0.001;
double dx=0.0003;
void main(void)
{ // main
clrscr();
cout.precision(5);
gotoxy(1,1);cout<<"wetdog4.cpp  Find exact solutions for the wet dog mode...";
gotoxy(2,4);cout<<" (The order can be 1, 2, 3, 4 or 5) ";
gotoxy(2,2);cout<<"What is the order n? ";cin>>nn;
gotoxy(2,4);cout<<"                                    ";
if (nn<0) return;
double X,y1,z1;
ang_delta=3.1415926*(39.0/180.0);
LB1:
gotoxy(1,3);cout<<"X?                        (enter X=-1 to quit)        ";
gotoxy(1,4);cout<<" (Enter a (positive) guess for frequency beta of the mode.)     ";
gotoxy(1,3);cout<<"X? ";cin>>X;
gotoxy(1,4);cout<<"                                                           ";
gotoxy(1,3);cout<<"X = "<<X<<"                                          ";
if (X<=0) return;
y1=0;
z1=0;

gotoxy(1,5);cout<<"Step 1: X="<<X<<" Y="<<y1<<" Z="<<z1<<"    ";

double ax2,ax1,ay2,ay1,az2,az1;
double bx2,bx1,by2,by1,bz2,bz1;
double cx2,cx1,cy2,cy1,cz2,cz1;
showpartialresults=1; EvalWetDog_u(X,y1,z1); showpartialresults=0;
EvalWetDog_u(X-0.5*dx,y1,z1); ax1=wdX; ay1=wdY; az1=wdZ;
EvalWetDog_u(X+0.5*dx,y1,z1); ax2=wdX; ay2=wdY; az2=wdZ;
EvalWetDog_u(X,y1-0.5*dx,z1); bx1=wdX; by1=wdY; bz1=wdZ;
EvalWetDog_u(X,y1+0.5*dx,z1); bx2=wdX; by2=wdY; bz2=wdZ;
EvalWetDog_u(X,y1,z1-0.5*dx); cx1=wdX; cy1=wdY; cz1=wdZ;
EvalWetDog_u(X,y1,z1+0.5*dx); cx2=wdX; cy2=wdY; cz2=wdZ;
divm[0][0]=(ax2-ax1)/dx;
divm[1][0]=(ay2-ay1)/dx;
divm[2][0]=(az2-az1)/dx;
divm[0][1]=(bx2-bx1)/dx;
divm[1][1]=(by2-by1)/dx;
divm[2][1]=(bz2-bz1)/dx;
divm[0][2]=(cx2-cx1)/dx;
divm[1][2]=(cy2-cy1)/dx;
divm[2][2]=(cz2-cz1)/dx;
int i1,i2;
gotoxy(40,11);cout<<"Derivative matrix:";
gotoxy(40,12);cout<<"------------------";
for (i1=0;i1<3;i1++)
for (i2=0;i2<3;i2++)
{ //
gotoxy(40+9*i1,13+i2);cout<<divm[i1][i2]<<"       ";
} //

// Compute and display strain matrix:
double biggeststraincomp=0;
gotoxy(10,17);cout<<"Strain Matrix: (xy element is zero for mode) ";
gotoxy(10,18);cout<<"-------------------------------------------- ";
for (i1=0;i1<3;i1++)
for (i2=0;i2<3;i2++)
{ //
emat[i1][i2]=0.5*(divm[i1][i2]+divm[i2][i1]);
if (fabs(emat[i1][i2])>biggeststraincomp) biggeststraincomp=fabs(emat[i1][i2]);
gotoxy(10+13*i1,19+i2);cout<<emat[i1][i2]<<"         ";
} //

goto LB1;
} // main

double LegenP(int n,double x)
{ //  LegenP - Legendre polynomials
double p;
if (n==0) return 1;
if (n==1) return x;
if (n==2) return 0.5*(3.0*x*x-1.0);
if (n==3) return 0.5*(5.0*x*x*x-3.0*x);
if (n==4) return 0.125*(35.0*x*x*x*x-30.0*x*x+3);
if (n==5) return 0.125*(63.0*x*x*x*x*x-70.0*x*x*x+15*x);
cout<<"* Legendre polynomial order not available *";exit(0);
} //  LegenP

double SphBess1j(int n,double x)
{ //  SphBess1j
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 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;
} //

void EvalWetDog_u(double X,double Y,double Z)
{ // EvalWetDog_u
double x2,y2,z2,theta2,r2,phi2,wdphi2,wdx2,wdy2,wdz2;

x2= cos(ang_delta)*X-sin(ang_delta)*Z;
y2= Y;
z2= sin(ang_delta)*X+cos(ang_delta)*Z;

if (showpartialresults==1)
{ //
gotoxy(1,6);cout<<"Step 2. Rotate by delta="<<180.0*ang_delta/3.1415926<<" degrees   ";
gotoxy(1,7);cout<<"x="<<x2<<" y="<<y2<<" z="<<z2<<"   ";
} //

r2=sqrt(x2*x2+y2*y2+z2*z2);
if (z2==0)
{ //
theta2=3.1415926*0.5;
} //
else
{ //
if (z2>0) theta2=atan(sqrt(x2*x2+y2*y2)/z2);
if (z2<0) theta2=3.1415926-atan(sqrt(x2*x2+y2*y2)/fabs(z2));
} //

phi2=0;
if (x2!=0)
{ //
if (x2>0) phi2=atan(y2/x2);
if (x2<0) phi2=3.1415926+atan(y2/x2);
} //
else
{ //
if (y2>0) phi2= 3.1415926*0.5;
if (y2<0) phi2=-3.1415926*0.5;
}//

if (showpartialresults==1)
{ // .
gotoxy(1,8);cout<<"Step 3: Convert to spherical coordinates ";
gotoxy(1,9);cout<<"r="<<r2<<" theta="<<180.0*theta2/3.1415926<<" deg  phi="<<180.0*phi2/3.1415926<<" deg   ";
} // .

double theta22,theta21;
double LP2,LP1;
LP2=LegenP(nn,cos(theta2+0.5*dtheta));
LP1=LegenP(nn,cos(theta2-0.5*dtheta));
wdphi2=SphBess1j(nn,r2)*(LP2-LP1)/dtheta;
if (showpartialresults==1)
{
gotoxy(1,10);cout<<"Step 4: u phi = "<<wdphi2<<"   ";
}
wdx2=-wdphi2*sin(phi2);
wdy2= wdphi2*cos(phi2);
wdz2=0;
if (showpartialresults==1)
{ //..
gotoxy(1,11);cout<<"Step 5: ux="<<wdx2<<" uy="<<wdy2<<" uz="<<wdz2<<"     ";
} //..

wdX = cos(ang_delta)*wdx2+sin(ang_delta)*wdz2;
wdY = wdy2;
wdZ =-sin(ang_delta)*wdx2+cos(ang_delta)*wdz2;

if (showpartialresults==1)
{ //..
gotoxy(1,12);cout<<"Step 6: uX="<<wdx2<<" uY="<<wdy2<<" uZ="<<wdz2<<"     ";
} //..

return;
} // EvalWetDog_u