#define VERN "ptrantwi.cpp"

// ptrantwi.cpp   Apr 7 04  curl sign error fixed

// plane3b2   Apr 7 04  noticed sign error in curl calculations,
//    meaning that results using this program are wrong.

//  plane3b.cpp  Mar 30 04
//    verifying I(L) expression, [eqTP2d] Eq. (47) in cross paper.
// 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);
double SphBess1j(int orderL,double x);
double SeriesSphBess1j(int orderL,double x);
double rPix(int orderL,int m,double x,double y,double z);
double rPiy(int orderL,int m,double x,double y,double z);
double rPiz(int orderL,int m,double x,double y,double z);
double CrPix(int orderL,int m,double x,double y,double z);
double CrPiy(int orderL,int m,double x,double y,double z);
double CrPiz(int orderL,int m,double x,double y,double z);
double C2rPix(int orderL,int m,double x,double y,double z);
double C2rPiy(int orderL,int m,double x,double y,double z);
double C2rPiz(int orderL,int m,double x,double y,double z);

double dPhidx(int n,int m,double x,double y,double z);
double dPhidy(int n,int m,double x,double y,double z);
double 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;

void main(void)
{ // main

clrscr();
gotoxy(1,1);cout<<" "<<VERN<<"  'Plane Transverse T-wave I'";
gotoxy(1,2);cout<<"               This verifies decomposition of x-polarized transverse";
gotoxy(1,3);cout<<"               plane wave into t-wave partial spherical waves";
gotoxy(1,4);cout<<"               This verifies the expression for Il in cross paper.";

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

double 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;
kL=1.1; kT = 1.5;
Rmax = 11.7/kT;
gotoxy(50,3);cout<<"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;

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<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%1000)==0)&&(count1>2.0))
{ //
//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<<"<1|1>="<<sum1/count1<<" ";
//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<<"I(L) = i^(L) (2L+1) / (L(L+1))";
complex cexactansw;
cexactansw=pow(complex(0.0,1.0),L1)*(2*L1+1)/(L1*(L1+1));
gotoxy(8,15);cout<<cexactansw<<"  ";

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


} //

//cf2x=complex(0,1.0)*exp(complex(0.0,kT*r2*cos(theta2)));  old
cf2x=exp(complex(0.0,kT*r2*cos(theta2))); // new
cf2y=0;
cf2z=0;

double f1x,f1y,f1z;

f1x = CrPix(L1,m1,x2,y2,z2);
f1y = CrPiy(L1,m1,x2,y2,z2);
f1z = CrPiz(L1,m1,x2,y2,z2);
//f1x = dPhidx(L1,m1,x2,y2,z2)/kL;
//f1y = dPhidy(L1,m1,x2,y2,z2)/kL;
//f1z = dPhidz(L1,m1,x2,y2,z2)/kL;

sum1 = sum1  +  f1x* f1x+f1y* f1y+f1z* f1z;
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;

randomize();

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

double LegenPabsm(int l,int m,double x)
{ //  LegenPabsm - Legendre polynomials
// Info on Legendre and associated Legendre polynomials at:
// http://mathworld.wolfram.com/LegendrePolynomial.html
// Warning: Not correct for m<0
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

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

/*
double LegenPabsm(int n,int m,double x)
{ //  LegenPabsm - Legendre polynomials
// Info on Legendre and associated Legendre polynomials at:
// http://mathworld.wolfram.com/LegendrePolynomial.html
if (m<0) m=-m;
if (m==0)
{ // m=0
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);
} // m=0
if (m==1)
{ // m=1
if (n==1) return -sqrt(1-x*x);
if (n==2) return -3*x*sqrt(1-x*x);
if (n==3) return 1.5*(1.0-5*x*x)*sqrt(1-x*x);
if (n==4) return (2.5)*x*(3-7*x*x)*sqrt(1-x*x);
} // m=1
if (m==2)
{ // m=2
if (n==2) return 3*(1-x*x);
if (n==3) return 15*x*(1-x*x);
if (n==4) return (7.5)*(7*x*x-1)*(1-x*x);
} //
if (m==3)
{ // m=3
if (n==3) return -15*(1-x*x)*sqrt(1-x*x);
if (n==4) return -105*x*(1-x*x)*sqrt(1-x*x);
} // m=3
if (m==4)
{ // m=4
if (n==4) return 105*(1-x*x)*(1-x*x);
} // m=4
cout<<"* Legendre polynomial order (l="<<n<<", m="<<m<<") not available *";exit(0);
} //  LegenPabsm
*/

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

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);
if ((n>=3)&&(x>12))
{ // use large x limit (see md30.htm) (non-singular expressions, wrong for small x)
if ((n%4)==0) return (pow(x,n)/(1+pow(x,n)))*sin(x)/x;
if ((n%4)==1) return (pow(x,n)/(1+pow(x,n)))*((sin(x)/(x*x))-(cos(x)/x));
if ((n%4)==2) return -(pow(x,n)/(1+pow(x,n)))*sin(x)/x;
if ((n%4)==3) return (pow(x,n)/(1+pow(x,n)))*((cos(x)/x)-(sin(x)/(x*x)));
} // use large x limit
return SeriesSphBess1j(n,x); // use series expression (valid x<14)
} //  SphBess1j

