#define VERN "varpr28j.cpp"
// [ Use with varpr27e.cpp (makes clm4) and varpr24g.cpp (testing) ]

#define ClmFileName "emp3d.clm"
//  28j  look for and trap pow domain errors
//  28i
//  28h speed up clm, etc.  works
//  28g June 30 correct dr.  Works, but is slow.
//  28f Jun 30 - reads clm4 file (with varying dr step)
//  28d
// varpr28c Jun 28  use new_PLm which is generalized to any order L
//  - special demonstration version for p05a.clm
//  28b Jun 27  change clmq and clmc() not to use L and m
//     - imag(E) was displayed erroneously
//  varpr28a  Jun 27
//  23z June 26  try to improve method (3) - no progress
//     got rid of "old exact"
//  23y - method 3 for dipole moment, based on equation [12.5]; test
//  23x  Jun 20  corrected exact ellipsoid dipole moment formula
//               ( but new formula is only a conjecture )
//             - introduce eps0 constant
//  23w Jun 18  using corrected formulas for dipole moment with
//               factor of 3/(2*eps8out+1)
//  23v Jun 18  still trying...
// varpr23u  Jun 18  try to resolve discrepancy
//     Now find discrepancy between d2aLm peak and daLm cusp.
//  23t  Jun 16  Charge density integration to calc. dipole moment
//  23s  June 16  Clean up dipole moment calculation
//  23r  Jun 15  Calculate dipole moment by integrating P.  Works
//    okay only if I introduce a reasonable cutoff in radius.
//    Something wrong happens at very high radius.
//  23q  Jun 15  calculate again to see if we get right outside field
//  23p  Jun 14  exact ellipsoid expression: prolate and oblate
//   - introduce fact() extended factorial such that 0! = 1
//  23n  wLm(r)  improved and simplified dr5 method; works well
//  23m  Jun 13  ...make it faster?
//       -  numerical singularity catching
//  23k  Jun 12  reads clm3 format file  - change mag820 to mag920
//        - compare with linearization of exact ellipsoid solution
//  23j
//  23i  Jun 12  - I realize auto step size has a possible weakness
//                 since it may not adjust to a very sudden change
//              - Noticed H() and K() are much faster now that divide
//                by zero problem is gone.
//              - tried increasing LMAX and HowManyaLm, as well as
//                setting m=1
//  - improved auto step size to handle the "crisis region"
// varpr23h  Jun 11 - HowManyaLm=2  ceL and cem()
//      - variable integration step size based on limiting rate
//        of change of variables.  works well
//      - checked for zero denominators in K()
//      - founds missing abs()
//      - fixed factored integration for zero integrand case
// 23g  use look up tables for H and K values.  works fine.
// 23f  still has strange K1v output sometimes; gone now
// 23e  added missing factor of 0.488603 (part of S10)
//      - put in new PLm going up to L=8
//      - warning nodivis=50 temporarily only!
//      - I get wrong results and strange printouts for K1v!
// 23c  reads "clm2" format with complex values
//      ( .clm file is generated by varper20.cpp )
//  23b  - cleaned up
//    ( keep in mind that permittivity will eventually be complex )
//  23  Jun 8  simple working version, finds E inside object
//     - seems to work OK
//  22  Jun 8  very crude integration
// varper21 Jun 7   read in .clm file into arrays
//      clm() and qlm() smoothed functions using arrays
//     also have smoothed derivative functions dclm and dqlm()
// 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 Eincz 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 <dos.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 L0 (theta,phi) = Y L 0                       for m = 0
// 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)

// [Note: Be cautious about these, since sign conventions differ!]
// 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

// [Note: Be cautious about these, since sign conventions differ!]
// 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 new_PLm(int L,int m,double x);
double fact(int i3);
complex YLm(int L,int m,double theta,double phi);
double Pi,eps0;
double pow5(double x,double y);
double pow7(double x,double y,int from);
double SLm(int L,int m,double theta,double phi);

#define LMAX 15
double plmc[LMAX+1];
double factlut[LMAX+LMAX+1];
complex aLm[LMAX+1][LMAX+LMAX+1];
complex daLm[LMAX+1][LMAX+LMAX+1];
complex d2aLm[LMAX+1][LMAX+LMAX+1];
complex wLm[LMAX+1][LMAX+LMAX+1];
complex dwLm[LMAX+1][LMAX+LMAX+1];
complex d2wLm[LMAX+1][LMAX+LMAX+1];
complex dLm[LMAX+1][LMAX+LMAX+1];

#define HowManyaLm 5
int     HowManycLm=3; // PPP

complex Fmat[HowManyaLm][HowManyaLm+1];

char readclmline(void);
FILE *clmfin;

#define MAXNCLM 795
int nclm;
double clmr[MAXNCLM];
int i9alasttime=0;
complex clm(int ipc,double r);
complex dclm(int ipc,double r);
complex qlm(int ipc,double r);
complex dqlm(int ipc,double r);
double H(int L1,int m1,int L2,int m2,int L3,int m3);
double K(int L1,int m1,int L2,int m2,int L3,int m3);
double dSLm_theta(int L,int im,double theta,double phi);
double dSLm_phi(int L,int im,double theta,double phi);
complex dYLm_theta(int L,int m,double theta,double phi);
double dPLm_theta(int L,int im,double theta);

#define KLUTMAX 302
int klutln;
int klutl1[KLUTMAX];
int klutm1[KLUTMAX];
int klutl2[KLUTMAX];
int klutm2[KLUTMAX];
int klutl3[KLUTMAX];
int klutm3[KLUTMAX];
double klutva[KLUTMAX];
int hlutln;
int hlutl1[KLUTMAX];
int hlutm1[KLUTMAX];
int hlutl2[KLUTMAX];
int hlutm2[KLUTMAX];
int hlutl3[KLUTMAX];
int hlutm3[KLUTMAX];
double hlutva[KLUTMAX];
int ceL(int icase);
int cem(int icase);

#define MAXNOPOTCOEF 5
int pcL[MAXNOPOTCOEF],pcm[MAXNOPOTCOEF];
float fqr[MAXNOPOTCOEF],fqi[MAXNOPOTCOEF];
float fcr[MAXNOPOTCOEF],fci[MAXNOPOTCOEF];
complex clmq[MAXNOPOTCOEF][MAXNCLM];
complex clmc[MAXNOPOTCOEF][MAXNCLM];

int iwhichaxisforE;
int iwhichaxisforp;
long int count0to0=0;

