#define VERN "sphmode4.cpp"

// sphmode4.cpp  Jul 8   Corrected TA normalization
//   Take old values of TA and multiply by (0.53301)^2 to get new.
// sphmode2 Jul 8 fixed error in determ. ratio
//  Jul 8 minor change
// sphmode.cpp  Jul 7
// spheroi.cpp  Jul 7 look for non-initialization error
//                  ( I didn't find any )
//  4
//  3  Jul 6  - add spheroidal modes
//   2  Jul 6 - generalized to include toroidal modes
//     - (Big puzzle that XZ component of stress does not vanish!)
// hollsph1.cpp  Jul 5   start to work on modes of hollow sphere
//  - SphBess() now 1st kind or 2nd kind
//  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 <graphics.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <dos.h>

void makescr(void);

double LegenP(int n,double x);
double SphBess(int kind,int n,double x);
double SeriesSphBess(int kind,int n,double x);
double rPix(int n,double x,double y,double z);
double rPiy(int n,double x,double y,double z);
double rPiz(int n,double x,double y,double z);
double CrPix(int n,double x,double y,double z);
double CrPiy(int n,double x,double y,double z);
double CrPiz(int n,double x,double y,double z);
double C2rPix(int n,double x,double y,double z);
double C2rPiy(int n,double x,double y,double z);
double C2rPiz(int n,double x,double y,double z);
double THETAx(int n,double x,double y,double z);
double THETAy(int n,double x,double y,double z);
double THETAz(int n,double x,double y,double z);
double factorial(double arg);
void Eval_u(double x1,double y1,double z1);
double angle_delta,uX,uY,uZ;
double divm[3][3];
double toremat[3][3];
double sphemat[3][3];
double emat[3][3];

int modetype,besselkind=-1,nn;
double Pi;

// Properties of material:
//double lambdaLame1=0.58385;
double lambdaLame1;
double ShearModG=0.38415;
//double nuPoisson=0.301575;
double nuPoisson=0.301575;
double Cs,Cl,ks,kl;
double dx=0.0001;
double alpha,beta,ObjRad=1.0; // ObjRad=radius of sphere

