int ExactEllipsoid=1;
#define VERN "varpr27e.cpp"
//  [ use with varpr28j.cpp ]
#define OutFileName "emp3d.clm"

// varpr27e Jun 29 - separate function  - exploit symmetry  - clm4 format
//  27d  Jun 29  vary dr step
//  27c  Jun 29  corrected true ellipse surface
//  27b  SPECIAL (WRONG) VERSION for true elliptical surface.
// varpr27a.cpp  Jun 27
// varpr20m  Jun 21  - supress near-zero values of cLm and qLm
// varpr20j  Jun 20  - higher order less singular interpolation
//  20i  Jun 16 minor changes
//  20h  Jun 13 put in singularity catches
//  20g  Jun 13  numerical singularity catching
//  20f  Jun 12 mag8 -> mag920  (changes normalization of this)
//   - put L and m for each mode   - new format clm3
//   - a strange minor problem is slight nonzero values of some
//     coefficients.  This causes problems in varpr23k.cpp
//  20e  Jun 12  - install upgraded PLm()
//  20d  Jun 9
//  20c   write "clm2" format (complex values)
//  20b  Jun 8  try and fail to have negative permittivity, because
//    of necessity to take log of permittivity.
// varper20  Jun 7   use better integration to make ".clm" file
//                ( assuming no phi dependence! )
// 19  Jun 7  smooth permittivity variation at surface
//      - linear interpolation  - cosine interpolation
//      - zero surface thickness means that 00 coefficient has infinite
//        derivative at one place.
//      - using linear permittivity interpolation, note that there is
//        greater singularity at the inner point than the outer point.
//        (This is only the case for L=2 m=0, not for other m values)
//      - write to type clm1 (.clm) file
// 18  Jun 6   constant dielectric ellipsoid
// 17    eval S derivative product integral  K( | ; )
// varper16  June 5    evaluate S lm product integrals  H( ; ; )
// varper15  - tried having (SPH,2,0) variation of permittivity, but it
//         doesn't work well.
//       - fixed error where E0z was applied to wrong L=1 component
//   14  - numerical derivatives for eps() [shifted sphere]
//       - doesn't work well
//   13  - solve system to get field coefficients at center.
//       - set up for r-dependent real-valued permittivity
//       - test for epsin<0;  A kind of instability appears!
//         Okay: -0.7 -0.8 -0.9 -0.95
//         okay: -1.1 -1.2 -1.3 -1.4 -1.6 -1.7 -1.8 -1.9 -2 -2.1 -2.2 -2.3
//         okay: -2.5  -2.9
//         NOT okay:   -0.6  -1.0  -1.5   -3.0
//   12  - omit constant V term, draw matrix on screen
//  varper11
//   - slow steps when eps is changing rapidly
//   - larger steps when r is larger
//  varper10  Jun 3
//    - some problems when m7 is not zero.
//    - could theta derivatives give problems if dtheta allows
//      theta to go beyond limits?
//    - fixed error where theta5 was out of bounds
//   9    integrates potential RSH coefficients outwards, but
//    gets inconsistent results due to random choice of angle points.
//  - also, it seems that this method cannot properly handle a
//    discontinuous jump in permittivity
//  - seems to work okay for eps(r)
//  - needs to be generalized for complex valued permittivity
//    and potential.
// varper8   solve matrix (works okay)
// varper7  set up problem as solution of matrix, to get second
//  derivative of RSH coefficients.
// varper6  June 2
// varper5.cpp June 1   Define real spherical harmonics and
//   check their real-valued-ness and orthonormality
// varper4   check orthonormality of YLm()   (works fine)
// varper3.cpp   - want to rule out domain error;  check carefully
//  - also checked against expressions for spherical harmonics
// varper2.cpp  June 1
//  June 1  define associated Legendre.  Enter website expressions for L<=4
//  and also implement derivative formula definition and check.  They agree
//  for positive m but not negative m.
// varper1.cpp  May 31, 2003
//   L = 1 m = -1 case does not agree

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

