// scp77c Dec 26  toroidal
// scp73c.cpp Dec 25  precise root finding (just like scp76c.cpp)
// scp73.cpp Dec 24  Fixed css[][] error!
//  scp70k7.cpp Dec 23 L=0 case   HAS ERROR [corrected in scp73]
//   k6g   L=2 Dec 22
//   k6f   S = beta(vl/(vt*3.14159)    real S agrees for L=2
//     L=1 works too, and results seem to agree with Saviot's.
//    e   tried L=1, but I also get numerical noise
//   I now think that the imaginary part is positive, not negative.
//   There seems to be a need to normalize the determinant in
//   some way, relative to the matrix element magnitude, I guess.
//    70k6d   betai too.   seems to work, but zero is awfully wide!
//  Also, I find that numerical noise appears when the imaginary
//  part gets larger, so I can't clearly see the zero.  But
//  all variables are double.
//    70k6c   abs det shown along real axis only.  seems to work
//  70k6b  Dec 22   kwvr kwvi    seems to work for L=2 spheroidal
// scp70k6  Dec 22  complex freq, determ... complex k
//   - first, switched back to two materials only.
//double thick3=0.150; // Ni0.75Ag0.25 Dmax 2.3 nm case
//double thick3=0.404; // for Ni0.50Ag0.50 Dmax 3 nm cas
// 70k5  Ag shell with Ni core
// 70k4  dec 16   special for Ag shell portales.pdf
//  Dmax=3nm  Ag0.5Ni1-0.5   partrad=1.096 nm   thickness=0.404
// 70k3  Dec 11    try to plot 3 displacement amplitudes
// scp70k.cpp Dec 10 orderL=2 spheroidal -- figures for md24.htm
//  j  incident longit
// scp70i Dec 9 2002   Used to make figures for md23.htm
// scp70h.cpp  sees to work, three materials  L=1 only.
//  70g  Generalized linear equation solution.  wavenummax
//  70f
//  70e minor changes.  works.
//  70b   A,B --> B,C   A is now incoming amplitude. still works
//  70  new factor arrangement in order to support four variables
//   per extra layer.  Seems to work correctly now.
//  scp67q  CdS in GeO2
//  - inspired by Tanaka Onari and Arai Phys Rev B 47 1993 p 1237
double Btenfactor=1;
// scp67p Dec 8  works fine
//  67n  cleanup() now much faster
//  67m  corrected i234 range error in cleanup()
//  67j  cleanup, works  - most time is spent in double loop in "cleanup".
//  67i more cleanup - still works
//  67h multiple factors in multiplybypow()
//  67g  multiplyby( ,"B") divideby()
//  67f  try again on cleanup (67e has bug) - still works okay
//  67d  Dec 7   demonstrate energy conservation of outgoing vs ingoing
//   67c  error add 0
//   67b  works okay [I don't know why version 67 stopped working]
// scp67.cpp  Dec 6  Incident SPH can be LONG or TRAN
//    - here is an example of the L=1 spheroidal peak (translational)
//  66  Dec 6  For L=2, for nu=0.25, check that peaks are at right
//      location.  Incident wave is spheroidal transverse...
//  65  Do not allow overwriting of an function that is in use
//      (implemented for some functions only)
//  64  Dec 6  L=1  same material.  B not constant, A not zero.
//    If D=1 and B=2 then it should match.
//    I found the specific reason for the problem (double use of a var)
//    It now works correctly for L=1 and L=2.
//  63  increased number of terms.  Result is now very different.
//     - so I think the previous problem was partly due to
//       variables having too few terms, and the program not catching it.
//     SPH Transverse amplitude is now approximately constant, and
//     the SPH LONG amplitude is much smaller.  But this does not
//     seem exactly right.  The value of B is very close to 2.
//  62  ngetm(,)  nsetm(,,)   still works the same
// http://mathworld.wolfram.com/
// gives general formula:
// jn(x) = (-1)^n x^n (d / xdx)^n (sin(x)/x)
// In particular:  j0(x) = sin(x)/x

//  New scheme:
// nm[0][]   sin(theta)
// nm[1][]   cos(theta)
// nm[2 was 6][]   sin(kLp r)
// nm[3 was 7][]   A
// nm[4 was 8][]   cos(kLp r)
// nm[5][]   kLp (longitudinal speed of sound in particle)
// nm[6 was 11][]  B
// nm[7 was 3][]   r
// nm[8 was 10][]  sin(kTp r)
// nm[9 was 15][]  C
// nm[10 was 12][]  cos(kTp r)
// nm[11 was 9][]   kTp (transverse speed of sound in particle)
// nm[12 was 19][]  D
//    13  unused
// nm[14][]  sin(kLm r)
// nm[15 new][]         E
// nm[16][]  cos(kLm r)
// nm[17 was 13][]  kLm (longitudinal speed of sound in matrix)
// nm[18 new][]         F
//    19  unused
// nm[20 was 18][]  sin(kTm r)
// nm[21 new][]         G
// nm[22 was 20][]  cos(kTm r)
// nm[23 was 17][]  kTm (transverse speed of sound in matrix)
// nm[24 new][]         H

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
#include <dos.h>
#include <complex.h>
#define NFUNC 28
#define NTRM 4932
#define NMATER 2
#define NFACT 1+12*NMATER

#define ntrm0 822
#define ntrm1 822
#define ntrm2 822
#define ntrm3 822
#define ntrm4 822
#define ntrm5 822

double far cr[NTRM],ci[NTRM];
char far nm0[NFACT][ntrm0];
char far nm1[NFACT][ntrm1];
char far nm2[NFACT][ntrm2];
char far nm3[NFACT][ntrm3];
char far nm4[NFACT][ntrm4];
char far nm5[NFACT][ntrm4];

int far bt[NFUNC]; // base term
int far nt[NFUNC]; // number of terms
int far mnt[NFUNC]; // max. number of terms
void dispterm(int it2,int ic,int ir);
void ddr(int func2,int func1); //  Take Deriv of func1 with respect to r
void ddth(int func2,int func1); // Take Deriv of func1 with respect to theta
void ddph(int func2,int func1); // Take Deriv of func1 with respect to phi
int irow2;
int ic2;
void gradient(int ior2,int iot2,int iop2,int func);
void curl(int ior2,int iot2,int iop2,int ior1,int iot1,int iop1);
void divergence(int func2,int ior1,int iot1,int iop1);
void cleanup(int func);
void multiply(int func3,int func2,int func1);
void add(int func3,int func2,int func1);
void subtract(int func3,int func2,int func1);
void legendre(int func1,int ordern);
void sphbessel(int func1,int orderL,int imater,int ipolar,int wavetype);
void copyfunc(int func2,int func1);
void stress(int sigrr,int sigrth,int sigrph,int ur,int uth,int uph,int imater);
void multiplyby(int func,double r1);
void multiplyby(int func,complex cr1);
void multiplyby(int func,char *str1);
void divideby(int func,char *str1);
void multiplybypow(int func,char *str1,int ipow);
void settozero(int func);
int ngetm(int i1,int i2);
void nsetm(int i1,int i2,int i3);
double rho[NMATER];
double css[2][NMATER];
double betar,partrad_nm,betai;
double frequencyomegar,frequencyomegai;
complex frequencyomegac;
double kwvr[2][NMATER];
double kwvi[2][NMATER];
void evaluate_func(int func2,int func1,double partrad_nm);
//complex constantcoef(int func1);
int mosttermseveradd=0;
int mosttermseverddr=0;
int mosttermsevercle=0;
void par(char *ch1);
void mat(char *ch1);
void mt3(char *ch1);
complex coefof(int func1,char *str1);
complex coefof(int func1,char ch1);
long int tm1,tm2,tm3,tm4;
void makescr(void);
complex varval[12];
complex cdet(void);
complex y[5][5];
complex y3[5][5];
#define Ny 2

