#define VERN "scp86a.cpp"
// Core-Shell Mode Normalizer: (Factored integral method)
// scp86a.cpp  Apr 1  generalize modetype, orderm in file format
// scp84c.cpp  Mar 29 - Two stage radial integral (2nd time done)
//     - seems fast but accurate to four digits
//  84b8  Mar 29 improve integral by using center value
// scp84b7.cpp  works well, factored integral
//  84b6
//  84b2   factored integral  (works well)
//  84b1   Monte Carlo r-theta integral (rather than MC xyz)
// scp84b.cpp  Mar 27
// scp84a8  Mar 26  corrected error in multiplybypow()
// scp84a7.cpp Mar 25  change def of A-F: A old = A/kLm etc.
// scp84a6.cpp
//     - sph bess 2nd kind sign error
//     - exchanged B and C, following Saviot
//     - input is .3 file, output is .4 file
// scp84a3.cpp    Used to make 2: file from 1: file  81a.txt -> 81a2.txt
// scp84a2.cpp  Mar 23
// scp84b4.cpp
//   b3
//   b2
// scp82b.cpp  Mar 23   find <uj|uj> for unnormalized modes  (works)
// scp82a4.cpp  Mar 22 03  Calculate total energy for modes
//          using Monte Carlo integration.
//   a3
// scp82a2  Reads in modes in .txt file calculated by scp81p4.cpp
//          and displays on screen as a test.
// 82a.cpp   read modes in from txt file
//  81p4  eval includes theta, write to file.
// scp81p2  Show that stress is continuous at boundary.
// scp81p  Mar 22  Show rr and r-theta stress near matrix surface.
//     (both seem to go to zero)
//  81n3   Plot r and theta components for eigenfunctions to check
//         for continuity  (seems fine) - detar needs to be small enough
// nm[0][]   sin(theta)
// nm[1][]   cos(theta)
// nm[2][]   sin(kLp r)
// nm[3][]   A
// nm[4][]   cos(kLp r)
// nm[5][]   kLp (longitudinal speed of sound in particle)
// nm[6][]   B
// nm[7][]   r
// nm[8][]   sin(kTp r)
// nm[9][]   C
// nm[10][]  cos(kTp r)
// nm[11][]  kTp (transverse speed of sound in particle)
// nm[12][]  D
//    13  unused
// nm[14][]  sin(kLm r)
// nm[15][]  E
// nm[16][]  cos(kLm r)
// nm[17][]  kLm (longitudinal speed of sound in matrix)
// nm[18][]  F
//    19  unused
// nm[20][]  sin(kTm r)
// nm[21][]  G
// nm[22][]  cos(kTm r)
// nm[23][]  kTm (transverse speed of sound in matrix)
// nm[24][]  H

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <graphics.h>
#include <conio.h>
#include <dos.h>
#include <time.h>
#include <string.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

#define CLP css[0][0]
#define CTP css[1][0]
#define CLM css[0][1]
#define CTM css[1][1]

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
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 etar,partrad_nm,matrad_nm,etai;
double frequencyomegar,frequencyomegai;
complex frequencyomegac;
double kwvr[2][NMATER];
double kwvi[2][NMATER];
void evaluate_func(int func2,int func1,double rad_nm,double theta);
int mosttermseveradd=0;
int mosttermseverddr=0;
int mosttermsevercle=0;
complex coefof(int func1,char *str1);
complex coefof(int func1,char ch1);
long int tm1,tm2,tm3,tm4;
void makescr(void);
FILE *indatfile,*outdatfile;
int outdatflag=0;
int outdatechoflag=1;
char indatfilename[40];
char outdatfilename[40];
char timestring[40];
char readinline(void);
long int licountit;
double sumuurhotot;
double sumrhotot;
double r_nm,ur,utheta,uphi;
double theta6;
double ur9(double r_nm,double theta6);
double utheta9(double r_nm,double theta6);
double uphi9(double r_nm,double theta6);
float nuzzz,Azzz,Bzzz,Czzz,Dzzz,Ezzz,Fzzz;
double rho9(double r_nm);
double h19(double r_nm,double theta6);
double h29(double r_nm,double theta6);
double h39(double r_nm,double theta6);
double h49(double r_nm,double theta6);
double Pi;
int orderL,orderm,modetype;
// modetype  1=spheroidal  2=torsional