//...."real spherical harmonics" which are denoted as Slm(theta,phi).
// One nice way of doing this is given at:
// http://www.uniovi.es/~quimica.fisica/qcq/harmonics/harmonic.html
// where they say
// S Lm (theta,phi) = 0.7071  ( Y L,m + Y L,-m)   for m > 0
// S 00 (theta,phi) = Y00
// S Lm (theta,phi) = 0.7071i ( Y L,m - Y L,-m)   for m < 0
// ( Wolfram calls something similar to these
//  "tesseral harmonics", but normalization is not made clear. )
// [ But this approach needs to be adapted with a factor of (-1)^m ]

// http://mathworld.wolfram.com/SphericalHarmonic.html
// Y L^m (th,ph) = sqrt( (2L+1)/(4pi) (L-m)!/(L+m)! )
//                   P L ^m (cos(th)) exp(i m phi)

// Y 0 0 = (1/2) 1/sqrt(pi)
// Y 1-1 = (1/2) sqrt(3/2pi) sin(theta) exp(-i phi)
// Y 1 0 = (1/2) sqrt(3/pi) cos(theta)
// Y 1 1 = (-1/2) sqrt(3/2pi) sin(theta) exp(i phi)
// Y 2-2 = (1/4) sqrt(15/2pi) sin(theta)^2 exp(-2iphi)
// Y 2-1 = (1/2) sqrt(15/2pi) sin(theta)cos(theta) exp(-i phi)
// Y 2 0 = (1/4) sqrt(5/pi) (3cos^2-1)
// Y 2 1 = (-1/2) sqrt(15/2pi) sin() cos() exp(i phi)
// Y 2 2 = (1/4) sqrt(15/2pi) sin()^2 exp(2 i phi)
// Y 3-3 = (1/8) sqrt(35/pi) sin()^3 exp(-3 i phi)
// Y 3-2 = (1/4) sqrt(105/2pi) sin()^2 cos() exp(-2iphi)
// Y 3-1 = (1/8) sqrt(21/pi) sin() (5cos()^2-1) exp(-iphi)
// Y 3 0 = (1/4) sqrt(7/pi) (5cos()^3 - 3 cos() )

// http://mathworld.wolfram.com/LegendrePolynomial.html
// Associated Legendre Functions:
// P L ^m (x) =   [ ?? double check later ]
//  (-1)^m / (2^L L!) (1-x^2)^(m/2) d(L+m)/dx(L+m) (x^2-1)^L

// P 0 0 = 1
// P 1 0 = x
// P 1 1 = (-1) sqrt(1-x^2)
// P 2 0 = (1/2) (3 x^2 - 1)
// P 2 1 = (-3x) sqrt(1-x^2)
// P 2 2 = 3 (1-x^2)

// Strangely, plus m and minus m are not the same!
// P L ^-m (x) = (-1)^m (L-m)!/(L+m)! P L ^m (x)
// Abramowitz and Stegun (1972) notation is that
// m can be a subscript or superscript.
// P L _m (x) = (-1)^m P L ^m (x)         (eq. 36)

// http://functions.wolfram.com/Polynomials/LegendreP2/
double PLm(int L,int m,double x);
double fact(int i3);
complex YLm(int L,int m,double theta,double phi);
double Pi;
double pow5(double x,double y);
double SLm(int L,int m,double theta,double phi);

double eval_eps(double r5,double theta5,double phi5);
double eval_depsdr(double r5,double theta5,double phi5);
double eval_depsdtheta(double r5,double theta5,double phi5);
double eval_depsdphi(double r5,double theta5,double phi5);

#define MAXNOPOTCOEF 6
int pcL[MAXNOPOTCOEF],pcm[MAXNOPOTCOEF];
int nopotcoef_infile;

#define LMAX 10
double plmc_[LMAX+1];
double factlut[LMAX+LMAX+1];
void calculatecLmqLm(double r8,int thetasteps,int symflag);

int Lshake,mshake;
complex cLm[MAXNOPOTCOEF];
complex qLm[MAXNOPOTCOEF];
complex old_cLm[MAXNOPOTCOEF];
complex old_qLm[MAXNOPOTCOEF];
double r8avg,mag920,surfthick;
double ella,ellb;
complex epsr8in,epsr8out;