void main(void)
{ // main
Pi = 4.0*atan(1.0);
eps0 = 8.854187817e-12;
int i43;
factlut[0]=1.0;
for (i43=1;i43<=2*LMAX;i43++)
{ //,, i43
factlut[i43]=factlut[i43-1]*i43;
} //,, i43
clrscr();
clmfin=fopen(ClmFileName,"rb");
int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();

char chEaxis;
gotoxy(2,1);cout<<VERN;
gotoxy(3,3);cout<<"Which axis to apply external E field along? (press x, y or z) ";chEaxis=getch();
iwhichaxisforE=-2;
if (chEaxis=='x') iwhichaxisforE=1;
if (chEaxis=='y') iwhichaxisforE=-1;
if (chEaxis=='z') iwhichaxisforE=0;
if (iwhichaxisforE<-1) {cout<<"* no such axis *";exit(0);}
gotoxy(3,3);cout<<"                                                                           ";

int i49;
klutln=0;
hlutln=0;
for (i49=0;i49<KLUTMAX;i49++)
{ //..i49
klutl1[i49]=-99;
klutm1[i49]=-99;
klutl2[i49]=-99;
klutm2[i49]=-99;
klutl3[i49]=-99;
klutm3[i49]=-99;
hlutl1[i49]=-99;
hlutm1[i49]=-99;
hlutl2[i49]=-99;
hlutm2[i49]=-99;
hlutl3[i49]=-99;
hlutm3[i49]=-99;
} //..i49

randomize();
gotoxy(2,1);cout<<VERN<<"  read clm file and calculate potential and dipole moment.";
char ch1,ch2,ch3,ch4;
char ch11;
fscanf(clmfin,"%c%c%c%c",&ch1,&ch2,&ch3,&ch4);
if ((ch1!='c')||(ch2!='l')||(ch3!='m')||(ch4!='4'))
{ //..
cout<<"* wrong file type: "<<ch1<<ch2<<ch3<<ch4<<" *\n";
cout<<"* .clm file must be type clm4 *\n";
cout<<"* the file you selected is type "<<ch1<<ch2<<ch3<<ch4<<" *";
return;
} //..
ch11=readclmline();
float fr8avg,feps8inr,feps8outr,feps8ini,feps8outi,fsurfthick,fmag920;
double r8,theta8,phi8,r8avg,r8surf,mag920,dr8;
double eps8,r82b,surfthick;
complex eps8in,eps8out;

fscanf(clmfin,"%f",&fr8avg); readclmline();
gotoxy(50,10);cout<<"r8avg = "<<fr8avg<<" ";
r8avg=fr8avg;
fscanf(clmfin,"%f%f",&feps8inr,&feps8ini); readclmline();
eps8in=complex(feps8inr,feps8ini);
gotoxy(50,11);cout<<"eps8in =  "<<eps8in<<" ";
fscanf(clmfin,"%f%f",&feps8outr,&feps8outi); readclmline();
eps8out=complex(feps8outr,feps8outi);
gotoxy(50,12);cout<<"eps8out = "<<eps8out<<" ";
cout.precision(1);
complex Esphere;
double Eincmag=10000.0;
double Eincx = 0;    // UUU (External electric field magnitude)
double Eincy = 0;    // UUU (External electric field magnitude)
double Eincz = 0; // UUU (External electric field magnitude)
double Eincmu[3];
complex pmu[3];
complex pmu2[3];
complex pmu3[3];
if (iwhichaxisforE==0)  Eincz = Eincmag;
if (iwhichaxisforE==1)  Eincx = Eincmag;
if (iwhichaxisforE==-1) Eincy = Eincmag; // ? not clear about sign here
Eincmu[0]=0;
Eincmu[1]=0;
Eincmu[2]=0;
Eincmu[1+iwhichaxisforE]=Eincmag;
if (iwhichaxisforE==-1) Eincmu[1+iwhichaxisforE]=-Eincmag;
Esphere = Eincmag*3.0/((eps8in/eps8out)+2.0); // exact solution for sphere
int iEr,iEi;
iEr=real(Esphere);
iEi=imag(Esphere);
gotoxy(50,17);
if (iEi>0)  cout<<"Exact Sphere: "<<iEr<<"+"<<iEi<<"i V/m ";
if (iEi<0)  cout<<"Exact Sphere: "<<iEr<<iEi<<"i V/m ";
if (iEi==0) cout<<"Exact Sphere: "<<iEr<<" V/m ";
fscanf(clmfin,"%f",&fsurfthick); readclmline();
cout.precision(3);
gotoxy(50,13);cout<<"surfthick = "<<fsurfthick<<" ";
fscanf(clmfin,"%f",&fmag920); readclmline();
gotoxy(50,14);cout<<"mag920 = "<<fmag920<<" ";
mag920 = fmag920;

double ellipsea,ellipseb,ellipsec,ellipsoid_e2,ellipsoid_linnz,ellipsoid_nz,ellipsoid_nx,ellipsoid_ny;
ellipsea=r8avg+mag920*SLm(2,0,0,0);
ellipseb=r8avg+mag920*SLm(2,0,Pi*0.5,0);
ellipsec=r8avg+mag920*SLm(2,0,Pi*0.5,Pi*0.5);
cout.precision(3);
gotoxy(50,18);cout<<"a="<<ellipsea<<" b="<<ellipseb<<" c="<<ellipsec<<"  ";
ellipsoid_e2 = 1 - ellipseb*ellipseb/(ellipsea*ellipsea);
///ellipsoid_linnz = (1.0/3.0) - (2.0/15.0)*ellipsoid_e2; // linearized
///gotoxy(50,21);cout<<"e2 = "<<ellipsoid_e2<<" lin nz="<<ellipsoid_linnz<<" ";
complex Eellipsoid,Ecenter;
///Eellipsoid = Eincmag / (1.0 + ((eps8in/eps8out)-1.0)*ellipsoid_linnz);
///gotoxy(50,22);cout<<"linearized E (z) = "<<Eellipsoid<<"   ";
ellipsoid_nz=1.0/3.0;
if ((ellipsea>ellipseb)&&(ellipseb==ellipsec))
{ // prolate: a > b = c
// Eq. (4.32) on page 25 of Landau, Lifschitz and Pitaevskii
// "Electrodynamics of Continuous Media" Pergamon 1984
double e=sqrt(1.0-ellipseb*ellipseb/(ellipsea*ellipsea));
ellipsoid_nz=((1-e*e)/(2*e*e*e))*(log((1+e)/(1-e))-2*e);
} // prolate: a > b = c
if ((ellipsea<ellipseb)&&(ellipseb==ellipsec))
{ // oblate: a < b = c
// Eq. (4.34) on page 25 of Landau, Lifschitz and Pitaevskii
// "Electrodynamics of Continuous Media" Pergamon 1984
double e=sqrt(ellipseb*ellipseb/(ellipsea*ellipsea)-1.0);
ellipsoid_nz=((1+e*e)/(e*e*e))*(e-atan(e));
} // oblate: a < b = c
ellipsoid_nx=0.5*(1-ellipsoid_nz);
ellipsoid_ny=0.5*(1-ellipsoid_nz);
double nmu[3];
nmu[1-1]=ellipsoid_ny;
nmu[1+0]=ellipsoid_nz;
nmu[1+1]=ellipsoid_nx;
cout.precision(4);
gotoxy(45,24);cout<<"exact nz="<<ellipsoid_nz<<" nx="<<ellipsoid_nx<<"  ";
if (iwhichaxisforE==0)
{ //..
Eellipsoid = Eincmu[1+iwhichaxisforE] / (1.0 + ((eps8in/eps8out)-1.0)*nmu[1+iwhichaxisforE]);
int i5r,i5i;
i5r=real(Eellipsoid);
i5i=imag(Eellipsoid);
gotoxy(45,25);
if (i5i>0)  cout<<"inside E paral (z) = "<<i5r<<"+"<<i5i<<"i V/m  ";
if (i5i==0) cout<<"inside E paral (z) = "<<i5r<<" V/m  ";
if (i5i<0)  cout<<"inside E paral (z) = "<<i5r<<i5i<<"i V/m  ";
} //..
if (iwhichaxisforE==1)
{ //..
Eellipsoid = Eincmu[1+iwhichaxisforE] / (1.0 + ((eps8in/eps8out)-1.0)*nmu[1+iwhichaxisforE]);
int i5r,i5i;
i5r=real(Eellipsoid);
i5i=imag(Eellipsoid);
gotoxy(45,25);
if (i5i>0)  cout<<"inside E perp (x) = "<<i5r<<"+"<<i5i<<"i V/m  ";
if (i5i==0) cout<<"inside E perp (x) = "<<i5r<<" V/m  ";
if (i5i<0)  cout<<"inside E perp (x) = "<<i5r<<i5i<<"i V/m  ";
} //..
if (iwhichaxisforE==-1)
{ //..
Eellipsoid = Eincmu[1+iwhichaxisforE] / (1.0 + ((eps8in/eps8out)-1.0)*nmu[1+iwhichaxisforE]);
gotoxy(45,25);cout<<"inside E perp (y) = "<<-Eellipsoid<<"   ";
} //..

complex exact_p;
// See equation 8.10 on page 41 of Landau, Lifshitz and Pitaevskii 1984
// ( But I adapted it to handle the permittivity of the matrix )

// New version: June 20th: (this is only a conjecture!)
// This is equation [8.2] in www.geocities.com/cofrest/md52.htm
exact_p = (eps8in-eps8out)/(eps8out+(eps8in-eps8out)*nmu[1+iwhichaxisforE]);
exact_p*= eps0*(4*Pi/3.0)*ellipsea*ellipseb*ellipsec*Eincmu[1+iwhichaxisforE];

if (iwhichaxisforE==-1) exact_p=-exact_p;

cout.precision(3);
gotoxy(1,15);cout<<"Exact p  "<<exact_p*1e9<<" nC-m";

int nopotcoef_infile;
fscanf(clmfin,"%d",&nopotcoef_infile);
gotoxy(55,1);cout<<"nopotcoef_infile "<<nopotcoef_infile<<" ";
if (nopotcoef_infile>MAXNOPOTCOEF)
{ //
gotoxy(1,22);cout<<"* Error 9303: nopotcoef_infile too big      *";
gotoxy(1,23);cout<<"* (.clm file is too long)                   *";
gotoxy(1,24);cout<<"* Please increase MAXNOPOTCOEF in varpr28   *";return;
gotoxy(1,25);cout<<"* or else decrease no. pot. coef in varpr27 *";return;
} //
int ipc;
for (ipc=0;ipc<nopotcoef_infile;ipc++)
{ //..
fscanf(clmfin,"%d",&pcL[ipc]);
fscanf(clmfin,"%d",&pcm[ipc]);
gotoxy(55,2+ipc);cout<<pcL[ipc]<<" "<<pcm[ipc]<<" ";
} //..

if (HowManycLm>=HowManyaLm-1)
{ //..
gotoxy(1,25);cout<<"## Warning: HowManyaLm is too small ##";
} //..
gotoxy(55,7);
if (HowManycLm<nopotcoef_infile) cout<<"Only "<<HowManycLm<<" will be used";
if (HowManycLm==nopotcoef_infile) cout<<"All "<<HowManycLm<<" will be used";
if ((HowManycLm>nopotcoef_infile)||(HowManycLm>MAXNOPOTCOEF))
{ //..
cout<<"* Error HowManycLm > nopotcoef_infile *";return;
} //..
ch11=readclmline();

LB323:
ch11=readclmline();
if (feof(clmfin)) {cout<<"* file ended *";return;}
if (ch11!='=') goto LB323;
int iclmline=0;
LB324:
fscanf(clmfin,"%c",&ch1);
if (ch1!='*') goto LB326;
float fr;
fscanf(clmfin,"%f",&fr);
for (ipc=0;ipc<nopotcoef_infile;ipc++)
{ //.. ipc
fscanf(clmfin,"%f",&fqr[ipc]);
fscanf(clmfin,"%f",&fqi[ipc]);
} //.. ipc
for (ipc=0;ipc<nopotcoef_infile;ipc++)
{ //.. ipc
fscanf(clmfin,"%f",&fcr[ipc]);
fscanf(clmfin,"%f",&fci[ipc]);
} //.. ipc
if (nopotcoef_infile>MAXNOPOTCOEF) {cout<<"* Error 5932 *";exit(0);}
cout.precision(6);

clmr[iclmline]=fr;
if (iclmline<MAXNCLM)
{ //.
for (ipc=0;ipc<nopotcoef_infile;ipc++)
{ //.. ipc
clmq[ipc][iclmline]=complex(fqr[ipc],fqi[ipc]);
} //.. ipc
for (ipc=0;ipc<nopotcoef_infile;ipc++)
{ //.. ipc
clmc[ipc][iclmline]=complex(fcr[ipc],fci[ipc]);
} //.. ipc
if (nopotcoef_infile>MAXNOPOTCOEF) {cout<<"** Error 5055 **";exit(0);}
} //.
ch11=readclmline();
iclmline++;
if (iclmline>=MAXNCLM) {cout<<"* Error: Too many clm lines *";return;}
goto LB324;
LB326:
nclm =iclmline;
double rinneredge,routeredge;
rinneredge=0.0;
int i18=0;
LB328:
if ((clmr[i18+1]-clmr[i18])>0.002) {i18++; goto LB328;}
rinneredge=clmr[i18];
i18=nclm-1;
LB330:
if ((clmr[i18]-clmr[i18-1])>0.002) {i18--; goto LB330;}
routeredge=clmr[i18];
gotoxy(30,2);cout<<"clm lines: "<<nclm<<"  ";
gotoxy(30,3);cout<<"r inner = "<<rinneredge<<" ";
gotoxy(30,4);cout<<"r outer = "<<routeredge<<" ";

double r5maxforp=routeredge; // PPP (can't be so big or p becomes inaccurate)
gotoxy(30,5);cout<<"r5 max for p "<<r5maxforp<<" ";
if (r5maxforp<routeredge)
{ //
char beepc=7;
cout<<beepc;
gotoxy(1,1);cout<<"### WARNING: r5maxforp is too low ###";
gotoxy(1,25);cout<<"### WARNING: r5maxforp is too low ###";
} //

gotoxy(2,9);cout<<"Input CLM File: "<<ClmFileName<<" ";

double r81,sfx33=320.0*2.0,sfy33=82;
double rmid=1.032;
int by33=240,bx33=320;
setcolor(8);
line(bx33+sfx33*(rinneredge-rmid),0,bx33+sfx33*(rinneredge-rmid),479);
line(bx33+sfx33*(routeredge-rmid),0,bx33+sfx33*(routeredge-rmid),479);


for (r81=rmid-320/sfx33;r81<rmid+320/sfx33;r81=r81+0.013/sfx33)
{ // r81
//putpixel(bx33+sfx33*(r81-rmid),by33,7);

/*
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*1*real(clm(0,r81)),14-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*1*real(clm(0,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*1*imag(clm(0,r81)),14-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*1*imag(clm(0,r81)),1);
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*2*real(clm(1,r81)),13-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*2*real(clm(1,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*3*real(clm(2,r81)),12-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*3*real(clm(2,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*4*real(clm(3,r81)),11-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*4*real(clm(3,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-1*sfy33*1.5*5*real(clm(4,r81)),10-8);
putpixel(bx33+sfx33*(r81-rmid),2+by33-1*sfy33*1.5*5*real(clm(4,r81)),7);
*/

/*
putpixel(bx33+sfx33*(r81-rmid),  by33-0.007*sfy33*real(dclm(0,r81)),14-8);
putpixel(bx33+sfx33*(r81-rmid),  by33-0.007*sfy33*real(dclm(1,r81)),13-8);
putpixel(bx33+sfx33*(r81-rmid),  by33-0.007*sfy33*real(dclm(2,r81)),12-8);
putpixel(bx33+sfx33*(r81-rmid),  by33-0.007*sfy33*real(dclm(3,r81)),11-8);
putpixel(bx33+sfx33*(r81-rmid),  by33-0.007*sfy33*real(dclm(4,r81)),10-8);
*/

/*
putpixel(bx33+sfx33*(r81-rmid),  by33-sfy33*real(clm(2,0,r81)),13);
putpixel(bx33+sfx33*(r81-rmid),2+by33-sfy33*real(clm(2,0,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-sfy33*real(clm(4,0,r81)),12);
putpixel(bx33+sfx33*(r81-rmid),2+by33-sfy33*real(clm(4,0,r81)),7);
putpixel(bx33+sfx33*(r81-rmid),  by33-sfy33*imag(clm(2,0,r81)),13);
putpixel(bx33+sfx33*(r81-rmid),2+by33-sfy33*imag(clm(2,0,r81)),6);
putpixel(bx33+sfx33*(r81-rmid),  by33-sfy33*imag(clm(4,0,r81)),12);
putpixel(bx33+sfx33*(r81-rmid),2+by33-sfy33*imag(clm(4,0,r81)),6);
*/
} // r81

//getch();
//cleardevice();

/*
double r23; complex v1,v2;
LB939:
gotoxy(1,21);cout<<"r23?                                      ";
gotoxy(1,21);cout<<"r23? ";cin>>r23;
if (r23<0) return;
v1=clm(0,r23);
v2=qlm(0,r23);
gotoxy(1,23);cout<<"clm = "<<v1<<"          ";
gotoxy(1,24);cout<<"qlm = "<<v2<<"          ";
goto LB939;
*/

// temporarily

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;

double r5,theta5,phi5,dr5,r50,mindrever,rmindr;
mindrever=1e10;

double r5max;

int icase7,L7,m7;
icase7=-1;
LBnextcol:
icase7++;
if (icase7<HowManyaLm) {L7=ceL(icase7);m7=cem(icase7);}

r50 = 0.02; // PPP starting radius

if (icase7<HowManyaLm)
{ //...
int icase6;
for (icase6=0;icase6<HowManyaLm;icase6++)
{ // icase6
int L6,m6;
L6=ceL(icase6);m6=cem(icase6);
if (L6>LMAX)
{ //L6>LMAX
gotoxy(1,24);cout<<"* Error LMAX too small 1234  LMAX="<<LMAX<<" *";
gotoxy(1,25);cout<<"* HowManyaLm="<<HowManyaLm<<" *";
exit(0);
} //L6>LMAX
aLm[L6][L6+m6] = 0;
daLm[L6][L6+m6] = 0;
wLm[L6][L6+m6] = 0;
dwLm[L6][L6+m6] = 0;
} // icase6
} //...

if (icase7==HowManyaLm)
{ //...
int icase6;
for (icase6=0;icase6<HowManyaLm;icase6++)
{ // icase6
int L6,m6;
L6=ceL(icase6); m6=cem(icase6);
wLm[L6][L6+m6] =  dLm[L6][L6+m6];
dwLm[L6][L6+m6] = 0;
aLm[L6][L6+m6] =  wLm[L6][L6+m6]*pow5(r50,L6);
daLm[L6][L6+m6] = wLm[L6][L6+m6]*L6*pow5(r50,L6-1);
} // icase6
pmu[1-1]=0;
pmu[1+0]=0;
pmu[1+1]=0;
pmu2[1-1]=0;
pmu2[1+0]=0;
pmu2[1+1]=0;
pmu3[1-1]=0;
pmu3[1+0]=0;
pmu3[1+1]=0;
} //...

r5max=100.0; // PPP  30 made no difference; couple of checks
gotoxy(60,4);cout<<"r max = "<<r5max<<" ";

// Prefactor should be called "d":

if (icase7<HowManyaLm)
{ //..
wLm[L7][L7+m7]  = 1.0;
dwLm[L7][L7+m7] = 0.0;
aLm[L7][L7+m7]  = wLm[L7][L7+m7]*pow5(r50,L7);
daLm[L7][L7+m7] = wLm[L7][L7+m7]*L7*pow5(r50,L7-1);
} //..

r5 = r50;
gotoxy(60,5);cout<<"r0 = "<<r5<<" ";

int irow2,icol,ic2,ir2,ic3;

licount=0;
int ans1,ans2,ans3; // count actual number of steps
ans1=0; ans2=0; ans3=0;
LBnextr:
licount++;

// Calculate d2aLm here:
int icase1;
for (icase1=0;icase1<HowManyaLm;icase1++)
{ // icase1
int L1,m1; L1=ceL(icase1); m1=cem(icase1);
d2aLm[L1][L1+m1]= 0;
d2aLm[L1][L1+m1]+= 2*r5*daLm[L1][L1+m1];
d2aLm[L1][L1+m1]-= L1*(L1+1)*aLm[L1][L1+m1];
for (int icase2=0;icase2<HowManyaLm;icase2++)
{ // icase2
int L2,m2; L2=ceL(icase2); m2=cem(icase2);
for (int i3case3=0;i3case3<HowManycLm;i3case3++)
{ // i3case3
int L3,m3; L3=pcL[i3case3]; m3=pcm[i3case3];
d2aLm[L1][L1+m1]+= r5*r5*daLm[L2][L2+m2]*dclm(i3case3,r5)*H(L1,m1,L2,m2,L3,m3);
d2aLm[L1][L1+m1]+= aLm[L2][L2+m2]*clm(i3case3,r5)*K(L1,m1,L2,m2,L3,m3);
} // i3case3
} // icase2
d2aLm[L1][L1+m1]=d2aLm[L1][L1+m1]/(-r5*r5);
} // icase1

int dns1,dns2,dns3;
// Note: dns1 does not have to be large if there is no
// structure inside.
// Expected to have error less than 1%
//dns1=50; dns2=500; dns3=500; // 1% accuracy
//dns1=1000; dns2=3050; dns3=3000; // expect 0.15% accuracy
dns1=3000; dns2=10000; dns3=10000; //
//dns1=3000; dns2=10000; dns3=10000; // 0.1% accuracy
//dns1=4000; dns2=10000; dns3=10000; // 0.1% accuracy
//dns1=2500; dns2=25000; dns3=25000; // 0.02% accuracy
//if (random(1000)==0)
//{ //
//gotoxy(1,23);cout<<dns1<<" "<<dns2<<" "<<dns3<<"   ";
//gotoxy(1,24);cout<<ans1<<" "<<ans2<<" "<<ans3<<"   ";
//} //
if (r5<rinneredge)
{ // nanoparticle interior
dr5=rinneredge/dns1;
if (r5+dr5>rinneredge) dr5=rinneredge-r5+0.001*(routeredge-rinneredge);
ans1++;
} // nanoparticle interior
else if ((r5>=rinneredge)&&(r5<=routeredge))
{ //
dr5=(routeredge-rinneredge)/dns2;
double drclm=clmr[i9alasttime+1]-clmr[i9alasttime];
if (drclm<dr5) dr5=drclm; // Jul 1 try 0.5 - made no difference
ans2++;
} //
else if (r5>routeredge)
{ // outside nanoparticle in matrix
dr5=routeredge*r5*r5/dns3;
ans3++;
} // outside nanoparticle in matrix

int icase8,L8,m8;
for (icase8=0;icase8<HowManyaLm;icase8++)
{ // icase8
L8=ceL(icase8); m8=cem(icase8);
wLm[L8][L8+m8]=aLm[L8][L8+m8]*pow5(r5,(-L8));
dwLm[L8][L8+m8]=daLm[L8][L8+m8]*pow5(r5,(-L8)) + aLm[L8][L8+m8]*(-L8)*pow5(r5,(-L8)-1);
d2wLm[L8][L8+m8]=d2aLm[L8][L8+m8]*pow5(r5,(-L8))
+2*daLm[L8][L8+m8]*(-L8)*pow5(r5,(-L8)-1)
+aLm[L8][L8+m8]*(-L8)*((-L8)-1)*pow5(r5,(-L8)-2);
} // icase8

for (icase8=0;icase8<HowManyaLm;icase8++)
{ // icase8
L8=ceL(icase8); m8=cem(icase8);
 wLm[L8][L8+m8] = wLm[L8][L8+m8]+dr5*dwLm[L8][L8+m8]+0.5*dr5*dr5*d2wLm[L8][L8+m8];
dwLm[L8][L8+m8] =dwLm[L8][L8+m8]+dr5*d2wLm[L8][L8+m8];
} // icase8

/*
// June 18th:
if (icase7==HowManyaLm)
{ //icase7
int ipx;
ipx=bx33+sfx33*(r5-1.0);
complex a1=daLm[1][1+0]*dqlm(0,0,r5);
putpixel(ipx,by33-0.000022*real(a1),10);
} //icase7
*/

r5=r5+dr5;

for (icase8=0;icase8<HowManyaLm;icase8++)
{ // icase8
L8=ceL(icase8); m8=cem(icase8);
aLm[L8][L8+m8]=wLm[L8][L8+m8]*pow5(r5,L8);
daLm[L8][L8+m8]=dwLm[L8][L8+m8]*pow5(r5,L8) + wLm[L8][L8+m8]*L8*pow5(r5,L8-1);
d2aLm[L8][L8+m8]=d2wLm[L8][L8+m8]*pow5(r5,L8)
+ 2*dwLm[L8][L8+m8]*L8*pow5(r5,L8-1)
+ wLm[L8][L8+m8]*L8*(L8-1)*pow5(r5,L8-2);
} // icase8


if (icase7==HowManyaLm)
{ // integrate dipole moment - Method 1 (polarization)
complex dpmu,curly,round;
int case22,L22,m22; int case23,L23,m23;
for (iwhichaxisforp=-1;iwhichaxisforp<=1;iwhichaxisforp++)
{ // iwhichaxisforp
// Evaluate contents of curly bracket of equation [...]
curly=0.0;
for (case22=0;case22<HowManycLm;case22++)
{ // case22
L22=pcL[case22];m22=pcm[case22];
for (case23=0;case23<HowManyaLm;case23++)
{ // case23
L23=ceL(case23);m23=cem(case23);
complex fac1,fac2,fac3;
fac1=sqrt(4.0*Pi/3.0);
fac2=-qlm(case22,r5);
if (L22==0) fac2=fac2+sqrt(4.0*Pi);
fac3=daLm[L23][L23+m23]*H(L23,m23,1,iwhichaxisforp,L22,m22);
fac3=fac3+(1.0/r5)*aLm[L23][L23+m23]*K(L22,m22,1,iwhichaxisforp,L23,m23);
curly=curly+fac1*fac2*fac3;
// In other words, add this term on just once during summation:
} // case23
} // case22
round=curly+4*Pi*(1.0-eps8out)*Eincmu[1+iwhichaxisforp];
dpmu=round*(3.0/(2*eps8out+1.0))*complex(eps0,0)*r5*r5*dr5;
if (r5<r5maxforp) pmu[1+iwhichaxisforp]=pmu[1+iwhichaxisforp]+dpmu;
} // iwhichaxisforp
putpixel(bx33+sfx33*(r5-1),475-fabs(abs(pmu[1+iwhichaxisforE]))*2.51e8,15);
putpixel(bx33+sfx33*(r5-1),475-fabs(abs(exact_p))*2.51e8,8);
} // integrate dipole moment - Method 1 (polarization)

complex curly,fac1,fac2,fac3;
if (icase7==HowManyaLm)
{ // integrate dipole moment - Method 2 (bound charge density)
int case22,L22,m22; int case23,L23,m23;
for (iwhichaxisforp=-1;iwhichaxisforp<=1;iwhichaxisforp++)
{ // iwhichaxisforp
curly=0;
for (case22=0;case22<HowManycLm;case22++)
{ // case22
L22=pcL[case22];m22=pcm[case22];
for (case23=0;case23<HowManyaLm;case23++)
{ // case23
L23=ceL(case23);m23=cem(case23);
fac1=qlm(case22,r5);
if (L22==0) fac1=fac1-sqrt(4*Pi);
fac2=r5*r5*d2aLm[L23][L23+m23]+2*r5*daLm[L23][L23+m23]-L23*(L23+1)*aLm[L23][L23+m23];
fac3=H(1,iwhichaxisforp,L23,m23,L22,m22);
curly=curly+fac1*fac2*fac3;
curly=curly+r5*r5*dqlm(case22,r5)*daLm[L23][L23+m23]*H(1,iwhichaxisforp,L23,m23,L22,m22);
curly=curly+qlm(case22,r5)*aLm[L23][L23+m23]*K(1,iwhichaxisforp,L23,m23,L22,m22);
} // case23
} // case22
if (r5<r5maxforp) pmu2[1+iwhichaxisforp]+=curly*r5*dr5*sqrt(4*Pi/3.0)*eps0;
} // iwhichaxisforp
putpixel(bx33+sfx33*(r5-1),1+475-fabs(abs(pmu2[1+iwhichaxisforE]))*2.51e8,11);
putpixel(bx33+sfx33*(r5-1),475-fabs(abs(pmu2[1+iwhichaxisforE]))*2.51e8,13);
} // integrate dipole moment - Method 2

if (icase7==HowManyaLm)
{ // integrate dipole moment - Method (3) (potential)
int case22,L22,m22; int case23,L23,m23;
for (iwhichaxisforp=-1;iwhichaxisforp<=1;iwhichaxisforp++)
{ // iwhichaxisforp
complex aavg;
aavg = aLm[1][1+iwhichaxisforp];
curly=r5*r5*r5*d2aLm[1][1+iwhichaxisforp]+2*r5*r5*daLm[1][1+iwhichaxisforp]-2*r5*aavg;
if (r5<r5maxforp) pmu3[1+iwhichaxisforp]-=curly*dr5*sqrt(4*Pi/3.0)*eps0;
} // iwhichaxisforp
putpixel(bx33+sfx33*(r5-1),475-fabs(abs(pmu3[1+iwhichaxisforE]))*2.51e8,14);
complex c53;
c53=pmu3[1+0]*1e9;
if (abs(c53)<0.01) c53=0;
if ((licount%30)==0){gotoxy(1,17);cout<<"(3) pz = "<<c53<<" nC-m ";cout.precision(2);cout<<100*abs((pmu3[1+0]-exact_p)/exact_p)<<"% ";}
} // integrate dipole moment - Method (3)

int icase4; int iL4,im4;
for (icase4=0;icase4<HowManyaLm;icase4++)
{ // icase4
iL4 = ceL(icase4); im4 = cem(icase4);
complex a1;
if (3+icase4<=25)
{ // ..25
a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
if (icase7<HowManyaLm) Fmat[icase4][icase7] = a1;

if (icase7==HowManyaLm)
{ //...
int ipx,iL4,im4;
ipx=bx33+sfx33*(r5-1.0);
putpixel(ipx, 2+by33,8);
putpixel(ipx,-2+by33,8);
putpixel(ipx, 1+by33,8);
putpixel(ipx,-1+by33,8);
putpixel(ipx,   by33,8);
iL4=1;im4=0; a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
putpixel(ipx,by33+0.00871*real(a1),14);
iL4=3;im4=0; a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
putpixel(ipx,by33+.5871*real(a1),13);
iL4=5;im4=0; a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
putpixel(ipx,by33+.5871*real(a1),12);
iL4=7;im4=0; a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
putpixel(ipx,by33+.5871*real(a1),11);
} //...

if ((licount%130)==0)
{ // ...licount%10)==0
iL4=ceL(icase4); im4=cem(icase4);
a1=aLm[iL4][iL4+im4]/pow5(r5,iL4);
if (abs(a1)<0.001) a1=complex(0.0,0);
gotoxy(7+8*icase7,4+icase4);cout<<" "<<real(a1)<<"        ";
} // ...

if (((licount%30)==0)&&(icase7==0))
{ // ...licount%10)==0
int ipx;
ipx=bx33+sfx33*(r5-1.0);
gotoxy(2,4+icase4);cout<<iL4<<" "<<im4<<" ";
//putpixel(ipx,by33-87.1*real(a1),14-icase4);
//if ((ipx%3)==0) putpixel(ipx,2+by33-87.1*real(a1),14);
//putpixel(ipx,by33-87.1*imag(a1),14-icase4);
//if ((ipx%3)==0) putpixel(ipx,2+by33-87.1*imag(a1),13);
//if (abs(a1)<0.001) a1=complex(0.0,0);
//gotoxy(7+8*icase7,4+icase4);cout<<real(a1)<<"        ";

//a1=daLm[iL4][iL4+im4]/(pow5(r5,iL4-1)*iL4);
//putpixel(ipx,by33-37.1*real(a1),14-icase4);
//if ((ipx%3)==0) putpixel(ipx,2+by33-37.1*real(a1),11);
//putpixel(ipx,by33-37.1*imag(a1),14-icase4);
//if ((ipx%3)==0) putpixel(ipx,2+by33-37.1*imag(a1),11);
} // ...licount%10)==0


//putpixel(ipx,1+by33-87.1*real(a1),1+icase4);
//putpixel(ipx,1+by33-87.1*imag(a1),3+icase4);
//putpixel(ipx,2+by33-87.1*real(a1),1+icase4);
//putpixel(ipx,2+by33-87.1*imag(a1),3+icase4);

} // ..25
} // icase4

if ((licount%100)==0)
{ // ...licount%10)==0
cout.precision(3);
gotoxy(60,10);cout<<"r = "<<r5<<"  ";
gotoxy(2,4-1);cout<<"L m  eLm ";
} // ...licount%10)==0

if (r5>r5max)
{ // done column
if (icase7==HowManyaLm-1) goto LBinvertc;
if (icase7==HowManyaLm) goto LBthatsall;
goto LBnextcol;
} // done column

if (!kbhit()) goto LBnextr;
gotoxy(1,25);cout<<"* interrupted since you hit a key *";
getch();
return;

LBinvertc:

for (irow2=0;irow2<HowManyaLm;irow2++)
{ // irow2
Fmat[irow2][HowManyaLm]=0;
} // irow2
Fmat[0][HowManyaLm]=Eincmu[1+iwhichaxisforE]/(-0.488603);

// Display matrix 'Fmat':
cout.precision(3);
int icol2;
for (icol2=0;icol2<=HowManyaLm;icol2++)
{ // icol2
for (irow2=0;irow2<HowManyaLm;irow2++)
{ // irow2
//gotoxy(1+16*icol2,11+irow2);cout<<Fmat[irow2][icol2]<<"   ";
} // irow2
} // icol2

// Solve the linear system 'Fmat' of equations:
for (ic2=0;ic2<HowManyaLm;ic2++)
{ // ic2
complex r2;
r2=Fmat[ic2][ic2];
if (r2==0) r2=1;
for (ic3=0;ic3<=HowManyaLm;ic3++)
{ // ic3
Fmat[ic2][ic3]=Fmat[ic2][ic3]/r2;
} // ic3
for (ir2=0;ir2<HowManyaLm;ir2++)
{ //..
if (ir2!=ic2)
{ //.
r2=Fmat[ir2][ic2]/Fmat[ic2][ic2];
for (ic3=0;ic3<=HowManyaLm;ic3++)
{ // ic3
Fmat[ir2][ic3]=Fmat[ir2][ic3]-Fmat[ic2][ic3]*r2;
} // ic3
} //.
} //..
} // ic2

// Obtain dLm (origin RSH coefficients) from solved system:
for (irow2=0;irow2<HowManyaLm;irow2++)
{ // irow2
int L2,m2;
L2=ceL(irow2);m2=cem(irow2);
dLm[L2][L2+m2]=Fmat[irow2][HowManyaLm];
} // irow2

/*
// Display the solved system:
cout.precision(3);
for (icol2=HowManyaLm;icol2<=HowManyaLm;icol2++)
{ // icol2
for (irow2=0;irow2<HowManyaLm;irow2++)
{ // irow2
complex a1;
a1=Fmat[irow2][icol2];
if (abs(a1)<0.001) a1=complex(0,0);
gotoxy(2+6*icol2,17+irow2);cout<<real(a1)<<"   ";
} // irow2
} // icol2
*/

gotoxy(2,22);
Ecenter=-Fmat[0][HowManyaLm]*0.488603;
iEr = real(Ecenter);
iEi = imag(Ecenter);
if (iwhichaxisforE==0)  { cout<<"center Ez = ";}
if (iwhichaxisforE==1)  { cout<<"center Ex = ";}
if (iwhichaxisforE==-1) { cout<<"center Ey = ";}
if (iEi>0) cout<<iEr<<"+"<<iEi<<"i V/m";
if (iEi<0) cout<<iEr<<iEi<<"i V/m";
if (iEi==0) cout<<iEr<<" V/m";
cout.precision(2);
cout<<" ["<<100.0*abs((Ecenter-Eellipsoid)/Eellipsoid)<<" %]";

// The right-most column are the dLm values.
goto LBnextcol;

LBthatsall:
cout.precision(3);
complex c53;
c53=pmu3[1+1]*1e9;
if (abs(c53)<0.01) c53=0;
gotoxy(1,16);cout<<"(3) px = "<<c53<<" nC-m ";
//c53=pmu3[1+0]*1e9;
//if (abs(c53)<0.01) c53=0;
//gotoxy(1,17);cout<<"(3) pz = "<<c53<<" nC-m ";
c53=pmu2[1+0]*1e9;
if (abs(c53)<0.01) c53=0;
gotoxy(1,18);cout<<"(2) pz = "<<c53<<" nC-m ";
c53=pmu[1+1]*1e9;
if (abs(c53)<0.01) c53=0;
gotoxy(1,19);cout<<"(1) px = "<<c53<<" nC-m ";
c53= pmu[1+0]*1e9;
if (abs(c53)<0.01) c53=0;
gotoxy(1,20);cout<<"(1) pz = "<<c53<<" nC-m ";
gotoxy(1,25);cout<<"finished";
getch();
return;
} // main