void main(void)
{ // main
clrscr();
gotoxy(1,1);cout<<"Find spheroidal modes:";
gotoxy(2,4);cout<<"(0 = 0.301575 fcc value) ";
gotoxy(2,2);cout<<"nu? ";cin>>nuPoisson; // Poisson ratio="<<nuPoisson<<" ";
if (nuPoisson==0) nuPoisson=0.301575;
gotoxy(2,4);cout<<"                         ";
lambdaLame1=ShearModG/((1.0/(2.0*nuPoisson))-1.0);
Cs=sqrt(ShearModG); // transverse wave speed
Cl=sqrt(lambdaLame1+2.0*ShearModG); // longit. wave speed
//nuPoisson=0.5*(1.0-2.0*(Cs/Cl)*(Cs/Cl))/(1.0-(Cs/Cl)*(Cs/Cl));
gotoxy(55,2);cout<<"lambda="<<lambdaLame1<<" ";
gotoxy(55,3);cout<<"Cs="<<Cs<<" ";
gotoxy(55,4);cout<<"Cl="<<Cl<<" ";
gotoxy(55,5);cout<<"Poisson="<<nuPoisson<<" ";
Pi = 4.0*atan(1.0);
cout.precision(5);
//int gdriver=VGA; int gmode=VGAHI;
//initgraph(&gdriver,&gmode,"c:\\tc\\bgi"); cleardevice();
//struct palettetype pal; getpalette(&pal);
int ic;
double rinbyR;
//gotoxy(1,1);cout<<VERN<<" Find solutions for the wet dog mode for hollow sphere:";
//gotoxy(40,2);cout<<"(inner to outer radius ratio)";
//gotoxy(1,2);cout<<"Wbat is rin/R? ";cin>>rinbyR;
//if (rinbyR<=0) return;
//gotoxy(40,2);cout<<"                             ";
//gotoxy(1,5);cout<<"(1=toroidal 2=spheroidal 3=wet dog)";
//gotoxy(1,3);cout<<"Mode type? ";cin>>modetype;
//if (modetype<=0) return;
//gotoxy(1,5);cout<<"                                   ";
gotoxy(10,3);cout<<"What is the order n? ";cin>>nn;
if (nn<0) return;

double x1,y1,z1;
double ax2,ax1,ay2,ay1,az2,az1;
double bx2,bx1,by2,by1,bz2,bz1;
double cx2,cx1,cy2,cy1,cz2,cz1;
double vala1,vala2,valb1,valb2;
double valXZa1,valXZa2,valXZb1,valXZb2;
double valXXa1,valXXa2,valXXb1,valXXb2;
int i1,i2,countroots;
double sphstressxx;
double torstressxx;
double divu,deter15,deter150,deter,TAfactor15,dbeta;

LB1:
deter150=0;
countroots=0;
gotoxy(1,4);cout<<"beta?                        (enter x=-1 to quit)        ";
gotoxy(1,4);cout<<"beta? ";cin>>beta;
gotoxy(1,4);cout<<"beta = "<<beta<<"                                          ";
if (beta<0) return;
gotoxy(20,4);cout<<"db? ";cin>>dbeta;
if (dbeta<0) return;
gotoxy(5,5);cout<<"( was "<<beta<<" )   ";
x1=1.0;
y1=0;
z1=0;
//gotoxy(1,4);cout<<"1: ("<<x1<<" "<<y1<<" "<<z1<<")  ";

double exx0;
//besselkind=1;
//Eval_u(x1,y1,z1);
//gotoxy(1,10);cout<<"wd (1) u= ("<<uX<<", "<<uY<<", "<<uZ<<")       ";
int iad=-1;

gotoxy(1,11);cout<<"Toroidal (zero divergence part)";
gotoxy(1,17);cout<<"Spheroidal (zero curl part)";
LB2:
beta=beta+dbeta;
gotoxy(65,1);cout<<beta<<" ";
alpha=(Cs/Cl)*beta; // 0.53301 = Cs/Cl = 0.61980 / 1.16282
gotoxy(48,1);cout<<"al "<<alpha<<" ";
// ks r= beta r/ObjRad; ?
ks=beta/ObjRad;
kl=alpha/ObjRad;

iad=-1;
for (angle_delta=0.0;angle_delta<0.7854;angle_delta+=0.2618)
{ //
iad++;
cout.precision(2);
gotoxy(5+17*iad,9);cout<<180.0*angle_delta/Pi<<"ø  ";
cout.precision(5);
//--------------------------------------------------------
modetype=1; // "toroidal (i.e. zero divergence part)"

besselkind=1;
Eval_u(x1+0.5*dx,y1,z1);ax2=uX; ay2=uY; az2=uZ;
Eval_u(x1-0.5*dx,y1,z1);ax1=uX; ay1=uY; az1=uZ;
Eval_u(x1,y1+0.5*dx,z1);bx2=uX; by2=uY; bz2=uZ;
Eval_u(x1,y1-0.5*dx,z1);bx1=uX; by1=uY; bz1=uZ;
Eval_u(x1,y1,z1+0.5*dx);cx2=uX; cy2=uY; cz2=uZ;
Eval_u(x1,y1,z1-0.5*dx);cx1=uX; cy1=uY; cz1=uZ;
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;
for (i1=0;i1<3;i1++) for (i2=0;i2<3;i2++)
toremat[i1][i2]=0.5*(divm[i1][i2]+divm[i2][i1]);
// (Divergence is zero for "toroidal" part)
divu=toremat[0][0]+toremat[1][1]+toremat[2][2];
torstressxx=lambdaLame1*divu+2.0*ShearModG*toremat[0][0];
gotoxy(1,12);cout<<"eXX";
gotoxy(1,13);cout<<"eXY";
gotoxy(1,14);cout<<"eXZ";
gotoxy(1,15);cout<<"sXX";
gotoxy(5+17*iad,12);cout<<toremat[0][0]<<"    ";
gotoxy(5+17*iad,13);cout<<toremat[0][1]<<"    ";
gotoxy(5+17*iad,14);cout<<toremat[0][2]<<"    ";
gotoxy(5+17*iad,15);cout<<torstressxx<<"   ";
//----------------------------------------------
modetype=2; // "spheroidal" i.e. zero curl part

besselkind=1;
Eval_u(1.0+0.5*dx,y1,z1);ax2=uX; ay2=uY; az2=uZ;
Eval_u(1.0-0.5*dx,y1,z1);ax1=uX; ay1=uY; az1=uZ;
Eval_u(1.0,y1+0.5*dx,z1);bx2=uX; by2=uY; bz2=uZ;
Eval_u(1.0,y1-0.5*dx,z1);bx1=uX; by1=uY; bz1=uZ;
Eval_u(1.0,y1,z1+0.5*dx);cx2=uX; cy2=uY; cz2=uZ;
Eval_u(1.0,y1,z1-0.5*dx);cx1=uX; cy1=uY; cz1=uZ;
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;
for (i1=0;i1<3;i1++) for (i2=0;i2<3;i2++)
sphemat[i1][i2]=0.5*(divm[i1][i2]+divm[i2][i1]);
divu=sphemat[0][0]+sphemat[1][1]+sphemat[2][2];
sphstressxx=lambdaLame1*divu+2.0*ShearModG*sphemat[0][0];
gotoxy(1,18);cout<<"eXX";
gotoxy(1,19);cout<<"eXY";
gotoxy(1,20);cout<<"eXZ";
gotoxy(1,21);cout<<"div";
gotoxy(1,22);cout<<"sXX";
gotoxy(1,25);cout<<"TA";
gotoxy(5+17*iad,18);cout<<sphemat[0][0]<<"    ";
gotoxy(5+17*iad,19);cout<<sphemat[0][1]<<"    ";
gotoxy(5+17*iad,20);cout<<sphemat[0][2]<<"    ";
gotoxy(5+17*iad,21);cout<<divu<<"   ";
gotoxy(5+17*iad,22);cout<<sphstressxx<<"   ";
// Found error in following line July 8:   (fixed here)
deter=sphstressxx*toremat[0][2] - torstressxx*sphemat[0][2];
if (iad==2) deter15=deter;
gotoxy(5+17*iad,24);cout<<deter<<"    ";
gotoxy(5+17*iad,25); cout<<"             "; gotoxy(5+17*iad,25);
double TAfactor=999;
if (torstressxx!=0) TAfactor=sphstressxx/torstressxx;
if (torstressxx!=0) {cout<<TAfactor;}
if (iad==2) TAfactor15=TAfactor;

} // angle_delta

if ((deter15*deter150<0)&&(countroots<20))
{ //
countroots++;
gotoxy(60,1+countroots);cout<<" "<<countroots<<": "<<beta<<" "<<TAfactor15<<"  ";
gotoxy(60,2+countroots);cout<<"------------------";
} //
deter150=deter15;
if (dbeta<=0) goto LB1;
if (!kbhit()) goto LB2;
getch();
goto LB1;



gotoxy(1,25);cout<<"Press any key to exit..";
getch();
return;



} // 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 SphBess(int kind,int n,double x)
{ //  SphBess   Spherical Bessel functions
if ((kind<1)||(kind>2))
{cout<<"** kind="<<kind<<" Press a key.";getch();exit(0);}
if (kind==1)
{ // first kind
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);
} // first kind
if (kind==2)
{ // 2nd kind
if (n==0) return -cos(x)/x;
if (n==1) return -((cos(x)/(x*x))+(sin(x)/x));
if (n==2) return -(
                   (
                    (3.0/(x*x*x))-(1.0/x)
                   )*cos(x)
                 + (3.0/(x*x))*sin(x)
                  );
} // 2nd kind
return SeriesSphBess(kind,n,x); // use series expression (valid x<14)
} //  SphBess