void main(void)
{ // main
Pi = 4.0*atan(1.0);
int i43;
factlut[0]=1.0;
for (i43=1;i43<=2*LMAX;i43++)
{ //,, i43
factlut[i43]=factlut[i43-1]*i43;
} //,, i43

clrscr();
int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();

FILE *fout;
fout=fopen(OutFileName,"wb");

randomize();
gotoxy(2,1);cout<<VERN<<"  Generate .clm file (clm4 format - varying dr) ";
if (ExactEllipsoid==1) {gotoxy(2,2);cout<<"Exact Ellipsoid";}
fprintf(fout,"%s%s%s%s%c%c","clm4 ",VERN,"  ",OutFileName,13,10);

gotoxy(2,8);cout<<"Output File: "<<OutFileName;

double x3,y3,z3,theta3,phi3,r3,r32,rxy;
long int licount;
int iL1,im1,iL2,im2,iL3,im3;
double S1,S2a,S3a,S2b,S3b,sum1,sum2,sum3;
int irow3;
int nocoef,L1,m1,L2,m2,L3,m3;
long int lidr8,maxlidr8;
maxlidr8=pow(2,20); // dictated by machine long int overflow limit
gotoxy(2,4);cout<<"maxlidr8 = "<<maxlidr8<<" ";
double r8,minpossdr8,maxalloweddr8;
double x8,y8,z8,r82;
double r8min,r8max;

minpossdr8=0.000001; // related to 6 digits in file format

r8avg=1.0; // UUU
maxalloweddr8=r8avg*0.05; // PPP
epsr8in = complex(4,0); // UUU
epsr8out= complex(1,0);  // UUU
surfthick=0.001; // PPP / UUU
gotoxy(2,9);cout<<"mag920? ";cin>>mag920;
if (mag920<=-1) return;
Lshake=2; // UUU
mshake=0; // UUU
ella=r8avg+mag920*SLm(Lshake,mshake,0,0);
ellb=r8avg+mag920*SLm(Lshake,mshake,0.5*Pi,0);


pcL[0] = 0; pcm[0] = 0;  // UUU
pcL[1] = 2; pcm[1] = 0;  // UUU
pcL[2] = 4; pcm[2] = 0;  // UUU
pcL[3] = 6; pcm[3] = 0;  // UUU
pcL[4] = 8; pcm[4] = 0;  // UUU
nopotcoef_infile=5;      // UUU
if (nopotcoef_infile>MAXNOPOTCOEF) {cout<<"** Error 3425 **";exit(0);}

gotoxy(50,10);cout<<"eps in = "<<epsr8in<<" ";
gotoxy(50,11);cout<<"eps out = "<<epsr8out<<" ";
gotoxy(50,12);cout<<"surf thick = "<<surfthick<<" ";
fprintf(fout,"%12.6f%s%c%c",r8avg," = r8avg (m)",13,10);
fprintf(fout,"%12.6f%12.6f%s%c%c",real(epsr8in),imag(epsr8in)," = epsr8in (rel)",13,10);
fprintf(fout,"%12.6f%12.6f%s%c%c",real(epsr8out),imag(epsr8out)," = epsr8out (rel)",13,10);
gotoxy(50,15);cout<<"mag920 "<<mag920<<" ";
fprintf(fout,"%12.6f%s%c%c",surfthick," = surf. thickness (m)",13,10);
fprintf(fout,"%12.6f%3d%3d%s%c%c",mag920,Lshake,mshake," = mag SLm (m)",13,10);

fprintf(fout,"%4d",nopotcoef_infile);
int i1;for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
fprintf(fout,"%4d%4d",pcL[i1],pcm[i1]);
} // i1
fprintf(fout,"%s%c%c"," n: L m ",13,10);