double 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*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return x*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return x*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi x-component

double 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*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return y*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return y*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi y-component

double 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*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta);
if (m>0)  return z*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  return z*SphBess1j(n,kT*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
} //
else return 0;
} //  r Pi z-component

double CrPix(int n,int m,double x,double y,double z)
{ //  CurlrPix
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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 =  (rPiz(n,m,x,y2,z)-rPiz(n,m,x,y1,z))/dx;
ans -= (rPiy(n,m,x,y,z2)-rPiy(n,m,x,y,z1))/dx;
return ans;
} //  CurlrPix

double CrPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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,m,x,y,z2)-rPix(n,m,x,y,z1))/dx;
ans -= (rPiz(n,m,x2,y,z)-rPiz(n,m,x1,y,z))/dx;
return ans;
} //  CurlrPiy

double CrPiz(int n,int m,double x,double y,double z)
{ //  CrPiz
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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,m,x2,y,z)-rPiy(n,m,x1,y,z))/dx;
ans -= (rPix(n,m,x,y2,z)-rPix(n,m,x,y1,z))/dx;
return ans;
} //  CrPiz

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

double C2rPix(int n,int m,double x,double y,double z)
{ //  curl of curl of rPi - x component
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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,m,x,y2,z)-CrPiz(n,m,x,y1,z))/dx;
ans -= (CrPiy(n,m,x,y,z2)-CrPiy(n,m,x,y,z1))/dx;
return ans;
} //  curl of curl of rPi - x component

double C2rPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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,m,x,y,z2)-CrPix(n,m,x,y,z1))/dx;
ans -= (CrPiz(n,m,x2,y,z)-CrPiz(n,m,x1,y,z))/dx;
return ans;
} //  CurlrPiy

double C2rPiz(int n,int m,double x,double y,double z)
{ //  C2rPiz
// April 7 04: curl sign error fixed
double x1,x2,y1,y2,z1,z2,dx=0.02,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 =  (CrPiy(n,m,x2,y,z)-CrPiy(n,m,x1,y,z))/dx;
ans -= (CrPix(n,m,x,y2,z)-CrPix(n,m,x,y1,z))/dx;
return ans;
} //  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;
} //

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 dPhidx(int orderL,int m,double x,double y,double z)
{ // x derivative of dPhid
double x1,x2,dx=0.001,ans,dPhid2=0,dPhid1=0;
double r,costheta;
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=SphBess1j(orderL,kL*r)*LegenPabsm(orderL,m,costheta);
if (m>0)  dPhid2=SphBess1j(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBess1j(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=SphBess1j(orderL,kL*r)*LegenPabsm(orderL,m,costheta);
if (m>0)  dPhid1=SphBess1j(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBess1j(orderL,kL*r)*LegenPabsm(orderL,m,costheta)*sin(m*phi);
}
ans=(dPhid2-dPhid1)/dx;
return ans;
} // x derivative of dPhid

double dPhidy(int n,int m,double x,double y,double z)
{ // y derivative of dPhid
double y1,y2,dx=0.001,ans,dPhid2=0,dPhid1=0;
double r,costheta;
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=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid2=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBess1j(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=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid1=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
ans=(dPhid2-dPhid1)/dx;
return ans;
} // y derivative of dPhid

double dPhidz(int n,int m,double x,double y,double z)
{ // z derivative of dPhid
double z1,z2,dx=0.001,ans,dPhid2=0,dPhid1=0;
double r,costheta;
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=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid2=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid2=SphBess1j(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=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta);
if (m>0)  dPhid1=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*cos(m*phi);
if (m<0)  dPhid1=SphBess1j(n,kL*r)*LegenPabsm(n,m,costheta)*sin(m*phi);
}
ans=(dPhid2-dPhid1)/dx;
return ans;
} // z derivative of dPhid