#define VERN "fluxp0.cpp"

// WARNING:  Curls have sign error

// fluxp0.cpp  Apr 4 04  flux due to p-wave
//  plane4b.cpp  Mar 30   s-wave
// plane4   Mar 27
// plane3.cpp  Mar 27   transverse plane wave --> t-wave partial wave
// plane2.cpp  Mar 27
// plane1.cpp  Mar 27, 2004
//  pert14  Jan 19 2003 diagonal element table with speeds
//  pert13  Jan 18  make html table
//  Functions were taken from c:\dice\d4\curldiv.cpp

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

double LegenPabsm(int orderL,int orbm,double x);
complex SphBessh1(int orderL,double x);  // Hankel 1st kind (outgoing)
complex rPix(int orderL,int m,double x,double y,double z);
complex rPiy(int orderL,int m,double x,double y,double z);
complex rPiz(int orderL,int m,double x,double y,double z);
complex CrPix(int orderL,int m,double x,double y,double z);
complex CrPiy(int orderL,int m,double x,double y,double z);
complex CrPiz(int orderL,int m,double x,double y,double z);
complex C2rPix(int orderL,int m,double x,double y,double z);
complex C2rPiy(int orderL,int m,double x,double y,double z);
complex C2rPiz(int orderL,int m,double x,double y,double z);

complex dPhidx(int n,int m,double x,double y,double z);
complex dPhidy(int n,int m,double x,double y,double z);
complex dPhidz(int n,int m,double x,double y,double z);
double factorial(double arg);

double C(int i,int j,int k,int l);
double ekl(int k,int l,int ipol,int orderL,double x,double y,double z);
double ekl_j(int k,int l,int j,int ipol,int orderL,int m,double x,double y,double z);

double kL,kT;
double vLm,omega,rhom; // frequency

void main(void)
{ // main

clrscr();
gotoxy(1,1);cout<<" "<<VERN<<"  This find total power P";
gotoxy(1,2);cout<<"             into p-wave scattered wave";

double x1,y1;
int L1,m1;
gotoxy(3,3);cout<<"L ? ";cin>>L1;
if (L1<0) return;
m1=0;
gotoxy(3,3);cout<<"L = "<<L1<<" ";
gotoxy(3,4);cout<<"m = "<<m1<<" ";

double Rmin,Rmax,x2,y2,z2,r22,r2,theta2,xy2,phi2,x3,y3,z3;
double sum1=0.0,count1=0.0,sum2,sum3,sum2r,sum2i;
complex csum2;
rhom=1.5;
kL=1.2; kT = 0.7;
vLm=1.6; // longitudinal speed of sound
omega = vLm * kL;
Rmax = 15.7/kT;
Rmin = 15.1/kT;
gotoxy(50,3);cout<<"Rmin = "<<Rmin<<"  Rmax = "<<Rmax<<"  ";


complex cf2x,cf2y,cf2z;

double f1,f2,f3;
gotoxy(50,4);cout<<"kL = "<<kL<<"  ";
gotoxy(50,5);cout<<"kT = "<<kT<<"  ";

sum1=0.0;count1=0.0;sum2=0.0;sum3=0.0;
csum2=0.0;
sum2r=0.0;
sum2i=0.0;

randomize();

long int li1=0; // iteration counter
LB1:
li1++;
if ((li1%100)==0) {gotoxy(5,5);cout<<li1<<" ";}
LB2:
x2 = Rmax * (2.0*0.001*random(1000)-1.0);
y2 = Rmax * (2.0*0.001*random(1000)-1.0);
z2 = Rmax * (2.0*0.001*random(1000)-1.0);
r22=x2*x2+y2*y2+z2*z2;
if (r22>Rmax*Rmax) goto LB2;
if (r22<Rmin*Rmin) goto LB2;
if (r22<0.1) goto LB2;
r2 = sqrt(r22);
xy2 = sqrt(x2*x2+y2*y2);
theta2 = (0.5*3.1415926)-atan(z2/xy2);
if (x2==0.0)
{ // x2 zero
if (y2>0) phi2 =  0.5*3.1415926;
if (y2<0) phi2 = -0.5*3.1415926;
} // x2 zero
else
{ //..
phi2 = atan(y2/x2);
if (x2<0) phi2=phi2+3.1415926;
} //..

z3 = r2 * cos(theta2);
x3 = r2 * sin(theta2) * cos(phi2);
y3 = r2 * sin(theta2) * sin(phi2);

if (((li1%100)==0)&&(count1>2.0))
{ // 1000

cout.precision(3);
//gotoxy(40,4);cout<<"x2 = "<<x2<<"  "<<x3<<"        ";
//gotoxy(40,5);cout<<"y2 = "<<y2<<"  "<<y3<<"        ";
//gotoxy(40,6);cout<<"z2 = "<<z2<<"  "<<z3<<"        ";
//gotoxy(40,7);cout<<"diff = "<<(x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)<<"        ";
gotoxy(8,8);cout<<"avg of P = "<<sum1/count1<<" watts";
float Pavg =sum1/count1;
gotoxy(8,10);cout<<"P*(kL*kL)*(2L+1)/(2*pi*omega*omega*vLm*rhom) = "<<Pavg*(kL*kL)*(2*L1+1)/(2.0*3.1415926*omega*omega*vLm*rhom)<<"  ";
//gotoxy(8,9);cout<<"<1|2>="<<csum2/count1<<" ";

//complex cnumeransw;
//cnumeransw =csum2/sum1;
// gotoxy(8,10);cout<<"<1|2>/<1|1> = "<<cnumeransw<<"     ";

//gotoxy(8,13);cout<<"Exact answer is supposed to be:";
//gotoxy(8,14);cout<<"H(L) = i^(L+1) (2L+1) / (L(L+1))";
//complex cexactansw;
//cexactansw=pow(complex(0.0,1.0),L1+1)*(2*L1+1)/(L1*(L1+1));
//gotoxy(8,15);cout<<cexactansw<<"  ";

//gotoxy(8,16);cout<<"% err: "<<100*abs(cexactansw-cnumeransw)/abs(cnumeransw)<<" % ";

} // 1000

complex cf1x,cf1y,cf1z;


double vx,vy,vz,S,P;
// &&&&&

// Find displacement field (complex)

cf1x = dPhidx(L1,m1,x2,y2,z2)/kL;
cf1y = dPhidy(L1,m1,x2,y2,z2)/kL;
cf1z = dPhidz(L1,m1,x2,y2,z2)/kL;

vx = 0.7071*omega*abs(cf1x); // rms speed along x
vy = 0.7071*omega*abs(cf1y); // rms speed along y
vz = 0.7071*omega*abs(cf1z); // rms speed along z

S = 0.5 * rhom * (vx*vx + vy*vy + vz*vz) * vLm;

S = S * 2.0; // total enrgy is double of kinetic energy alone

P = S * 4.0 * 3.1415926 * ((Rmin+Rmax)*0.5)* ((Rmin+Rmax)*0.5);

sum1 = sum1 + P;
//sum2r = sum2r + real(f1x*cf2x+f1y*cf2y+f1z*cf2z);
//sum2i = sum2i + imag(f1x*cf2x+f1y*cf2y+f1z*cf2z);
//csum2 = complex(sum2r,sum2i);

count1=count1+1.0;

if (!kbhit()) goto LB1;


gotoxy(1,25);cout<<"All done...";
getch();
return;
} // main