double new_PLm(int L,int m,double x)
{ //   new_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=pow7(2.0,L,3);
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:
if (am>0) PLm1 = PLm1 * pow7(s,am,6);

} // 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;
} //  new_PLm   associated Legendre function


double fact(int i3)
{ //  fact
if (i3<0) { cout<<"* fact negative factorial *"; exit(0); }
double fact1;
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 (m superscript)
complex YLm1;
YLm1 = (2.0*L+1.0)/(4.0*Pi);
YLm1=YLm1 * fact(L-m) / fact(L+m);
YLm1=sqrt(YLm1);
YLm1=YLm1 * new_PLm(L,m,cos(theta));
YLm1=YLm1 * exp(complex(0.0,m*phi));
return YLm1;
} //   YLm (m superscript)

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<<"** 7 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 = pow7(x,y,9);
return pow51;
} //   pow5

double pow7(double x,double y,int from)
{ //   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<<"* pow7(): 0 to negative power y = "<<y<<" from"<<from<<"*";exit(0);}
count0to0++;
gotoxy(1,25);cout<<"## pow7: 0^0 given as 0 from "<<from<<" ("<<count0to0<<" times) ##";
return 0.0;
} //
else if (x<0.0)
{ //
if (y!=ceil(y)) {cout<<"* pow7(): negative to noninteger: "<<x<<"^"<<y<<" from"<<from<<"*";exit(0);}
} //
pow51 = pow(x,y);
return pow51;
} //   pow7

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)
{ //.
SLm1 = real(p707*(YLm(L,m,theta,phi)+pow7(-1.0,m,12)*YLm(L,-m,theta,phi)));
} //.
else if (m<0)
{ //..
SLm1 = real(p707*complex(0,1.0)*(YLm(L,m,theta,phi)-pow7(-1.0,m,15)*YLm(L,-m,theta,phi)));
} //..
return SLm1;
} // SLm   Real Spherical Harmonics