void main(void)
{ // main
clrscr();
Pi=4.0*atan(1.0);
randomize();
gotoxy(1,1);cout<<VERN;
gotoxy(1,2);cout<<"Core-Shell Mode Normalizer: (Factored integral method)";
gotoxy(1,6);cout<<"( .3 extension added automatically  'q' to quit ) ";

gotoxy(1,4);
cout<<"Input data file name?  ";
cin>>indatfilename;

gotoxy(1,6);cout<<"                                                  ";
if (indatfilename[0]=='q') return;
int i1;
i1=strlen(indatfilename);
if ((i1>2)&&(indatfilename[i1-2]=='.')) indatfilename[i1-2]=0;
if ((i1>4)&&(indatfilename[i1-4]=='.')) indatfilename[i1-4]=0;
strcpy(outdatfilename,indatfilename);
i1=strlen(indatfilename);
indatfilename[i1+0]='.';
indatfilename[i1+1]='3';
indatfilename[i1+2]=0;
gotoxy(1,4);cout<<"Input from: "<<indatfilename<<"                  ";

char ch1,ch2,ch3,ch4,ch5;
indatfile=fopen(indatfilename,"rb");
fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
if ((ch1!='3')||(ch2!=':'))
{ // file type check
gotoxy(1,6);cout<<"* Error. File is of wrong type *";
gotoxy(1,7);cout<<"* ( Must start with '3:' ) * ";
gotoxy(1,8);cout<<"* Input file name was: "<<indatfilename<<" ";
if (ch1<' ') ch1='?';
if (ch2<' ') ch2='?';
if (ch3<' ') ch3='?';
gotoxy(1,9);cout<<"* First three characters were: "<<ch1<<ch2<<ch3<<" ";
gotoxy(1,15);cout<<"Press any key to exit. ";
getch();
return;
} // file type check

i1=strlen(outdatfilename);
outdatfilename[i1]='.';
outdatfilename[i1+1]='4';
outdatfilename[i1+2]=0;
gotoxy(1,5);cout<<"Output to: "<<outdatfilename<<"                  ";
outdatfile=fopen(outdatfilename,"wb");
fprintf(outdatfile,"%s%s%s","4: ",VERN," ");
time_t thetime;
time(&thetime);
strcpy(timestring,ctime(&thetime));
timestring[strlen(timestring)-1]=0;
fprintf(outdatfile,"%s%s",timestring," ");


ch5=readinline(); // Top line

ch5=readinline(); // Particle: Avg Si
float r5; // Must be f_l_o_a_t, not d_o_u_b_l_e

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);partrad_nm=r5;
fprintf(outdatfile,"%s%9.1f","Rp ",r5);ch5=readinline();

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);rho[0]=r5;
fprintf(outdatfile,"%s%9.1f","rho",r5);ch5=readinline();
gotoxy(1,15);cout<<"rhoP "<<rho[0]<<" kg/m^3";

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);css[0][0]=r5;
fprintf(outdatfile,"%s%9.1f","CLP",r5);ch5=readinline();

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);css[1][0]=r5;
fprintf(outdatfile,"%s%9.1f","CTP",r5);ch5=readinline();

ch5=readinline();  // Matrix SiO2

fscanf(indatfile,"%c%c",&ch1,&ch2);
fscanf(indatfile,"%f",&r5);matrad_nm=r5;
fprintf(outdatfile,"%s%9.1f","Rm ",r5);ch5=readinline();

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);rho[1]=r5;
fprintf(outdatfile,"%s%9.1f","rho",r5);ch5=readinline();
gotoxy(1,18);cout<<"rhoM "<<rho[1]<<" kg/m^3";

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);css[0][1]=r5;
fprintf(outdatfile,"%s%9.1f","CLM",r5);ch5=readinline();
gotoxy(1,19);cout<<"CLM "<<r5<<" m/s";

fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
fscanf(indatfile,"%f",&r5);css[1][1]=r5;
fprintf(outdatfile,"%s%9.1f","CTM",r5);ch5=readinline();
gotoxy(1,20);cout<<"CTM "<<r5<<" m/s";

modetype=-1;
fscanf(indatfile,"%c%c%c",&ch1,&ch2,&ch3);
if ((ch1=='S')&&(ch2=='p')&&(ch3=='h')) modetype=1;
if ((ch1=='T')&&(ch2=='o')&&(ch3=='r')) modetype=2;
if (modetype<1) {cout<<"* Error Unknown mode type: "<<modetype<<" "<<ch1<<ch2<<ch3<<" ";getch();return;}
if (modetype>2) {cout<<"* Error Unknown mode type: "<<modetype<<" "<<ch1<<ch2<<ch3<<" ";getch();return;}
gotoxy(1,11);
if (modetype==1) cout<<"Sph ";
if (modetype==2) cout<<"Tor ";
if (modetype==1) fprintf(outdatfile,"%s","Spheroidal  Order l = ");
if (modetype==2) fprintf(outdatfile,"%s","Torsional   Order l = ");
fscanf(indatfile,"%c%c%c%c%c",&ch1,&ch2,&ch3,&ch4,&ch5);
fscanf(indatfile,"%c%c%c%c%c",&ch1,&ch2,&ch3,&ch4,&ch5);
fscanf(indatfile,"%c%c%c%c%c",&ch1,&ch2,&ch3,&ch4,&ch5);
fscanf(indatfile,"%c%c%c%c%c",&ch1,&ch2,&ch3,&ch4,&ch5);
gotoxy(5,11);cout<<"Order L = "<<ch5<<" ";
orderL=ch5-'0';
orderm=0;
fprintf(outdatfile,"%c",ch5);
//if (orderL!=2) {cout<<"* Err L not 2 *";return;}
ch5=readinline(); // rest of line: Spheroidal  order L = 2