fprintf(fout,"%s%c%c","[  eps rel = sum qlm Slm  b = log(eps/eps0)= sum clm Slm  ]",13,10);
fprintf(fout,"%s%c%c","=========r==========q[0]==iiiiiiiiiiii=======q[1]==iiiiiiiiiiii=======q40==iiiiiiiiiiii=======c[0]==iiiiiiiiiiii=======c[1]==iiiiiiiiiiii=======c40==iiiiiiiiiiii xx",13,10);

//gotoxy(2,11);cout<<"dr8 = "<<dr8<<" m";

//r8min=1.0-0.5*1.05*(surfthick+1.5*fabs(mag920));
//r8max=1.0+0.5*1.05*(surfthick+1.5*fabs(mag920));
r8min=0.1;
r8max=1.5;

double sfx11=320*16.0;
int bx11=320;

for (r8=r8min;r8<r8max;r8=r8+0.05)
{ // r8
calculatecLmqLm(r8,600,1);
putpixel(bx11+sfx11*(r8-r8avg),240-15*abs(qLm[0]),14);
putpixel(bx11+sfx11*(r8-r8avg),240-30*abs(qLm[1]),13);
putpixel(bx11+sfx11*(r8-r8avg),240-40*abs(qLm[2]),12);
putpixel(bx11+sfx11*(r8-r8avg),240-50*abs(qLm[3]),11);
} // r8

r8=r8min;
if (r8<0) return;
gotoxy(2,12);cout<<"r8 min = "<<r8min<<" ";
gotoxy(2,13);cout<<"r8 max = "<<r8max<<" ";

double maxchange;
maxchange=0.03; // PPP
gotoxy(2,14);cout<<"Max change "<<maxchange<<" ";

int counter1;
counter1=0;
long int lifastcalls=0;
double olddr8;
lidr8=8192;  // must be a power of 2
LB_nextr:
counter1++;
gotoxy(5,5);cout<<"r8 = "<<r8<<"    ";
gotoxy(5,6);cout<<"nclm = "<<counter1<<" steps ";
gotoxy(5,7);cout<<"[ "<<lifastcalls<<" test steps ]";
if (kbhit()) return;

calculatecLmqLm(r8,3000,1);

putpixel(bx11+sfx11*(r8-r8avg),240-15*abs(qLm[0]),14);
putpixel(bx11+sfx11*(r8-r8avg),240-30*abs(qLm[1]),13);
putpixel(bx11+sfx11*(r8-r8avg),240-40*abs(qLm[2]),12);
putpixel(bx11+sfx11*(r8-r8avg),240-50*abs(qLm[3]),11);

// June 21 03:
for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
if (abs(qLm[i1])<10*minpossdr8) qLm[i1]=complex(0,0);
if (abs(cLm[i1])<10*minpossdr8) cLm[i1]=complex(0,0);
} // i1

if (r8<r8max)
{ //.. write to file

fprintf(fout,"%c%12.6f",'*',r8);

for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
fprintf(fout,"%12.6f%12.6f",real(qLm[i1]),imag(qLm[i1]));
} // i1

for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
fprintf(fout,"%12.6f%12.6f",real(cLm[i1]),imag(cLm[i1]));
} // i1

fprintf(fout,"%s%c%c"," xx ",13,10);

for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
old_qLm[i1]=qLm[i1];
old_cLm[i1]=cLm[i1];
} // i1

double possnewr8,possdr8,maxdiff,diff8;

if (((lidr8*minpossdr8<0.5*maxalloweddr8)||(r8>r8avg))&&(lidr8<maxlidr8))
 lidr8=lidr8*2;
olddr8=lidr8*minpossdr8;
gotoxy(60,24);cout<<"r8="<<r8<<" ";
LB532:
olddr8=lidr8*minpossdr8;
putpixel(bx11+sfx11*(r8-r8avg),475-4.0*log(lidr8*1.0)/log(2.0),11);
possnewr8=r8+olddr8;

lifastcalls++;
calculatecLmqLm(possnewr8,600,1);

int worstone;
maxdiff=0;
for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
diff8 = abs(old_qLm[i1]-qLm[i1]);
if (diff8>maxdiff) {maxdiff=diff8;worstone=i1;}
diff8 = abs(old_cLm[i1]-cLm[i1]);
if (diff8>maxdiff) {maxdiff=diff8;worstone=100+i1;}
} // i1