void main(void)
{ // main
clrscr();
randomize();

int bx0=24,by0=480-40;
int i77;
double betamax,wavenummax=171*1.6;
double sfx2=(300-bx0-0.5*bx0)/wavenummax;
double sfx3;

double cnms=1.0e-9*2.998e8; // speed of light in nm/s

//rho[0]=4.82; css[0][0]=4.25; css[1][0]=1.86; par("CdS"); // CdS
//rho[0]=5.81; css[0][0]=3.57; css[1][0]=1.54; par("CdSe"); // CdSe
// Saviot et al. for CdSe: vl=3.7e5 cm/s vt=1.54e5 cm/s
//rho[0]=10; css[0][0]=3.00; css[1][0]=1.732; par("nu=0.25 mat'l.");
//rho[0]=10; css[0][0]=3.65; css[1][0]=1.66; par("Ag");
//rho[0]=0.125 ; css[0][0]=1; css[1][0]=0.3; par("very soft mat'l.");
//rho[0]=8.908; css[0][0]=5.72; css[1][0]=3.00; par("Ni"); // nu=0.31
//rho[0]=5.33; css[0][0]=5.25; css[1][0]=3.35; par("Ge");
//rho[0]=2.33; css[0][0]=9.0; css[1][0]=5.2; par("Si"); // rough avg.
rho[0]=2.34; css[0][0]=8.43; css[1][0]=5.84; par("Si"); // Lucien Saviot's preferred values
//rho[2]=2.33; css[0][2]=9.0; css[1][2]=5.2; mt3("Si"); // rough avg.
//rho[0]=2.2 ; css[0][0]=5.95; css[1][0]=3.76; par("SiO2");
//rho[2]=10.28; css[0][2]=3.65; css[1][2]=1.66; mt3("Ag"); // nu=0.37

//rho[1]=3.6 ; css[0][1]=3.43; css[1][1]=2.14; mat("GeO2"); // GeO2
//rho[1]=2.2 ; css[0][1]=5.95; css[1][1]=3.76; mat("SiO2");
rho[1]=2.2 ; css[0][1]=5.72; css[1][1]=3.75; mat("SiO2"); // Lucien Saviot's preferred values
//rho[2]=2.2 ; css[0][2]=5.95; css[1][2]=3.76; mt3("SiO2");
//rho[2]=2.2 ; css[0][2]=5.95; css[1][2]=3.76; mt3("SiO2");
//rho[1]=2.33; css[0][1]=9.0; css[1][1]=5.2; mat("Si"); // rough avg.
//rho[1]=0.25 ; css[0][1]=0.5; css[1][1]=0.15; mat("imag. mat'l.");
//rho[1]=0.125 ; css[0][1]=1; css[1][1]=0.3; mat("very soft mat'l.");
//rho[2]=0.125 ; css[0][2]=1; css[1][2]=0.3; mt3("very soft mat'l.");
//rho[1]=0.5 ; css[0][1]=1; css[1][1]=0.3; mat("imag. mat'l.");
//rho[2]=0.5 ; css[0][2]=1; css[1][2]=0.3; mt3("imag. mat'l.");
//rho[1]=1 ; css[0][1]=2; css[1][1]=0.6; mat("imag. mat'l.");
//rho[2]=1 ; css[0][2]=2; css[1][2]=0.6; mt3("imag. mat'l.");
//rho[1]=1.5 ; css[0][1]=3; css[1][1]=0.9; mat("imag. mat'l.");
//rho[2]=1.5 ; css[0][2]=3; css[1][2]=0.9; mt3("imag. mat'l.");
//rho[1]=2 ; css[0][1]=4; css[1][1]=1.2; mat("imag. mat'l.");
//rho[2]=2 ; css[0][2]=4; css[1][2]=1.2; mt3("imag. mat'l.");
//rho[2]=1 ; css[0][2]=2; css[1][2]=0.7; mat("imag. mat'l.");
//rho[2]=1.1 ; css[0][2]=3; css[1][2]=1.9; mt3("soft mat'l.");
//rho[1]=2 ; css[0][1]=4; css[1][1]=1.3; mat("imag. mat'l.");
//rho[1]=2.5 ; css[0][1]=5; css[1][1]=1.6; mat("imag. mat'l.");
//rho[1]=7 ; css[0][1]=5; css[1][1]=2.5; mat("imag. mat'l.");
//rho[1]=10; css[0][1]=30.000; css[1][1]=17.32; mat("nu=0.25 mat'l.");
//rho[2]=10; css[0][2]=30.00; css[1][2]=17.32; mat("nu=0.25 mat'l.");
//rho[1]=2.2; css[0][1]=5.4; css[1][1]=3.332; mat("borosil glass ???");
// There are many kinds of borosilicate glass.

int imater;
for (imater=0;imater<NMATER;imater++)
{ //
rho[imater]=rho[imater]*1000; // convert to kg/m^3
css[0][imater]=css[0][imater]*1000; // convert to m/s
css[1][imater]=css[1][imater]*1000;
} //

partrad_nm=1.75; // Particle radius (in nanometers)

//partrad_nm=1.096; // Particle radius (in nanometers) 50 50 cas

gotoxy(5,14);cout<<rho[0]*0.001<<" g/cc cl="<<css[0][0]<<" ct="<<css[1][0]<<" m/s";
//gotoxy(5,16);cout<<rho[2]*0.001<<" g/cc cl="<<css[0][2]<<" ct="<<css[1][2]<<" m/s";
gotoxy(5,18);cout<<rho[1]*0.001<<" g/cc cl="<<css[0][1]<<" ct="<<css[1][1]<<" m/s";

// Function is a series of terms.
int func;
for (func=0;func<NFUNC;func++)
{ // ..
bt[func]=-1; nt[func]=-1; mnt[func]=0;
// nt[]=-1 signifies function not in use.
} // ..

int it1;int ifac;int ifu;
bt[0]=0;mnt[0]=176; if (bt[0]+mnt[0]>=NTRM) {cout<<"* Err 1110";exit(0);}
for (ifu=1;ifu<NFUNC;ifu++)
{ //.
bt[ifu]=bt[ifu-1]+mnt[ifu-1];
mnt[ifu]=176;
if (bt[ifu]+mnt[ifu]>=NTRM) {cout<<"* Err 111 >NTRM "<<ifu;getch();exit(0);}
} //.

// Real: (standing)
#define First_kind  12000
#define Second_kind 12001
// Complex: (moving)
#define In_going    12002
#define Out_going   12003

// Material types:
#define PARTIC 7000
#define MATRIX 7001
#define MATER3 7002

// Polarizations:
#define LONGIT 5000
#define TRANSV 5001

int orderL;
gotoxy(1,1);cout<<"order L? "; cin>>orderL;
if (orderL<=0) return;
gotoxy(5+14+2,19);cout<<"Toroidal l = "<<orderL<<" ";


// Inside particle:
// (1 2 3) = displacment field inside particle

sphbessel(4,orderL,PARTIC,TRANSV,First_kind);
legendre(5,orderL);multiply(6,4,5);nt[4]=-1;nt[5]=-1;
multiplyby(6,"r B"); settozero(7);settozero(8);
curl(1,2,3,6,7,8); nt[6]=-1; nt[7]=-1; nt[8]=-1;

// (1 2 3) = displacment field inside particle

sphbessel(8,orderL,MATRIX,TRANSV,Out_going);
multiplyby(8,"C");
legendre(9,orderL); multiply(10,8,9); nt[8]=-1;
nt[9]=-1; multiplyby(10,"r"); settozero(8); settozero(9);
curl(5,6,7,10,8,9); nt[8]=-1; nt[9]=-1; nt[10]=-1;

// 1 2 3 = displacement field inside
// 5 6 7 = displacement field outside

stress(8,9,10,1,2,3,PARTIC);
stress(11,12,13,5,6,7,MATRIX);

//  8  9 10 = stress field inside particle
// 11 12 13 = stress field in glass matrix

subtract(20,3,7); // phi continuity

subtract(21,13,10); // r-phi force balance





double reS,imS;
gotoxy(40,5);cout<<"Initial guess:";
gotoxy(40,6);cout<<"Re(S)? "; cin>>reS;
gotoxy(40,7);cout<<"Im(S)? "; cin>>imS;

double bestabsdet1=1e20,bestbetar,bestbetai;
bestbetar=reS*3.1415926*css[0][0]/css[1][0];
bestbetai=imS*3.1415926*css[0][0]/css[1][0];

long int countit=20-2;
LB555:
countit++;
if ((countit%10)==0)
{ //
gotoxy(1,10);cout<<countit<<" ";
betar=bestbetar+0.0001*(random(2000)-1000);
betai=bestbetai+0.0001*(random(2000)-1000);
} //
if ((countit%10)==1)
{ //
betar=bestbetar+0.00001*(random(2000)-1000);
betai=bestbetai+0.00001*(random(2000)-1000);
} //
if ((countit%10)==2)
{ //
betar=bestbetar+0.000001*(random(2000)-1000);
betai=bestbetai+0.000001*(random(2000)-1000);
} //
if ((countit%10)==3)
{ //
betar=bestbetar+0.0000001*(random(2000)-1000);
betai=bestbetai+0.0000001*(random(2000)-1000);
} //
if ((countit%10)==4)
{ //
betar=bestbetar+0.00000001*(random(2000)-1000);
betai=bestbetai+0.00000001*(random(2000)-1000);
} //
if ((countit%10)==5)
{ //
betar=bestbetar+0.000000001*(random(2000)-1000);
betai=bestbetai+0.000000001*(random(2000)-1000);
} //
if ((countit%10)==6)
{ //
betar=bestbetar+0.0000000001*(random(2000)-1000);
betai=bestbetai+0.0000000001*(random(2000)-1000);
} //
if ((countit%10)==7)
{ //
betar=bestbetar+0.00000000001*(random(2000)-1000);
betai=bestbetai+0.00000000001*(random(2000)-1000);
} //
if ((countit%10)==8)
{ //
betar=bestbetar+0.000000000001*(random(2000)-1000);
betai=bestbetai+0.000000000001*(random(2000)-1000);
} //
if ((countit%10)==9)
{ //
betar=bestbetar+0.0000000000001*(random(2000)-1000);
betai=bestbetai+0.0000000000001*(random(2000)-1000);
} //


gotoxy(3,2);cout<<"S imag = "<<betai*css[1][0]/(css[0][0]*3.14159)<<"  ";
gotoxy(3,3);cout<<"beta i = "<<betai<<"  ";

gotoxy(3,5);cout<<"omega imag = "<<betai*css[1][0]/(partrad_nm)<<" ns^-1";
if (betai!=0)
{
gotoxy(3,6);cout<<"tau = "<<(partrad_nm)/(betai*css[1][0])<<" ns";
}
//beta = omega * R / CT
//omega = beta * CT / R;
//wavenum = 0.01*beta*css[1][0]/(0.30*2.0*3.14159*partrad_nm);

gotoxy(25,13);cout<<2*partrad_nm<<" nm diam";
frequencyomegar=betar*css[1][0]/partrad_nm;
frequencyomegai=betai*css[1][0]/partrad_nm;

for (imater=0;imater<NMATER;imater++)
{ //.
kwvr[0][imater]=frequencyomegar/css[0][imater];
kwvr[1][imater]=frequencyomegar/css[1][imater];
kwvi[0][imater]=frequencyomegai/css[0][imater];
kwvi[1][imater]=frequencyomegai/css[1][imater];
} //.

nt[8]=-1;
nt[9]=-1;
nt[10]=-1;
nt[11]=-1;
evaluate_func(0+8,20,partrad_nm);
evaluate_func(1+8,21,partrad_nm);
//evaluate_func(2+8,22,partrad_nm);
//evaluate_func(3+8,23,partrad_nm);

cleanup(0+8);
cleanup(1+8);

y[0][0]=coefof(0+8,'B');
y[1][0]=coefof(0+8,'C');

y[0][1]=coefof(1+8,'B');
y[1][1]=coefof(1+8,'C');

complex det1;
double absdet1;

det1=cdet()*exp(-1*betai)/(betar); // good for L=0
absdet1=abs(det1);

if (absdet1<bestabsdet1)
{ //
bestabsdet1=absdet1;
gotoxy(50,9);cout<<"best abs det1: ";
gotoxy(50,10);cout<<bestabsdet1<<"  ";
bestbetar=betar;
bestbetai=betai;
cout.precision(4);
gotoxy(50,18);cout<<"best re S  "<<bestbetar*css[1][0]/(css[0][0]*3.1415926)<<"     ";
gotoxy(50,19);cout<<"best im S  "<<bestbetai*css[1][0]/(css[0][0]*3.1415926)<<"     ";
double wavenum;
wavenum = 0.01*betar*css[1][0]/(cnms*2.0*3.1415926*partrad_nm);
cout.precision(2);
gotoxy(50,21);cout<<"wavenum "<<wavenum<<" cm^-1  ";
cout.precision(4);
gotoxy(50,23);cout<<"tau = "<<
1000.0*(partrad_nm)/(betai*css[1][0])<<" ps  ";
} //

if (!kbhit()) goto LB555;
getch();
return;
} // main