LB100:
ch5=readinline();
if ((!kbhit())&&(ch5!='=')) goto LB100;
cout<<"\n";
outdatechoflag=0;

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

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

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++)
{ //. 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);}
} //. ifu

// 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

if (outdatflag) {fprintf(indatfile,"%s%d%c%c","Spheroidal  Order l = ",orderL,13,10);}

// (1,2,3) = displacement field in nanoparticle E F
// (4,5,6) = displacement field in matrix A B C D
// A = SPH LONG 1st kind matrix
// B = SPH TRAN 1st kind matrix
// C = SPH LONG 2nd kind matrix
// D = SPH TRAN 2nd kind matrix
// E = SPH LONG 1st kind particle
// F = SPH TRAN 1st kind particle

if ((modetype==1)&&(orderL==0))
{ // spheroidal order L=0
sphbessel(7,orderL,MATRIX,LONGIT,First_kind);
multiplyby(7,"A"); divideby(7,"kLm");
sphbessel(8,orderL,MATRIX,LONGIT,Second_kind);
multiplyby(8,"B"); divideby(8,"kLm");
add(9,7,8); nt[7]=-1; nt[8]=-1;
legendre(7,orderL);
multiply(8,7,9); nt[7]=-1; nt[9]=-1;
gradient(4,5,6,8); nt[8]=-1; // (4,5,6) = displacement field in matrix
// C: SPH LONG 1st kind particle
sphbessel(9,orderL,PARTIC,LONGIT,First_kind);
multiplyby(9,"C"); divideby(9,"kLp");
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1; nt[9]=-1;
gradient(1,2,3,7); nt[7]=-1; // (1,2,3) = displacement field in particle
} // spheroidal order L=0

if (modetype==2)
{ // torsional order L
sphbessel(7,orderL,MATRIX,TRANSV,First_kind);
multiplyby(7,"A r"); //divideby(7,"kTm");
sphbessel(8,orderL,MATRIX,TRANSV,Second_kind);
multiplyby(8,"B r"); //divideby(8,"kTm");
add(9,7,8); nt[7]=-1; nt[8]=-1;
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1;nt[9]=-1;
settozero(8); settozero(9);
curl(4,5,6,7,8,9); nt[7]=-1; nt[8]=-1; nt[9]=-1;
sphbessel(9,orderL,PARTIC,TRANSV,First_kind);
multiplyby(9,"C r"); //divideby(9,"kTm");
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1;nt[9]=-1;
settozero(8); settozero(9);
curl(1,2,3,7,8,9); nt[7]=-1; nt[8]=-1; nt[9]=-1;
} // torsional order L

if ((modetype==1)&&(orderL>=1))
{ // spheroidal order L>=1
sphbessel(7,orderL,MATRIX,LONGIT,First_kind);
multiplyby(7,"A");
divideby(7,"kLm");
sphbessel(8,orderL,MATRIX,LONGIT,Second_kind);
multiplyby(8,"C");
divideby(8,"kLm");
add(9,7,8); nt[7]=-1; nt[8]=-1;
legendre(7,orderL);
multiply(8,7,9); nt[7]=-1; nt[9]=-1;
gradient(4,5,6,8); // (4,5,6) = displacement field in matrix
nt[8]=-1; // release

sphbessel(7,orderL,MATRIX,TRANSV,First_kind);
multiplyby(7,"B r");
divideby(7,"kTm");
sphbessel(8,orderL,MATRIX,TRANSV,Second_kind);
multiplyby(8,"D r");
divideby(8,"kTm");
add(9,7,8); nt[7]=-1; nt[8]=-1;
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1;nt[9]=-1;
settozero(8); settozero(9);
curl(10,11,12,7,8,9); nt[7]=-1; nt[8]=-1; nt[9]=-1;
curl(7,8,9,10,11,12); nt[10]=-1; nt[11]=-1; nt[12]=-1;
add(4,4,7); nt[7]=-1;
add(5,5,8); nt[8]=-1;
add(6,6,9); nt[9]=-1;

// E: SPH LONG 1st kind particle
sphbessel(9,orderL,PARTIC,LONGIT,First_kind);
multiplyby(9,"E");
divideby(9,"kLp");
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1; nt[9]=-1;
gradient(1,2,3,7); // (1,2,3) = displacement field in particle
nt[7]=-1; // release

// F: SPH TRAN 1st kind particle
sphbessel(9,orderL,PARTIC,TRANSV,First_kind);
multiplyby(9,"F r");
divideby(9,"kTp");
legendre(8,orderL);
multiply(7,8,9); nt[8]=-1;nt[9]=-1;
settozero(8); settozero(9);
curl(10,11,12,7,8,9); nt[7]=-1; nt[8]=-1; nt[9]=-1;
curl(7,8,9,10,11,12); nt[10]=-1; nt[11]=-1; nt[12]=-1;
add(1,1,7); nt[7]=-1;
add(2,2,8); nt[8]=-1;
add(3,3,9); nt[9]=-1;
} // spheroidal order L>=1