char readclmline(void)
{ // readclmline
int i3=0;
unsigned char ch1,ch1f;
ch1f=255;
LB111:
i3++;
if (i3>320) {gotoxy(1,10);cout<<"** Error: readclmline line too long **";exit(0);}
fscanf(clmfin,"%c",&ch1);
if (feof(clmfin)) {gotoxy(1,10);cout<<"* readclmline end of file *";exit(0);}
if (ch1f==255) ch1f=ch1;
if (ch1!=10) goto LB111;
return ch1f;
} // readclmline

complex clm(int ipc,double r)
{ //    clm
int i9a,i9b;
complex clm1=0;
double dr;
i9a=i9alasttime;
i9b=i9alasttime+1;

if (r<clmr[i9a])
{ // r outside low now
if (r<=clmr[0]) return clmc[ipc][0];
LB395:
if (r<clmr[i9a]) {i9a--;goto LB395;}
i9b=i9a+1;
i9alasttime=i9a;
} // r outside low now
else if (r>clmr[i9b])
{ // r outside high now
if (r>=clmr[nclm-1]) return clmc[ipc][nclm-1];
LB396:
if (r>clmr[i9b]) {i9b++;goto LB396;}
i9a=i9b-1;
i9alasttime=i9a;
} // r outside high now
/*
if (i9alasttime<0) {cout<<"* error 9834 *";exit(0);}
if (i9alasttime>nclm-2) {cout<<"* error 9835 *";exit(0);}
if (i9b!=i9a+1) {cout<<"* Error 9885 r="<<r<<" i9a="<<i9a<<" i9b="<<i9b<<"  *";exit(0);}
if (i9a<0) {cout<<"* error 9887 *";exit(0);}
if (i9b>nclm-1) {cout<<"* error 9888 *";exit(0);}
*/
dr=clmr[i9b]-clmr[i9a];
return ((clmr[i9b]-r)/dr)*clmc[ipc][i9a]+((r-clmr[i9a])/dr)*clmc[ipc][i9b];
} //   clm