void ddr(int func2,int func1)
{ // ddr
if (func2==func1) {cout<<"* Err ddr 0 ";getch();exit(0);}
if (func2<0)      {cout<<"* Err ddr 1 ";getch();exit(0);}
if (func2>=NFUNC) {cout<<"* Err ddr 2 ";getch();exit(0);}
if (func1<0)      {cout<<"* Err ddr 3 ";getch();exit(0);}
if (func1>=NFUNC) {cout<<"* Err ddr 4 ";getch();exit(0);}
if (mnt[func1]<=0){cout<<"* Err ddr 5 mnt<0 ";getch();exit(0);}
if (mnt[func2]<=0){cout<<"* Err ddr 6 mnt<0";getch();exit(0);}
if (nt[func1]<0)  {cout<<"* Err ddr 7 nt<0";getch();exit(0);}
// Take derivative with respect to r:
int it2;
int it1;
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //.. it1
int i234;
for (i234=0;i234<12*NMATER;i234=i234+6)
{ // i234
// Take derivative with respect to sin(r) term
if (it2>=bt[func2]+mnt[func2])
                  {cout<<"* Err ddr 8 >mnt";getch();exit(0);}
if (it2>=NTRM)    {cout<<"* Err ddr 9 >NTRM";getch();exit(0);}
cr[it2]=cr[it1]*ngetm(i234+2,it1);
ci[it2]=ci[it1]*ngetm(i234+2,it1);
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(i234+5,it2,ngetm(i234+5,it1)+1); // k
nsetm(i234+2,it2,ngetm(i234+2,it1)-1); // sin(k r)
nsetm(i234+4,it2,ngetm(i234+4,it1)+1); // cos(k r)
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
// Take derivative with respect to cos(r) term
if (it2>=bt[func2]+mnt[func2])
                  {cout<<"* Err ddr 10 >mnt";getch();exit(0);}
if (it2>=NTRM)    {cout<<"* Err ddr 11 >NTRM";getch();exit(0);}
cr[it2]=-cr[it1]*ngetm(i234+4,it1);
ci[it2]=-ci[it1]*ngetm(i234+4,it1);
for (ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(i234+5,it2,ngetm(i234+5,it1)+1); // k
nsetm(i234+2,it2,ngetm(i234+2,it1)+1); // sin(k r)
nsetm(i234+4,it2,ngetm(i234+4,it1)-1); // cos(k r)
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
} // i234
// Take derivative with respect to r^n3 term
if (it2>=bt[func2]+mnt[func2]) {cout<<"ddr error 12 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"ddr error 13 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1]*ngetm(7,it1); //
ci[it2]=ci[it1]*ngetm(7,it1); //
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(7,it2,ngetm(7,it1)-1); // r //
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
} //.. it1
nt[func2]=it2-bt[func2];
if (nt[func2]>mnt[func2]) {cout<<"* Err ddr 14 nt>mnt";getch();exit(0);}
if (nt[func2]>mosttermseverddr) {mosttermseverddr=nt[func2];}
return;
} // ddr