int whichmode;
whichmode=-1;
LBnextmode:
whichmode++;
int i15,i16,i17m;
fscanf(indatfile,"%d%d%d",&i15,&i16,&i17m);
if (i15!=modetype)
{ //..
cout<<"* Error i15 != modetype *";return;
} //..
if (i16!=orderL)
{ //.
cout<<"* Error i16 != orderL *";return;
} //.
if (i17m!=orderm)
{ //.
cout<<"* Error i17m != orderm *";return;
} //.
fscanf(indatfile,"%f",&nuzzz);
if (feof(indatfile))
{ //..
fclose(outdatfile);
return;
} //..
fscanf(indatfile,"%f",&Azzz);
fscanf(indatfile,"%f",&Bzzz);
fscanf(indatfile,"%f",&Czzz);
fscanf(indatfile,"%f",&Dzzz);
fscanf(indatfile,"%f",&Ezzz);
fscanf(indatfile,"%f",&Fzzz);
cout.precision(3);
gotoxy(1,7+whichmode%18);
cout<<""<<whichmode<<". "<<i15<<" "<<i16<<" "<<nuzzz<<" "<<Azzz<<" "<<Bzzz<<"  ";
ch5=readinline();

gotoxy(40,1);cout<<"mode = "<<whichmode;
gotoxy(40,2);cout<<"nu = "<<nuzzz<<" cm^-1";
frequencyomegar=2.0*Pi*100.0*c_nms*nuzzz;
frequencyomegai=0.0;

gotoxy(40,4);cout<<"omega = "<<frequencyomegar<<" rad/ns ";

int imater;
// Wavevectors are in nm^-1 ?
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];
} //.

gotoxy(40,5);cout<<"kLP = "<<kwvr[0][0]<<" nm^-1 ";
gotoxy(40,6);cout<<"kTP = "<<kwvr[1][0]<<" nm^-1 ";
gotoxy(40,7);cout<<"kLM = "<<kwvr[0][1]<<" nm^-1 ";
gotoxy(40,8);cout<<"kTM = "<<kwvr[1][1]<<" nm^-1 ";

double int1,int1a,int1b,int1c,int2,int2a,int2b,int2c;
double int3,int3a,int3b,int3c,int4,int4a,int4b,int4c;
double r2,theta2,r0,theta0,dr01,dr02,dtheta0;
double maxh;
{ //..
double hval,theta21,r21;
maxh=0;

for (int i5=0;i5<25;i5++)
{ // i5
r21=matrad_nm*0.0001*(1+random(10000));
theta21  = Pi*0.0001*random(10000);
hval=fabs(h19(r21,theta21));
if (hval>maxh) { r2=r21; theta2=theta21; maxh=hval; }
} // i5
} //..
if (kbhit()) {getch();return;}

dr01=matrad_nm/600.0;  // small part of integral
dr02=(matrad_nm-partrad_nm)/1200.0; // major part
dtheta0=Pi/600.0;  // no sharp features

int1a=0;
for (r0=0.5*dr01;r0<partrad_nm;r0=r0+dr01)
{ // r0
double r7; r7=h19(r0,theta2); int1a+=dr01*r7;
} // r0
for (r0=partrad_nm+0.5*dr02;r0<matrad_nm;r0=r0+dr02)
{ // r0
double r7; r7=h19(r0,theta2); int1a+=dr02*r7;
} // r0

int1b=0;
for (theta0=0.5*dtheta0;theta0<Pi;theta0=theta0+dtheta0)
{ // r0
double r7; r7=h19(r2,theta0); int1b+=dtheta0*r7;
} // r0
int1c=h19(r2,theta2);
//if (int1c==0.0) {cout<<"* Err int1c=0 *";return;}
int1=0;
if (int1c!=0) int1=int1a*int1b/int1c;
if (kbhit()) {getch();return;}

{ //..
double hval,theta21,r21;
maxh=0;
for (int i5=0;i5<25;i5++)
{ // i5
r21=matrad_nm*0.0001*(1+random(10000));
theta21=Pi*0.0001*random(10000);
hval=fabs(h29(r21,theta21));
if (hval>maxh) {r2=r21; theta2=theta21;maxh=hval;}
} // i5
} //..

int2a=0;
for (r0=0.5*dr01;r0<partrad_nm;r0=r0+dr01)
{ // r0
int2a+=dr01*h29(r0,theta2);
} // r0
for (r0=partrad_nm+0.5*dr02;r0<matrad_nm;r0=r0+dr02)
{ // r0
int2a+=dr02*h29(r0,theta2);
} // r0

