// scp78.cpp Jan 1 2003  A=1/f (from md30.htm)
double thick3=0.03;
// 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
// 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 <stdio.h>
#include <conio.h>
#include <dos.h>
#include <time.h>
#include <graphics.h>
#include <complex.h>
#define NFUNC 28
#define NTRM 4932
#define NMATER 3
#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 beta,frequencyomega,partrad_nm;
double kwv[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];

void main(void)
{ // main
int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();

setcolor(8);rectangle(0,479-301,301,479);

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

setcolor(6);
line(bx0,479-301,bx0,by0);
/*
for (i77=1;i77<=wavenummax;i77=i77+1)
{ //
setcolor(8);
line(bx0+sfx2*i77,479-301,bx0+sfx2*i77,by0);
} //
*/
for (i77=10;i77<=wavenummax;i77=i77+10)
{ //
setcolor(8);
line(bx0+sfx2*i77,479-301,bx0+sfx2*i77,by0);
} //
setcolor(6);
line(bx0,by0,bx0+wavenummax*sfx2,by0);


//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]=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[1]=3.6 ; css[0][1]=3.43; css[1][1]=2.14; mat("GeO2"); // GeO2
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.72; css[1][2]=3.75; mat("SiO2"); // Lucien Saviot's preferred values
//rho[1]=2.2 ; css[0][1]=5.95; css[1][1]=3.76; mat("SiO2");
//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]=2; 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=0.5*3.5; // Particle radius (in nanometers)
double cnms=1.0e-9*2.998e8; // speed of light in nm/s
betamax=100.0*wavenummax*(cnms*2.0*3.1415926*partrad_nm)/css[1][0];

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";

// Calculate Poisson ratios:
double nup,num,nu3;
//nup=0.5*(1.0-2.0*(css[1][0]/css[0][0])*(css[1][0]/css[0][0]))/(1.0-(css[1][0]/css[0][0])*(css[1][0]/css[0][0]));
//gotoxy(60,6);cout<<"nu p = "<<nup<<" ";
//num=0.5*(1.0-2.0*(css[1][1]/css[0][1])*(css[1][1]/css[0][1]))/(1.0-(css[1][1]/css[0][1])*(css[1][1]/css[0][1]));
//gotoxy(60,7);cout<<"nu m = "<<num<<" ";
//nu3=0.5*(1.0-2.0*(ct3/cl3)*(ct3/cl3))/(1.0-(ct3/cl3)*(ct3/cl3));
//gotoxy(60,3);cout<<"nu 3 = "<<nu3<<" ";

// 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;
orderL=2;
gotoxy(5+14+2,19);cout<<"Spheroid. l = "<<orderL<<" ";

/*
ic2=1; irow2=4;
func=1;
gotoxy(ic2,irow2-1);cout<<" ";
gotoxy(ic2,irow2);  cout<<" ";
for (it1=bt[func];it1<bt[func]+nt[func];it1++)
{ //.. it1
dispterm(it1,ic2,irow2); irow2=irow2+2;
} //..
if (nt[func]==0) {gotoxy(ic2,irow2);cout<<"0";irow2=irow2+2;}
if (nt[func]<0)  {gotoxy(ic2,irow2);cout<<"?????";irow2=irow2+2;}
*/

sphbessel(2,orderL,PARTIC,LONGIT,First_kind);
legendre(1,orderL);multiply(0,1,2);nt[1]=-1;nt[2]=-1;
multiplyby(0,"B");gradient(1,2,3,0);nt[0]=-1;

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

add(1,1,6); add(2,2,7);
add(3,3,8); nt[6]=-1;nt[7]=-1;nt[8]=-1;

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

int incidenttransverse=0;
int incidentlongit=1-incidenttransverse;

sphbessel(6,orderL,MATRIX,LONGIT,Out_going);
multiplyby(6,"D");
if (incidentlongit==1)
{ //
gotoxy(8,20);cout<<"Incident: SPH LONGIT A = 1/f";
sphbessel(7,orderL,MATRIX,LONGIT,In_going);
multiplyby(7,"A"); add(6,6,7); nt[7]=-1;
} //
legendre(5,orderL);
multiply(4,6,5); nt[5]=-1; nt[6]=-1; gradient(5,6,7,4);
sphbessel(8,orderL,MATRIX,TRANSV,Out_going);
multiplyby(8,"E");
if (incidenttransverse==1)
{ //..
gotoxy(8,20);cout<<"Incident: SPH TRANSV A = 1/f";
sphbessel(11,orderL,MATRIX,TRANSV,In_going);
multiplyby(11,"A"); add(8,8,11); nt[11]=-1;
} //..
legendre(9,orderL); multiply(10,8,9); nt[8]=-1;
nt[9]=-1; multiplyby(10,"r"); settozero(8); settozero(9);
curl(11,12,13,10,8,9); nt[8]=-1; nt[9]=-1; nt[10]=-1;
curl(8,9,10,11,12,13); nt[11]=-1; nt[12]=-1; nt[13]=-1;
add(5,5,8); add(6,6,9); add(7,7,10); 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

sphbessel(14,orderL,MATER3,LONGIT,Out_going);
multiplyby(14,"F");
sphbessel(15,orderL,MATER3,LONGIT,In_going);
multiplyby(15,"G");
add(14,14,15);
nt[15]=-1;
legendre(15,orderL);
multiply(17,15,14); nt[15]=-1; nt[14]=-1;
gradient(14,15,16,17); nt[17]=-1;

sphbessel(18,orderL,MATER3,TRANSV,Out_going);
multiplyby(18,"H");
sphbessel(19,orderL,MATER3,TRANSV,In_going);
multiplyby(19,"I");
add(19,19,18); nt[18]=-1;
multiplyby(19,"r");
legendre(18,orderL);
multiply(17,18,19);
settozero(18);
settozero(19);
curl(20,21,22,17,18,19);
nt[17]=-1;
nt[18]=-1;
nt[19]=-1;
curl(17,18,19,20,21,22);
nt[20]=-1;
nt[21]=-1;
nt[22]=-1;
add(14,14,17);
add(15,15,18);
add(16,16,19);
// 14,15,16 is displacement in material #3
nt[17]=-1;
nt[18]=-1;
nt[19]=-1;
stress(17,18,19,14,15,16,MATER3);

subtract(20,1,14); // r continuity
subtract(21,5,14); // r continuity

nt[5]=-1;
nt[14]=-1;

subtract(22,2,15); // theta continuity
subtract(23,15,6); // theta continuity

nt[6]=-1;
nt[15]=-1;

subtract(24,8,17); // rr force balance
subtract(25,17,11); // rr force balance
nt[8]=-1;
nt[11]=-1;
nt[17]=-1;


//  8  9  10 = stress field inside particle
// 17 18  19  stress field inside material #3
// 11 12  13 = stress field in glass matrix

//  1  2   3 = displacement field inside
// 14 15  16 is displacement in material #3
//  5  6   7 = displacement field outside

subtract(26,9,18); // r-theta force balance
subtract(27,18,12); // r-theta force balance
nt[9]=-1;
nt[12]=-1;
nt[18]=-1;

beta=0.02;

//gotoxy(20,10);cout<<"add: "<<mosttermseveradd<<" terms";
//gotoxy(20,11);cout<<"ddr: "<<mosttermseverddr<<" terms";
//gotoxy(20,12);cout<<"cle: "<<mosttermsevercle<<" terms";

setcolor(15);
settextstyle(2,0,6);
outtextxy(bx0+0*sfx2-3, 440,"0");
//outtextxy(bx0+5*sfx2-9, 440,"5");
outtextxy(bx0+10*sfx2-9,440,"10");
//outtextxy(bx0+15*sfx2-9,440,"15");
outtextxy(bx0+20*sfx2-9,440,"20");
//outtextxy(bx0+25*sfx2-9,440,"25");
outtextxy(bx0+30*sfx2-9,440,"30");
outtextxy(bx0+40*sfx2-9,440,"40");
outtextxy(bx0+50*sfx2-9,440,"50");
outtextxy(bx0+60*sfx2-9,440,"60");
outtextxy(bx0+70*sfx2-9,440,"70");
settextstyle(2,0,6);
outtextxy(bx0+50+2,459,"Raman shift (cm  )");
settextstyle(2,0,5);
outtextxy(bx0+195+1,459-5,"-1");
settextstyle(2,1,6);
//setcolor(13);
//outtextxy(-1,190,"abs(B) (long)");
setcolor(15);
outtextxy(-1,340,"abs(displ.):");
setcolor(14);
outtextxy(-1,185+90,"radial");
setcolor(12);
outtextxy(-1,130+90,"theta");
setcolor(10);
outtextxy(-1,90+90,"phi");

tm1=clock(); tm4=0;

float old1=0,old2=0,old3=0,old4=0,oldx=-1;
LB555:
//beta = omega * R / CT
//omega = beta * CT / R;

if (beta<betamax)
{ //
beta=beta+0.03;
} //
else
{ //..
makescr();
return;
} //..

//gotoxy(5,13);cout<<partrad_nm<<" nm radius";
gotoxy(25,13);cout<<2*partrad_nm<<" nm diam";
frequencyomega=beta*css[1][0]/partrad_nm;
gotoxy(25,15);cout<<thick3<<" nm thick";

for (imater=0;imater<NMATER;imater++)
{ //.
kwv[0][imater]=frequencyomega/css[0][imater];
kwv[1][imater]=frequencyomega/css[1][imater];
} //.

evaluate_func(0+8,20,partrad_nm);
evaluate_func(1+8,22,partrad_nm);
evaluate_func(2+8,24,partrad_nm);
evaluate_func(3+8,26,partrad_nm);
evaluate_func(4+8,21,(partrad_nm+thick3));
evaluate_func(5+8,23,(partrad_nm+thick3));
evaluate_func(6+8,25,(partrad_nm+thick3));
evaluate_func(7+8,27,(partrad_nm+thick3));

evaluate_func(4,1,partrad_nm);
cleanup(4);
evaluate_func(5,2,partrad_nm);
cleanup(5);
evaluate_func(6,3,partrad_nm);
cleanup(6);

// Most of the time is spent here, in cleanup:
cleanup(0+8);
cleanup(1+8);
cleanup(2+8);
cleanup(3+8);
cleanup(4+8);
cleanup(5+8);
cleanup(6+8);
cleanup(7+8);

complex crb4,c3,c4,c5;
complex c31;
int iv1,iv2,iv3;

// Eliminate variables ...E,D,C
int topeqno=7+8,noeqused=8;
for (iv1='C'+noeqused-2;iv1>='C';iv1--)
{ // iv1
iv2=topeqno+'C'-(iv1+1); //iv2 = equation to eliminate
for (iv3=iv2+1;iv3<=topeqno;iv3++)
{ // iv3
crb4=coefof(iv2,iv1);
c5=coefof(iv3,iv1);
if (crb4==0) {cout<<"* Err crb4==0 1 *";getch();exit(0);}
crb4=-c5/crb4;
multiplyby(iv2,crb4);
add(iv3,iv3,iv2);
} // iv3
} // iv1

// Find values of all variables assuming A=1.
//varval[0]=1; // "A"   OLD VERSIONS
if (incidenttransverse==1)
{ // Transverse spheroidal mode: (from md30.htm)
gotoxy(60,16);cout<<"Using tran sph";
float L=orderL,w=frequencyomega,pi=3.1415926;
float dw=1.0,kBT=1.0;
float Ctm=css[1][1],rhom=rho[1];
varval[0]=sqrt(dw*kBT*(2*L+1)/(L*(L+1)*pi*pi*rhom*Ctm*w*w)); // "A"
} //
if (incidenttransverse==0)
{ // Longitudinal spheroidal mode (from md30.htm)
gotoxy(60,15);cout<<"Using long sph";
float L=orderL,w=frequencyomega,pi=3.1415926;
float dw=1.0; // width of frequency bin
float kBT=1.0; // kB times temperature T
float Clm=css[0][1],rhom=rho[1];
varval[0]=sqrt(dw*kBT*(2*L+1)/(pi*pi*rhom*Clm*w*w)); // "A"
} //

for (iv1=1;iv1<=noeqused;iv1++)
{ // iv1
iv3=topeqno+1-iv1;
varval[iv1]=0;
for (iv2=0;iv2<iv1;iv2++)
{ //. iv2
varval[iv1]+=varval[iv2]*coefof(iv3,'A'+iv2);
} //. iv2
varval[iv1]*=-1;
c31=coefof(iv3,'A'+iv1);
if (c31==0) {cout<<"* Err c31==0 1 *";getch();exit(0);}
varval[iv1]=varval[iv1]/c31;
} // iv1

double wavenum;
wavenum = 0.01*beta*css[1][0]/(0.2998*2.0*3.1415926*partrad_nm);

double absrda,absthda,absphda;
{
complex c40,c41,c42,c43,rda,thda,phda;
rda=coefof(4,'B')*varval[1]+coefof(4,'C')*varval[2];
absrda=abs(rda);
thda=coefof(5,'B')*varval[1]+coefof(5,'C')*varval[2];
absthda=abs(thda);
phda=coefof(6,'B')*varval[1]+coefof(6,'C')*varval[2];
absphda=abs(phda);
}

float sfy=10;

float new0,new1,new2,new3,new4;
new1=by0-500000000.*absrda*sfy;
new2=by0-500000000.*absthda*sfy;
new3=by0-500000000.*absphda*sfy;
if (oldx>0)
{ //
setcolor(14);line(oldx,old1,bx0+sfx2*wavenum,new1);
setcolor(12);line(oldx,old2,bx0+sfx2*wavenum,new2);
setcolor(10);line(oldx,old3,bx0+sfx2*wavenum,new3);
} //

old3=new3;
old2=new2;
old1=new1;
oldx=bx0+sfx2*wavenum;

//putpixel(bx0+sfx*beta,479-r67c*sfy,14);
//putpixel(bx0+sfx*beta,479-r67d*sfy,15);
//putpixel(bx0+sfx*beta,479-(r67c*r67c+1.8*r67d*r67d)*sfy,13); // for L=2
//putpixel(bx0+sfx*beta,479-(r67c*r67c+0.7*r67d*r67d)*sfy,13); // for L=1
if (!kbhit()) goto LB555;
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<<"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 (mnt[func2]<nt[func1])
 {cout<<"* Err evalobj 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)
{ //
r2=kwv[ipol][imat];
r2=pow(r2,ngetm(i7+5,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(i7+2,it1)!=0)
{ //
r2=sin(kwv[ipol][imat]*partrad_nm);
r2=pow(r2,ngetm(i7+2,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(i7+4,it1)!=0)
{ //
r2=cos(kwv[ipol][imat]*partrad_nm);
r2=pow(r2,ngetm(i7+4,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
} // 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)
/*
complex constantcoef(int func1)
{ //    constantcoef
complex cr3=0.0;
int i7,flag1;
for (int it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //..
flag1=1;
for (i7=3;i7<NFACT;i7=i7+3) if (ngetm(i7,it1)!=0) flag1=0;
if (flag1==1) cr3=complex(cr[it1],ci[it1]);
} //..
return cr3;
} //    constantcoef
*/
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

void makescr(void)
{ // makescr
gotoxy(60,1);cout<<"Make .scr ";
gotoxy(60,2);cout<<"file? ";
gotoxy(60,3);cout<<"(y/n/q)  ";
char ch678;
ch678=getch();
if (ch678=='q') exit(0);
if (ch678!='y') return;
char scrfilename[40];
FILE *scrfout;
gotoxy(60,4);cout<<"Filename";
gotoxy(55,8);cout<<".scr extension is needed";
gotoxy(60,5);cout<<"? ";cin>>scrfilename;

scrfout=fopen(scrfilename,"wb");
unsigned char ch1b,ch2b,ch3b,ch4b,ch5,ch6;
int ipybot=478;
int ipytop=ipybot-300;
int scrhgt=ipybot-ipytop;
int scrx0=1;
int scry0=ipytop;
setcolor(15);
int scrwid=300;
line(scrx0,scry0-1,scrwid+scrwid-1,scry0-1);
line(scrx0,scry0+scrhgt,scrwid+scrwid-1,scry0+scrhgt);
ch3b=scrwid/256;
ch1b=scrwid%256;
ch4b=scrhgt/256;
ch2b=scrhgt%256;
fprintf(scrfout,"%c%c%c%c%c%c%c%c%c%c%c%c%c",'s','c','r',ch3b,ch1b,ch4b,ch2b,'0','0','0','0',13,10);
int ipixcol,ix,iy;
gotoxy(65,6);cout<<"saving ";
for (iy=scry0;iy<scry0+scrhgt;iy++)
{ // iy
gotoxy(65,7);cout<<iy;
for (ix=scrx0;ix<scrx0+scrwid;ix++)
{ //
ipixcol=getpixel(ix,iy);
//putpixel(ix-100,iy-100,ipixcol);
char ch1;
ch1=ipixcol+'A';
fprintf(scrfout,"%c",ch1);
} //
fprintf(scrfout,"%c%c",13,10);
} // iy
fclose(scrfout);
return;
} // makescr