void ddth(int func2,int func1)
{ // ddth
if (func2==func1)  {cout<<"* Err ddth 0 ";getch();exit(0);}
if (func1<0)       {cout<<"* Err ddth 1 func1";getch();exit(0);}
if (func1>=NFUNC)  {cout<<"* Err ddth 2 func1";getch();exit(0);}
if (func2<0)       {cout<<"* Err ddth 3 func2";getch();exit(0);}
if (func2>=NFUNC)  {cout<<"* Err ddth 4 func2";getch();exit(0);}
if (mnt[func1]<=0) {cout<<"* Err ddth 5 mnt<0";getch();exit(0);}
if (mnt[func2]<=0) {cout<<"* Err ddth 6 mnt<0";getch();exit(0);}
if (nt[func1]<0)   {cout<<"* Err ddth 7 nt< 0";getch();exit(0);}
int it1,it2;
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //.. it1
// Take derivative with respect to sin(th) term
if (it2>=bt[func2]+mnt[func2]) {cout<<"ddth error 1 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"ddth error 4 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1]*ngetm(0,it1);
ci[it2]=ci[it1]*ngetm(0,it1);
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(0,it2,ngetm(0,it1)-1); // sin(theta)
nsetm(1,it2,ngetm(1,it1)+1); // cos(theta)
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
// derivative with respect to cos(theta) term
if (it2>=bt[func2]+mnt[func2]) {cout<<"ddth error 2 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"ddth error 4 >NTRM Press a key";getch();exit(0);}
cr[it2]=-cr[it1]*ngetm(1,it1);
ci[it2]=-ci[it1]*ngetm(1,it1);
for (ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(0,it2,ngetm(0,it1)+1); // sin(theta)
nsetm(1,it2,ngetm(1,it1)-1); // cos(theta)
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
} //.. it1
nt[func2]=it2-bt[func2];
if (nt[func2]>mnt[func2]) {cout<<"* Err ddth nt>mnt";getch();exit(0);}
return;
} // ddth

void ddph(int func2,int func1)
{ // ddph
if (func2==func1) {cout<<"* Err ddph 0 ";getch();exit(0);}
if (func2<0)      {cout<<"* Err ddph 1 ";getch();exit(0);}
if (func2>=NFUNC) {cout<<"* Err ddph 2 ";getch();exit(0);}
if (func1<0)      {cout<<"* Err ddph 3 ";getch();exit(0);}
if (func1>=NFUNC) {cout<<"* Err ddph 4 ";getch();exit(0);}
if (mnt[func1]<=0){cout<<"* Err ddph 5 mnt<0 ";getch();exit(0);}
if (mnt[func2]<=0){cout<<"* Err ddph 6 mnt<0";getch();exit(0);}
if (nt[func1]<0)  {cout<<"* Err ddph 7 nt<0";getch();exit(0);}
nt[func2]=0; return;
} // ddph

void gradient(int ior2,int iot2,int iop2,int func)
{ // gradient
if (nt[func]<0)   {cout<<"* Err grad 0 nt<0";getch();exit(0);}
if (nt[ior2]!=-1) {cout<<"* Err grad 1 var "<<ior2 <<" in use. ";getch();exit(0);}
if (nt[iot2]!=-1) {cout<<"* Err grad 2 var "<<iot2 <<" in use. ";getch();exit(0);}
if (nt[iop2]!=-1) {cout<<"* Err grad 3 var "<<iop2 <<" in use. ";getch();exit(0);}
ddr(ior2,func);
ddth(iot2,func);
divideby(iot2,"r");
ddph(iop2,func);
divideby(iop2,"r sin(theta)");
return;
} // gradient

void divergence(int func2,int ior1,int iot1,int iop1)
{ // divergence
int ifaux;
if (nt[func2]!=-1) {cout<<"* Err diverg 0 var "<<func2 <<" in use. ";getch();exit(0);}
ifaux=NFUNC-1; // auxilliary function (temp variable)
if (nt[ifaux]!=-1) {cout<<"* Err diverg 1";getch();exit(0);}
if (ifaux==func2)  {cout<<"* Err diverg 2";getch();exit(0);}
multiplyby(ior1,"r*r");
ddr(ifaux,ior1);
divideby(ior1,"r*r");
divideby(ifaux,"r*r");
copyfunc(func2,ifaux);
multiplyby(iot1,"sin(theta)");
ddth(ifaux,iot1);
divideby(iot1,"sin(theta)");
divideby(ifaux,"r sin(theta)");
add(func2,func2,ifaux);
ddph(ifaux,iop1);
divideby(ifaux,"r sin(theta)");
add(func2,func2,ifaux);
nt[ifaux]=-1;
return;
} // divergence

void cleanup(int func)
{ // cleanup

if (nt[func]>mosttermsevercle) {mosttermsevercle=nt[func];}
int it1,it2,it3,i234,i345;
int flag1,flagsame;

for (i234=0;i234<12*NMATER;i234=i234+6)
{ // i234
//  Look for reducible factors of cos(r)^2, cos(kLp r), cos(kTp r):
flag1=0;
doagain1:
for (it1=bt[func];it1<bt[func]+nt[func];it1++)
{ //.. it1
if (ngetm(i234+4,it1)>=2)
{ //.
flag1=1;
nsetm(i234+4,it1,ngetm(i234+4,it1)-2); // divide by cos(r)^2
it2=bt[func]+nt[func];nt[func]=nt[func]+1;
if (it2>=bt[func]+mnt[func]) {cout<<"* Err clean 2 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM)               {cout<<"* Err clean 3 >NTRM Press a key";getch();exit(0);}
cr[it2]=-1*cr[it1]; ci[it2]=-1*ci[it1];
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(i234+2,it2,ngetm(i234+2,it1)+2); // multiply by sin(r)^2
} //.
} // it1
if (flag1==1) goto doagain1;
} // i234

// Next, look for reducible factors of cos(theta)^2
flag1=0;
doagain2:
for (it1=bt[func];it1<bt[func]+nt[func];it1++)
{ //.. it1
if (ngetm(1,it1)>=2)
{ //.
nsetm(1,it1,ngetm(1,it1)-2); // divide by cos(theta)^2
it2=bt[func]+nt[func]; nt[func]=nt[func]+1;
if (it2>=bt[func]+mnt[func]) {cout<<"* Err clean 4 >mnt";getch();exit(0);}
if (it2>=NTRM)               {cout<<"* Err clean 5 >NTRM";getch();exit(0);}
cr[it2]=-1*cr[it1]; ci[it2]=-1*ci[it1];
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
nsetm(0,it2,ngetm(0,it1)+2); // sin(theta)
} //.
} //.. it1
if (flag1==1) goto doagain2;

int ifac;
// Now add up terms of equal powers in all factors:
//  (thus eliminating extra terms)
it1=bt[func];
LBnextit1:
if (it1>=bt[func]+nt[func]-1) goto done3;
it2=it1+1;
LBnextit2:
flagsame=1;
// About half of the time is spent in this following short loop:
for (ifac=0;ifac<NFACT;ifac++)
{ //
if (ngetm(ifac,it1)!=ngetm(ifac,it2)) {flagsame=0;break;}
} //
if (flagsame==1)
{ //. found two terms
cr[it1]=cr[it1]+cr[it2];
ci[it1]=ci[it1]+ci[it2];
cr[it2]=0; ci[it2]=0;
it3=bt[func]+nt[func]-1; // last term in function
cr[it2]=cr[it3]; ci[it2]=ci[it3];
for (ifac=0;ifac<NFACT;ifac++)
{ //.
nsetm(ifac,it2,ngetm(ifac,it3));
} //.
nt[func]=nt[func]-1;
it3=bt[func]+nt[func]-1; // last term in function
if ((cr[it1]==0)&&(ci[it1]==0))
{ //..
cr[it1]=cr[it3]; ci[it1]=ci[it3];
for (ifac=0;ifac<NFACT;ifac++)
{ //.
nsetm(ifac,it1,ngetm(ifac,it3));
} //.
nt[func]=nt[func]-1;
goto LBnextit1;
} //..
if (it2<=bt[func]+nt[func]-1) {goto LBnextit2;}
} //. found two terms
if (it2+1<=bt[func]+nt[func]-1) {it2++;goto LBnextit2;}
if (it1+1<=bt[func]+nt[func]-2) {it1++;goto LBnextit1;}
done3:
return;
} // cleanup
//------------------------------------------
void multiply(int func3,int func2,int func1)
{ // multiply
if (nt[func3]!=-1){cout<<"* Err mult 0 var "<<func3<<" in use. ";getch();exit(0);}
if ((func3==func1)||(func3==func2))
                  {cout<<"* Err mult 1  Press a key";getch();exit(0);}
if (func3<0)      {cout<<"* Err mult 2 ";getch();exit(0);}
if (func3>=NFUNC) {cout<<"* Err mult 3 ";getch();exit(0);}
if (func2<0)      {cout<<"* Err mult 4 ";getch();exit(0);}
if (func2>=NFUNC) {cout<<"* Err mult 5 ";getch();exit(0);}
if (func1<0)      {cout<<"* Err mult 6 ";getch();exit(0);}
if (func1>=NFUNC) {cout<<"* Err mult 7 ";getch();exit(0);}
if (mnt[func1]<=0){cout<<"* Err mult 8 mnt<0 ";getch();exit(0);}
if (mnt[func2]<=0){cout<<"* Err mult 9 mnt<0";getch();exit(0);}
if (mnt[func3]<=0){cout<<"* Err mult 10 mnt<0";getch();exit(0);}
int it1,it2,it3,ifac;
it3=bt[func3];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ // it1
for (it2=bt[func2];it2<bt[func2]+nt[func2];it2++)
{ // it2
if (it3>=bt[func3]+mnt[func3]) {cout<<"abc error 2 >mnt Press a key";getch();exit(0);}
if (it3>=NTRM) {cout<<"jkl error 4 >NTRM Press a key";getch();exit(0);}
cr[it3]=cr[it1]*cr[it2]-ci[it1]*ci[it2];
ci[it3]=cr[it1]*ci[it2]+ci[it1]*cr[it2];
for (ifac=0;ifac<NFACT;ifac++)
{ //.
nsetm(ifac,it3,ngetm(ifac,it1)+ngetm(ifac,it2));
} //.
it3++;
} // it2
} // it1
nt[func3]=it3-bt[func3];
if (nt[func3]>mnt[func3]) {cout<<"Err mult nt>mnt";getch();exit(0);}
cleanup(func3);
return;
} // multiply

void legendre(int func1,int ordern)
{ // legendre
if (nt[func1]!=-1) {cout<<"* Err legend 0 var "<<func1<<" in use. ";getch();exit(0);}
if (ordern<0)      {cout<<"* Err legend 5 ordern<0 ";getch();exit(0);}
if (func1<0)       {cout<<"* Err legend 6 ";getch();exit(0);}
if (func1>=NFUNC)  {cout<<"* Err legend 7 ";getch();exit(0);}
if (mnt[func1]<=0) {cout<<"* Err legend 8 mnt<0 ";getch();exit(0);}
int it1,ifac;
it1=bt[func1];
for (it1=bt[func1];it1<bt[func1]+(ordern/2)+1;it1++)
{ //
for (ifac=0;ifac<NFACT;ifac++)
{ //.
nsetm(ifac,it1,0);
} //.
nsetm(1,it1,ordern-2*(it1-bt[func1]));
ci[it1] = 0;
} //
if (ordern==0) cr[bt[func1]+0] = 1;
if (ordern==1) cr[bt[func1]+0] = 1;
if (ordern==2){cr[bt[func1]+0] = 1.5;cr[bt[func1]+1] = -0.5;}
if (ordern==3){cr[bt[func1]+0] = 2.5;cr[bt[func1]+1] = -1.5;}
if (ordern==4)
{ //.
cr[bt[func1]+0]=35.0/8.0;cr[bt[func1]+1]=-30.0/8.0;cr[bt[func1]+2]=3.0/8.0;
} //.
if (ordern==5)
{ //.
cr[bt[func1]+0]=63.0/8.0;cr[bt[func1]+1]=-70.0/8.0;cr[bt[func1]+2]=15.0/8.0;
} //.
if (ordern>5)  {cout<<"* Err lege 1 order L too high ";getch();exit(0);}
nt[func1]=(ordern/2)+1;
return;
} // legendre

void copyfunc(int func2,int func1)
{ // copyfunc
if (func2<0)      {cout<<"* Err copy 4 ";getch();exit(0);}
if (func2>=NFUNC) {cout<<"* Err copy 5 ";getch();exit(0);}
if (func1<0)      {cout<<"* Err copy 6 ";getch();exit(0);}
if (func1>=NFUNC) {cout<<"* Err copy 7 ";getch();exit(0);}
if (nt[func1]<0 ) {cout<<"* Err copy 7a";getch();exit(0);}
if (mnt[func2]<nt[func1]) {cout<<"Err copyfunc 10 mnt<nt";getch();exit(0);}
nt[func2]=nt[func1];
int it1,ifac,it2;
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //..
cr[it2]=cr[it1]; ci[it2]=ci[it1];
for (ifac=0;ifac<NFACT;ifac++) nsetm(ifac,it2,ngetm(ifac,it1));
it2++;
} //..
return;
} // copyfunc