int2b=0;
for (theta0=0.5*dtheta0;theta0<Pi;theta0=theta0+dtheta0)
{ // r0
double r7; r7=h29(r2,theta0); int2b+=dtheta0*r7;
} // r0
int2c=h29(r2,theta2);
int2=0;
if (int2c!=0) int2=int2a*int2b/int2c;
if (kbhit()) {getch();return;}

//---------
{ //..
double hval,theta21,r21;
maxh=0;
for (int i5=0;i5<25;i5++)
{ // i5
r21=matrad_nm*0.0001*(1+random(10000));
theta21=Pi*0.0001*random(10000);
hval=fabs(h39(r21,theta21));
if (hval>maxh) {r2=r21; theta2=theta21;maxh=hval;}
} // i5
} //..

int3a=0;
for (r0=0.5*dr01;r0<partrad_nm;r0=r0+dr01)
{ // r0
int3a+=dr01*h39(r0,theta2);
} // r0
for (r0=partrad_nm+0.5*dr02;r0<matrad_nm;r0=r0+dr02)
{ // r0
int3a+=dr02*h39(r0,theta2);
} // r0

int3b=0;
for (theta0=0.5*dtheta0;theta0<Pi;theta0=theta0+dtheta0)
{ // r0
double r7; r7=h39(r2,theta0); int3b+=dtheta0*r7;
} // r0
int3c=h39(r2,theta2);
int3=0;
if (int3c!=0) int3=int3a*int3b/int3c;
if (kbhit()) {getch();return;}

{ //..
double hval,theta21,r21;
maxh=0;
for (int i5=0;i5<15;i5++)
{ // i5
r21=matrad_nm*0.0001*(1+random(10000));
theta21=Pi*0.0001*random(10000);
hval=fabs(h49(r21,theta21));
if (hval>maxh) {r2=r21; theta2=theta21;maxh=hval;}
} // i5
} //..
int4a=0;
for (r0=0.5*dr01;r0<partrad_nm;r0=r0+dr01)
{ // r0
int4a+=dr01*h49(r0,theta2);
} // r0
for (r0=partrad_nm+0.5*dr02;r0<matrad_nm;r0=r0+dr02)
{ // r0
int4a+=dr02*h49(r0,theta2);
} // r0

int4b=0;
for (theta0=0.5*dtheta0;theta0<Pi;theta0=theta0+dtheta0)
{ // r0
int4b+=dtheta0*h49(r2,theta0);
} // r0
int4c=h49(r2,theta2);
int4=int4a*int4b/int4c;

cout.precision(5);
gotoxy(40,11);cout<<"int 1 "<<int1<<"   ";
gotoxy(40,12);cout<<"int 2 "<<int2<<"   ";
gotoxy(40,13);cout<<"int 3 "<<int3<<"   ";
gotoxy(40,14);cout<<"int 4 "<<int4<<"   ";

gotoxy(40,16);cout<<"int = "<<(int1+int2+int3)/int4<<" ";

double r33;
r33=sqrt(int4/(int1+int2+int3));
fprintf(outdatfile,"%1d%2d%3d",i15,i16,i17m);
fprintf(outdatfile,"%10.4f",nuzzz);
fprintf(outdatfile,"%10.4f",Azzz*r33);
fprintf(outdatfile,"%10.4f",Bzzz*r33);
fprintf(outdatfile,"%10.4f",Czzz*r33);
fprintf(outdatfile,"%10.4f",Dzzz*r33);
if (fabs(r33*Ezzz)<999.0)
 fprintf(outdatfile,"%10.4f",Ezzz*r33);
else
 fprintf(outdatfile,"%10.2e",Ezzz*r33);
if (fabs(r33*Fzzz)<999.0)
 fprintf(outdatfile,"%10.4f",Fzzz*r33);
else
 fprintf(outdatfile,"%10.2e",Fzzz*r33);
fprintf(outdatfile,"%c%c",13,10);
if (!kbhit()) goto LBnextmode;

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