if ((maxdiff>maxchange)&&(lidr8>1)) {lidr8=lidr8/2; goto LB532;}

gotoxy(55,22);cout<<"worst="<<worstone<<"  ";
gotoxy(55,23);cout<<"lidr8="<<lidr8<<"  ";

r8=possnewr8;
putpixel(bx11+sfx11*(r8-r8avg),475,7);
putpixel(bx11+sfx11*(r8-r8avg),0+475-4.0*log(lidr8*1.0)/log(2.0),14);
putpixel(bx11+sfx11*(r8-r8avg),1+475-4.0*log(lidr8*1.0)/log(2.0),14);
putpixel(bx11+sfx11*(r8-r8avg),2+475-4.0*log(lidr8*1.0)/log(2.0),14);
putpixel(bx11+sfx11*(r8-r8avg),3+475-4.0*log(lidr8*1.0)/log(2.0),14);

goto LB_nextr;
} //.. write to file

fprintf(fout,"%s%c%c","============",13,10);
fclose(fout);
gotoxy(1,25);cout<<"* done *";getch();return;

} // main

double PLm(int L,int m,double x)
{ //   PLm   associated Legendre function   ( superscript m in Wolfram and ciks )
// Updated: June 10, 2003 up to order L=8
// Reference: http://ciks.cbt.nist.gov/~garbocz/paper134/node10.html
// "Appendix A: List of associated Legendre functions"
// Note that their sign convention differs from that of Wolfram.
// I don't know which is more common.
double PLm1=-1234.5,s,x2,x3,x4,x5,x6,x7,x8,s2,s3,s4,s5,s6,s7,s8;
int am,floorhalfL,ik,ik2; double d3;
x2=x*x;
s = sqrt(1-x2); // ciks.cbt.nist.gov notation
s2=s*s;
am=abs(m);
if (am>L) {cout<<"* Error: PLm |m|>L *";exit(0);}
// I think these are following the Abramowitz and Stegun definition,
// and include the Condon-Shortley phase.
// (Abramowitz and Stegun use subscript m for the other case)
if (L>=9)
{ // L >= 9

floorhalfL = L/2; // intentional integer division
for (ik=0;ik<=L;ik++) plmc_[ik]=0;
for (ik=0;ik<=floorhalfL;ik++)
{ // ik
if ((ik%2)==1)
 plmc_[L-2*ik] = -1;
else
 plmc_[L-2*ik] = 1;
plmc_[L-2*ik]*= factlut[2*L-2*ik];
plmc_[L-2*ik]/= factlut[L-ik];
plmc_[L-2*ik]/= factlut[ik];
plmc_[L-2*ik]/= factlut[L-2*ik];
} // ik
d3=pow(2.0,L);
for (ik=0;ik<=L;ik++) plmc_[ik]/=d3;

// Take derivative am times:
for (ik=0;ik<am;ik++)
{ // ik
plmc_[0]=0;
int minik2;
minik2=((am-ik-1)/2)*2;
if (minik2<0) minik2=0;
minik2=minik2+(1+L-ik)%2;
for (ik2=minik2;ik2<L-ik;ik2=ik2+2)
{ // ik2
plmc_[ik2]=plmc_[ik2+1]*(ik2+1);
plmc_[ik2+1]=0;
} // ik2
plmc_[L-ik]=0;
} // ik

// Evaluate polynomial:
PLm1 = plmc_[L-am];
for (ik=L-am;ik>=2;ik=ik-2)
{ //,,
PLm1 = x2*PLm1 + plmc_[ik-2];
} //,,
if (((L-am)%2)==1) PLm1=PLm1*x;

// Multiply by prefactor:
PLm1 = PLm1 * pow(s,am);

} // L >= 9
else
{ // L < 9
if (L>=3)
{ // >=3
x3=x2*x; s3=s2*s;
x4=x3*x; s4=s3*s;
if (L>=5)
{ // >=5
x5=x4*x; s5=s4*s;
x6=x5*x; s6=s5*s;
if (L>=7)
{ // >=7
x7=x6*x; s7=s6*s;
x8=x7*x; s8=s7*s;
if (L>=9)
{ // >=9
cout<<"* Error L>=9 here 5324 * ";exit(0);
} // >=9
} // >=7
} // >=5
} // >=3
// Note: In ciks.cbt.nist.gov webpage, s = sqrt(1-x*x)

if (L==0) PLm1=1;
else if (L==1)
{ // L=1
if      (am==0) PLm1=x;
else if (am==1) PLm1=s;
} // L=1
else if (L==2)
{ // L=2
if      (am==0) PLm1=0.5*(3*x2-1.0);
else if (am==1) PLm1=3.0*x*s;
else if (am==2) PLm1=3.0*s2;
} // L=2
else if (L==3)
{ // L=3
if      (am==0) PLm1=(1.0/2.0)*(5*x3-3.0*x);
else if (am==1) PLm1=(3.0/2.0)*(5.0*x2-1.0)*s;
else if (am==2) PLm1=15.0*x*s2;
else if (am==3) PLm1=15.0*s3;
} // L=3
else if (L==4)
{ // L=4
if      (am==0) PLm1=(1.0/8.)*(35*x4-30*x2+3.0);
else if (am==1) PLm1=2.5*(7*x3-3.0*x)*s;
else if (am==2) PLm1=7.5*(7*x2-1)*s2;
else if (am==3) PLm1=105.0*x*s3;
else if (am==4) PLm1=105.0*s4;
} // L=4
else if (L==5)
{ //     L= 5
if      (am==0) PLm1=(1.0/8.0)*(63*x5-70*x3+15*x);
else if (am==1) PLm1=(15.0/8.0)*(21*x4-14*x2+1)*s;
else if (am==2) PLm1=(105.0/2.0)*(3*x3-x)*s2;
else if (am==3) PLm1=(105.0/2.0)*(9*x2-1)*s3;
else if (am==4) PLm1=945.0*x*s4;
else if (am==5) PLm1=945.0*s5;
} //     L= 5
else if (L==6)
{ //     L= 6
if      (am==0) PLm1=(1.0/16.0)*(231*x6-315*x4+105*x2-5);
else if (am==1) PLm1=(21.0/8.0)*(33*x5-30*x3+5*x)*s;
else if (am==2) PLm1=(105.0/8.0)*(33*x4-18*x2+1)*s2;
else if (am==3) PLm1=(315.0/2.0)*(11*x3-3*x)*s3;
else if (am==4) PLm1=(945.0/2.0)*(11*x2-1.0)*s4;
else if (am==5) PLm1=10395.0*x*s5;
else if (am==6) PLm1=10395.0*s6;
} //     L= 6
else if (L==7)
{ //     L= 7
if      (am==0) PLm1=(1.0/16.0)*(429*x7-693*x5+315*x3-35*x);
else if (am==1) PLm1=(7.0/16.0)*(429*x6-495*x4+135*x2-5)*s;
else if (am==2) PLm1=(63.0/8.0)*(143*x5-110*x3+15*x)*s2;
else if (am==3) PLm1=(315.0/8.0)*(143*x4-66*x2+3)*s3;
else if (am==4) PLm1=(3465.0/2.0)*(13*x3-3*x)*s4;
else if (am==5) PLm1=(10395.0/2.0)*(13*x2-1)*s5;
else if (am==6) PLm1=135135.0*x*s6;
else if (am==7) PLm1=135135.0*s7;
} //     L= 7
else if (L==8)
{ //     L= 8
if      (am==0) PLm1=(1.0/128.0)*(6435*x8-12012*x6+6930*x4-1260*x2+35);
else if (am==1) PLm1=(9.0/16.0)*x*(715*x6-1001*x4+385*x2-35)*s;
else if (am==2) PLm1=(315.0/16.0)*(143*x6-143*x4+33*x2-1)*s2;
else if (am==3) PLm1=(3465.0/8.0)*x*(39*x4-26*x2+3)*s3;
else if (am==4) PLm1=(10395.0/8.0)*(65*x4-26*x2+1)*s4;
else if (am==5) PLm1=(135135.0/2.0)*(5*x3-x)*s5;
else if (am==6) PLm1=(135135.0/2.0)*(15*x2-1)*s6;
else if (am==7) PLm1=2027025.0*x*s7;
else if (am==8) PLm1=2027025.0*s8;
} //     L= 8
} // L < 9

if (m<0)
{ // m<0
PLm1 = PLm1*factlut[L-am]/factlut[L+am];
if ((am%2)==1) PLm1*=-1;
} // m<0
// [ Strangely, PLm for plus m and minus m are not the same! ]
// P L ^-m (x) = (-1)^m (L-m)!/(L+m)! P L ^m (x)
// eq. (42) on http://mathworld.wolfram.com/LegendrePolynomial.html
// Also eq (30) in http://ciks.cbt.nist.gov/~garbocz/paper134/node10.html
// I am assuming that zero factorial is 1 in this.
return PLm1;
} //  PLm   associated Legendre function