complex qlm(int ipc,double r)
{ //    qlm
int i9a,i9b;
complex clm1=0;
double dr;
i9a=i9alasttime;
i9b=i9alasttime+1;

if (r<clmr[i9a])
{ // r outside low now
if (r<=clmr[0]) return clmc[ipc][0];
LB395:
if (r<clmr[i9a]) {i9a--;goto LB395;}
i9b=i9a+1;
i9alasttime=i9a;
} // r outside low now
else if (r>clmr[i9b])
{ // r outside high now
if (r>=clmr[nclm-1]) return clmc[ipc][nclm-1];
LB396:
if (r>clmr[i9b]) {i9b++;goto LB396;}
i9a=i9b-1;
i9alasttime=i9a;
} // r outside high now
/*
if (i9alasttime<0) {cout<<"* error 1834 *";exit(0);}
if (i9alasttime>nclm-2) {cout<<"* error 1835 *";exit(0);}
if (i9b!=i9a+1) {cout<<"* Error 1885 r="<<r<<" i9a="<<i9a<<" i9b="<<i9b<<"  *";exit(0);}
if (i9a<0) {cout<<"* error 1887 *";exit(0);}
if (i9b>nclm-1) {cout<<"* error 1888 *";exit(0);}
*/
dr=clmr[i9b]-clmr[i9a];
return ((clmr[i9b]-r)/dr)*clmq[ipc][i9a]+((r-clmr[i9a])/dr)*clmq[ipc][i9b];
} //   qlm