void dispterm(int it2,int ic,int ir)
{ // dispterm
if (it2<0)     {cout<<"* Error dispte 0 ";exit(0);}
if (it2>=NTRM) {cout<<"* Error dispte 1  ";exit(0);}
if (ic<1)      {cout<<"* Error dispte 2  ";exit(0);}
if (ic>70)     {cout<<"* Error dispte 3  ";exit(0);}
if (ir<2)      {cout<<"* Error dispte 4 - disp. row < 2 ";exit(0);}
//if (ir>25) {cout<<"Error 32 irow off screen";getch();exit(0);}
if (ir>25) {ir=ir-24;ic=ic+40;}
// Display a single term:
//ir= // row on screen
//ic= // starting column
gotoxy(ic,ir);
if ((cr[it2]!=0)&&(ci[it2]!=0)) {cout<<"(";ic=ic+1;}
if (cr[it2]<0) {cout<<"-";ic=ic+1;}
if (
    ((cr[it2]!=0)&&(ci[it2]!=0)) ||
    ((cr[it2]!=1)&&(cr[it2]!=-1)&&(ci[it2]==0)) ||
    ((cr[it2]==0)&&(ci[it2]==0))
   )
{ // print real part
double r1=fabs(cr[it2]);
if (r1<0.001)
{ //
cout<<"0";ic=ic+1;
} //
else if ((r1<10)&&(ceil(r1)==r1))
{ //
cout<<r1;ic=ic+1;
} //
else if ((r1<100)&&(ceil(r1)==r1))
{ //
cout<<r1;ic=ic+2;
} //
else
{ //
  cout<<r1;ic=ic+5;
} //
} // print real part
if ((cr[it2]!=0)&&(fabs(ci[it2])>0.001)) {gotoxy(ic,ir);cout<<"+";ic=ic+1;}
if (ci[it2]<-0.001) {gotoxy(ic,ir);cout<<"-";ic=ic+1;}
if ((ci[it2]!=0)&&(ci[it2]!=1)&&(ci[it2]!=-1))
{ // print imag part
double r1=fabs(ci[it2]);
if (r1<0.001)
{ //
cout<<" ";ic=ic+1;
} //
else if ((r1<10)&&(ceil(r1)==r1))
{ //
cout<<r1;ic=ic+1;
} //
else if ((r1<100)&&(ceil(r1)==r1))
{ //
cout<<r1;ic=ic+2;
} //
else
{ //
cout<<r1;ic=ic+5;
} //
} // print imag part
if (fabs(ci[it2])>0.001) {gotoxy(ic,ir);cout<<"i";ic=ic+1;}
if ((cr[it2]!=0)&&(ci[it2]!=0))
{
cout<<")";ic=ic+1;
}
if ((cr[it2]!=1)||(ci[it2]!=0)) {gotoxy(ic,ir);cout<<" ";ic=ic+1;}

int anyfactors=0;
int ifac;
for (ifac=0;ifac<NFACT;ifac++)
{ //
if (ngetm(ifac,it2)!=0) anyfactors=1;
} //
if ((cr[it2]==1)&&(ci[it2]==0)&&(anyfactors==0))
{ //
cout<<"1";ic=ic+1;
} //
for (int ivar2=0;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(5+6*ivar2,it2)!=0)
{ // kLp
gotoxy(ic,ir);
if (ivar2==0) cout<<"kLp  ";
if (ivar2==1) cout<<"kTp  ";
if (ivar2==2) cout<<"kLm  ";
if (ivar2==3) cout<<"kTm  ";
if (ivar2==4) cout<<"kL3  ";
if (ivar2==5) cout<<"kT3  ";
if (ivar2>5) {cout<<"* Err dispte 3553 ";getch();exit(0);}
ic=ic+3;
if (ngetm(5+6*ivar2,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(5+6*ivar2,it2)<<"  ";ic=ic+1;
} //..
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} // kLp
} // ivar2

for (ivar2=0;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(2+6*ivar2,it2)!=0)
{ //..
gotoxy(ic,ir);
if (ivar2==-1){cout<<"sin(r)  ";ic=ic+6;}
if (ivar2==0) {cout<<"sin(kLp r)  ";ic=ic+10;}
if (ivar2==1) {cout<<"sin(kTp r)  ";ic=ic+10;}
if (ivar2==2) {cout<<"sin(kLm r)  ";ic=ic+10;}
if (ivar2==3) {cout<<"sin(kTm r)  ";ic=ic+10;}
if (ivar2==4) {cout<<"sin(kL3 r)  ";ic=ic+10;}
if (ivar2==5) {cout<<"sin(kT3 r)  ";ic=ic+10;}
if (ivar2>5) {cout<<"* Err dispte 3554 ";getch();exit(0);}
gotoxy(ic,ir-1);if (ngetm(2+6*ivar2,it2)!=1) {cout<<ngetm(2+6*ivar2,it2)<<" ";ic=ic+1;}
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} //..
} // ivar2

for (ivar2=0;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(4+6*ivar2,it2)!=0)
{ //..
gotoxy(ic,ir);
if (ivar2==-1) {cout<<"cos(r)  ";ic=ic+6;}
if (ivar2==0) {cout<<"cos(kLp r)  ";ic=ic+10;}
if (ivar2==1) {cout<<"cos(kTp r)  ";ic=ic+10;}
if (ivar2==2) {cout<<"cos(kLm r)  ";ic=ic+10;}
if (ivar2==3) {cout<<"cos(kTm r)  ";ic=ic+10;}
if (ivar2==4) {cout<<"cos(kL3 r)  ";ic=ic+10;}
if (ivar2==5) {cout<<"cos(kT3 r)  ";ic=ic+10;}
if (ivar2>5) {cout<<"* Err dispte 3555 ";getch();exit(0);}
gotoxy(ic,ir-1);if (ngetm(4+6*ivar2,it2)!=1) {cout<<ngetm(4+6*ivar2,it2)<<" ";ic=ic+1;}
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} //..
} // ivar2