double SeriesSphBess(int kind,int n,double x)
{ // infinite series spherical bessel function of 1st kind
if (kind==2) {cout<<"* 2nd kind not available. Press any key.";getch();exit(0);}
if (kind!=1) {cout<<"* No such kind. Press any key."<<kind;getch();exit(0);}
// 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.)
// Caution: This is common normalization and not the normalization
// used in Z. Ye 2000.
if (x>14) {cout<<"* x ="<<x<<" too large for series sph bess func. Press any key.";getch();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 rPix(int n,double x,double y,double z)
{ //  rPi x-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
costheta=z/r;
return x*SphBess(besselkind,n,ks*r)*LegenP(n,costheta);
} //
else return 0;
} //  r Pi x-component

double rPiy(int n,double x,double y,double z)
{ //  r Pi y-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
costheta=z/r;
return y*SphBess(besselkind,n,ks*r)*LegenP(n,costheta);
} //
else return 0;
} //  r Pi y-component

double rPiz(int n,double x,double y,double z)
{ //  r Pi z-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //.
costheta=z/r;
return z*SphBess(besselkind,n,ks*r)*LegenP(n,costheta);
} //.
else return 0;
} //  r Picz-component

double CrPix(int n,double x,double y,double z)
{ //  CurlrPix
double ans;
ans = (rPiz(n,x,y+0.5*dx,z)-rPiz(n,x,y-0.5*dx,z))/dx;
ans-= (rPiy(n,x,y,z+0.5*dx)-rPiy(n,x,y,z-0.5*dx))/dx;
return ans;
} //  CurlrPix