double fact(int i3)
{ //   fact
double fact1;
if (i3<0) {cout<<"* Error factorial of "<<i3<<" *";exit(0);}
fact1=1.0;
int i4;
for (i4=2;i4<=i3;i4++)
{ // i4
fact1=fact1*i4;
} // i4
return fact1;
} //   fact

complex YLm(int L,int m,double theta,double phi)
{ //    YLm (superscript m)
complex YLm1;
YLm1=(2.0*L+1.0)/(4.0*Pi);
YLm1=YLm1 * fact(L-m) / fact(L+m);
YLm1=sqrt(YLm1);
YLm1=YLm1 * PLm(L,m,cos(theta));
YLm1=YLm1 * exp(complex(0.0,m*phi));
return YLm1;
} //   YLm (superscript m)

double pow5(double x,double y)
{ //   pow5
// June 13th: upgraded to give error for for 0 to negative power
double pow51;
if (x==0.0)
{ //
// June 13th:
if (y<0) {cout<<"* pow5(): 0 to negative power y = "<<y<<" *";exit(0);}
cout<<"** Warning 0^0 : given as 0 **";
return 0.0;
} //
else if (x<0.0)
{ //
if (y!=ceil(y)) {cout<<"* pow5(): negative to noninteger: "<<x<<"^"<<y<<" *";exit(0);}
} //
pow51 = pow(x,y);
return pow51;
} //   pow5