if (ngetm(7,it2)!=0) //
{ //
gotoxy(ic,ir);cout<<"r  ";ic=ic+1;
if (ngetm(7,it2)!=1) //
{ //..
gotoxy(ic,ir-1);cout<<ngetm(7,it2)<<" ";ic=ic+1; //
} //..
ic=ic+1;
} //
if (ngetm(0,it2)!=0)
{ //
gotoxy(ic,ir);cout<<"sin(th)  ";ic=ic+7;
if (ngetm(0,it2)!=1)
{ //.
gotoxy(ic,ir-1);cout<<ngetm(0,it2)<<" ";ic=ic+1;
} //.
ic=ic+1;
} //
if (ngetm(1,it2)!=0)
{ //
gotoxy(ic,ir);cout<<"cos(th)  ";ic=ic+7;
if (ngetm(1,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(1,it2)<<" ";ic=ic+1;
} //..
ic=ic+1;
} //

// Display coefficient variables: A, B, C, D...
int ivar; char sch1[4];
for (ivar=0;ivar<2*NMATER;ivar++)
{ // ivar
sch1[0]='A'+ivar; sch1[1]=' ';sch1[2]=' ';sch1[3]=0;
if (ngetm(3+3*ivar,it2)!=0)
{ // A
gotoxy(ic,ir);cout<<sch1;ic=ic+1;
if (ngetm(3+3*ivar,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(3+3*ivar,it2)<<" ";ic=ic+1;
} //..
ic=ic+1;
} // A
} // ivar

return;
} // dispterm

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<<"* mult error 11 >mnt Press a key";getch();exit(0);}
if (it3>=NTRM) {cout<<"* mult error 12 >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)
// Mar 25 03 - add kLp etc.
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;}
if ((str1[ic1]=='r')&&(flag1c==1)) {i7=7; ic1b=ic1+2;}
if ((str1[ic1]=='k')&&(str1[ic1+1]=='L')&&(str1[ic1+2]=='p')&&(flag1c==1))
{ i7=5;ic1b=ic1+4; }
if ((str1[ic1]=='k')&&(str1[ic1+1]=='T')&&(str1[ic1+2]=='p')&&(flag1c==1))
{ i7=11;ic1b=ic1+4; }
if ((str1[ic1]=='k')&&(str1[ic1+1]=='L')&&(str1[ic1+2]=='m')&&(flag1c==1))
{ i7=17;ic1b=ic1+4; }
if ((str1[ic1]=='k')&&(str1[ic1+1]=='T')&&(str1[ic1+2]=='m')&&(flag1c==1))
{ i7=23;ic1b=ic1+4; }
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
// Mar 24 03 - corrected 2nd kind sign error
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
cr[it1]=-1; ci[it1]=0; // Mar 24 '03
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 rad_nm,double theta)
{ // 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;
// Step 1: Copy terms from func1 to func2
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;
// Only factors A, B, C... are kept after evaluation
if ((ifac>=3)&&(((ifac-3)%3)==0)) flag3=0;
if (flag3==1) nsetm(ifac,it2,0);
} //..
it1++;
} //. it2
int ipol,imat,i7;
// Step 2: !!!!!
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //. it1
if (ngetm(0,it1)!=0)
{ // sin(theta)
cr[it2]=cr[it2]*pow(sin(theta),ngetm(0,it1));
ci[it2]=ci[it2]*pow(sin(theta),ngetm(0,it1));
} // sin(theta)
if (ngetm(1,it1)!=0)
{ // cos(theta)
cr[it2]=cr[it2]*pow(cos(theta),ngetm(1,it1));
ci[it2]=ci[it2]*pow(cos(theta),ngetm(1,it1));
} // cos(theta)
if (ngetm(7,it1)!=0)
{ // r
r1=rad_nm;
if ((r1==0.0)&&(ngetm(7,it1)==0)) {cout<<"* err pow 1234 *";exit(0);}
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]);
if ((cr2==0.0)&&(ngetm(i7+5,it1)==0)) {cout<<"* err pow 1235 *";exit(0);}
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]*rad_nm,kwvi[ipol][imat]*rad_nm);
cr2=sin(cr1);
if ((cr2==0.0)&&(ngetm(i7+2,it1)==0)) {cout<<"* err pow 1236 *";exit(0);}
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]*rad_nm,kwvi[ipol][imat]*rad_nm);
cr2=cos(cr1);
if ((cr2==0.0)&&(ngetm(i7+4,it1)==0)) {cout<<"* err pow 1237 *";exit(0);}
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;
return;
} // nsetm

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