double CrPiy(int n,double x,double y,double z)
{ //  CurlrPiy
double x1,x2,y1,y2,z1,z2,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans = (rPix(n,x,y,z2)-rPix(n,x,y,z1))/dx;
ans -=  (rPiz(n,x2,y,z)-rPiz(n,x1,y,z))/dx;
return ans;
} //  CurlrPiy

double CrPiz(int n,double x,double y,double z)
{ //  CrPiz
double x1,x2,y1,y2,z1,z2,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans = (rPiy(n,x2,y,z)-rPiy(n,x1,y,z))/dx;
ans -=  (rPix(n,x,y2,z)-rPix(n,x,y1,z))/dx;
return ans;
} //  CrPiz

//--------------------

double C2rPix(int n,double x,double y,double z)
{ //  curl of curl of rPi - x component
double x1,x2,y1,y2,z1,z2,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans = (CrPiz(n,x,y2,z)-CrPiz(n,x,y1,z))/dx;
ans -=  (CrPiy(n,x,y,z2)-CrPiy(n,x,y,z1))/dx;
return ans;
} //  curl of curl of rPi - x component

double C2rPiy(int n,double x,double y,double z)
{ //  CurlrPiy
double x1,x2,y1,y2,z1,z2,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (CrPix(n,x,y,z2)-CrPix(n,x,y,z1))/dx;
ans -= (CrPiz(n,x2,y,z)-CrPiz(n,x1,y,z))/dx;
return ans;
} //  CurlrPiy

double C2rPiz(int n,double x,double y,double z)
{ //   C2rPiz
double ans;
ans =  (CrPiy(n,x+0.5*dx,y,z)-CrPiy(n,x-0.5*dx,y,z))/dx;
ans -= (CrPix(n,x,y+0.5*dx,z)-CrPix(n,x,y-0.5*dx,z))/dx;
return ans;
} //   C2rPiz

double factorial(double arg)
{ //   factorial
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;
} //   factorial