double SLm(int L,int m,double theta,double phi)
{ // SLm    Real Spherical Harmonics
//One nice way of doing this is given at:
//http://www.uniovi.es/~quimica.fisica/qcq/harmonics/harmonic.html
//where they say
// S Lm (theta,phi) = 0.7071  ( Y L,m + Y L,-m)   for m > 0
// S L0 (theta,phi) = Y L0                        for m = 0
// S Lm (theta,phi) = 0.7071i ( Y L,m - Y L,-m)   for m < 0
//  [ But apparently this doesn't work with the kind of spherical
//    harmonics that have m as a superscript.  But maybe it would
//    work the spherical harmonics with m as a subscript. ]
double p707; p707 = 1.0/sqrt(2.0);
double SLm1;
if (m==0)
{ //,
SLm1 = real(YLm(L,m,theta,phi));
} //,
else if (m>0)
{ //.    m>0
SLm1 = real(p707*(YLm(L,m,theta,phi)+pow(-1.0,m)*YLm(L,-m,theta,phi)));
} //.    m>0
else if (m<0)
{ //..   m<0
SLm1 = real(p707*complex(0,1.0)*(YLm(L,m,theta,phi)-pow(-1.0,m)*YLm(L,-m,theta,phi)));
} //..   m<0
return SLm1;
} // SLm   Real Spherical Harmonics

void calculatecLmqLm(double r8,int thetasteps,int symflag)
{ // calculate cLm qLm for radius r
int i1;
double theta8,phi8,r8surf;
double dtheta8=Pi/thetasteps;
double r81,r82b;
complex epsr8;
double x,y,x2;

for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
cLm[i1]=0;qLm[i1]=0;
} // i1