complex dclm(int ipc,double r)
{ //    dclm
int i9a,i9b;
complex clm1=0;
double dr;
i9a=i9alasttime;
i9b=i9alasttime+1;

if (r<clmr[i9a])
{ // r outside low now
if (r<=clmr[0]) return 0;
LB395:
if (r<clmr[i9a]) {i9a--;goto LB395;}
i9b=i9a+1;
i9alasttime=i9a;
} // r outside low now
else if (r>clmr[i9b])
{ // r outside high now
if (r>=clmr[nclm-1]) return 0;
LB396:
if (r>clmr[i9b]) {i9b++;goto LB396;}
i9a=i9b-1;
i9alasttime=i9a;
} // r outside high now
/*
if (i9alasttime<0) {cout<<"* error 29834 *";exit(0);}
if (i9alasttime>nclm-2) {cout<<"* error 29835 *";exit(0);}
if (i9b!=i9a+1) {cout<<"* Error 29885 r="<<r<<" i9a="<<i9a<<" i9b="<<i9b<<"  *";exit(0);}
if (i9a<0) {cout<<"* error 29887 *";exit(0);}
if (i9b>nclm-1) {cout<<"* error 29888 *";exit(0);}
*/
return (clmc[ipc][i9b]-clmc[ipc][i9a])/(clmr[i9b]-clmr[i9a]);
} //   dclm