double THETAx(int n,double x,double y,double z)
{ // x derivative of THETA
double x1,x2,ans,THETA2=0,THETA1=0;
double r,costheta;
x2=x+0.5*dx; r=sqrt(x2*x2+y*y+z*z);
if (r!=0) { costheta=z/r; THETA2=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
x1=x-0.5*dx; r=sqrt(x1*x1+y*y+z*z);
if (r!=0) { costheta=z/r; THETA1=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // x derivative of THETA

double THETAy(int n,double x,double y,double z)
{ // y derivative of THETA
double y1,y2,ans,THETA2=0,THETA1=0;
double r,costheta;
y2=y+0.5*dx; r=sqrt(x*x+y2*y2+z*z);
if (r!=0) { costheta=z/r; THETA2=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
y1=y-0.5*dx; r=sqrt(x*x+y1*y1+z*z);
if (r!=0) { costheta=z/r; THETA1=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // y derivative of THETA

double THETAz(int n,double x,double y,double z)
{ // z derivative of THETA
double z1,z2,ans,THETA2=0,THETA1=0; double r,costheta;
z2=z+0.5*dx; r=sqrt(x*x+y*y+z2*z2);
if (r!=0) { costheta=z2/r; THETA2=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
z1=z-0.5*dx; r=sqrt(x*x+y*y+z1*z1);
if (r!=0) { costheta=z1/r; THETA1=SphBess(besselkind,n,kl*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // z derivative of THETA

void makescr(void)
{ // makescr
gotoxy(70,5);cout<<"Make .scr ";
gotoxy(70,6);cout<<"file? ";
gotoxy(70,7);cout<<"(y/n/q)  ";
char ch678;
ch678=getch();
if (ch678=='q') exit(0);
if (ch678!='y') return;
char scrfilename[40];
FILE *scrfout;
gotoxy(65,8);cout<<"Filename";
gotoxy(65,9);cout<<"? ";cin>>scrfilename;

scrfout=fopen(scrfilename,"wb");
unsigned char ch1b,ch2b,ch3b,ch4b,ch5,ch6;
int scrwid=640;
int ipybot=400;
int ipytop=300;
int scrhgt=ipybot-ipytop;
int scrx0=0;
int scry0=ipytop;
setcolor(15);
line(scrx0,scry0-1,scrwid+scrwid-1,scry0-1);
line(scrx0,scry0+scrhgt,scrwid+scrwid-1,scry0+scrhgt);
ch3b=scrwid/256;
ch1b=scrwid%256;
ch4b=scrhgt/256;
ch2b=scrhgt%256;
fprintf(scrfout,"%c%c%c%c%c%c%c%c%c%c%c%c%c",'s','c','r',ch3b,ch1b,ch4b,ch2b,'0','0','0','0',13,10);
int ipixcol,ix,iy;
gotoxy(65,10);cout<<"saving ";
for (iy=scry0;iy<scry0+scrhgt;iy++)
{ // iy
gotoxy(65,11);cout<<iy;
for (ix=scrx0;ix<scrx0+scrwid;ix++)
{ //
ipixcol=getpixel(ix,iy);
//putpixel(ix-100,iy-100,ipixcol);
char ch1;
ch1=ipixcol+'A';
fprintf(scrfout,"%c",ch1);
} //
fprintf(scrfout,"%c%c",13,10);
} // iy
fclose(scrfout);
return;
} //

void Eval_u(double X,double Y,double Z)
{ // Eval_u (evaluate displacement field (uX,uY,uZ) at point (XYZ) )
double x2,y2,z2,theta2,r2,phi2;
double wdx2,wdy2,wdz2,wdr2,wdtheta2,wdphi2;

// Rotate from (X,Y,Z) to (x2,y2,z2) coordinates
x2= cos(angle_delta)*X-sin(angle_delta)*Z;
y2= Y;
z2= sin(angle_delta)*X+cos(angle_delta)*Z;

//gotoxy(1,5);cout<<"2: ("<<x2<<" "<<y2<<" "<<z2<<")  ";
// Convert from (x2 y2 z2) to (r2 theta2 phi2)

//gotoxy(1,6);cout<<"r = "<<r2<<"   theta = "<<180.0*theta2/Pi<<" deg  phi= "<<180.0*phi2/Pi<<" deg   ";

if (modetype==1)
{ // 1 toroidal mode
wdx2=C2rPix(nn,x2,y2,z2);
wdy2=C2rPiy(nn,x2,y2,z2);
wdz2=C2rPiz(nn,x2,y2,z2);
} // 1
if (modetype==2)
{ // 2 spheroidal
wdx2=THETAx(nn,x2,y2,z2);
wdy2=THETAy(nn,x2,y2,z2);
wdz2=THETAz(nn,x2,y2,z2);
} // 2 spheroidal
if (modetype==3)
{ // 3 wet dog mode
double LP2,LP1;
r2=sqrt(x2*x2+y2*y2+z2*z2);
if (z2==0)
{ //.
theta2=Pi*0.5;
} //.
else
{ //
if (z2>0)
 theta2=atan(sqrt(x2*x2+y2*y2)/z2);
else
 theta2=Pi-atan(sqrt(x2*x2+y2*y2)/fabs(z2));
} //

phi2=0;
if (x2!=0)
{ //
if (x2>0)
 phi2=atan(y2/x2);
else
 phi2=Pi+atan(y2/x2);
} //
else
{ //.
if (y2>0)
 phi2 =  Pi*0.5;
else
 phi2 = -Pi*0.5;
}//.
LP2=LegenP(nn,cos(theta2+0.5*dx));
LP1=LegenP(nn,cos(theta2-0.5*dx));
wdphi2=SphBess(besselkind,nn,r2)*(LP2-LP1)/dx;
wdtheta2=0;
wdr2=0;
wdx2=wdr2*sin(theta2)*cos(phi2)+(-wdphi2*sin(phi2))+cos(theta2)*cos(phi2)*wdtheta2;
wdy2=wdr2*sin(theta2)*sin(phi2)+( wdphi2*cos(phi2))+cos(theta2)*sin(phi2)*wdtheta2;
wdz2=wdr2*cos(theta2)-sin(theta2)*wdtheta2;
} // 3 wet dog mode
//gotoxy(1,7);cout<<"u phi = "<<wdphi2<<" ";
//gotoxy(1,8);cout<<"wd (2) u = "<<wdx2<<" "<<wdy2<<", "<<wdz2<<")      ";
// Rotate back from (x2,y2,z2) system to (X,Y,Z) coordinates.
uX = cos(angle_delta)*wdx2+sin(angle_delta)*wdz2;
uY = wdy2;
uZ =-sin(angle_delta)*wdx2+cos(angle_delta)*wdz2;
//gotoxy(40,10);cout<<uX<<" "<<uY<<" "<<uZ<<"  ";
return;
} // Eval_u

/*
void func1(double r1,double *r2)
{ //
(*r2) = r1+2.0;
return;
} //
*/