phi8=0.5; // value does not matter here
double maxtheta8;
maxtheta8=Pi;
if (symflag==1) maxtheta8=0.5*Pi;
for (theta8=0.5*dtheta8;theta8<maxtheta8;theta8=theta8+dtheta8)
{ // theta8

if (ExactEllipsoid!=1)
{ // normal case
r8surf = r8avg + mag920*SLm(Lshake,mshake,theta8,phi8);
} // normal case
else
{ // exact elipsoid
double tpar;
if (theta8<0.25*Pi)
{ //.
tpar=atan(ella*tan(theta8)/ellb);
} //.
else if (theta8<0.75*Pi)
{ //..
tpar=0.5*Pi+(ellb/ella)*tan(theta8-0.5*Pi);
} //..
else
{ //...
tpar=Pi-atan(ella*tan(Pi-theta8)/ellb);
} //...
r8surf=sqrt(ella*ella*cos(tpar)*cos(tpar)+ellb*ellb*sin(tpar)*sin(tpar));
} // exact elipsoid

r81=r8surf-0.5*surfthick;
r82b=r8surf+0.5*surfthick;

if (r8<=r81)
 epsr8 = epsr8in;
else if (r8>r82b)
 epsr8 = epsr8out;
else
{ //..
// y = x simplest version
// y = a x + b x^3    dy = a + 3 b x^2
// a + b = 1          a + 3b = 0    b = (-1/3)a   a + (-1/3) a = 1  a = 1.5 b = -0.5
// y = a x + b x^3 + c x^5    dy=a+3bx^2+5cx^4   d2y=6bx+20cx^3
// a + b + c = 1    a + (-2/3)a + (-2/3)(-6/20)a = 1
// a + 3b + 5c = 0   a + 3b -1.5b = 0  a + 1.5b=0  b = -2/3a
// 6b + 20c = 0     c = -6/20 b     a - 2/3a + (1/5)a = 1   a = 1.875
// a = 1.875  b = -1.25  c = 0.375
// -----------------------------
// y =   a x + b x^3 + c x^5 +  dx^7
// dy =  a   + 3bx^2 + 5cx^4 + 7dx^6
// d2y = 6 b x + 20 c x^3 + 42dx^5
// d3y = 6b + 60cx^2 + 210 dx^3
// a +  b +  c +  d = 1
// a + 3b + 5c + 7d = 0
// 6b + 20c + 42d = 0
// 6b + 60c + 210d = 0
//      40c + 168d = 0      d = -5/21 c
// 6b + 20c - 10c = 0  6b + 10c = 0     c = -3/5 b   d = 1/7 b
// a + 3b - 3b + b = 0  b = -a
// a -a + 3/5 a - 1/7a = 1
//  a = 2.1875  b = -2.1875 c = 1.3125 d = -0.3125

x=2.0*(r8-r8surf)/surfthick;
x2=x*x;
// y=x
// y=1.5*x-0.5*x*x*x;
// y=1.875*x-1.25*x*x*x+0.375*x*x*x*x*x;
//y=2.1875*x-2.1875*x*x*x+1.3125*x*x*x*x*x-0.3125*x*x*x*x*x*x*x;
y=x*(2.1875-x2*(2.1875-1.3125*x2+0.3125*x2*x2));
// Note that x ranges from -1 to 1 in the transition region
epsr8 = 0.5*(epsr8in+epsr8out) + 0.5*(epsr8out-epsr8in)*y;
} //..

for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
qLm[i1]+=    epsr8   * SLm(pcL[i1],pcm[i1],theta8,phi8)*2*Pi*sin(theta8)*dtheta8;
cLm[i1]+= log(epsr8) * SLm(pcL[i1],pcm[i1],theta8,phi8)*2*Pi*sin(theta8)*dtheta8;
} // i1

} // theta8

if (symflag==1)
{ // special case of symmetry
for (i1=0;i1<nopotcoef_infile;i1++)
{ // i1
qLm[i1]*=2;
cLm[i1]*=2;
} // i1
} // special case of symmetry

return;
} // calculate cLm qLm