void curl(int ior2,int iot2,int iop2,int ior1,int iot1,int iop1)
{ // curl
if (nt[ior2]!=-1) {cout<<"* Err curl 0 var "<<ior2 <<" in use. ";getch();exit(0);}
if (nt[iot2]!=-1) {cout<<"* Err curl 1 var "<<iot2 <<" in use. ";getch();exit(0);}
if (nt[iop2]!=-1) {cout<<"* Err curl 2 var "<<iop2 <<" in use. ";getch();exit(0);}
// This assumes that there is no phi dependence of any
// of the components.  One aux. function is needed.
int ifaux1; ifaux1=NFUNC-1;
if (nt[ifaux1]!=-1) {cout<<"* Err curl 3 var in use";getch();exit(0);}
if (ifaux1==ior2) {cout<<"* Err curl 4 var confl";getch();exit(0);}
if (ifaux1==iot2) {cout<<"* Err curl 5 var confl";getch();exit(0);}
if (ifaux1==iop2) {cout<<"* Err curl 6 var confl";getch();exit(0);}
//-----------------------
multiplyby(iop1,"sin(theta)");
ddth(ior2,iop1);
divideby(iop1,"sin(theta)");
divideby(ior2,"r sin(theta)");
//------------------------------
multiplyby(iop1,"r");
ddr(iot2,iop1);
divideby(iop1,"r");
multiplyby(iot2,-1.0);
divideby(iot2,"r");
//------------------------------
multiplyby(iot1,"r");
ddr(iop2,iot1);
divideby(iot1,"r");
ddth(ifaux1,ior1);
multiplyby(ifaux1,-1.0);
add(iop2,iop2,ifaux1);
divideby(iop2,"r");
nt[ifaux1]=-1; // release aux var
return;
} // curl