complex dqlm(int ipc,double r)
{ //    dqlm
int i9a,i9b;
complex clm1=0;
double dr;
i9a=i9alasttime;
i9b=i9alasttime+1;

if (r<clmr[i9a])
{ // r outside low now
if (r<=clmr[0]) return 0;
LB395:
if (r<clmr[i9a]) {i9a--;goto LB395;}
i9b=i9a+1;
i9alasttime=i9a;
} // r outside low now
else if (r>clmr[i9b])
{ // r outside high now
if (r>=clmr[nclm-1]) return 0;
LB396:
if (r>clmr[i9b]) {i9b++;goto LB396;}
i9a=i9b-1;
i9alasttime=i9a;
} // r outside high now
/*
if (i9alasttime<0) {cout<<"* error 39834 *";exit(0);}
if (i9alasttime>nclm-2) {cout<<"* error 39835 *";exit(0);}
if (i9b!=i9a+1) {cout<<"* Error 39885 r="<<r<<" i9a="<<i9a<<" i9b="<<i9b<<"  *";exit(0);}
if (i9a<0) {cout<<"* error 39887 *";exit(0);}
if (i9b>nclm-1) {cout<<"* error 39888 *";exit(0);}
*/
return (clmq[ipc][i9b]-clmq[ipc][i9a])/(clmr[i9b]-clmr[i9a]);
} //   dqlm


double H(int L1,int m1,int L2,int m2,int L3,int m3)
{ //   H(L1 m1;L2 m2;L3 m3)
double H1v=-1;
if ((L1==L2)&&(m1==m2)&&(L3==0)&&(m3==0)) return 1.0/sqrt(4.0*Pi);
if (((L1!=L2)||(m1!=m2))&&(L3==0)&&(m3==0)) return 0.0;
int i56;

for (i56=0;i56<hlutln;i56++)
{ // i56
int flag5=1;
if (L1!=hlutl1[i56]) flag5=0;
if (m1!=hlutm1[i56]) flag5=0;
if (L2!=hlutl2[i56]) flag5=0;
if (m2!=hlutm2[i56]) flag5=0;
if (L3!=hlutl3[i56]) flag5=0;
if (m3!=hlutm3[i56]) flag5=0;
if (flag5==1) return hlutva[i56];
} // i56

// Here is where we put in the efficient integration algorithm:
double theta0,phi0,dtheta,dphi,maxval,newval,theta1,phi1,int_th,int_phi,Hval,nodivis;
nodivis=300.0; // PPP Numerical experiments show very good convergence
//very quickly as the number of divisions increases.
maxval=0.0;

for (i56=0;i56<20;i56++)
{ // i56
theta1=Pi*0.0001*random(10000);
phi1=2.0*Pi*0.0001*random(10000);
newval=SLm(L1,m1,theta1,phi1)*SLm(L2,m2,theta1,phi1)*SLm(L3,m3,theta1,phi1)*sin(theta1);
if (fabs(newval)>maxval) {maxval=fabs(newval);theta0=theta1;phi0=phi1;}
} // i56

Hval=0.0;
if (maxval!=0)
{ // maxval!=0
int_th=0;
dtheta= Pi/nodivis;
dphi= 2.0*Pi/nodivis;
for (theta1=0.5*dtheta;theta1<Pi;theta1=theta1+dtheta)
{ //
int_th=int_th+dtheta*
SLm(L1,m1,theta1,phi0)*SLm(L2,m2,theta1,phi0)*SLm(L3,m3,theta1,phi0)*sin(theta1);
} //
int_phi=0;
for (phi1=0.5*dphi;phi1<2.0*Pi;phi1=phi1+dphi)
{ //
int_phi=int_phi+dphi*
SLm(L1,m1,theta0,phi1)*SLm(L2,m2,theta0,phi1)*SLm(L3,m3,theta0,phi1)*sin(theta0);
} //
double den1;
den1=SLm(L1,m1,theta0,phi0)*SLm(L2,m2,theta0,phi0)*SLm(L3,m3,theta0,phi0)*sin(theta0);
if (den1!=0)
{
Hval=int_th*int_phi/den1;
}
else
{ //
cout<<"* Error H den1=0 abcd theta0="<<theta0<<" *";exit(0);
} //
} // maxval!=0

if (hlutln<KLUTMAX)
{ // add to look up table
hlutl1[hlutln]=L1;
hlutm1[hlutln]=m1;
hlutl2[hlutln]=L2;
hlutm2[hlutln]=m2;
hlutl3[hlutln]=L3;
hlutm3[hlutln]=m3;
hlutva[hlutln]=Hval;
hlutln++;
gotoxy(40,23);cout<<"hlutln "<<hlutln<<" ";
} // add to look up table
else if (random(10)==0)
{ //..
gotoxy(40,24);cout<<"* H look up table full *";
gotoxy(40,25);cout<<"* max table size = "<<KLUTMAX<<" *";
} //..
return Hval;
} //   H(L1 m1;L2 m2;L3 m3)

