//  pert12  increase Ro and Rm
//  pert11.cpp  Jan 17   be careful of k r < 12 !
//  pert10.cpp  Jan 16  much faster
//  pert9   general associated Legendre functions
//  pert8.cpp  use weighted Monte Carlo Jan 16
//  pert7c.cpp  Jan 16   try to avoid large fluctuations
//  pert7   Jan 16   allow l>2
//  pert6.cpp  Jan 15
//  pert4.cpp  Jan 14  orbm  |p l m>
//  Note that |1 2 2> different from |1 2 -2>.
//  pert3b Jan 14  tried larger dx and solved roundoff problem.
//  pert3
//  pert2  Jan 14   evaluate random matrix elements <p1 l1 | p2 l2>
//    functions are apparently orthonormalized
// pert1.cpp  Jan 13, 2002   perturbation approach to anisotropicity
//   - calculation of matrix element <i|j>
//  Functions were taken from c:\dice\d4\curldiv.cpp

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

double LegenP(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 ksguess,klguess;
double THETAx(int n,int m,double x,double y,double z);
double THETAy(int n,int m,double x,double y,double z);
double THETAz(int n,int m,double x,double y,double z);
double factorial(double arg);
double u(int iaxis,int ipol,int orderL,int m,double x,double y,double z);
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 C11,C12,C44;
double densityrho=2.33; // Si silicon

void main(void)
{ // main
randomize();
clrscr();
C11=166.0; // Si at 300 K   //www.ioffe.rssi.ru
C12=64.0;
C44=80.0;
// C44 = 0.5*(C11-C12)  if isotropic
// 2*C44 = C11- C12     if isotropic
// C12 = C11 - 2 * C44; if isotropic
// C12 = C11 - 2 * C44; // make it isotropic keeping [100] speeds constant!
double Clong,Ctrans,omega;
int iClong,iCtrans;
iClong= 1000.0*sqrt(C11/densityrho);
iCtrans=1000.0*sqrt(C44/densityrho);
cout.precision(1);
gotoxy(1,25);cout<<"CL[100]= "<<iClong<<" m/s  CT[100]= "<<iCtrans<<" m/s ";
ksguess=0.9853;
klguess=0.9851;
omega = ksguess*Ctrans;
//gotoxy(40,24);cout<<"omegaż = "<<omega*omega<<" if transv ";
omega = klguess*Clong;
//gotoxy(40,25);cout<<"omegaż = "<<omega*omega<<" if longit ";

double a1byk2=0;
gotoxy(1,1);cout<<"Starting...";
int ipol1,ipol2;
int orderL1,orderL2,m1,m2;
double x,y,z;
long int countit;
double f1,f2,sum1,sum2,sum3,Ro=40,Rm=300,f2biggest;
int countrowmin=2;
int countrow=countrowmin-1,Lmax=6;
LB0:
countrow++;
if (countrow>25) countrow=countrowmin;
ipol1=random(3)+1; ipol1=3;
ipol2=random(3)+1; //ipol2=ipol1;
ipol2=3;
orderL1=random(Lmax+1); if (ipol1>1) orderL1=1+random(Lmax);
orderL1=2;
orderL2=random(Lmax+1); if (ipol2>1) orderL2=1+random(Lmax);
//orderL2=10;
orderL2=orderL2+(orderL2%2); // make L2 even
m1=random(2*orderL1+1)-orderL1; m1=0;
m2=random(2*orderL2+1)-orderL2; //m2=m1;
m2=m2+(m2%2); // make m2 even
//m2=m1;
//if (random(3)==0) m2=m1-4;
//if (random(3)==0) m2=m1+4;
if (m2>orderL2) m2=m1;
if (m2<-orderL2) m2=m1;
//if (m2==m1) m2=random(2*orderL2+1)-orderL2;
//if (m2==m1) m2=random(2*orderL2+1)-orderL2;
//if (m2==m1) m2=random(2*orderL2+1)-orderL2;
//ipol2=1; orderL2=4; m2=4;
//ipol2=ipol1; orderL2=orderL1; m2=m1;
//ipol1=1; orderL1=6; m1=0;
//ipol2=1, orderL2=6; m2=0;
gotoxy(20,countrow);cout<<"< "<<ipol1<<" k "<<orderL1<<" "<<m1<<" | H' | "<<ipol2<<" k "<<orderL2<<" "<<m2<<" >  ";
countit=0;
sum1=0;
sum2=0;
sum3=0;
f2biggest=0;
long int countitmax=10000;
gotoxy(1,6);cout<<countitmax<<" iter.";
LB1:
countit++;
if (countit>countitmax)
{//.
if (fabs(a1byk2)<1.00) countrow--;
goto LB0;
}//.
if ((countit%100)==0)
{ //,,
cout.precision(2);
gotoxy(1,4);cout<<countit<<" ";
double a1;
//a1=(2.0/3.1415926)*sum1*Rm*Rm/(countit);
//if (fabs(a1)<0.002) a1=0;
//gotoxy(52,countrow);cout<<a1<<"   ";
a1=sum1/sqrt(sum2*sum3);
if (fabs(a1)<0.002) a1=0;
if (ipol1==1) {gotoxy(50,countrow);a1byk2=a1/(klguess*klguess);cout<<a1/(klguess*klguess)<<" kż   ";}
if (ipol1==2) {gotoxy(50,countrow);a1byk2=a1/(ksguess*ksguess);cout<<a1/(ksguess*ksguess)<<" kż   ";}
if (ipol1==3) {gotoxy(50,countrow);a1byk2=a1/(ksguess*ksguess);cout<<a1/(ksguess*ksguess)<<" kż   ";}
int ispeed;
if ((ipol1==ipol2)&&(orderL1==orderL2)&&(m1==m2))
{ //...
if (ipol1==1) {gotoxy(65,countrow);ispeed=1000.0*sqrt(a1)/klguess;cout<<ispeed<<" m/s ";}
if (ipol1==2) {gotoxy(65,countrow);ispeed=1000.0*sqrt(a1)/ksguess;cout<<ispeed<<" m/s ";}
if (ipol1==3) {gotoxy(65,countrow);ispeed=1000.0*sqrt(a1)/ksguess;cout<<ispeed<<" m/s ";}
} //...
} //,,
int X=1,Y=2,Z=3,SPHL=1,TORT=2,SPHT=3,iax;
double x5,y5,z5,r5,r6;
LB2:
x5=0.0001*(random(20000)-10000);
y5=0.0001*(random(20000)-10000);
z5=0.0001*(random(20000)-10000);
r5=sqrt(x5*x5+y5*y5+z5*z5);
if (r5<0.001*0.001) goto LB2;
x5=x5/r5;
y5=y5/r5;
z5=z5/r5;
r6=Ro+(Rm-Ro)*0.0001*random(10000);
x=x5*r6;
y=y5*r6;
z=z5*r6;
//if (x*x+y*y+z*z>Rm*Rm) goto LB2;
//if (x*x+y*y+z*z<12.1*12.1) goto LB2;
//        sigma i j = C i j k l  e k l
// f iax =  sigma iax j,j = C iax j k l {(d/dxj)  e k l}
for (iax=1;iax<=3;iax++)
{ // iax
ksguess=1.2363;klguess=1.2361;
f1=u(iax,ipol1,orderL1,m1,x,y,z);
ksguess=1.2363;klguess=1.2361;
int ij,ik,il;
f2=0.0;
for (ij=1;ij<=3;ij++)
{ // ij
for (ik=1;ik<=3;ik++)
{ // ik
for (il=ik;il<=3;il++)
{ // il
if (C(iax,ij,ik,il)!=0)
{ //..
double dc1;
dc1=(-1.0/densityrho)*C(iax,ij,ik,il)*ekl_j(ik,il,ij,ipol2,orderL2,m2,x,y,z);
if (ik!=il) dc1=dc1*2.0;
f2+=dc1;
} //..
} // il
} // ik
} // ij
//if (fabs(f2)>f2biggest) {f2biggest=f2;gotoxy(1,countrow);cout<<f2biggest<<" "<<x<<" "<<y<<" "<<z<<" ";}
sum1=sum1+r6*r6*f1*f2;
f2=u(iax,ipol2,orderL2,m2,x,y,z);
sum2=sum2+r6*r6*f1*f1;
sum3=sum3+r6*r6*f2*f2;
} // iax
if (!kbhit()) goto LB1;


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

double u(int iaxis,int ipol,int orderL,int m,double x,double y,double z)
{ // u
double rf1;
if (iaxis<1) {exit(0);}
if (iaxis>3) {exit(0);}
if (ipol<1) {exit(0);}
if (ipol>3) {exit(0);}
if (ipol==1)
{ // 1 spheroidal
rf1=sqrt(2*orderL+1.0);
if (iaxis==1) return rf1*THETAx(orderL,m,x,y,z);
if (iaxis==2) return rf1*THETAy(orderL,m,x,y,z);
if (iaxis==3) return rf1*THETAz(orderL,m,x,y,z);
} // 1 spheroidal longit.
if (ipol==2)
{ // 2 toroidal (trans.)
rf1=sqrt((2*orderL+1.0)/(orderL*(orderL+1.0)))*ksguess;
if (iaxis==1) return rf1*CrPix(orderL,m,x,y,z);
if (iaxis==2) return rf1*CrPiy(orderL,m,x,y,z);
if (iaxis==3) return rf1*CrPiz(orderL,m,x,y,z);
} // 2 toroidal (trans.)
if (ipol==3)
{ // 3 spheroidal transv.
rf1=sqrt((2*orderL+1.0)/(orderL*(orderL+1.0)));
if (iaxis==1) return rf1*C2rPix(orderL,m,x,y,z);
if (iaxis==2) return rf1*C2rPiy(orderL,m,x,y,z);
if (iaxis==3) return rf1*C2rPiz(orderL,m,x,y,z);
} // 3 spheroidal transv.
exit(0);
} // u

double LegenP(int l,int m,double x)
{ //  LegenP - Legendre polynomials
// 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;
} //  LegenP
//-----------------------------
/*
double LegenP(int n,int m,double x)
{ //  LegenP - 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);
} //  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);
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,ksguess*r)*LegenP(n,m,costheta);
if (m>0)  return x*SphBess1j(n,ksguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  return x*SphBess1j(n,ksguess*r)*LegenP(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,ksguess*r)*LegenP(n,m,costheta);
if (m>0)  return y*SphBess1j(n,ksguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  return y*SphBess1j(n,ksguess*r)*LegenP(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,ksguess*r)*LegenP(n,m,costheta);
if (m>0)  return z*SphBess1j(n,ksguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  return z*SphBess1j(n,ksguess*r)*LegenP(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
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,x,y,z2)-rPiy(n,m,x,y,z1))/dx;
ans -= (rPiz(n,m,x,y2,z)-rPiz(n,m,x,y1,z))/dx;
return ans;
} //  CurlrPix

double CrPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
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,x2,y,z)-rPiz(n,m,x1,y,z))/dx;
ans -= (rPix(n,m,x,y,z2)-rPix(n,m,x,y,z1))/dx;
return ans;
} //  CurlrPiy

double CrPiz(int n,int m,double x,double y,double z)
{ //  CrPiz
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,y2,z)-rPix(n,m,x,y1,z))/dx;
ans -= (rPiy(n,m,x2,y,z)-rPiy(n,m,x1,y,z))/dx;
return ans;
} //  CrPiz

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

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

double C2rPiy(int n,int m,double x,double y,double z)
{ //  CurlrPiy
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,x2,y,z)-CrPiz(n,m,x1,y,z))/dx;
ans -= (CrPix(n,m,x,y,z2)-CrPix(n,m,x,y,z1))/dx;
return ans;
} //  CurlrPiy

double C2rPiz(int n,int m,double x,double y,double z)
{ //  C2rPiz
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,y2,z)-CrPix(n,m,x,y1,z))/dx;
ans -= (CrPiy(n,m,x2,y,z)-CrPiy(n,m,x1,y,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 THETAx(int n,int m,double x,double y,double z)
{ // x derivative of THETA
double x1,x2,dx=0.001,ans,THETA2=0,THETA1=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) THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA2=SphBess1j(n,klguess*r)*LegenP(n,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) THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*sin(m*phi);
}
ans=(THETA2-THETA1)/dx;
return ans;
} // x derivative of THETA

double THETAy(int n,int m,double x,double y,double z)
{ // y derivative of THETA
double y1,y2,dx=0.001,ans,THETA2=0,THETA1=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) THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA2=SphBess1j(n,klguess*r)*LegenP(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) THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*sin(m*phi);
}
ans=(THETA2-THETA1)/dx;
return ans;
} // y derivative of THETA

double THETAz(int n,int m,double x,double y,double z)
{ // z derivative of THETA
double z1,z2,dx=0.001,ans,THETA2=0,THETA1=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) THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA2=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA2=SphBess1j(n,klguess*r)*LegenP(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) THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta);
if (m>0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*cos(m*phi);
if (m<0)  THETA1=SphBess1j(n,klguess*r)*LegenP(n,m,costheta)*sin(m*phi);
}
ans=(THETA2-THETA1)/dx;
return ans;
} // z derivative of THETA

double C(int i,int j,int k,int l)
{ // C
int n1=-1,n2=-1;
if ((i==1)&&(j==1)) n1=1;
if ((i==2)&&(j==2)) n1=2;
if ((i==3)&&(j==3)) n1=3;
if ((i==2)&&(j==3)) n1=4;
if ((i==3)&&(j==2)) n1=4;
if ((i==1)&&(j==3)) n1=5;
if ((i==3)&&(j==1)) n1=5;
if ((i==1)&&(j==2)) n1=6;
if ((i==2)&&(j==1)) n1=6;
if ((k==1)&&(l==1)) n2=1;
if ((k==2)&&(l==2)) n2=2;
if ((k==3)&&(l==3)) n2=3;
if ((k==2)&&(l==3)) n2=4;
if ((k==3)&&(l==2)) n2=4;
if ((k==1)&&(l==3)) n2=5;
if ((k==3)&&(l==1)) n2=5;
if ((k==1)&&(l==2)) n2=6;
if ((k==2)&&(l==1)) n2=6;
if (n2==-1) exit(0);
if ((n1==1)&&(n2==1)) return C11;
if ((n1==2)&&(n2==2)) return C11;
if ((n1==3)&&(n2==3)) return C11;
if ((n1==1)&&(n2==2)) return C12;
if ((n1==2)&&(n2==1)) return C12;
if ((n1==1)&&(n2==3)) return C12;
if ((n1==3)&&(n2==1)) return C12;
if ((n1==2)&&(n2==3)) return C12;
if ((n1==3)&&(n2==2)) return C12;
if ((n1==4)&&(n2==4)) return C44;
if ((n1==5)&&(n2==5)) return C44;
if ((n1==6)&&(n2==6)) return C44;
return 0.0;
} // C

//double u(int iaxis,int ipol,int orderL,double x,double y,double z);

double ekl(int k,int l,int ipol,int orderL,int m,double x,double y,double z)
{ // strain tensor (Landau & Lifshitz def'n)
double dukdxl,dx=0.001,duldxk;
if (l==1) dukdxl=(u(k,ipol,orderL,m,x+dx,y,z)-u(k,ipol,orderL,m,x-dx,y,z))/(2*dx);
if (l==2) dukdxl=(u(k,ipol,orderL,m,x,y+dx,z)-u(k,ipol,orderL,m,x,y-dx,z))/(2*dx);
if (l==3) dukdxl=(u(k,ipol,orderL,m,x,y,z+dx)-u(k,ipol,orderL,m,x,y,z-dx))/(2*dx);
duldxk=dukdxl;
if (k!=l)
{ //
if (k==1) duldxk=(u(l,ipol,orderL,m,x+dx,y,z)-u(l,ipol,orderL,m,x-dx,y,z))/(2*dx);
if (k==2) duldxk=(u(l,ipol,orderL,m,x,y+dx,z)-u(l,ipol,orderL,m,x,y-dx,z))/(2*dx);
if (k==3) duldxk=(u(l,ipol,orderL,m,x,y,z+dx)-u(l,ipol,orderL,m,x,y,z-dx))/(2*dx);
} //
return 0.5*(dukdxl+duldxk);
} //

double ekl_j(int k,int l,int j,int ipol,int orderL,int m,double x,double y,double z)
{ //  ekl_j
if (j<1) {cout<<"ekl_j error 1 *";exit(0);}
if (j>3) {cout<<"ekl_j error 2 *";exit(0);}
if (ipol<1) {cout<<"ekl_j error 3 *";exit(0);}
if (ipol>3) {cout<<"ekl_j error 4 *";exit(0);}
double a1=-12345.0,dx=0.001;
if (j==1) a1=(ekl(k,l,ipol,orderL,m,x+dx,y,z)-ekl(k,l,ipol,orderL,m,x-dx,y,z))/(2*dx);
if (j==2) a1=(ekl(k,l,ipol,orderL,m,x,y+dx,z)-ekl(k,l,ipol,orderL,m,x,y-dx,z))/(2*dx);
if (j==3) a1=(ekl(k,l,ipol,orderL,m,x,y,z+dx)-ekl(k,l,ipol,orderL,m,x,y,z-dx))/(2*dx);
if (a1!=-12345.0) return a1;
cout<<"ekl_j error";exit(0);
} //  ekl_j