void add(int func3,int func2,int func1)
{ // add
// It is okay for any of these three functions to be the same.
if (nt[func1]+nt[func2]>mnt[func3])
                  {
gotoxy(1,20);cout<<cout<<"* Err add 0 ";
gotoxy(1,21);cout<<"mnt["<<func3<<"]="<<mnt[func3]<<" (too small)";
gotoxy(1,22);cout<<"nt["<<func1<<"]="<<nt[func1]<<" ";
gotoxy(1,23);cout<<"nt["<<func2<<"]="<<nt[func2]<<" ";
getch();exit(0);
}
if (func1<0)      {cout<<"* Err add 1 ";getch();exit(0);}
if (func1>=NFUNC) {cout<<"* Err add 2 ";getch();exit(0);}
if (func2<0)      {cout<<"* Err add 3 ";getch();exit(0);}
if (func2>=NFUNC) {cout<<"* Err add 4 ";getch();exit(0);}
if (func3<0)      {cout<<"* Err add 5 ";getch();exit(0);}
if (func3>=NFUNC) {cout<<"* Err add 6 ";getch();exit(0);}
if (nt[func1]<0 ) {cout<<"* Err add 7a";getch();exit(0);}
if (nt[func2]<0 ) {cout<<"* Err add 7b";getch();exit(0);}
int i3,it2,it3,orignt1,orignt2;
if (func2==func3) { i3=func2; func2=func1; func1=i3; }
orignt1=nt[func1]; orignt2=nt[func2];
copyfunc(func3,func1);
nt[func3]=orignt1+orignt2;
if (nt[func3]>mnt[func3]) {cout<<"* Err add 9";getch();exit(0);}
if (nt[func3]>mosttermseveradd) {mosttermseveradd=nt[func3];}
it3=bt[func3]+orignt1;
for (it2=bt[func2];it2<bt[func2]+orignt2;it2++)
{ // it2
if (it3>=NTRM)    {cout<<"* Err add 8 >NTRM";getch();exit(0);}
cr[it3]=cr[it2]; ci[it3]=ci[it2];
for (int ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it3,ngetm(ifac,it2));
} //
it3=it3+1;
} // it2
cleanup(func3);
return;
} // add

void stress(int sigrr,int sigrth,int sigrph,int ur,int uth,int uph,int imater)
{ // strain
if (nt[sigrr]!=-1)   {cout<<"* Err stress 0 var "<<sigrr<<" in use. ";getch();exit(0);}
if (nt[sigrth]!=-1)  {cout<<"* Err stress 1 var "<<sigrth<<" in use. ";getch();exit(0);}
if (nt[sigrph]!=-1)  {cout<<"* Err stress 2 var "<<sigrph<<" in use. ";getch();exit(0);}
double lambdalame;
double shearmodmu;
imater=imater-PARTIC;
if (imater<0)        {cout<<"* Err stress 4 var";getch();exit(0);}
if (imater>=NMATER)  {cout<<"* Err stress 4 var";getch();exit(0);}
shearmodmu=rho[imater]*css[1][imater]*css[1][imater];
lambdalame=rho[imater]*css[0][imater]*css[0][imater]-2*shearmodmu;

// Reference http://www.rochester.edu/College/ME/gans/ME445/stress.pdf
// c:/cofrest/stress.pdf
int err,ethth,ephph,lamdilat,ifaux1;
ifaux1=NFUNC-1; err=NFUNC-2; ethth=NFUNC-3; ephph=NFUNC-4; lamdilat=NFUNC-5;
if (nt[ifaux1]!=-1) {cout<<"* err stress 0 var in use.";getch();exit(0);}
if (nt[err]!=-1)    {cout<<"* err stress 1 var in use.";getch();exit(0);}
if (nt[ethth]!=-1)  {cout<<"* err stress 2 var in use.";getch();exit(0);}
if (nt[ephph]!=-1)  {cout<<"* err stress 3 var in use.";getch();exit(0);}
if (nt[lamdilat]!=-1) {cout<<"* err stress 4 var in use.";getch();exit(0);}
if (lamdilat==sigrr)  {cout<<"* err stress 5 var conflic.";getch();exit(0);}
if (lamdilat==sigrth) {cout<<"* err stress 6 var conflic.";getch();exit(0);}
if (lamdilat==sigrph) {cout<<"* err stress 7 var conflic.";getch();exit(0);}
if (ephph==sigrr)   {cout<<"* err stress 8 var conflic.";getch();exit(0);}
if (ephph==sigrth)  {cout<<"* err stress 9 var conflic.";getch();exit(0);}
if (ephph==sigrph)  {cout<<"* err stress 10 var conflic.";getch();exit(0);}
if (ethth==sigrr)   {cout<<"* err stress 11 var conflic.";getch();exit(0);}
if (ethth==sigrth)  {cout<<"* err stress 12 var conflic.";getch();exit(0);}
if (ethth==sigrph)  {cout<<"* err stress 13 var conflic.";getch();exit(0);}
if (err==sigrr)     {cout<<"* err stress 14 var conflic.";getch();exit(0);}
if (err==sigrth)    {cout<<"* err stress 15 var conflic.";getch();exit(0);}
if (err==sigrph)    {cout<<"* err stress 16 var conflic.";getch();exit(0);}
if (ifaux1==sigrr)  {cout<<"* err stress 17 var conflic.";getch();exit(0);}
if (ifaux1==sigrth) {cout<<"* err stress 18 var conflic.";getch();exit(0);}
if (ifaux1==sigrph) {cout<<"* err stress 19 var conflic.";getch();exit(0);}
//- rr strain ----------------
ddr(err,ur);
//- theta theta strain ----------------
ddth(ethth,uth);
add(ethth,ethth,ur);
divideby(ethth,"r");
//- phi phi strain ----------------
copyfunc(ephph,uth);
divideby(ephph,"sin(theta)");
multiplyby(ephph,"cos(theta)");
add(ephph,ephph,ur);
divideby(ephph,"r");
//-----------------------
ddr(sigrph,uph);
copyfunc(ifaux1,uph);
multiplyby(ifaux1,-1.0);
divideby(ifaux1,"r");
add(sigrph,sigrph,ifaux1);
multiplyby(sigrph,0.5);
//---------------------
ddth(sigrth,ur);
divideby(sigrth,"r");
ddr(ifaux1,uth);
add(sigrth,sigrth,ifaux1);
copyfunc(ifaux1,uth);
multiplyby(ifaux1,-1.0);
divideby(ifaux1,"r");
add(sigrth,sigrth,ifaux1);
multiplyby(sigrth,0.5);
//- Divergence -----------------------
add(lamdilat,err,ethth);
add(lamdilat,lamdilat,ephph);
//- Find stress:
multiplyby(lamdilat,lambdalame);
copyfunc(sigrr,err);
multiplyby(sigrr,2.0*shearmodmu);
add(sigrr,sigrr,lamdilat);
multiplyby(sigrph,2.0*shearmodmu);
multiplyby(sigrth,2.0*shearmodmu);
// Release auxilliary variables:
nt[ifaux1]=-1; nt[err]=-1; nt[ethth]=-1; nt[ephph]=-1; nt[lamdilat]=-1;
return;
} // strain

void settozero(int func)
{ // settozero
nt[func]=0; return;
} // settozero

void multiplyby(int func,double r1)
{ // multiplyby (real)
if (r1==0.0) {settozero(func);return;}
for (int it1=bt[func];it1<bt[func]+nt[func];it1++)
 {cr[it1]=r1*cr[it1];ci[it1]=r1*ci[it1];}
return;
} // multiplyby (real)

void multiplyby(int func,char *str1)
{ // multiplyby (string)
multiplybypow(func,str1,1); return;
} // multiplyby (string)