double LegenPabsm(int l,int m,double x)
{ //  LegenPabsm - Legendre polynomials
// WARNING: Not valid for negative m
// Info on Legendre and associated Legendre polynomials at:
// http://mathworld.wolfram.com/LegendrePolynomial.html
if (l>4) goto Lpast4;
if (m<0) m=-m;
if (m==0)
{ // m=0
if (l==0) return 1;
if (l==1) return x;
if (l==2) return 0.5*(3.0*x*x-1.0);
if (l==3) return 0.5*(5.0*x*x*x-3.0*x);
if (l==4) return 0.125*(35.0*x*x*x*x-30.0*x*x+3);
if (l==5) return 0.125*(63.0*x*x*x*x*x-70.0*x*x*x+15*x);
} // m=0
if (m==1)
{ // m=1
if (l==1) return -sqrt(1-x*x);
if (l==2) return -3*x*sqrt(1-x*x);
if (l==3) return 1.5*(1.0-5*x*x)*sqrt(1-x*x);
if (l==4) return (2.5)*x*(3-7*x*x)*sqrt(1-x*x);
} // m=1
if (m==2)
{ // m=2
if (l==2) return 3*(1-x*x);
if (l==3) return 15*x*(1-x*x);
if (l==4) return (7.5)*(7*x*x-1)*(1-x*x);
} //
if (m==3)
{ // m=3
if (l==3) return -15*(1-x*x)*sqrt(1-x*x);
if (l==4) return -105*x*(1-x*x)*sqrt(1-x*x);
} // m=3
if (m==4)
{ // m=4
if (l==4) return 105*(1-x*x)*(1-x*x);
} // m=4

Lpast4:
if (m<0) m=-m;
double total=0.0,factor;
int k,kmax;
kmax=l/2;
for (k=0;k<=kmax;k++)
{ // k
factor=0.0;
if (l-2*k-m>=0)
{ //.
factor=1.0;
factor=factor*pow(-1.0,k);
factor=factor*factorial(2*l-2*k);
factor=factor/factorial(k);
factor=factor/factorial(l-k);
factor=factor/factorial(l-(2*k+m));
// Must handle exception pow(0,0)...
if ((x!=0)||(l!=(2*k+m))) factor=factor*pow(x,l-(2*k+m));
} //.
total=total+factor;
} // k
total=total/pow(2.0,l);
total=total*pow(-1.0,m);
total=total*pow(fabs(1.0-x*x),m/2.0);
return total;
} //  LegenPabsm