char readinline(void)
{ // readinline
//cout<<"\n";
char ch1,chf;
int countchar=0;
LB333:
countchar++;
fscanf(indatfile,"%c",&ch1);
if (outdatechoflag==1) fprintf(outdatfile,"%c",ch1);
if (feof(indatfile)) return 'x';
if (countchar==1) chf=ch1;
//if (ch1>=32) cout<<ch1;
if ((countchar<80)&&(ch1!=10)) goto LB333;
//fscanf(indatfile,"%c",&ch1);
return chf;
} // readinline

double ur9(double r_nm,double theta6)
{ // ur9
double ur;
if (modetype==2) return 0.0;

if ((modetype==1)&&(orderL>=1))
{ // SPH L>=1
if (r_nm<=partrad_nm)
{ // r_nm
evaluate_func(13,1,r_nm,theta6);
cleanup(13);
complex ce,cf;
ce=coefof(13,'E');
cf=coefof(13,'F');
ur=real(ce)*Ezzz+real(cf)*Fzzz;
nt[13]=-1;
} // r_nm
else
{ // r_nm
evaluate_func(13,4,r_nm,theta6);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
ur=real(ca)*Azzz
  +real(cb)*Bzzz
  +real(cc)*Czzz
  +real(cd)*Dzzz;
nt[13]=-1;
} // r_nm
} // SPH L>=1

if ((modetype==1)&&(orderL==0))
{ // SPH L=0
if (r_nm<=partrad_nm)
{ // r_nm
evaluate_func(13,1,r_nm,theta6);
cleanup(13);
complex cc,cf;
cc=coefof(13,'C');
ur=real(cc)*Czzz;
nt[13]=-1;
} // r_nm
else
{ // r_nm
evaluate_func(13,4,r_nm,theta6);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
ur=real(ca)*Azzz
  +real(cb)*Bzzz;
nt[13]=-1;
} // r_nm
} // SPH L=0

return ur;
} // ur9

double utheta9(double r_nm,double theta6)
{ //   utheta9
double utheta;
if (modetype==2) return 0.0;
if ((modetype==1)&&(orderL==0)) return 0.0;
if (r_nm<=partrad_nm)
{ // r_nm
evaluate_func(13,2,r_nm,theta6);
cleanup(13);
complex ce,cf;
ce=coefof(13,'E');
cf=coefof(13,'F');
utheta=real(ce)*Ezzz
      +real(cf)*Fzzz;
nt[13]=-1;
} // r_nm
else
{ // r_nm
evaluate_func(13,5,r_nm,theta6);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
utheta=real(ca)*Azzz
      +real(cb)*Bzzz
      +real(cc)*Czzz
      +real(cd)*Dzzz;
nt[13]=-1;
} // r_nm
return utheta;
} // utheta9

double uphi9(double r_nm,double theta6)
{ //   uphi9
double uphi;
if (modetype!=2) return 0.0;

if (r_nm<=partrad_nm)
{ // r_nm
evaluate_func(13,3,r_nm,theta6);
cleanup(13);
complex cc;
cc=coefof(13,'C');
uphi=real(cc)*Czzz;
nt[13]=-1;
} // r_nm
else
{ // r_nm
evaluate_func(13,6,r_nm,theta6);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
uphi=real(ca)*Azzz
    +real(cb)*Bzzz;
nt[13]=-1;
} // r_nm
return uphi;
} // uphi9

double rho9(double r_nm)
{ //   rho9
double rho6_kgm3=rho[0]; // inside nanoparticle
if (r_nm>partrad_nm) rho6_kgm3=rho[1]; // in matrix
return rho6_kgm3;
} //   rho9

double h19(double r_nm,double theta6)
{ // h19
double h1,ur,rho_kgm3;
ur= ur9(r_nm,theta6);
rho_kgm3=rho9(r_nm);
h1= ur*ur*rho_kgm3*2*Pi*sin(theta6)*r_nm*r_nm;
return h1;
} // h19

double h29(double r_nm,double theta6)
{ //   h29
double h2,utheta,rho_kgm3;
utheta= utheta9(r_nm,theta6);
rho_kgm3=rho9(r_nm);
h2= utheta*utheta*rho_kgm3*2.0*Pi*sin(theta6)*r_nm*r_nm;
return h2;
} //   h29

double h39(double r_nm,double theta6)
{ //   h39
double h3,uphi,rho_kgm3;
uphi= uphi9(r_nm,theta6);
rho_kgm3=rho9(r_nm);
h3= uphi*uphi*rho_kgm3*2.0*Pi*sin(theta6)*r_nm*r_nm;
return h3;
} //   h39

double h49(double r_nm,double theta6)
{ //   h49
double h4;
h4=2.0*Pi*sin(theta6)*r_nm*r_nm*rho9(r_nm);
return h4;
} //   h49