void multiplybypow(int func,char *str1,int ipow)
{ // multiplybypow (string)
if (nt[func]<0) {cout<<"Err mult pow 0 func="<<func<<" ";getch();exit(0);}
int i7,ic1,ic1b,flag1c,flagtheta;
ic1=0;
LBnc:
flag1c=0;
if (str1[ic1+1]==' ') flag1c=1;
if (str1[ic1+1]==0)   flag1c=1;
if (str1[ic1+1]=='*') flag1c=1;
flag1c=1;
i7=-1;
ic1b=ic1+1;
if ((str1[ic1]>='A')&&(str1[ic1]<='L')&&(flag1c==1))
   {i7=3+3*(str1[ic1]-'A');ic1b=ic1+2;}
// {i7=3+3*(str1[ic1]-'A');ic1b=ic1+2;} // ^^^^
if ((str1[ic1]=='r')&&(flag1c==1)) {i7=7; ic1b=ic1+2;}
flagtheta=1;
if (str1[ic1+3]!='(') flagtheta=0;
if (str1[ic1+4]!='t') flagtheta=0;
if (str1[ic1+5]!='h') flagtheta=0;
if (str1[ic1+6]!='e') flagtheta=0;
if (str1[ic1+7]!='t') flagtheta=0;
if (str1[ic1+8]!='a') flagtheta=0;
if (str1[ic1+9]!=')') flagtheta=0;
if (flagtheta==1)
{ // (theta)
if ((str1[ic1+0]=='s')&&(str1[ic1+1]=='i')&&(str1[ic1+2]=='n'))
 {i7=0; ic1b=ic1+11;}
if ((str1[ic1+0]=='c')&&(str1[ic1+1]=='o')&&(str1[ic1+2]=='s'))
 {i7=1; ic1b=ic1+11;}
} // theta
if (i7==-1)    {cout<<"* err multiplyby(s) 1 '"<<str1<<"' ";getch();exit(0); }
if (i7>=NFACT) {cout<<"* err multiplyby(s) 2 '"<<str1<<"'>=NFACT";getch();exit(0);}
for (int it1=bt[func];it1<bt[func]+nt[func];it1++)
 nsetm(i7,it1,ngetm(i7,it1)+ipow);
ic1=ic1b;
if ((str1[ic1-2]==0)||(str1[ic1-1]==0)||(str1[ic1]==0)) return;
if (ic1<80) goto LBnc;
} // multiplybypow (string)

void divideby(int func,char *str1)
{ // divby (string)
multiplybypow(func,str1,-1); return;
} // divby (string)

void subtract(int func3,int func2,int func1)
{ // subtract
if (nt[func3]!=-1) {cout<<"* Err subtra 0 var "<<func3<<" in use. ";getch();exit(0);}
if (func1<0)       {cout<<"* Err subtract 1 ";getch();exit(0);}
if (func1>=NFUNC)  {cout<<"* Err subtract 2 ";getch();exit(0);}
if (func2<0)       {cout<<"* Err subtract 3 ";getch();exit(0);}
if (func2>=NFUNC)  {cout<<"* Err subtract 4 ";getch();exit(0);}
if (func3<0)       {cout<<"* Err subtract 5 ";getch();exit(0);}
if (func3>=NFUNC)  {cout<<"* Err subtract 6 ";getch();exit(0);}
if (nt[func1]<0)   {cout<<"* Err subtract 7a";getch();exit(0);}
if (nt[func2]<0)   {cout<<"* Err subtract 7b";getch();exit(0);}
// It is okay for any of these three functions to be the same.
int i3,it2,it3,signflag=-1,orignt1,orignt2;
int ifac;
if (func2==func3) { i3=func2; func2=func1; func1=i3; signflag*=-1; }
orignt1=nt[func1]; orignt2=nt[func2];
nt[func3]=orignt1+orignt2;
if (nt[func3]>mnt[func3]) {cout<<"* Err subtract 7c";getch();exit(0);}
it3=bt[func3];
for (it2=bt[func1];it2<bt[func1]+orignt1;it2++)
{ // it2
if (it3>=NTRM)   {cout<<"subtract error 8 >NTRM Press a key";getch();exit(0);}
cr[it3]=cr[it2]*signflag; ci[it3]=ci[it2]*signflag;
for (ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it3,ngetm(ifac,it2));
} //
it3=it3+1;
} //
signflag=signflag*(-1);
it3=bt[func3]+orignt1;
for (it2=bt[func2];it2<bt[func2]+orignt2;it2++)
{ // it2
if (it3>=bt[func3]+mnt[func3])
               {cout<<"* Err subtract 9  >mnt ";getch();exit(0);}
if (it3>=NTRM) {cout<<"* Err subtract 10 >NTRM ";getch();exit(0);}
cr[it3]=cr[it2]*signflag;
ci[it3]=ci[it2]*signflag;
for (ifac=0;ifac<NFACT;ifac++)
{ //
nsetm(ifac,it3,ngetm(ifac,it2));
} //
it3=it3+1;
} // it2
cleanup(func3);
return;
} // subtract

void sphbessel(int func1,int orderL,int imater,int ipolar,int wavetype)
{ // sphbessel
if (nt[func1]!=-1)    {cout<<"* Err sphbes 0 var "<<func1<<" in use. ";getch();exit(0);}
imater=imater-7000;
ipolar=ipolar-5000;
int kfactor;
kfactor=1+2*imater+ipolar;
if (imater==-1) kfactor=0;
wavetype=wavetype-12000;
if (wavetype<0)       {cout<<"* Err sphbes wavetype 1 ";getch();exit(0);}
if (wavetype>=4)      {cout<<"* Err sphbes wavetype 2 ";getch();exit(0);}
if (kfactor<=0)       {cout<<"* Err sphbes kfactor <=0 ";getch();exit(0);}
if (kfactor>2*NMATER) {cout<<"* Err sphbes kfactor>4 ";getch();exit(0);}
if (orderL<0)         {cout<<"* Err sphbes 5 ";getch();exit(0);}
if (func1<0)          {cout<<"* Err sphbes 6 ";getch();exit(0);}
if (func1>=NFUNC)     {cout<<"* Err sphbes 7 ";getch();exit(0);}
if (mnt[func1]<=0)    {cout<<"* Err sphbes 8 mnt<0 ";getch();exit(0);}
int it1,ifac,ifaux,it2,iL;
ifaux=NFUNC-1;
if (nt[ifaux]!=-1)    {cout<<"* Err sphbess 9";getch();exit(0);}
it1=bt[ifaux];
for (ifac=0;ifac<NFACT;ifac++) nsetm(ifac, it1, 0);
cr[it1]=1; ci[it1]=0;
nsetm(7,it1,-1); // r
if ((wavetype==0)||(wavetype==2)||(wavetype==3))
{ // 0 = 1st
nsetm(2+6*(kfactor-1),it1,1);  // sin(r)
nsetm(5+6*(kfactor-1),it1,-1); // k
} // 0 = 1st
if (wavetype==1)
{ // 1 = 2nd
nsetm(4+6*(kfactor-1),it1,1); // cos(r)
nsetm(5+6*(kfactor-1),it1,-1); // k
} // 1 = 2nd
nt[ifaux]=1;
if ((wavetype==2)||(wavetype==3))
{ // 3 = outgoing 2 = ingoing
it1=bt[ifaux]+1;
for (ifac=0;ifac<NFACT;ifac++) nsetm(ifac,it1,0);
cr[it1]=0;
if (wavetype==3) ci[it1]=1; // outgoing
if (wavetype==2) ci[it1]=-1; // ingoing
nsetm(7, it1, -1); // r
nsetm(4+6*(kfactor-1),it1,1); // cos(r)
nsetm(5+6*(kfactor-1),it1,-1); // k
nt[ifaux]=2;
} // 3 = outgoing  2 = ingoing
for (iL=0;iL<orderL;iL++)
{ //. iL
ddr(func1,ifaux);
it1=bt[func1];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //..
nsetm(7,it1,ngetm(7,it1)-1); // r //
nsetm(5+6*(kfactor-1),it1,ngetm(5+6*(kfactor-1),it1)-2); // k factor
} //..
copyfunc(ifaux,func1);
} //. iL
for (iL=0;iL<orderL;iL++)
{ //. iL
it1=bt[ifaux];
for (it1=bt[ifaux];it1<bt[ifaux]+nt[ifaux];it1++)
{ //..
cr[it1]=cr[it1]*(-1); ci[it1]=ci[it1]*(-1);
nsetm(7,it1,ngetm(7,it1)+1); // r //
nsetm(5+6*(kfactor-1),it1,ngetm(5+6*(kfactor-1),it1)+1); // k factor
} //..
} //. iL
copyfunc(func1,ifaux);
nt[ifaux]=-1;
return;
} // sphbessel