complex SphBessh1(int n,double x)
{ //    SphBessh1
complex eix = exp(complex(0.0,x)), c1 = complex(1.0,0), ci = complex(0,1.0);
if (n==0) return -ci*eix/x;
if (n==1) return -eix*( (c1/x) + (ci/(x*x)) );
if (n==2) return  eix*( (ci/x) + (-3*c1/(x*x)) + (-3*ci/(x*x*x)) );
cout<<"* SphBessh1: Sorry. Cannot do it for l>2. 1234 *";
exit(0);
} //  SphBess1j

complex rPix(int n,int m,double x,double y,double z)
{ //  rPi x-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926;
if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) return x*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return x*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return x*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi x-component

complex rPiy(int n,int m,double x,double y,double z)
{ //  r Pi y-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) return y*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return y*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return y*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi y-component

complex rPiz(int n,int m,double x,double y,double z)
{ //  r Pi z-component
double r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) return z*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return z*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return z*SphBessh1(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi z-component

complex CrPix(int n,int m,double x,double y,double z)
{ //  CurlrPix
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (rPiy(n,m,x,y,z2)-rPiy(n,m,x,y,z1))/dx;
cans -= (rPiz(n,m,x,y2,z)-rPiz(n,m,x,y1,z))/dx;
return cans;
} //  CurlrPix

complex CrPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (rPiz(n,m,x2,y,z)-rPiz(n,m,x1,y,z))/dx;
cans -= (rPix(n,m,x,y,z2)-rPix(n,m,x,y,z1))/dx;
return cans;
} //  CurlrPiy

complex CrPiz(int n,int m,double x,double y,double z)
{ //  CrPiz
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (rPix(n,m,x,y2,z)-rPix(n,m,x,y1,z))/dx;
cans -= (rPiy(n,m,x2,y,z)-rPiy(n,m,x1,y,z))/dx;
return cans;
} //  CrPiz

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

complex  C2rPix(int n,int m,double x,double y,double z)
{ //  curl of curl of rPi - x component
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (CrPiy(n,m,x,y,z2)-CrPiy(n,m,x,y,z1))/dx;
cans -= (CrPiz(n,m,x,y2,z)-CrPiz(n,m,x,y1,z))/dx;
return cans;
} //  curl of curl of rPi - x component

complex C2rPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (CrPiz(n,m,x2,y,z)-CrPiz(n,m,x1,y,z))/dx;
cans -= (CrPix(n,m,x,y,z2)-CrPix(n,m,x,y,z1))/dx;
return cans;
} //  CurlrPiy

complex  C2rPiz(int n,int m,double x,double y,double z)
{ //  C2rPiz
// NOTE: Sign error, here and below.
double x1,x2,y1,y2,z1,z2,dx=0.02,ans;
complex cans;
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;
cans =  (CrPix(n,m,x,y2,z)-CrPix(n,m,x,y1,z))/dx;
cans -= (CrPiy(n,m,x2,y,z)-CrPiy(n,m,x1,y,z))/dx;
return cans;
} //  C2rPiz

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;
} //

complex dPhidx(int orderL,int m,double x,double y,double z)
{ // x derivative of Phi
double x1,x2,dx=0.001,ans; complex dPhid2=0,dPhid1=0;
double r,costheta;
complex cans;
x2=x+0.5*dx; r=sqrt(x2*x2+y*y+z*z);
if (r!=0)
{ //
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) dPhid2=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta);
if (m>0)  dPhid2=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*sin(m*phi);
} //
x1=x-0.5*dx; r=sqrt(x1*x1+y*y+z*z);
if (r!=0)
{
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) dPhid1=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta);
if (m>0)  dPhid1=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBessh1(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*sin(m*phi);
}
cans=(dPhid2-dPhid1)/dx;
return cans;
} // x derivative of Phi

complex dPhidy(int n,int m,double x,double y,double z)
{ // y derivative of Phi
double y1,y2,dx=0.001,ans; complex dPhid2=0,dPhid1=0;
double r,costheta;
complex cans;
y2=y+0.5*dx; r=sqrt(x*x+y2*y2+z*z);
if (r!=0)
{
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
y1=y-0.5*dx; r=sqrt(x*x+y1*y1+z*z);
if (r!=0)
{
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926; if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z/r;
if (m==0) dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
cans=(dPhid2-dPhid1)/dx;
return cans;
} // y derivative of Phi

complex dPhidz(int n,int m,double x,double y,double z)
{ // z derivative of Phi
double z1,z2,dx=0.001; complex dPhid2=0,dPhid1=0;
double r,costheta;
complex cans;
z2=z+0.5*dx; r=sqrt(x*x+y*y+z2*z2);
if (r!=0)
{
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926;
if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
 costheta=z2/r;
if (m==0) dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
z1=z-0.5*dx; r=sqrt(x*x+y*y+z1*z1);
if (r!=0)
{
double phi=0.5*3.1415926; if (y<0) phi=-0.5*3.1415926;
if (x!=0) phi=atan(y/x); if (x<0) phi=phi+3.1415926;
costheta=z1/r;
if (m==0) dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBessh1(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
cans=(dPhid2-dPhid1)/dx;
return cans;
} // z derivative of Phi