double K(int L1,int m1,int L2,int m2,int L3,int m3)
{ //   K(L1 m1|L2 m2;L3 m3)
double theta0,phi0,dtheta,dphi,maxval,newval,theta1,phi1,int_th,int_phi;
double Kval1,Kval2,Kval,nodivis;
double den1;
int i56;
if ((L3==0)&&(m3==0)) return 0.0;

// Next check in look up table of already calculated values:
for (i56=0;i56<klutln;i56++)
{ // i56
int flag5=1;
if (L1!=klutl1[i56]) flag5=0;
if (m1!=klutm1[i56]) flag5=0;
if (L2!=klutl2[i56]) flag5=0;
if (m2!=klutm2[i56]) flag5=0;
if (L3!=klutl3[i56]) flag5=0;
if (m3!=klutm3[i56]) flag5=0;
if (flag5==1) return klutva[i56];
} // i56

nodivis=300.0; // PPP Numerical experiments show the convergence is good

maxval=0.0;
for (i56=0;i56<20;i56++)
{ // i56
theta1=Pi*0.0001*(1+random(10000-1)); // excludes endpoints
phi1=2.0*Pi*0.0001*random(10000);
newval=SLm(L1,m1,theta1,phi1)*dSLm_theta(L2,m2,theta1,phi1)*dSLm_theta(L3,m3,theta1,phi1)*sin(theta1);
if (fabs(newval)>maxval) {maxval=fabs(newval);theta0=theta1;phi0=phi1;}
} // i56
Kval1=0.0;
if (maxval!=0)
{ // maxval!=0
int_th=0;
dtheta= Pi/nodivis;
dphi= 2.0*Pi/nodivis;
for (theta1=0.5*dtheta;theta1<Pi;theta1=theta1+dtheta)
{ //
int_th=int_th+dtheta*
SLm(L1,m1,theta1,phi0)*dSLm_theta(L2,m2,theta1,phi0)*dSLm_theta(L3,m3,theta1,phi0)*sin(theta1);
} //
int_phi=0;
for (phi1=0.5*dphi;phi1<2.0*Pi;phi1=phi1+dphi)
{ //
int_phi=int_phi+dphi*
SLm(L1,m1,theta0,phi1)*dSLm_theta(L2,m2,theta0,phi1)*dSLm_theta(L3,m3,theta0,phi1)*sin(theta0);
} //
den1=SLm(L1,m1,theta0,phi0)*dSLm_theta(L2,m2,theta0,phi0)*dSLm_theta(L3,m3,theta0,phi0)*sin(theta0);
if (den1!=0)
{
Kval1=int_th*int_phi/den1;
}
else
{ //
cout<<"* Error K den1=0 abcd theta0="<<theta0<<" *";exit(0);
} //
} // maxval!=0

maxval=0.0;
for (i56=0;i56<20;i56++)
{ // i56
theta1=Pi*0.0001*(1+random(10000-1)); // excludes endpoints
phi1=2.0*Pi*0.0001*random(10000);
den1=sin(theta1);
if (den1!=0)
{ //
newval=SLm(L1,m1,theta1,phi1)*dSLm_phi(L2,m2,theta1,phi1)*dSLm_phi(L3,m3,theta1,phi1)/den1;
} //
else
{ //
cout<<"* Error K den1=0 a2 theta0="<<theta0<<" *";exit(0);
} //
if (fabs(newval)>maxval) {maxval=fabs(newval);theta0=theta1;phi0=phi1;}
} // i56

Kval2=0.0;
if (maxval!=0)
{ // maxval!=0
int_th=0;
dtheta= Pi/nodivis;
dphi= 2.0*Pi/nodivis;
for (theta1=0.5*dtheta;theta1<Pi;theta1=theta1+dtheta)
{ //
den1=sin(theta1);
if (den1!=0)
{
int_th=int_th+dtheta*
SLm(L1,m1,theta1,phi0)*dSLm_phi(L2,m2,theta1,phi0)*dSLm_phi(L3,m3,theta1,phi0)/den1;
}
else
{ //
cout<<"* Error K den1=0 b3 theta1="<<theta1<<" *";exit(0);
} //
} //

int_phi=0;
for (phi1=0.5*dphi;phi1<2.0*Pi;phi1=phi1+dphi)
{ //
den1=sin(theta0);
if (den1!=0)
{
int_phi=int_phi+dphi*
SLm(L1,m1,theta0,phi1)*dSLm_phi(L2,m2,theta0,phi1)*dSLm_phi(L3,m3,theta0,phi1)/den1;
}
else
{ //
cout<<"* Error K den1=0 c4 theta0="<<theta0<<" *";exit(0);
} //
} //
double z1,z2,z3;
z1=int_th*int_phi;
den1=SLm(L1,m1,theta0,phi0)*dSLm_phi(L2,m2,theta0,phi0)*dSLm_phi(L3,m3,theta0,phi0);
if (den1!=0)
{ //
Kval2=z1*sin(theta0)/den1;
} //
else
{ //.
cout<<"* Error K den1=0 d5 *";
cout<<"* L1="<<L1<<" m1="<<m1<<" L2="<<L2<<" m2="<<m2<<" L3="<<L3<<" m3="<<m3<<" *";
cout<<"* theta0="<<theta0<<" phi0="<<phi0<<" *\n";
cout<<"* z1="<<z1<<" *\n";
cout<<"* sin(theta0)="<<sin(theta0)<<" *\n";
cout<<"* maxval= "<<maxval<<" *";
exit(0);
} //.
} // maxval!=0

Kval=Kval1+Kval2;

if (klutln<KLUTMAX)
{ // add to look up table
klutl1[klutln]=L1;
klutm1[klutln]=m1;
klutl2[klutln]=L2;
klutm2[klutln]=m2;
klutl3[klutln]=L3;
klutm3[klutln]=m3;
klutva[klutln]=Kval;
klutln++;
gotoxy(40,24);cout<<"klutln "<<klutln<<" ";
} // add to look up table
else if (random(20)==0)
{ //..
gotoxy(40,24);cout<<"* K look up table full *";
gotoxy(40,25);cout<<"* max table size = "<<KLUTMAX<<" *";
} //..

return Kval;
} //   K(L1 m1|L2 m2;L3 m3)

double dSLm_theta(int L,int im,double theta,double phi)
{ //   dSLm_theta    Real Spherical Harmonics
double SLm1,p707; p707 = 1.0/sqrt(2.0);
if (im==0)
{ //,
SLm1 = real(dYLm_theta(L,im,theta,phi));
} //,
else if (im>0)
{ //.
SLm1 = real(p707*(dYLm_theta(L,im,theta,phi)+pow7(-1.0,im,18)*dYLm_theta(L,-im,theta,phi)));
} //.
else if (im<0)
{ //..
SLm1 = real(p707*complex(0,1.0)*(dYLm_theta(L,im,theta,phi)-pow7(-1.0,im,21)*dYLm_theta(L,-im,theta,phi)));
} //..
return SLm1;
} //  dSLm_theta   Real Spherical Harmonics

double dSLm_phi(int L,int im,double theta,double phi)
{ //   dSLm_phi    Real Spherical Harmonics
double SLm1,p707; p707 = 1.0/sqrt(2.0);
if (im==0)
{ //,
SLm1 = real(complex(0,im)*YLm(L,im,theta,phi));
} //,
else if (im>0)
{ //.
SLm1 = real(p707*(complex(0,im)*YLm(L,im,theta,phi)+pow7(-1.0,im,24)*complex(0,-im)*YLm(L,-im,theta,phi)));
} //.
else if (im<0)
{ //..
SLm1 = real(p707*complex(0,1.0)*(complex(0,im)*YLm(L,im,theta,phi)-pow7(-1.0,im,27)*complex(0,-im)*YLm(L,-im,theta,phi)));
} //..
return SLm1;
} //  dSLm_phi   Real Spherical Harmonics

complex dYLm_theta(int L,int m,double theta,double phi)
{ //    dYLm_theta
// June 13th: upgraded for singular cases: ( L=0, catch L=|m| )
complex YLm1;
if (L==0) return 0.0;
YLm1 = (2.0*L+1.0)/(4.0*Pi);
YLm1=YLm1 * fact(L-m) / fact(L+m);
YLm1=sqrt(YLm1);
YLm1=YLm1 * dPLm_theta(L,m,theta);
YLm1=YLm1 * exp(complex(0.0,m*phi));
return YLm1;
} //   dYLm_theta

double dPLm_theta(int L,int im,double theta)
{ //   dPLm_theta
// For |m|<L this is based on equation (40) in
// http://mathworld.wolfram.com/LegendrePolynomial.html
double dP1,mu;
mu = cos(theta);
if ((fabs(theta)<=1e-6)||(fabs(Pi-theta)<1e-6)) return 0;
if (L==0) return 0;
if (abs(im)<L)
{ //.. |m|<L
dP1 = (L*mu*new_PLm(L,im,mu)-(L+im)*new_PLm(L-1,im,mu))/sqrt(1.0-mu*mu);
} //.. |m|<L
else
{ //.. |m|=L cases
if (im==L)
 dP1=( fact(2*L-1) / (pow7(2.0,L-1,30)*fact(L-1)) )*L*pow5(sin(theta),L-1)*cos(theta);
if (im==-L)
 dP1=(pow7(-0.5,L,33)/fact(L))*L*pow5(sin(theta),L-1)*cos(theta);
} //.. |m|=L cases
return dP1;
} //   dPLm_theta

int ceL(int icase)
{ //ceL
if (icase==0) return 1;
if (icase==1) return 3;
if (icase==2) return 5;
if (icase==3) return 7;
if (icase==4) return 9;
if (icase==5) return 11;
if (icase==6) return 13;
if (icase==7) return 15;
if (icase==8) return 17;
if (icase==9) return 19;
if (icase==10) return 21;
cout<<"* Error ceL *";exit(0);
} //ceL

int cem(int icase)
{ //cem
if (icase==0) return iwhichaxisforE;
if (icase==1) return iwhichaxisforE;
if (icase==2) return iwhichaxisforE;
if (icase==3) return iwhichaxisforE;
if (icase==4) return iwhichaxisforE;
if (icase==5) return iwhichaxisforE;
if (icase==6) return iwhichaxisforE;
if (icase==7) return iwhichaxisforE;
if (icase==8) return iwhichaxisforE;
if (icase==9) return iwhichaxisforE;
if (icase==10) return iwhichaxisforE;
cout<<"* Error cem *";exit(0);
} //cem