void evaluate_func(int func2,int func1,double partrad_nm)
{ // evaluate_func
if (nt[func2]!=-1) {cout<<"* Err evalfunc 0: func "<<func2<<" in use.";getch();exit(0);}
if (mnt[func2]<nt[func1])
 {cout<<"* Err evalfunc 1 mnt<nt";getch();exit(0);}
double r1,r2,r3;
int it1,it2,ifac;
it1=bt[func1];
nt[func2]=nt[func1];
for (it2=bt[func2];it2<bt[func2]+nt[func2];it2++)
{ //. it2
cr[it2]=cr[it1]; ci[it2]=ci[it1];
int flag3;
for (ifac=0;ifac<NFACT;ifac++)
{ //..
nsetm(ifac,it2,ngetm(ifac,it1));
flag3=1;
if ((ifac>=3)&&(((ifac-3)%3)==0)) flag3=0;
if (flag3==1) nsetm(ifac,it2,0);
} //..
it1++;
} //. it2
int ipol,imat,i7;
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //. it1
if (ngetm(7,it1)!=0) //
{ // r
r1=partrad_nm; r2=pow(r1,ngetm(7,it1)); //
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} // r
for (ipol=0;ipol<2;ipol++)
{ // ipol
for (imat=0;imat<NMATER;imat++)
{ // imat
i7=6*ipol+12*imat;
if (ngetm(i7+5,it1)!=0)
{ //
complex cr2,cr3;
cr2=complex(kwvr[ipol][imat],kwvi[ipol][imat]);
cr2=pow(cr2,ngetm(i7+5,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //
if (ngetm(i7+2,it1)!=0)
{ //
complex cr1,cr2,cr3;
cr1=complex(kwvr[ipol][imat]*partrad_nm,kwvi[ipol][imat]*partrad_nm);
cr2=sin(cr1);
cr2=pow(cr2,ngetm(i7+2,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //
if (ngetm(i7+4,it1)!=0)
{ //
complex cr1,cr2,cr3;
cr1=complex(kwvr[ipol][imat]*partrad_nm,kwvi[ipol][imat]*partrad_nm);
cr2=cos(cr1);
cr2=pow(cr2,ngetm(i7+4,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //
} // imat
} // ipol
it2++;
} //. it1
return;
} // evaluate_func

void multiplyby(int func,complex cr1)
{ // multiplyby (complex)
double ar1,ai1,r3;
ar1=real(cr1); ai1=imag(cr1);
if ((ar1==0.0)&&(ai1==0.0)) {settozero(func);return;}
for (int it1=bt[func];it1<bt[func]+nt[func];it1++)
{ //.
r3=ar1*cr[it1]-ai1*ci[it1];
ci[it1]=ar1*ci[it1]+ai1*cr[it1];
cr[it1]=r3;
} //.
return;
} // multiplyby (complex)

int ngetm(int i1,int i2)
{ // ngetm
if (i1 < 0  )  {cout<<"* Err ngetm 0";getch();exit(0);}
if (i1>=NFACT) {cout<<"* Err ngetm 1";getch();exit(0);}
if (i2 < 0  )  {cout<<"* Err ngetm 3";getch();exit(0);}
if (i2>=NTRM)  {cout<<"* Err ngetm 4";getch();exit(0);}
int i3;
if      ((i2>=0)          &&(i2<ntrm0))             i3=nm0[i1][i2-0];
else if ((i2>=ntrm0)      &&(i2<ntrm0+ntrm1))       i3=nm1[i1][i2-ntrm0];
else if ((i2>=ntrm0+ntrm1)&&(i2<ntrm0+ntrm1+ntrm2)) i3=nm2[i1][i2-(ntrm0+ntrm1)];
else if ((i2>=ntrm0+ntrm1+ntrm2)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3)) i3=nm3[i1][i2-(ntrm0+ntrm1+ntrm2)];
else if ((i2>=ntrm0+ntrm1+ntrm2+ntrm3)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)) i3=nm4[i1][i2-(ntrm0+ntrm1+ntrm2+ntrm3)];
else if ((i2>=ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3+ntrm4+ntrm5)) i3=nm5[i1][i2-(ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)];
return i3;
} // ngetm

void nsetm(int i1,int i2,int i3)
{ // nsetm
if (i1<0)      {cout<<"* Err nsetm 0   ";getch();exit(0);}
if (i2<0)      {cout<<"* Err nsetm 1   ";getch();exit(0);}
if (i1>=NFACT) {cout<<"* Err nsetm 2   ";getch();exit(0);}
if (i2>=NTRM)  {cout<<"* Err nsetm 3   ";getch();exit(0);}
if      ((i2>=0)          &&(i2<ntrm0))             nm0[i1][i2-0]=i3;
else if ((i2>=ntrm0)      &&(i2<ntrm0+ntrm1))       nm1[i1][i2-ntrm0]=i3;
else if ((i2>=ntrm0+ntrm1)&&(i2<ntrm0+ntrm1+ntrm2)) nm2[i1][i2-(ntrm0+ntrm1)]=i3;
else if ((i2>=ntrm0+ntrm1+ntrm2)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3)) nm3[i1][i2-(ntrm0+ntrm1+ntrm2)]=i3;
else if ((i2>=ntrm0+ntrm1+ntrm2+ntrm3)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)) nm4[i1][i2-(ntrm0+ntrm1+ntrm2+ntrm3)]=i3;
else if ((i2>=ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)&&(i2<ntrm0+ntrm1+ntrm2+ntrm3+ntrm4+ntrm5)) nm5[i1][i2-(ntrm0+ntrm1+ntrm2+ntrm3+ntrm4)]=i3;

//if      ((i2>=0)          &&(i2<ntrm0))             nm0[i1][i2]=i3;
//else if ((i2>=ntrm0)      &&(i2<ntrm0+ntrm1))       nm1[i1][i2-ntrm0]=i3;
//else if ((i2>=ntrm0+ntrm1)&&(i2<ntrm0+ntrm1+ntrm2)) nm2[i1][i2-(ntrm0+ntrm1)]=i3;
return;
} // nsetm

void par(char *ch1)
{ // par
gotoxy(5,13);cout<<"Particle: "<<ch1; return;
} // par

void mat(char *ch1)
{ // mat
gotoxy(5,17);cout<<"Matrix: "<<ch1; return;
} // mat

void mt3(char *ch1)
{ // mat
gotoxy(5,15);cout<<"Mat 3: "<<ch1; return;
} // mat

complex coefof(int func1,char ch1)
{ //
complex c1;
char str1[3];
str1[0]=ch1;
str1[1]=' ';
c1=coefof(func1,str1);
return c1;
} //
complex coefof(int func1,char *str1)
{ //    coefof
int i7=3+3*(str1[0]-'A');
complex cr3=0.0;
for (int it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //..
if (ngetm(i7,it1)!=0) cr3=complex(cr[it1],ci[it1]);
} //..
return cr3;
} //    coefof

complex cdet(void)
{ // det   take determinant of matrix y[][] (y is unchanged]
//  y3[][] is temporary variable
complex det2;
int i1,j1,k1;
for (i1=0;i1<Ny;i1++)
for (j1=0;j1<Ny;j1++)
 y3[i1][j1]=y[i1][j1];

det2=y3[0][0];
for (i1=0;i1<=Ny-2;i1++)
{ // i1
if (y3[i1][i1]==0.0) {cout<<"* Err determ y["<<i1<<","<<i1<<"] is zero.";getch();exit(0);}
for (j1=i1+1;j1<=Ny-1;j1++)
{ // j1.
for (k1=i1+1;k1<=Ny-1;k1++)
{ // k1..
y3[j1][k1]=y3[j1][k1]-y3[j1][i1]*y3[i1][k1]/y3[i1][i1];
} // k1..
} // j1.
det2=det2*y3[i1+1][i1+1];
} // i1
return det2;
} // det   take determinant of matrix y[][] (y is unchanged]